function periphOutput = takanen2013periphery(insig,fs,outputPlot)
%TAKANEN2013PERIPHERY Process the input through the model of periphery         
%   Usage: periphOutput = takanen2013periphery(insig,fs,outputPlot);
%          periphOutput = takanen2013periphery(insig,fs);
%
%   Input parameters:
%        insig      : binaural input signal to be processed. Optionally,
%                     the output of the cochlear model by Verhulst et. al.
%                     2012 can be used as well
%        fs         : sampling rate
%        outputPlot : boolean value that defines whether the periphery
%                     model output at different frequency bands are plotted
%
%   Output parameters:
%        periphOutput : Structure consisting of the following elements
%
%                       `periphOutput.left`
%                          Left ear "where" stream output
% 
%                       `periphOutput.right`  
%                           Right ear "where" stream output
%
%                       `periphOutput.fc`
%                           Characteristic frequencies
%
%                       `periphOutput.ventralLeft`
%                           Left hemisphere "what" stream output 
%
%                       `periphOutput.ventralLeft`
%                           Right hemisphere "what" stream output
%
%   This function processes the binaural input signal through the model of
%   periphery presented by Takanen, Santala, Pulkki 2013, the model which
%   consists of a nonlinear time-domain model of cochlea and model of
%   cochlear nucleus. The processing contains the following steps:
%
%   1) the binaural input signal is processed with the nonlinear time-
%      domain model of cochlea by Verhulst et. al. 2013 to obtain the
%      velocity of basilar membrane movement at different positions
%
%   2) the obtained velocity information is half-wave rectified
%
%   3) the half-waves are replaced with Gaussian pulses centered around
%      the local maxima of the half-waves
%
%   4) the frequency-dependent delays of the cochlea model are
%      compensated
%
%   See also: takanen2013, verhulst2012

%   References: takanen2013a pulkki2009 verhulst2012
%
%   AUTHOR: Marko Takanen, Olli Santala, Ville Pulkki
%
%   COPYRIGHT (C) 2013 Aalto University
%                      School of Electrical Engineering
%                      Department of Signal Processing and Acoustics
%                      Espoo, Finland

%% ------ Check the input arguments ---------------------------------------
if nargin<3
    outputPlot =0;
end

%specify the characteristic frequencies of the model
fc = erbtofreq(4:39);
spl = 60; % sound pressure level of the input signal 
%% ------ Process the binaural input with the cochlear model --------------

if isstruct(insig)
    % input to the periphery model can be also a output of the Verhulst
    % cochlear model
    cochlear = insig;
else
    norm_factor=max(abs(insig(:,1)));     %first channel as reference for rms normalization
    insig(:,1)=insig(:,1)./norm_factor;
    insig(:,2)=insig(:,2)./norm_factor;
    [V,Y,E,CF] = verhulst2012(insig,fs,fc,[spl spl]);
    cochlear.velocityLeft=V(:,:,1);
    cochlear.velocityRight=V(:,:,2);
end

%load a vector of frequency specific delays computed for the Verhulst model
%outputs at fc
load takanen2013_cochleardelays -mat
cochlear.delaysV = velocitydelays;

%use the velocity of the basilar membrane movement as input to CN model
processed.left = cochlear.velocityLeft;processed.right = cochlear.velocityRight;
[nVals nBands] = size(processed.left);

%% ------ Half-wave rectification -----------------------------------------
processed.left= processed.left.*(processed.left >0);
processed.right= processed.right.*(processed.right >0);
periphOutput.ventralLeft = processed.right;
periphOutput.ventralRight = processed.left;

