MAP

% clear global -regexp OMEParams DRNLParams IHC_cilia_RPParams IHCpreSynapseParams

% clear global -regexp  AN_IHCsynapseParams MacGregorParams MacGregorMultiParams

% clear global -regexp dt ANdt  savedBFlist saveAN_spikesOrProbability saveMAPparamsName...

%     savedInputSignal OMEextEarPressure TMoutput OMEoutput ARattenuation ...

%     DRNLoutput IHC_cilia_output IHCrestingCiliaCond IHCrestingV...

%     IHCoutput ANprobRateOutput ANoutput savePavailable ANtauCas  ...

%     CNtauGk CNoutput  ICoutput ICmembraneOutput ICfiberTypeRates ...

%     MOCattenuation

global dt ANdt  savedBFlist saveAN_spikesOrProbability saveMAPparamsName...

    savedInputSignal OMEextEarPressure TMoutput OMEoutput ARattenuation ...

    DRNLoutput IHC_cilia_output IHCrestingCiliaCond IHCrestingV...

    IHCoutput ANprobRateOutput ANoutput savePavailable ANtauCas  ...

    CNtauGk CNoutput  ICoutput ICmembraneOutput ICfiberTypeRates ...

    MOCattenuation 

%   global ANdt ICoutput ANprobRateOutput

%  'spike' matrices using a lower sample rate (see ANdt).

%  They are represented with dt=1/sampleRate (not ANdt)

ANspeedUpFactor=AN_IHCsynapseParams.ANspeedUpFactor;  % e.g.5 times

% AN (spikes) is computed at a lower sample rate when spikes required

%  so it has a reduced segment length (see 'ANspeeUpFactor' above)

ANdt=dt*ANspeedUpFactor;

OMEhighPassHighCutOff=OMEParams.OMEstapesLPcutoff;

CtBM=10e-9*10^(DRNLParams.CtBMdB/20);

CtS=CtBM/DRNLa;

lengthAbsRefractory= round(AN_refractory_period/ANdt);

AN_ydt= repmat(AN_IHCsynapseParams.y*ANdt, nANfibers,1);

AN_ldt= repmat(AN_IHCsynapseParams.l*ANdt, nANfibers,1);

AN_xdt= repmat(AN_IHCsynapseParams.x*ANdt, nANfibers,1);

AN_rdt= repmat(AN_IHCsynapseParams.r*ANdt, nANfibers,1);

t=ANdt:ANdt:5*CNspikeToCurrentTau;

% figure(98), plot(t,CNalphaFunction)

% epochsPerSpike=round(ICspikeWidth/ANdt);

t=ANdt:ANdt:3*ICspikeToCurrentTau;

        abs_x = abs(nonlinOutput);

        y(abs_x<CtS)   = DRNLa * nonlinOutput(abs_x<CtS);

        y(abs_x>=CtS)  = sign(nonlinOutput(abs_x>=CtS)) * DRNLa * CtS .* ...

            exp(   DRNLc * log(  abs_x(abs_x>=CtS)/CtS  )   );

        nonlinOutput=y;

        % % Boltzmann compression function

        % y=(nonlinOutput * DRNLa*100000);

        % holdY=y;

        % y=abs(y);

        % s=10; u=0.0;

        % x=1./(1+exp(-(y-u)/s))-0.5;

        % nonlinOutput=sign(nonlinOutput).*x/10000;

        % if segmentStartPTR==10*segmentLength+1

        %     figure(90)

        %     plot(holdY,'b'), hold on

        %     plot(nonlinOutput, 'r'), hold off

        %     ylim([-1e-5 1e-5])

        %     pause(1)

        % end

        %       second filter removes distortion products

               % we need a running account of cumProb=II(1-p(t))

            %% Estimate efferent effect:  0 < ARattenuation > 1

            %  acoustic reflex

                % Nick's new code

%                 for idx=1:nBFs

%                     [smoothedRates, MOCprobBoundary{idx}] = ...

%                         filter(MOCfilt_b, MOCfilt_a, rates(idx,:), ...

%                         MOCprobBoundary{idx});

% %                     smoothedRates=smoothedRates-MOCrateThresholdProb;

% %                     smoothedRates(smoothedRates<0)=0;

% %                     x =  (1- smoothedRates* rateToAttenuationFactorProb); %ORIGINAL 

%  

%                     %NEW !!!

%                     x = -20*log10(  max(smoothedRates/MOCrateThresholdProb,1)  )*rateToAttenuationFactorProb; %dB attenuation

%                     x = 10.^(x/20);

%                     x = max(x,10^(-35/20));    

% %                     %ALSO - filter at the end - this will stop rapid attack

% %                     %and slow decay

% %                     [x, MOCprobBoundary{idx}] = ...

% %                         filter(MOCfilt_b, MOCfilt_a, x, ...

% %                         MOCprobBoundary{idx});

%                     

%                                         

%                     MOCattenuation(idx,segmentStartPTR:segmentEndPTR)= ...

%                         x;                                                            

%                 end

                releaseProb=vesicleReleaseRate(:,t)*ANdt;

            % plotInstructions.displaydt=ANdt;

            % input is from AN so ANdt is used throughout

                    % ANpsth has a bin width of ANdt

                    CNtimeSinceLastSpike=CNtimeSinceLastSpike-ANdt;

                        CN_Gk/CN_cap).*(CN_Ek-CN_E))*ANdt;

                    dGk=-CN_Gk*ANdt./tauGk + CN_b*s;

                    dTh=-(CN_Th-CN_Th0)*ANdt/CN_tauTh + CN_c*s;

                    CNtimeSinceLastSpike=CNtimeSinceLastSpike-ANdt;

                        (CN_Gk/CN_cap).*(CN_Ek-CN_E))*ANdt;

                    dGk=-CN_Gk*ANdt./tauGk +CN_b*s;

            % plotInstructions.displaydt=ANdt;

                            (IC_Gk/IC_cap).*(IC_Ek-IC_E))*ANdt;

                        dGk=-IC_Gk*ANdt/IC_tauGk +IC_b*s;

                            .*(IC_Ek-IC_E))*ANdt;

                        dGk=-IC_Gk*ANdt/IC_tauGk +IC_b*s;

                        dTh=-(IC_Th-Th0)*ANdt/IC_tauTh +s*IC_c;

                        /(nCellsPerTau*ANdt);

                % ARis based on LSR units. LSR channels are 1:nBF

                    ARAttSeg=mean(ICspikes(1:nBFs,:),1)/ANdt;

%                     ARAttSeg(ARAttSeg<0)=0;   % prevent negative strengths

                    % scale up to dt from ANdt

                    ARattenuation(ARattenuation<0)=0.001;

                rates=ICspikes(HSRbegins:end,:)/ANdt;

                    % expand timescale back to model dt from ANdt

                % plotInstructions.displaydt=ANdt;

function MAPrunner(MAPparamsName, AN_spikesOrProbability, ...

    signalCharacteristics, paramChanges, showMapOptions)

dbstop if error

restorePath=path;

addpath (['..' filesep 'MAP'],    ['..' filesep 'wavFileStore'], ...

    ['..' filesep 'utilities'])

%% #3 pure tone, harmonic sequence or speech file input

signalType= signalCharacteristics.type;

sampleRate= signalCharacteristics.sampleRate;

duration=signalCharacteristics.duration;                 % seconds

rampDuration=signalCharacteristics.rampDuration;              % raised cosine ramp (seconds)

beginSilence=signalCharacteristics.beginSilence;               

endSilence=signalCharacteristics.endSilence;                  

toneFrequency= signalCharacteristics.toneFrequency;            % or a pure tone (Hz)

leveldBSPL=signalCharacteristics.leveldBSPL;

BFlist=-1;

%% Generate stimuli