%% ------ Impulse generation ----------------------------------------------
%search of the local maximas of each half-wave separately for left and 
%right ear signals
for chanInd=1:nBands
    left = processed.left(:,chanInd);
    right = processed.right(:,chanInd);
    processed.left(:,chanInd) = zeros(nVals,1);
    processed.right(:,chanInd) = zeros(nVals,1);

    %start end end points for each half-wave
    startingPoints = strfind((left'>0),[0 1]);
    endpoints = [strfind((left'>0),[1 0])+1 length(left)];
    if (isempty(startingPoints) && ~isempty(endpoints))
        startingPoints = 1;
    else
        if startingPoints(1)>=endpoints(1)
            startingPoints = [1 startingPoints];
        end
    end
    
    rmsVals=zeros(size(startingPoints));
    locations = rmsVals;
    nsamples = endpoints(1:length(startingPoints))-startingPoints+1;
    %compute the rms values of each half-wave and position it at the local
    %maximum
    for locInd=1:length(startingPoints)
        rmsVals(locInd)= norm(left(startingPoints(locInd):endpoints(locInd)))/sqrt(nsamples(locInd));
        [unnecessary, locations(locInd)] = max(left(startingPoints(locInd):endpoints(locInd)));
    end
    processed.left(locations+startingPoints-1,chanInd) = rmsVals';

    %processing of the right channel
    
    %start end end points for each half-wave
    startingPoints = strfind((right'>0),[0 1]);
    endpoints = [strfind((right'>0),[1 0])+1 length(right)];
    if (isempty(startingPoints) && ~isempty(endpoints))
        startingPoints = 1;
    else
        if startingPoints(1)>=endpoints(1)
            startingPoints = [1 startingPoints];
        end
    end
    %compute the rms values of each half-wave and position it at the local
    %maximum
    rmsVals= zeros(size(startingPoints));locations = rmsVals;
    nsamples = endpoints(1:length(startingPoints))-startingPoints+1;
    for locInd=1:length(startingPoints)
        rmsVals(locInd)= norm(right(startingPoints(locInd):endpoints(locInd)))/sqrt(nsamples(locInd));
        [unnecessary, locations(locInd)] = max(right(startingPoints(locInd):endpoints(locInd)));
    end
    processed.right(locations+startingPoints-1,chanInd) = rmsVals';

end
%% ------ Convolution with Gaussian window-functions ----------------------
% the width of the Gaussian window in the peripheral hearing model depends
% on the center frequency
N = zeros(size(fc));
N(fc<800) = round((2./fc(fc<800))*fs);
indMid = find(((800<=fc).*(fc<=2800))==1);N(indMid) = round(0.0024*(0.6+0.4*fc(indMid)./800)*fs);
N(fc>2800) = round(0.0048*fs);
% the constant alpha is set to 20
alpha=20;

periphOutput.left = zeros(nVals,nBands);periphOutput.right = periphOutput.left;
for chanInd=1:nBands
    left = processed.left(:,chanInd);
    right = processed.right(:,chanInd);
    n = -N(chanInd)/2:1:N(chanInd)/2;
    winFunction = (exp(-0.5*(alpha*n/(N(chanInd)/2)).^2))';
    %convolution with the gaussian window function
    left = (1/sum(winFunction))*conv(left,winFunction,'same');
    right = (1/sum(winFunction))*conv(right,winFunction,'same');

    %% compensation for the cochlear model delays
    periphOutput.left(:,chanInd) = [left(cochlear.delaysV(chanInd)+1:end);zeros(cochlear.delaysV(chanInd),1)];
    periphOutput.right(:,chanInd) = [right(cochlear.delaysV(chanInd)+1:end);zeros(cochlear.delaysV(chanInd),1)];
    periphOutput.ventralLeft(:,chanInd) = [periphOutput.ventralLeft(cochlear.delaysV(chanInd)+1:end,chanInd);zeros(cochlear.delaysV(chanInd),1)];
    periphOutput.ventralRight(:,chanInd) = [periphOutput.ventralRight(cochlear.delaysV(chanInd)+1:end,chanInd);zeros(cochlear.delaysV(chanInd),1)];
end

periphOutput.fc = fc;
%% ------ Plot the periphery model output, if desired ---------------------
if outputPlot
    figure(90);
    g(1) = subplot(6,1,1);plot((0:length(insig(:,1))-1)./(fs/1000),insig(:,1));ylabel('Input');set(gca,'yTick',[]);
    g(2) = subplot(6,1,2);plot((0:size(periphOutput.left,1)-1)./(fs/1000),periphOutput.left(:,4));title([num2str(round(fc(4))) ' Hz']);
    g(3) = subplot(6,1,3);plot((0:size(periphOutput.left,1)-1)./(fs/1000),periphOutput.left(:,9));title([num2str(round(fc(9))) ' Hz']);
    g(4) = subplot(6,1,4);plot((0:size(periphOutput.left,1)-1)./(fs/1000),periphOutput.left(:,16));title([num2str(round(fc(16))) ' Hz']);
    g(5) = subplot(6,1,5);plot((0:size(periphOutput.left,1)-1)./(fs/1000),periphOutput.left(:,22));title([num2str(round(fc(22))) ' Hz']);
    g(6) = subplot(6,1,6);plot((0:size(periphOutput.left,1)-1)./(fs/1000),periphOutput.left(:,28));title([num2str(round(fc(28))) ' Hz']);
    xlabel('Time [ms]');
    linkaxes(g,'x');
    clear g
end