switch signalType

    case 'tones'

        % Create pure tone stimulus

        dt=1/sampleRate; % seconds

        time=dt: dt: duration;

        inputSignal=sum(sin(2*pi*toneFrequency'*time), 1);

        amp=10^(leveldBSPL/20)*28e-6;   % converts to Pascals (peak)

        inputSignal=amp*inputSignal;

        % apply ramps

        % catch rampTime error

        if rampDuration>0.5*duration, rampDuration=duration/2; end

        rampTime=dt:dt:rampDuration;

        ramp=[0.5*(1+cos(2*pi*rampTime/(2*rampDuration)+pi)) ...

            ones(1,length(time)-length(rampTime))];

        inputSignal=inputSignal.*ramp;

        ramp=fliplr(ramp);

        inputSignal=inputSignal.*ramp;

        % add silence

        intialSilence= zeros(1,round(beginSilence/dt));

        finalSilence= zeros(1,round(endSilence/dt));

        inputSignal= [intialSilence inputSignal finalSilence];

    case 'file'

        %% file input simple or mixed

        [inputSignal sampleRate]=wavread(fileName);

        dt=1/sampleRate;

        inputSignal=inputSignal(:,1);

        targetRMS=20e-6*10^(leveldBSPL/20);

        rms=(mean(inputSignal.^2))^0.5;

        amp=targetRMS/rms;

        inputSignal=inputSignal*amp;

        intialSilence= zeros(1,round(0.1/dt));

        finalSilence= zeros(1,round(0.2/dt));

        inputSignal= [intialSilence inputSignal' finalSilence];

end

%% run the model

MAP1_14(inputSignal, sampleRate, BFlist, ...

    MAPparamsName, AN_spikesOrProbability, paramChanges);

%% the model run is now complete. Now display the results

UTIL_showMAP(showMapOptions, paramChanges)

path(restorePath)

function [P, BFlist, sacf, boundaryValue] = ...

    filteredSACF(inputSignalMatrix, method, params)

% UTIL_filteredSACF computes within-channel, running autocorrelations (acfs)

%  and finds the sum across channels (sacf).

%  The SACF is smoothed to give the 'p function' (P).

%  	method.dt:			the signal sampling interval in seconds

%   method.segmentNo:

%   method.nonlinCF

%   

%   filteredSACFParams.lags: an array of lags to be computed (seconds)

%   filteredSACFParams.acfTau: time constant (sec) of the running acf calculations

%     if acfTau>1 it is assumed that Wiegrebe'sacf method

%   filteredSACFParams.Lambda: time constant  for smoothing thsacf to make P

%   filteredSACFParams.lagsProcedure identifies a strategy for omitting some lags.

%    Options are: 'useAllLags', 'omitShortLags', or 'useBernsteinLagWeights'

%   filteredSACFParams.usePressnitzer applies lower weights longer lags

%   parafilteredSACFParamsms.plotACFs (=1) creates movie of acf matrix (optional)

%

%  	P:	P function 	(time x lags), a low pass filtered version of sacf.

%  method: updated version of input method (to include lags used)

%   sacf:  	(time x lags)

boundaryValue=[];

% list of BFs must be repeated is many fiber types used

BFlist=method.nonlinCF;

nfibertypes=nChannels/length(BFlist);

BFlist=repmat(BFlist',2,1)';

dt=method.dt;

% adjust sample rate, if required

if isfield(params,'dt')

    [inputSignalMatrix dt]=UTIL_adjustDT(params.dt, method.dt, inputSignalMatrix);

    method.dt=dt;

end

% create acf movie

% disp(['lag procedure= ''' lagsProcedure ''''])

       % no action required lagWeights set above

    case 'useBernsteinLagWeights'

        lagWeights=BernsteinLags(BFlist, lags)';

        error ('acf: params.lagProcedure not recognised')

% Establish matrix of lag time constants 

else                      % use Wiegrebe method: tau= 2*lag (for example)

% P function lowpass filter decay (only one value needed)

P=zeros(nLags,inputLength+1);   % P must match segment length +1

sacf=zeros(nLags,inputLength);

    acf=zeros(nChannels,nLags);

    % boundaryValue picks up from a previous calculation

    acf=params.boundaryValue{1};

    P(: , 1)=params.boundaryValue{2}; % NB first value is last value of previous segment

    buffer=params.boundaryValue{3};

    % acf is a continuously changing channels x lags matrix

    % sacf is the vertical summary of acf ( a vector) and all values are kept and returned

    % P is the smoothed version of sacf and all values are kept and returned

    % E.g. if max(lags) is .04 then this is when the ACf will begin.

    %            AN                                       AN delayed                             weights               filtering

    acf= (repmat(inputSignalMatrix(:, timePointer), 1, nLags) .* inputSignalMatrix(:, timePointer-lagPointers)).*lagWeights *dt + acf.* acfDecayFactors;

    x=(mean(acf,1)./acfTaus)';

%     disp(num2str(x'))

    sacf(:,timeCounter)=x;

    P(:,timeCounter+1)=sacf(:,timeCounter)*(1-pDecayFactor)+P(:,timeCounter)*pDecayFactor;

    % plot at intervals of 200 points

    if plotACF && ~mod(timePointer,params.plotACFsInterval)

        %       mark cursor on chart to signal progress

        % this assumes that the user has already plotted

        % the signal in subplot(2,1,1) of figure (13)

        figure(13)

        hold on

        subplot(4,1,1)

        a =ylim;

        plot([time time], [a(1) a(1)+(a(2)-a(1))/4]) % current signal point marker

        %         plot ACFs one per channel

        subplot(2,1,2), cla

        cascadePlot(acf, lags, BFlist)

        xlim([min(lags) max(lags)])

        %         set(gca,'xscale','log')

        title(num2str(method.dt*timePointer))

        ylabel('BF'), xlabel('period (lag)')

        %         plot SACF

        subplot(4,1,2), hold off

        plot(lags,sacf(:,timeCounter)-min(sacf(:,timeCounter)))

        biggestSACF=max(biggestSACF, max(sacf(:,timeCounter)));

        if biggestSACF>0, ylim([0 biggestSACF]), end

        %         set(gca,'xscale','log')

        title('SACF')

P=P(:,1:end-1);  % correction for speed up above

    [a nTimePoints]=size(P);

    P=P.*PressnitzerWeights';

    sacf=sacf.*PressnitzerWeights';

method.acfLags=lags;

method.filteredSACFdt=dt;

boundaryValue{1}=acf(:,end-nLags+1:end);

boundaryValue{2}=P(:,end);

boundaryValue{3} = inputSignalMatrix(:, end-max(lagPointers)+1 : end);

method.displaydt=method.filteredSACFdt;

% if ~isfield(params, 'plotUnflteredSACF'), params.plotUnflteredSACF=0; end

% if method.plotGraphs

%     method.plotUnflteredSACF=params.plotUnflteredSACF;

%     if ~method.plotUnflteredSACF

%         method=filteredSACFPlot(P,method);

%     else

%         method=filteredSACFPlot(SACF,method);

%     end

% end

% ------------------------------------------------  plotting ACFs

function cascadePlot(toPlot, lags, BFs)

% % useful code

[nChannels nLags]=size(toPlot);

% cunning code to represent channels as parallel lines

[nRows nCols]=size(toPlot);

if nChannels>1

    % max(toPlot) defines the spacing between lines

    a=max(max(toPlot))*(0:nRows-1)';

    % a is the height to be added to each channel

    peakGain=10;

    % peakGain emphasises the peak height

    x=peakGain*toPlot+repmat(a,1,nCols);

    x=nRows*x/max(max(x));

else

    x=toPlot;                            % used when only the stimulus is returned

end

plot(lags, x','k')

ylim([0 nRows])

function  MAP1_14AP(inputSignal, sampleRate, BFlist, MAPparamsName, ...

    AN_spikesOrProbability, paramChanges)

% To test this function use test_MAP1_14 in this folder

%

% All arguments are mandatory.

%

%  BFlist is a vector of BFs but can be '-1' to allow MAPparams to choose

%  MAPparamsName='Normal';          % source of model parameters

%  AN_spikesOrProbability='spikes'; % or 'probability'

%  paramChanges is a cell array of strings that can be used to make last

%   minute parameter changes, e.g., to simulate OHC loss

%  e.g.  paramChanges{1}= 'DRNLParams.a=0;'; % disable OHCs

%  e.g.  paramchanges={};                    % no changes

% The model parameters are established in the MAPparams<***> file

%  and stored as global

restorePath=path;

addpath (['..' filesep 'parameterStore'])

CONVOLUTION_CHANGE_TEST = 0; %for debug

global OMEParams DRNLParams IHC_cilia_RPParams IHCpreSynapseParams

global AN_IHCsynapseParams MacGregorParams MacGregorMultiParams

% All of the results of this function are stored as global

global dt ANdt  savedBFlist saveAN_spikesOrProbability saveMAPparamsName...

    savedInputSignal OMEextEarPressure TMoutput OMEoutput ARattenuation ...

    DRNLoutput IHC_cilia_output IHCrestingCiliaCond IHCrestingV...

    IHCoutput ANprobRateOutput ANoutput savePavailable tauCas  ...

    CNoutput  ICoutput ICmembraneOutput ICfiberTypeRates MOCattenuation

% Normally only ICoutput(logical spike matrix) or ANprobRateOutput will be

% needed by the user; so the following will suffice

%   global ANdt ICoutput ANprobRateOutput

% Note that sampleRate has not changed from the original function call and

%  ANprobRateOutput is sampled at this rate

% However ANoutput, CNoutput and IC output are stored as logical

%  'spike' matrices using a lower sample rate (see ANdt).

% When AN_spikesOrProbability is set to probability,

%  no spike matrices are computed.

% When AN_spikesOrProbability is set to 'spikes',

%  no probability output is computed

% Efferent control variables are ARattenuation and MOCattenuation

%  These are scalars between 1 (no attenuation) and 0.

%  They are represented with dt=1/sampleRate (not ANdt)

%  They are computed using either AN probability rate output

%   or IC (spikes) output as approrpriate.

% AR is computed using across channel activity

% MOC is computed on a within-channel basis.

if nargin<1

    error(' MAP1_14 is not a script but a function that must be called')

end

if nargin<6

    paramChanges=[];

end

% Read parameters from MAPparams<***> file in 'parameterStore' folder

% Beware, 'BFlist=-1' is a legitimate argument for MAPparams<>

%  It means that the calling program allows MAPparams to specify the list

cmd=['method=MAPparams' MAPparamsName ...

    '(BFlist, sampleRate, 0, paramChanges);'];

eval(cmd);

BFlist=DRNLParams.nonlinCFs;

% save as global for later plotting if required

savedBFlist=BFlist;

saveAN_spikesOrProbability=AN_spikesOrProbability;

saveMAPparamsName=MAPparamsName;

dt=1/sampleRate;

duration=length(inputSignal)/sampleRate;

% segmentDuration is specified in parameter file (must be >efferent delay)

segmentDuration=method.segmentDuration;

segmentLength=round(segmentDuration/ dt);

segmentTime=dt*(1:segmentLength); % used in debugging plots

% all spiking activity is computed using longer epochs

ANspeedUpFactor=AN_IHCsynapseParams.ANspeedUpFactor;  % e.g.5 times

% inputSignal must be  row vector

[r c]=size(inputSignal);

if r>c, inputSignal=inputSignal'; end       % transpose

% ignore stereo signals

inputSignal=inputSignal(1,:);               % drop any second channel

savedInputSignal=inputSignal;

% Segment the signal

% The sgment length is given but the signal length must be adjusted to be a

% multiple of both the segment length and the reduced segmentlength

[nSignalRows signalLength]=size(inputSignal);

segmentLength=ceil(segmentLength/ANspeedUpFactor)*ANspeedUpFactor;

% Make the signal length a whole multiple of the segment length

nSignalSegments=ceil(signalLength/segmentLength);

padSize=nSignalSegments*segmentLength-signalLength;

pad=zeros(nSignalRows,padSize);

inputSignal=[inputSignal pad];

[ignore signalLength]=size(inputSignal);

% AN (spikes) is computed at a lower sample rate when spikes required

%  so it has a reduced segment length (see 'ANspeeUpFactor' above)

% AN CN and IC all use this sample interval

ANdt=dt*ANspeedUpFactor;

reducedSegmentLength=round(segmentLength/ANspeedUpFactor);

reducedSignalLength= round(signalLength/ANspeedUpFactor);

%% Initialise with respect to each stage before computing

%  by allocating memory,

%  by computing constants

%  by establishing easy to read variable names

% The computations are made in segments and boundary conditions must

%  be established and stored. These are found in variables with

%  'boundary' or 'bndry' in the name

%% OME ---

% external ear resonances

OMEexternalResonanceFilters=OMEParams.externalResonanceFilters;

[nOMEExtFilters c]=size(OMEexternalResonanceFilters);

% details of external (outer ear) resonances

OMEgaindBs=OMEexternalResonanceFilters(:,1);

OMEgainScalars=10.^(OMEgaindBs/20);

OMEfilterOrder=OMEexternalResonanceFilters(:,2);

OMElowerCutOff=OMEexternalResonanceFilters(:,3);

OMEupperCutOff=OMEexternalResonanceFilters(:,4);

% external resonance coefficients

ExtFilter_b=cell(nOMEExtFilters,1);

ExtFilter_a=cell(nOMEExtFilters,1);

for idx=1:nOMEExtFilters

    Nyquist=sampleRate/2;

    [b, a] = butter(OMEfilterOrder(idx), ...

        [OMElowerCutOff(idx) OMEupperCutOff(idx)]...

        /Nyquist);

    ExtFilter_b{idx}=b;

    ExtFilter_a{idx}=a;

end

OMEExtFilterBndry=cell(2,1);

OMEextEarPressure=zeros(1,signalLength); % pressure at tympanic membrane

% pressure to velocity conversion using smoothing filter (50 Hz cutoff)

tau=1/(2*pi*50);

a1=dt/tau-1; a0=1;

b0=1+ a1;

TMdisp_b=b0; TMdisp_a=[a0 a1];

% figure(9), freqz(TMdisp_b, TMdisp_a)

OME_TMdisplacementBndry=[];

% OME high pass (simulates poor low frequency stapes response)

OMEhighPassHighCutOff=OMEParams.OMEstapesLPcutoff;

Nyquist=sampleRate/2;

[stapesDisp_b,stapesDisp_a] = butter(1, OMEhighPassHighCutOff/Nyquist, 'high');

% figure(10), freqz(stapesDisp_b, stapesDisp_a)

OMEhighPassBndry=[];

% OMEampStapes might be reducdant (use OMEParams.stapesScalar)

stapesScalar= OMEParams.stapesScalar;

% Acoustic reflex

efferentDelayPts=round(OMEParams.ARdelay/dt);

% smoothing filter

a1=dt/OMEParams.ARtau-1; a0=1;

b0=1+ a1;

ARfilt_b=b0; ARfilt_a=[a0 a1];

ARattenuation=ones(1,signalLength);

ARrateThreshold=OMEParams.ARrateThreshold; % may not be used

ARrateToAttenuationFactor=OMEParams.rateToAttenuationFactor;

ARrateToAttenuationFactorProb=OMEParams.rateToAttenuationFactorProb;

ARboundary=[];

ARboundaryProb=0;

% save complete OME record (stapes displacement)

OMEoutput=zeros(1,signalLength);

TMoutput=zeros(1,signalLength);

%% BM ---

% BM is represented as a list of locations identified by BF

DRNL_BFs=BFlist;

nBFs= length(DRNL_BFs);

% DRNLchannelParameters=DRNLParams.channelParameters;

DRNLresponse= zeros(nBFs, segmentLength);

MOCrateToAttenuationFactor=DRNLParams.rateToAttenuationFactor;

rateToAttenuationFactorProb=DRNLParams.rateToAttenuationFactorProb;

MOCrateThresholdProb=DRNLParams.MOCrateThresholdProb;

% smoothing filter for MOC

a1=dt/DRNLParams.MOCtau-1; a0=1;

b0=1+ a1;

MOCfilt_b=b0; MOCfilt_a=[a0 a1];

% figure(9), freqz(stapesDisp_b, stapesDisp_a)

MOCboundary=cell(nBFs,1);

MOCprobBoundary=cell(nBFs,1);

MOCattSegment=zeros(nBFs,reducedSegmentLength);

MOCattenuation=ones(nBFs,signalLength);

if DRNLParams.a>0

    DRNLcompressionThreshold=10^((1/(1-DRNLParams.c))* ...

    log10(DRNLParams.b/DRNLParams.a));

else

    DRNLcompressionThreshold=inf;

end

DRNLlinearOrder= DRNLParams.linOrder;

DRNLnonlinearOrder= DRNLParams.nonlinOrder;

DRNLa=DRNLParams.a;

DRNLb=DRNLParams.b;

DRNLc=DRNLParams.c;

linGAIN=DRNLParams.g;

%

% gammatone filter coefficients for linear pathway

bw=DRNLParams.linBWs';

phi = 2 * pi * bw * dt;

cf=DRNLParams.linCFs';

theta = 2 * pi * cf * dt;

cos_theta = cos(theta);

sin_theta = sin(theta);

alpha = -exp(-phi).* cos_theta;

b0 = ones(nBFs,1);

b1 = 2 * alpha;

b2 = exp(-2 * phi);

z1 = (1 + alpha .* cos_theta) - (alpha .* sin_theta) * i;

z2 = (1 + b1 .* cos_theta) - (b1 .* sin_theta) * i;

z3 = (b2 .* cos(2 * theta)) - (b2 .* sin(2 * theta)) * i;

tf = (z2 + z3) ./ z1;

a0 = abs(tf);

a1 = alpha .* a0;

GTlin_a = [b0, b1, b2];

GTlin_b = [a0, a1];

GTlinOrder=DRNLlinearOrder;

GTlinBdry=cell(nBFs,GTlinOrder);

% nonlinear gammatone filter coefficients 

bw=DRNLParams.nlBWs';

phi = 2 * pi * bw * dt;

cf=DRNLParams.nonlinCFs';

theta = 2 * pi * cf * dt;

cos_theta = cos(theta);

sin_theta = sin(theta);

alpha = -exp(-phi).* cos_theta;

b0 = ones(nBFs,1);

b1 = 2 * alpha;

b2 = exp(-2 * phi);

z1 = (1 + alpha .* cos_theta) - (alpha .* sin_theta) * i;

z2 = (1 + b1 .* cos_theta) - (b1 .* sin_theta) * i;

z3 = (b2 .* cos(2 * theta)) - (b2 .* sin(2 * theta)) * i;

tf = (z2 + z3) ./ z1;

a0 = abs(tf);

a1 = alpha .* a0;

GTnonlin_a = [b0, b1, b2];

GTnonlin_b = [a0, a1];

GTnonlinOrder=DRNLnonlinearOrder;

GTnonlinBdry1=cell(nBFs, GTnonlinOrder);

GTnonlinBdry2=cell(nBFs, GTnonlinOrder);

% complete BM record (BM displacement)

DRNLoutput=zeros(nBFs, signalLength);

%% IHC ---

% IHC cilia activity and receptor potential

% viscous coupling between BM and stereocilia displacement

% Nyquist=sampleRate/2;

% IHCcutoff=1/(2*pi*IHC_cilia_RPParams.tc);

% [IHCciliaFilter_b,IHCciliaFilter_a]=...

%     butter(1, IHCcutoff/Nyquist, 'high');

a1=dt/IHC_cilia_RPParams.tc-1; a0=1;

b0=1+ a1;

% high pass (i.e. low pass reversed)

IHCciliaFilter_b=[a0 a1]; IHCciliaFilter_a=b0;

% figure(9), freqz(IHCciliaFilter_b, IHCciliaFilter_a)

IHCciliaBndry=cell(nBFs,1);

% IHC apical conductance (Boltzman function)

IHC_C= IHC_cilia_RPParams.C;

IHCu0= IHC_cilia_RPParams.u0;

IHCu1= IHC_cilia_RPParams.u1;

IHCs0= IHC_cilia_RPParams.s0;

IHCs1= IHC_cilia_RPParams.s1;

IHCGmax= IHC_cilia_RPParams.Gmax;

IHCGa= IHC_cilia_RPParams.Ga; % (leakage)

IHCGu0 = IHCGa+IHCGmax./(1+exp(IHCu0/IHCs0).*(1+exp(IHCu1/IHCs1)));

IHCrestingCiliaCond=IHCGu0;

% Receptor potential

IHC_Cab= IHC_cilia_RPParams.Cab;

IHC_Gk= IHC_cilia_RPParams.Gk;

IHC_Et= IHC_cilia_RPParams.Et;

IHC_Ek= IHC_cilia_RPParams.Ek;

IHC_Ekp= IHC_Ek+IHC_Et*IHC_cilia_RPParams.Rpc;

IHCrestingV= (IHC_Gk*IHC_Ekp+IHCGu0*IHC_Et)/(IHCGu0+IHC_Gk);

IHC_Vnow= IHCrestingV*ones(nBFs,1); % initial voltage

IHC_RP= zeros(nBFs,segmentLength);

% complete record of IHC receptor potential (V)

IHCciliaDisplacement= zeros(nBFs,segmentLength);

IHCoutput= zeros(nBFs,signalLength);

IHC_cilia_output= zeros(nBFs,signalLength);

%% pre-synapse ---

% Each BF is replicated using a different fiber type to make a 'channel'

% The number of channels is nBFs x nANfiberTypes

% Fiber types are specified in terms of tauCa

nANfiberTypes= length(IHCpreSynapseParams.tauCa);

tauCas= IHCpreSynapseParams.tauCa;

nANchannels= nANfiberTypes*nBFs;

synapticCa= zeros(nANchannels,segmentLength);

% Calcium control (more calcium, greater release rate)

ECa=IHCpreSynapseParams.ECa;

gamma=IHCpreSynapseParams.gamma;

beta=IHCpreSynapseParams.beta;

tauM=IHCpreSynapseParams.tauM;

mICa=zeros(nANchannels,segmentLength);

GmaxCa=IHCpreSynapseParams.GmaxCa;

synapse_z= IHCpreSynapseParams.z;

synapse_power=IHCpreSynapseParams.power;

% tauCa vector is established across channels to allow vectorization

%  (one tauCa per channel). Do not confuse with tauCas (one pre fiber type)

tauCa=repmat(tauCas, nBFs,1);

tauCa=reshape(tauCa, nANchannels, 1);

% presynapse startup values (vectors, length:nANchannels)

% proportion (0 - 1) of Ca channels open at IHCrestingV

mICaCurrent=((1+beta^-1 * exp(-gamma*IHCrestingV))^-1)...

    *ones(nBFs*nANfiberTypes,1);

% corresponding startup currents

ICaCurrent= (GmaxCa*mICaCurrent.^3) * (IHCrestingV-ECa);

CaCurrent= ICaCurrent.*tauCa;

% vesicle release rate at startup (one per channel)

% kt0 is used only at initialisation

kt0= -synapse_z * CaCurrent.^synapse_power;

%% AN ---

% each row of the AN matrices represents one AN fiber

% The results computed either for probabiities *or* for spikes (not both)

% Spikes are necessary if CN and IC are to be computed

nFibersPerChannel= AN_IHCsynapseParams.numFibers;

nANfibers= nANchannels*nFibersPerChannel;

AN_refractory_period= AN_IHCsynapseParams.refractory_period;

y=AN_IHCsynapseParams.y;

l=AN_IHCsynapseParams.l;

x=AN_IHCsynapseParams.x;

r=AN_IHCsynapseParams.r;

M=round(AN_IHCsynapseParams.M);

% probability            (NB initial 'P' on everything)

PAN_ydt = repmat(AN_IHCsynapseParams.y*dt, nANchannels,1);

PAN_ldt = repmat(AN_IHCsynapseParams.l*dt, nANchannels,1);

PAN_xdt = repmat(AN_IHCsynapseParams.x*dt, nANchannels,1);

PAN_rdt = repmat(AN_IHCsynapseParams.r*dt, nANchannels,1);

PAN_rdt_plus_ldt = PAN_rdt + PAN_ldt;

PAN_M=round(AN_IHCsynapseParams.M);

% compute starting values

Pcleft    = kt0* y* M ./ (y*(l+r)+ kt0* l);

Pavailable    = Pcleft*(l+r)./kt0;

Preprocess    = Pcleft*r/x; % canbe fractional

ANprobability=zeros(nANchannels,segmentLength);

ANprobRateOutput=zeros(nANchannels,signalLength);

lengthAbsRefractoryP= round(AN_refractory_period/dt);

% special variables for monitoring synaptic cleft (specialists only)

savePavailableSeg=zeros(nANchannels,segmentLength);

savePavailable=zeros(nANchannels,signalLength);

% spikes     % !  !  !    ! !        !   !  !

lengthAbsRefractory= round(AN_refractory_period/ANdt);

AN_ydt= repmat(AN_IHCsynapseParams.y*ANdt, nANfibers,1);

AN_ldt= repmat(AN_IHCsynapseParams.l*ANdt, nANfibers,1);

AN_xdt= repmat(AN_IHCsynapseParams.x*ANdt, nANfibers,1);

AN_rdt= repmat(AN_IHCsynapseParams.r*ANdt, nANfibers,1);

AN_rdt_plus_ldt= AN_rdt + AN_ldt;

AN_M= round(AN_IHCsynapseParams.M);

% kt0  is initial release rate

% Establish as a vector (length=channel x number of fibers)

kt0= repmat(kt0', nFibersPerChannel, 1);

kt0=reshape(kt0, nANfibers,1);

% starting values for reservoirs

AN_cleft    = kt0* y* M ./ (y*(l+r)+ kt0* l);

AN_available    = round(AN_cleft*(l+r)./kt0); %must be integer

AN_reprocess    = AN_cleft*r/x;

% output is in a logical array spikes = 1/ 0.

ANspikes= false(nANfibers,reducedSegmentLength);

ANoutput= false(nANfibers,reducedSignalLength);

%% CN (first brain stem nucleus - could be any subdivision of CN)

% Input to a CN neuorn is a random selection of AN fibers within a channel

%  The number of AN fibers used is ANfibersFanInToCN

% CNtauGk (Potassium time constant) determines the rate of firing of

%  the unit when driven hard by a DC input (not normally >350 sp/s)

% If there is more than one value, everything is replicated accordingly

ANavailableFibersPerChan=AN_IHCsynapseParams.numFibers;

ANfibersFanInToCN=MacGregorMultiParams.fibersPerNeuron;

CNtauGk=MacGregorMultiParams.tauGk; % row vector of CN types (by tauGk)

nCNtauGk=length(CNtauGk);

% the total number of 'channels' is now greater

nCNchannels=nANchannels*nCNtauGk;

nCNneuronsPerChannel=MacGregorMultiParams.nNeuronsPerBF;

tauGk=repmat(CNtauGk, nCNneuronsPerChannel,1);

tauGk=reshape(tauGk,nCNneuronsPerChannel*nCNtauGk,1);

% Now the number of neurons has been increased

nCNneurons=nCNneuronsPerChannel*nCNchannels;

CNmembranePotential=zeros(nCNneurons,reducedSegmentLength);

% establish which ANfibers (by name) feed into which CN nuerons

CNinputfiberLists=zeros(nANchannels*nCNneuronsPerChannel, ANfibersFanInToCN);

unitNo=1;

for ch=1:nANchannels

    % Each channel contains a number of units =length(listOfFanInValues)

    for idx=1:nCNneuronsPerChannel

        for idx2=1:nCNtauGk

            fibersUsed=(ch-1)*ANavailableFibersPerChan + ...

                ceil(rand(1,ANfibersFanInToCN)* ANavailableFibersPerChan);

            CNinputfiberLists(unitNo,:)=fibersUsed;

            unitNo=unitNo+1;

        end

    end

end

% input to CN units

AN_PSTH=zeros(nCNneurons,reducedSegmentLength);

% Generate CNalphaFunction function

%  by which spikes are converted to post-synaptic currents

CNdendriteLPfreq= MacGregorMultiParams.dendriteLPfreq;

CNcurrentPerSpike=MacGregorMultiParams.currentPerSpike;

CNspikeToCurrentTau=1/(2*pi*CNdendriteLPfreq);

t=ANdt:ANdt:5*CNspikeToCurrentTau;

CNalphaFunction= (1 / ...

    CNspikeToCurrentTau)*t.*exp(-t /CNspikeToCurrentTau);

CNalphaFunction=CNalphaFunction*CNcurrentPerSpike;

% figure(98), plot(t,CNalphaFunction)

% working memory for implementing convolution

CNcurrentTemp=...

    zeros(nCNneurons,reducedSegmentLength+length(CNalphaFunction)-1);

% trailing alphas are parts of humps carried forward to the next segment

CNtrailingAlphas=zeros(nCNneurons,length(CNalphaFunction));

CN_tauM=MacGregorMultiParams.tauM;

CN_tauTh=MacGregorMultiParams.tauTh;

CN_cap=MacGregorMultiParams.Cap;

CN_c=MacGregorMultiParams.c;

CN_b=MacGregorMultiParams.dGkSpike;

CN_Ek=MacGregorMultiParams.Ek;

CN_Eb= MacGregorMultiParams.Eb;

CN_Er=MacGregorMultiParams.Er;

CN_Th0= MacGregorMultiParams.Th0;

CN_E= zeros(nCNneurons,1);

CN_Gk= zeros(nCNneurons,1);

CN_Th= MacGregorMultiParams.Th0*ones(nCNneurons,1);

CN_Eb=CN_Eb.*ones(nCNneurons,1);

CN_Er=CN_Er.*ones(nCNneurons,1);

CNtimeSinceLastSpike=zeros(nCNneurons,1);

% tauGk is the main distinction between neurons

%  in fact they are all the same in the standard model

tauGk=repmat(tauGk,nANchannels,1);

CNoutput=false(nCNneurons,reducedSignalLength);

%% MacGregor (IC - second nucleus) --------

nICcells=nANchannels*nCNtauGk;  % one cell per channel

CN_PSTH=zeros(nICcells ,reducedSegmentLength);

ICspikeWidth=0.00015;   % this may need revisiting

epochsPerSpike=round(ICspikeWidth/ANdt);

if epochsPerSpike<1

    error(['MacGregorMulti: sample rate too low to support ' ...

        num2str(ICspikeWidth*1e6) '  microsec spikes']);

end

% short names

IC_tauM=MacGregorParams.tauM;

IC_tauGk=MacGregorParams.tauGk;

IC_tauTh=MacGregorParams.tauTh;

IC_cap=MacGregorParams.Cap;

IC_c=MacGregorParams.c;

IC_b=MacGregorParams.dGkSpike;

IC_Th0=MacGregorParams.Th0;

IC_Ek=MacGregorParams.Ek;

IC_Eb= MacGregorParams.Eb;

IC_Er=MacGregorParams.Er;

IC_E=zeros(nICcells,1);

IC_Gk=zeros(nICcells,1);

IC_Th=IC_Th0*ones(nICcells,1);

% Dendritic filtering, all spikes are replaced by CNalphaFunction functions

ICdendriteLPfreq= MacGregorParams.dendriteLPfreq;

ICcurrentPerSpike=MacGregorParams.currentPerSpike;

ICspikeToCurrentTau=1/(2*pi*ICdendriteLPfreq);

t=ANdt:ANdt:3*ICspikeToCurrentTau;

IC_CNalphaFunction= (ICcurrentPerSpike / ...

    ICspikeToCurrentTau)*t.*exp(-t / ICspikeToCurrentTau);

% figure(98), plot(t,IC_CNalphaFunction)

% working space for implementing alpha function

ICcurrentTemp=...

    zeros(nICcells,reducedSegmentLength+length(IC_CNalphaFunction)-1);

ICtrailingAlphas=zeros(nICcells, length(IC_CNalphaFunction));

ICfiberTypeRates=zeros(nANfiberTypes,reducedSignalLength);

ICoutput=false(nICcells,reducedSignalLength);

ICmembranePotential=zeros(nICcells,reducedSegmentLength);

ICmembraneOutput=zeros(nICcells,signalLength);

%% Main program %%  %%  %%  %%  %%  %%  %%  %%  %%  %%  %%  %%  %%  %%

%  Compute the entire model for each segment

segmentStartPTR=1;

reducedSegmentPTR=1; % when sampling rate is reduced

while segmentStartPTR<signalLength

    segmentEndPTR=segmentStartPTR+segmentLength-1;

    % shorter segments after speed up.

    shorterSegmentEndPTR=reducedSegmentPTR+reducedSegmentLength-1;

    inputPressureSegment=inputSignal...

        (:,segmentStartPTR:segmentStartPTR+segmentLength-1);

    % segment debugging plots

    % figure(98)

    % plot(segmentTime,inputPressureSegment), title('signalSegment')

    % OME ----------------------

    % OME Stage 1: external resonances. Add to inputSignal pressure wave

    y=inputPressureSegment;

    for n=1:nOMEExtFilters

        % any number of resonances can be used

        [x  OMEExtFilterBndry{n}] = ...

            filter(ExtFilter_b{n},ExtFilter_a{n},...

            inputPressureSegment, OMEExtFilterBndry{n});

        x= x* OMEgainScalars(n);

        % This is a parallel resonance so add it

        y=y+x;

    end

    inputPressureSegment=y;

    OMEextEarPressure(segmentStartPTR:segmentEndPTR)= inputPressureSegment;

    % OME stage 2: convert input pressure (velocity) to

    %  tympanic membrane(TM) displacement using low pass filter

    [TMdisplacementSegment  OME_TMdisplacementBndry] = ...

        filter(TMdisp_b,TMdisp_a,inputPressureSegment, ...

        OME_TMdisplacementBndry);

    % and save it

    TMoutput(segmentStartPTR:segmentEndPTR)= TMdisplacementSegment;

    % OME stage 3: middle ear high pass effect to simulate stapes inertia

    [stapesDisplacement  OMEhighPassBndry] = ...

        filter(stapesDisp_b,stapesDisp_a,TMdisplacementSegment, ...

        OMEhighPassBndry);

    % OME stage 4:  apply stapes scalar

    stapesDisplacement=stapesDisplacement*stapesScalar;

    % OME stage 5:    acoustic reflex stapes attenuation

    %  Attenuate the TM response using feedback from LSR fiber activity

    if segmentStartPTR>efferentDelayPts

        stapesDisplacement= stapesDisplacement.*...

            ARattenuation(segmentStartPTR-efferentDelayPts:...

            segmentEndPTR-efferentDelayPts);

    end

    % segment debugging plots

    % figure(98)

    % plot(segmentTime, stapesDisplacement), title ('stapesDisplacement')

    % and save

    OMEoutput(segmentStartPTR:segmentEndPTR)= stapesDisplacement;

    %% BM ------------------------------

    % Each location is computed separately

    for BFno=1:nBFs

        %            *linear* path

        linOutput = stapesDisplacement * linGAIN;  % linear gain

        for order = 1 : GTlinOrder

            [linOutput GTlinBdry{BFno,order}] = ...

                filter(GTlin_b(BFno,:), GTlin_a(BFno,:), linOutput, GTlinBdry{BFno,order});