MAP

function vectorStrength=testAN(targetFrequency,BFlist, levels, ...

    paramsName,paramChanges)

% testIHC used either for IHC I/O function ...

%  or receptive field (doReceptiveFields=1)

global IHC_VResp_VivoParams  IHC_cilia_RPParams IHCpreSynapseParams

global AN_IHCsynapseParams

    global ANoutput ANdt CNoutput ICoutput ICmembraneOutput tauCas

    global ARattenuation MOCattenuation

dbstop if error

restorePath=path;

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

    ['..' filesep 'parameterStore'],  ['..' filesep 'wavFileStore'],...

    ['..' filesep 'testPrograms'])

if nargin<5

    paramChanges=[];

end

if nargin<4

    paramsName='Normal';

end

if nargin<3

    levels=-10:10:80;

    % levels=80:10:90;

end

nLevels=length(levels);

toneDuration=.2;

rampDuration=0.002;

silenceDuration=.02;

localPSTHbinwidth=0.001;

% Use only the first frequency in the GUI targetFrequency box to defineBF

% targetFrequency=stimulusParameters.targetFrequency(1);

% BFlist=targetFrequency;

AN_HSRonset=zeros(nLevels,1);

AN_HSRsaturated=zeros(nLevels,1);

AN_LSRonset=zeros(nLevels,1);

AN_LSRsaturated=zeros(nLevels,1);

CNLSRrate=zeros(nLevels,1);

CNHSRsaturated=zeros(nLevels,1);

ICHSRsaturated=zeros(nLevels,1);

ICLSRsaturated=zeros(nLevels,1);

vectorStrength=zeros(nLevels,1);

AR=zeros(nLevels,1);

MOC=zeros(nLevels,1);

% ANoutput=zeros(200,200);

figure(15), clf

set(gcf,'position',[980   356   401   321])

figure(5), clf

set(gcf,'position', [980 34 400 295])

drawnow

%% guarantee that the sample rate is at least 10 times the frequency

sampleRate=50000;

while sampleRate< 10* targetFrequency

    sampleRate=sampleRate+10000;

end

%% adjust sample rate so that the pure tone stimulus has an integer

%% numver of epochs in a period

dt=1/sampleRate;

period=1/targetFrequency;

ANspeedUpFactor=5;  %anticipating MAP (needs clearing up)

ANdt=dt*ANspeedUpFactor;

ANdt=period/round(period/ANdt);

dt=ANdt/ANspeedUpFactor;

sampleRate=1/dt;

epochsPerPeriod=sampleRate*period;

%% main computational loop (vary level)    

levelNo=0;

for leveldB=levels

    levelNo=levelNo+1;

    fprintf('%4.0f\t', leveldB)

    amp=28e-6*10^(leveldB/20);

    time=dt:dt:toneDuration;

    rampTime=dt:dt:rampDuration;

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

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

    ramp=ramp.*fliplr(ramp);

    silence=zeros(1,round(silenceDuration/dt));

    % create signal (leveldB/ targetFrequency)

    inputSignal=amp*sin(2*pi*targetFrequency'*time);

    inputSignal= ramp.*inputSignal;

    inputSignal=[silence inputSignal];

    %% run the model

    AN_spikesOrProbability='spikes';

    showPlotsAndDetails=0;

    MAP1_14(inputSignal, 1/dt, BFlist, ...

        paramsName, AN_spikesOrProbability, paramChanges);

    nTaus=length(tauCas);

    %LSR (same as HSR if no LSR fibers present)

    [nANFibers nTimePoints]=size(ANoutput);

    numLSRfibers=nANFibers/nTaus;

    numHSRfibers=numLSRfibers;

    LSRspikes=ANoutput(1:numLSRfibers,:);

    PSTH=UTIL_PSTHmaker(LSRspikes, ANdt, localPSTHbinwidth);

    PSTHLSR=mean(PSTH,1)/localPSTHbinwidth;  % across fibers rates

    PSTHtime=localPSTHbinwidth:localPSTHbinwidth:...

        localPSTHbinwidth*length(PSTH);

    AN_LSRonset(levelNo)= max(PSTHLSR); % peak in 5 ms window

    AN_LSRsaturated(levelNo)= mean(PSTHLSR(round(length(PSTH)/2):end));

    % HSR

    HSRspikes= ANoutput(end- numHSRfibers+1:end, :);

    PSTH=UTIL_PSTHmaker(HSRspikes, ANdt, localPSTHbinwidth);

    PSTH=mean(PSTH,1)/localPSTHbinwidth; % sum across fibers (HSR only)

    AN_HSRonset(levelNo)= max(PSTH);

    AN_HSRsaturated(levelNo)= mean(PSTH(round(length(PSTH)/2): end));

    figure(5), subplot(2,2,2)

    hold off, bar(PSTHtime,PSTH, 'b')

    hold on,  bar(PSTHtime,PSTHLSR,'r')

    ylim([0 1000])

    xlim([0 length(PSTH)*localPSTHbinwidth])

    set(gcf,'name',[num2str(BFlist), ' Hz: ' num2str(leveldB) ' dB']);

    % AN - CV

    %  CV is computed 5 times. Use the middle one (3) as most typical

    cvANHSR=  UTIL_CV(HSRspikes, ANdt);

    % AN - vector strength

    PSTH=sum(HSRspikes);

    [PH, binTimes]=UTIL_periodHistogram...

        (PSTH, ANdt, targetFrequency);

    VS=UTIL_vectorStrength(PH);

    vectorStrength(levelNo)=VS;

    disp(['sat rate= ' num2str(AN_HSRsaturated(levelNo)) ...

        ';   phase-locking VS = ' num2str(VS)])

    title(['AN HSR: CV=' num2str(cvANHSR(3),'%5.2f') ...

        'VS=' num2str(VS,'%5.2f')])

    % CN - first-order neurons

    % CN LSR

    [nCNneurons c]=size(CNoutput);

    nLSRneurons=round(nCNneurons/nTaus);

    CNLSRspikes=CNoutput(1:nLSRneurons,:);

    PSTH=UTIL_PSTHmaker(CNLSRspikes, ANdt, localPSTHbinwidth);

    PSTH=sum(PSTH)/nLSRneurons;

    CNLSRrate(levelNo)=mean(PSTH(round(length(PSTH)/2):end))/localPSTHbinwidth;

    %CN HSR

    MacGregorMultiHSRspikes=...

        CNoutput(end-nLSRneurons:end,:);

    PSTH=UTIL_PSTHmaker(MacGregorMultiHSRspikes, ANdt, localPSTHbinwidth);

    PSTH=sum(PSTH)/nLSRneurons;

    PSTH=mean(PSTH,1)/localPSTHbinwidth; % sum across fibers (HSR only)

    CNHSRsaturated(levelNo)=mean(PSTH(length(PSTH)/2:end));

    figure(5), subplot(2,2,3)

    bar(PSTHtime,PSTH)

    ylim([0 1000])

    xlim([0 length(PSTH)*localPSTHbinwidth])

    cvMMHSR= UTIL_CV(MacGregorMultiHSRspikes, ANdt);

    title(['CN    CV= ' num2str(cvMMHSR(3),'%5.2f')])

    % IC LSR

    [nICneurons c]=size(ICoutput);

    nLSRneurons=round(nICneurons/nTaus);

    ICLSRspikes=ICoutput(1:nLSRneurons,:);

    PSTH=UTIL_PSTHmaker(ICLSRspikes, ANdt, localPSTHbinwidth);

    ICLSRsaturated(levelNo)=mean(PSTH(round(length(PSTH)/2):end))/localPSTHbinwidth;

    %IC HSR

    MacGregorMultiHSRspikes=...

        ICoutput(end-nLSRneurons:end,:);

    PSTH=UTIL_PSTHmaker(MacGregorMultiHSRspikes, ANdt, localPSTHbinwidth);

    PSTH=sum(PSTH)/nLSRneurons;

    PSTH=mean(PSTH,1)/localPSTHbinwidth; % sum across fibers (HSR only)

    ICHSRsaturated(levelNo)=mean(PSTH(length(PSTH)/2:end));

    AR(levelNo)=min(ARattenuation);

    MOC(levelNo)=min(MOCattenuation(length(MOCattenuation)/2:end));

    time=dt:dt:dt*size(ICmembraneOutput,2);

    figure(5), subplot(2,2,4)

    plot(time,ICmembraneOutput(2, 1:end),'k')

    ylim([-0.07 0])

    xlim([0 max(time)])

    title(['IC  ' num2str(leveldB,'%4.0f') 'dB'])

    drawnow

    figure(5), subplot(2,2,1)

    plot(20*log10(MOC), 'k'),

    title(' MOC'), ylabel('dB attenuation')

    ylim([-30 0])

end % level

figure(5), subplot(2,2,1)

plot(levels,20*log10(MOC), 'k'),

title(' MOC'), ylabel('dB attenuation')

ylim([-30 0])

xlim([0 max(levels)])

fprintf('\n')

toneDuration=2;

rampDuration=0.004;

silenceDuration=.02;

nRepeats=200;   % no. of AN fibers

fprintf('toneDuration  %6.3f\n', toneDuration)

fprintf(' %6.0f  AN fibers (repeats)\n', nRepeats)

fprintf('levels')

fprintf('%6.2f\t', levels)

fprintf('\n')

% ---------------------------------------------------- display parameters

figure(15), clf

nRows=2; nCols=2;

% AN rate - level ONSET functions

subplot(nRows,nCols,1)

plot(levels,AN_LSRonset,'ro'), hold on

plot(levels,AN_HSRonset,'ko'), hold off

ylim([0 1000]),  xlim([min(levels) max(levels)])

ttl=['tauCa= ' num2str(IHCpreSynapseParams.tauCa)];

title( ttl)

xlabel('level dB SPL'), ylabel('peak rate (sp/s)'), grid on

text(0, 800, 'AN onset', 'fontsize', 16)

% AN rate - level ADAPTED function

subplot(nRows,nCols,2)

plot(levels,AN_LSRsaturated, 'ro'), hold on

plot(levels,AN_HSRsaturated, 'ko'), hold off

ylim([0 400])

set(gca,'ytick',0:50:300)

xlim([min(levels) max(levels)])

set(gca,'xtick',[levels(1):20:levels(end)])

%     grid on

ttl=[   'spont=' num2str(mean(AN_HSRsaturated(1,:)),'%4.0f')...

    '  sat=' num2str(mean(AN_HSRsaturated(end,1)),'%4.0f')];

title( ttl)

xlabel('level dB SPL'), ylabel ('adapted rate (sp/s)')

text(0, 340, 'AN adapted', 'fontsize', 16), grid on

% CN rate - level ADAPTED function

subplot(nRows,nCols,3)

plot(levels,CNLSRrate, 'ro'), hold on

plot(levels,CNHSRsaturated, 'ko'), hold off

ylim([0 400])

set(gca,'ytick',0:50:300)

xlim([min(levels) max(levels)])

set(gca,'xtick',[levels(1):20:levels(end)])

%     grid on

ttl=[   'spont=' num2str(mean(CNHSRsaturated(1,:)),'%4.0f') '  sat=' ...

    num2str(mean(CNHSRsaturated(end,1)),'%4.0f')];

title( ttl)

xlabel('level dB SPL'), ylabel ('adapted rate (sp/s)')

text(0, 350, 'CN', 'fontsize', 16), grid on

% IC rate - level ADAPTED function

subplot(nRows,nCols,4)

plot(levels,ICLSRsaturated, 'ro'), hold on

plot(levels,ICHSRsaturated, 'ko'), hold off

ylim([0 400])

set(gca,'ytick',0:50:300)

xlim([min(levels) max(levels)])

set(gca,'xtick',[levels(1):20:levels(end)]), grid on

ttl=['spont=' num2str(mean(ICHSRsaturated(1,:)),'%4.0f') ...

    '  sat=' num2str(mean(ICHSRsaturated(end,1)),'%4.0f')];

title( ttl)

xlabel('level dB SPL'), ylabel ('adapted rate (sp/s)')

text(0, 350, 'IC', 'fontsize', 16)

set(gcf,'name',' AN CN IC rate/level')

peakVectorStrength=max(vectorStrength);

fprintf('\n')

disp('levels vectorStrength')

fprintf('%3.0f \t %6.4f \n', [levels; vectorStrength'])

fprintf('\n')

fprintf('Phase locking, max vector strength=\t %6.4f\n\n',...

    max(vectorStrength))

allData=[ levels'  AN_HSRonset AN_HSRsaturated...

    AN_LSRonset AN_LSRsaturated ...

    CNHSRsaturated CNLSRrate...

    ICHSRsaturated ICLSRsaturated];

fprintf('\n levels \tANHSR Onset \tANHSR adapted\tANLSR Onset \tANLSR adapted\tCNHSR\tCNLSR\tICHSR  \tICLSR \n');

UTIL_printTabTable(round(allData))

fprintf('VS (phase locking)= \t%6.4f\n\n',...

    max(vectorStrength))

UTIL_showStruct(IHC_cilia_RPParams, 'IHC_cilia_RPParams')

UTIL_showStruct(IHCpreSynapseParams, 'IHCpreSynapseParams')

UTIL_showStruct(AN_IHCsynapseParams, 'AN_IHCsynapseParams')

fprintf('\n')

disp('levels vectorStrength')

fprintf('%3.0f \t %6.4f \n', [levels; vectorStrength'])

fprintf('\n')

fprintf('Phase locking, max vector strength= \t%6.4f\n\n',...

    max(vectorStrength))

allData=[ levels'  AN_HSRonset AN_HSRsaturated...

    AN_LSRonset AN_LSRsaturated ...

    CNHSRsaturated CNLSRrate...

    ICHSRsaturated ICLSRsaturated];

fprintf('\n levels \tANHSR Onset \tANHSR adapted\tANLSR Onset \tANLSR adapted\tCNHSR\tCNLSR\tICHSR  \tICLSR \n');

UTIL_printTabTable(round(allData))

fprintf('VS (phase locking)= \t%6.4f\n\n',...

    max(vectorStrength))

path(restorePath)

function vectorStrength=testANprob(targetFrequency,BFlist, levels, ...

    paramsName, paramChanges)

% testIHC used either for IHC I/O function ...

%  or receptive field (doReceptiveFields=1)

global IHC_VResp_VivoParams  IHC_cilia_RPParams IHCpreSynapseParams

global AN_IHCsynapseParams

global ANprobRateOutput dt tauCas

global ARattenuation MOCattenuation

AN_spikesOrProbability='probability';

dbstop if error

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

    ['..' filesep 'parameterStore'],  ['..' filesep 'wavFileStore'],...

    ['..' filesep 'testPrograms'])

if nargin<5

    paramChanges=[];

end

if nargin<4

    paramsName='Normal';

end

if nargin<3

    levels=-10:10:80;

end

nLevels=length(levels);

toneDuration=.2;

rampDuration=0.002;

silenceDuration=.02;

localPSTHbinwidth=0.001;

% Use only the first frequency in the GUI targetFrequency box to defineBF

% targetFrequency=stimulusParameters.targetFrequency(1);

% BFlist=targetFrequency;

AN_HSRonset=zeros(nLevels,1);

AN_HSRsaturated=zeros(nLevels,1);

AN_LSRonset=zeros(nLevels,1);

AN_LSRsaturated=zeros(nLevels,1);

AR=zeros(nLevels,1);

MOC=zeros(nLevels,1);

figure(15), clf

set(gcf,'position',[607    17   368   321])

drawnow

%% guarantee that the sample rate is at least 10 times the frequency

sampleRate=50000;

while sampleRate< 10* targetFrequency

    sampleRate=sampleRate+10000;

end

%% adjust sample rate so that the pure tone stimulus has an integer

%% numver of epochs in a period

dt=1/sampleRate;

period=1/targetFrequency;

%% main computational loop (vary level)

levelNo=0;

for leveldB=levels

    levelNo=levelNo+1;

    fprintf('%4.0f\t', leveldB)

    amp=28e-6*10^(leveldB/20);

    time=dt:dt:toneDuration;

    rampTime=dt:dt:rampDuration;

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

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

    ramp=ramp.*fliplr(ramp);

    silence=zeros(1,round(silenceDuration/dt));

    % create signal (leveldB/ targetFrequency)

    inputSignal=amp*sin(2*pi*targetFrequency'*time);

    inputSignal= ramp.*inputSignal;

    inputSignal=[silence inputSignal];

    %% run the model

    showPlotsAndDetails=0;

    MAP1_14(inputSignal, 1/dt, BFlist, ...

        paramsName, AN_spikesOrProbability, paramChanges);

    nTaus=length(tauCas);

    %LSR (same as HSR if no LSR fibers present)

    [nANFibers nTimePoints]=size(ANprobRateOutput);

    numLSRfibers=1;

    numHSRfibers=numLSRfibers;

    LSRspikes=ANprobRateOutput(1:numLSRfibers,:);

    PSTH=UTIL_PSTHmaker(LSRspikes, dt, localPSTHbinwidth);

    PSTHLSR=PSTH/(localPSTHbinwidth/dt);  % across fibers rates

    PSTHtime=localPSTHbinwidth:localPSTHbinwidth:...

        localPSTHbinwidth*length(PSTH);

    AN_LSRonset(levelNo)= max(PSTHLSR); % peak in 5 ms window

    AN_LSRsaturated(levelNo)= mean(PSTHLSR(round(length(PSTH)/2):end));

    % HSR

    HSRspikes= ANprobRateOutput(end- numHSRfibers+1:end, :);

    PSTH=UTIL_PSTHmaker(HSRspikes, dt, localPSTHbinwidth);

    PSTH=PSTH/(localPSTHbinwidth/dt); % sum across fibers (HSR only)

    AN_HSRonset(levelNo)= max(PSTH);

    AN_HSRsaturated(levelNo)= mean(PSTH(round(length(PSTH)/2): end));

    figure(15), subplot(2,2,4)

    hold off, bar(PSTHtime,PSTH, 'b')

    hold on,  bar(PSTHtime,PSTHLSR,'r')

    ylim([0 1000])

    xlim([0 length(PSTH)*localPSTHbinwidth])

    set(gcf,'name',[num2str(BFlist), ' Hz: ' num2str(leveldB) ' dB']);

    AR(levelNo)=min(ARattenuation);

    MOC(levelNo)=min(MOCattenuation(length(MOCattenuation)/2:end));

    figure(15), subplot(2,2,3)

    plot(20*log10(MOC), 'k'),

    title(' MOC'), ylabel('dB attenuation')

    ylim([-30 0])

end % level

figure(15), subplot(2,2,3)

plot(levels,20*log10(MOC), 'k'),

title(' MOC'), ylabel('dB attenuation')

ylim([-30 0])

xlim([0 max(levels)])

fprintf('\n')

toneDuration=2;

rampDuration=0.004;

silenceDuration=.02;

nRepeats=200;   % no. of AN fibers

fprintf('toneDuration  %6.3f\n', toneDuration)

fprintf(' %6.0f  AN fibers (repeats)\n', nRepeats)

fprintf('levels')

fprintf('%6.2f\t', levels)

fprintf('\n')

% ---------------------------------------------------- display parameters

nRows=2; nCols=2;

% AN rate - level ONSET functions

subplot(nRows,nCols,1)

plot(levels,AN_LSRonset,'ro'), hold on

plot(levels,AN_HSRonset,'ko'), hold off

ylim([0 1000]),  xlim([min(levels) max(levels)])

ttl=['tauCa= ' num2str(IHCpreSynapseParams.tauCa)];

title( ttl)

xlabel('level dB SPL'), ylabel('peak rate (sp/s)'), grid on

text(0, 800, 'AN onset', 'fontsize', 16)

% AN rate - level ADAPTED function

subplot(nRows,nCols,2)

plot(levels,AN_LSRsaturated, 'ro'), hold on

plot(levels,AN_HSRsaturated, 'ko'), hold off

ylim([0 400])

set(gca,'ytick',0:50:300)

xlim([min(levels) max(levels)])

set(gca,'xtick',[levels(1):20:levels(end)])

%     grid on

ttl=[   'spont=' num2str(mean(AN_HSRsaturated(1,:)),'%4.0f')...

    '  sat=' num2str(mean(AN_HSRsaturated(end,1)),'%4.0f')];

title( ttl)

xlabel('level dB SPL'), ylabel ('adapted rate (sp/s)')

text(0, 340, 'AN adapted', 'fontsize', 16), grid on

allData=[ levels'  AN_HSRonset AN_HSRsaturated...

    AN_LSRonset AN_LSRsaturated ];

fprintf('\n levels \tANHSR Onset \tANHSR adapted\tANLSR Onset \tANLSR adapted\tCNHSR\tCNLSR\tICHSR  \tICLSR \n');

UTIL_printTabTable(round(allData))

UTIL_showStruct(IHC_cilia_RPParams, 'IHC_cilia_RPParams')

UTIL_showStruct(IHCpreSynapseParams, 'IHCpreSynapseParams')

UTIL_showStruct(AN_IHCsynapseParams, 'AN_IHCsynapseParams')

fprintf('\n')

disp('levels vectorStrength')

allData=[ levels'  AN_HSRonset AN_HSRsaturated...

    AN_LSRonset AN_LSRsaturated ];

fprintf('\n levels \tANHSR Onset \tANHSR adapted\tANLSR Onset \tANLSR adapted\tCNHSR\tCNLSR\tICHSR  \tICLSR \n');

UTIL_printTabTable(round(allData))

function testBM (BMlocations, paramsName,...

    relativeFrequencies, AN_spikesOrProbability, paramChanges)

% testBM generates input output functions for DRNL model for any number

% of locations.

% Computations are bast on a single channel model (channelBFs=BF)

% peak displacement (peakAmp) is measured.

%  if DRNLParams.useMOC is chosen, the full model is run (slow)

%  otherwise only DRNL is computed.

% Tuning curves are generated based on a range of frequencies reletove to

% the BF of the location.

%

global    DRNLParams

if nargin<5

    paramChanges=[];

end

if nargin<4

    AN_spikesOrProbability='probability';

end

savePath=path;

addpath (['..' filesep 'utilities'],['..' filesep 'MAP'])

levels=-10:10:90;   nLevels=length(levels);

% levels= 50;   nLevels=length(levels);

% refBMdisplacement is the displacement of the BM at threshold

% 1 nm disp at  threshold (9 kHz, Ruggero)

% ? adjust for frequency

refBMdisplacement= 1e-8; % adjusted for 10 nm at 1 kHz 

toneDuration=.200;

rampDuration=0.01;

silenceDuration=0.01;

sampleRate=30000;

dbstop if error

figure(3), clf

set(gcf,'position',[280   350   327   326])

set(gcf,'name','DRNL - BM')

pause(0.1)

finalSummary=[];

nBFs=length(BMlocations);

BFno=0; plotCount=0;

for BF=BMlocations

    BFno=BFno+1;

    plotCount=plotCount+nBFs;

    stimulusFrequencies=BF* relativeFrequencies;

    nFrequencies=length(stimulusFrequencies);

    peakAmpBM=zeros(nLevels,nFrequencies);

    peakAmpBMdB=NaN(nLevels,nFrequencies);

    peakEfferent=NaN(nLevels,nFrequencies);

    peakAREfferent=NaN(nLevels,nFrequencies);

    levelNo=0;

    for leveldB=levels

        disp(['level= ' num2str(leveldB)])

        levelNo=levelNo+1;

        freqNo=0;

        for frequency=stimulusFrequencies

            freqNo=freqNo+1;

            % Generate stimuli

            globalStimParams.FS=sampleRate;

            globalStimParams.overallDuration=...

                toneDuration+silenceDuration;  % s

            stim.phases='sin';

            stim.type='tone';

            stim.toneDuration=toneDuration;

            stim.frequencies=frequency;

            stim.amplitudesdB=leveldB;

            stim.beginSilence=silenceDuration;

            stim.rampOnDur=rampDuration;

            % no offset ramp

            stim.rampOffDur=rampDuration;

            doPlot=0;

            inputSignal=stimulusCreate(globalStimParams, stim, doPlot);

            inputSignal=inputSignal(:,1)';

            %% run the model

            MAPparamsName=paramsName;

            global DRNLoutput MOCattenuation ARattenuation

            MAP1_14(inputSignal, sampleRate, BF, ...

                MAPparamsName, AN_spikesOrProbability, paramChanges);

            DRNLresponse=DRNLoutput;

            peakAmp=max(max(...

                DRNLresponse(:, end-round(length(DRNLresponse)/2):end)));

            peakAmpBM(levelNo,freqNo)=peakAmp;

            peakAmpBMdB(levelNo,freqNo)=...

                20*log10(peakAmp/refBMdisplacement);

            peakEfferent(levelNo,freqNo)=min(min(MOCattenuation));

            peakAREfferent(levelNo,freqNo)=min(min(ARattenuation));

        end  % tone frequency

    end  % level

    %% analyses results and plot

    % BM I/O plot (top panel)

    figure(3)

    subplot(3,nBFs,BFno), cla

    plot(levels,peakAmpBMdB, 'linewidth',2)

    hold on, plot(levels, repmat(refBMdisplacement,1,length(levels)))

    hold off

    title(['BF=' num2str(BF,'%5.0f') ' - ' paramsName])

    xlabel('level')

    %     set(gca,'xtick',levels),  grid on

    if length(levels)>1,xlim([min(levels) max(levels)]), end

    ylabel(['dB re:' num2str(refBMdisplacement,'%6.1e') 'm'])

    ylim([-20 50])

    set(gca,'ytick',[-10 0 10 20 40])

    grid on

    %     legend({num2str(stimulusFrequencies')}, 'location', 'EastOutside')

    UTIL_printTabTable([levels' peakAmpBMdB], ...

        num2str([0 stimulusFrequencies]','%5.0f'), '%5.0f')

    finalSummary=[finalSummary peakAmpBMdB];

    % Tuning curve

    if length(relativeFrequencies)>2

        figure(3), subplot(3,nBFs, 2*nBFs+BFno)

        %         contour(stimulusFrequencies,levels,peakAmpBM,...

        %             [refBMdisplacement refBMdisplacement],'r')

        contour(stimulusFrequencies,levels,peakAmpBM,...

            refBMdisplacement.*[1 5 10 50 100])

        title(['tuning curve at ' num2str(refBMdisplacement) 'm']);

        ylabel('level (dB) at reference')

        xlim([100 10000])

        grid on

        hold on

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

    end

    % MOC contribution

    figure(3)

    subplot(3,nBFs,nBFs+BFno), cla

    plot(levels,20*log10(peakEfferent), 'linewidth',2)

    ylabel('MOC (dB attenuation)'), xlabel('level')

    title(['peak MOC: model= ' AN_spikesOrProbability])

    grid on

    if length(levels)>1, xlim([min(levels) max(levels)]), end

    % AR contribution

    hold on

    plot(levels,20*log10(peakAREfferent), 'r')

    hold off

end % best frequency

UTIL_showStructureSummary(DRNLParams, 'DRNLParams', 10)

    UTIL_printTabTable([levels' finalSummary], ...

        num2str([0 stimulusFrequencies]','%5.0f'), '%5.0f')

path(savePath);

function testFM(BFlist,paramsName,showPSTHs,paramChanges)

%   specify whether you want AN 'probability' or 'spikes'

%       spikes is more realistic but takes longer

%       refractory effect is included only for spikes

%

% specify the AN ANresponse bin width (normally 1 ms).

%      This can influence the measure of the AN onset rate based on the

%      largest bin in the ANresponse

%

% Demonstration is based on Harris and Dallos

global inputStimulusParams outerMiddleEarParams DRNLParams 

global IHC_VResp_VivoParams IHCpreSynapseParams  AN_IHCsynapseParams

dbstop if error

restorePath=path;

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

    ['..' filesep 'parameterStore'],  ['..' filesep 'wavFileStore'],...

    ['..' filesep 'testPrograms'])

if nargin<3

    paramChanges=[];

end

% masker and probe levels are relative to this threshold

thresholdAtCF=10; % dB SPL

thresholdAtCF=5; % dB SPL

if nargin<1, showPSTHs=1;end

sampleRate=50000;

% fetch BF from GUI: use only the first target frequency

maskerFrequency=BFlist;

maskerDuration=.1;

targetFrequency=maskerFrequency;

probeLeveldB=20+thresholdAtCF;	% H&D use 20 dB SL/ TMC uses 10 dB SL

probeDuration=0.008; % HD use 15 ms probe (fig 3).

probeDuration=0.015; % HD use 15 ms probe (fig 3).

rampDuration=.001;  % HD use 1 ms linear ramp

initialSilenceDuration=0.02;

finalSilenceDuration=0.05;      % useful for showing recovery

maskerLevels=[-80   10 20 30 40 60 ] + thresholdAtCF;

% maskerLevels=[-80   40 60 ] + thresholdAtCF;

dt=1/sampleRate;

figure(7), clf

set(gcf,'position',[613    36   360   247])

set(gcf,'name', ['forward masking: thresholdAtCF=' num2str(thresholdAtCF)])

if showPSTHs

    figure(8), clf

    set(gcf,'name', 'Harris and Dallos simulation')

    set(gcf,'position',[980    36   380   249])

end

% Plot Harris and Dallos result for comparison

gapDurations=[0.001	0.002	0.005	0.01	0.02	0.05	0.1	0.3];

HDmaskerLevels=[+10	+20	+30	+40	+60];

HDresponse=[

    1 1 1 1 1 1 1 1;

    0.8  	0.82	0.81	0.83	0.87	0.95	1	    1;

    0.48	0.5	    0.54	0.58	0.7	    0.85	0.95	1;

    0.3	    0.31	0.35	0.4	    0.5	    0.68	0.82	0.94;

    0.2 	0.27	0.27	0.29	0.42	0.64	0.75	0.92;

    0.15	0.17	0.18	0.23	0.3	     0.5	0.6	    0.82];

figure(7)

semilogx(gapDurations,HDresponse,'o'), hold on

legend(strvcat(num2str(maskerLevels')),'location','southeast')

title([ 'masker dB: ' num2str(HDmaskerLevels)])

%% Run the trials

maxProbeResponse=0;

levelNo=0;

resultsMatrix=zeros(length(maskerLevels), length(gapDurations));

summaryFiringRates=[];

% initial silence

time=dt: dt: initialSilenceDuration;

initialSilence=zeros(1,length(time));

% probe

time=dt: dt: probeDuration;

amp=28e-6*10^(probeLeveldB/20);

probe=amp*sin(2*pi.*targetFrequency*time);

% ramps

% catch rampTime error

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

rampTime=dt:dt:rampDuration;

% raised cosine ramp

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

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

%  onset ramp

probe=probe.*ramp;

%  offset ramp makes complete ramp for probe

ramp=fliplr(ramp);

% apply ramp mask to probe. Probe does not change below

probe=probe.*ramp;

% final silence

time=dt: dt: finalSilenceDuration;

finalSilence=zeros(1,length(time));

PSTHplotCount=0;

longestSignalDuration=initialSilenceDuration + maskerDuration +...

    max(gapDurations) + probeDuration + finalSilenceDuration ;

for maskerLeveldB=maskerLevels

    levelNo=levelNo+1;

    allDurations=[];

    allFiringRates=[];

    %masker

    time=dt: dt: maskerDuration;

    masker=sin(2*pi.*maskerFrequency*time);

    % masker ramps

    if rampDuration>0.5*maskerDuration

        % catch ramp duration error

        rampDuration=maskerDuration/2;

    end

    % NB masker ramp (not probe ramp)

    rampTime=dt:dt:rampDuration;

    % raised cosine ramp

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

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

    %  onset ramp

    masker=masker.*ramp;

    %  offset ramp

    ramp=fliplr(ramp);

    % apply ramp

    masker=masker.*ramp;

    amp=28e-6*10^(maskerLeveldB/20);

    maskerPa=amp*masker;

    for gapDuration=gapDurations

        time=dt: dt: gapDuration;

        gap=zeros(1,length(time));

        inputSignal=...

            [initialSilence maskerPa gap probe finalSilence];

        % **********************************  run MAP model

        global  ANprobRateOutput  tauCas  ...

    showPlotsAndDetails=0;

AN_spikesOrProbability='probability';

MAP1_14(inputSignal, 1/dt, targetFrequency, ...

    paramsName, AN_spikesOrProbability, paramChanges);

    [nFibers c]=size(ANprobRateOutput);

    nLSRfibers=nFibers/length(tauCas);

            ANresponse=ANprobRateOutput(end-nLSRfibers:end,:);

            ANresponse=sum(ANresponse)/nLSRfibers;

        % analyse results

        probeStart=initialSilenceDuration+maskerDuration+gapDuration;

        PSTHbinWidth=dt;

        responseDelay=0.005;

        probeTimes=probeStart+responseDelay:...

            PSTHbinWidth:probeStart+probeDuration+responseDelay;

        probeIDX=round(probeTimes/PSTHbinWidth);

        probePSTH=ANresponse(probeIDX);

        firingRate=mean(probePSTH);

        % NB this only works if you start with the lowest level masker

        maxProbeResponse=max([maxProbeResponse firingRate]);

        allDurations=[allDurations gapDuration];

        allFiringRates=[allFiringRates firingRate];

        %% show PSTHs

        if showPSTHs

            nLevels=length(maskerLevels);

            nDurations=length(gapDurations);

            figure(8)

            PSTHbinWidth=1e-3;

            PSTH=UTIL_PSTHmaker(ANresponse, dt, PSTHbinWidth);

            PSTH=PSTH*dt/PSTHbinWidth;

            PSTHplotCount=PSTHplotCount+1;

            subplot(nLevels,nDurations,PSTHplotCount)

            probeTime=PSTHbinWidth:PSTHbinWidth:...

                PSTHbinWidth*length(PSTH);

            bar(probeTime, PSTH)

            if strcmp(AN_spikesOrProbability, 'spikes')

                ylim([0 500])

            else

                ylim([0 500])

            end

            xlim([0 longestSignalDuration])

            grid on

            if PSTHplotCount< (nLevels-1) * nDurations+1

                set(gca,'xticklabel',[])

            end

            if ~isequal(mod(PSTHplotCount,nDurations),1)

                set(gca,'yticklabel',[])

            else

                ylabel([num2str(maskerLevels...

                    (round(PSTHplotCount/nDurations) +1))])

            end

            if PSTHplotCount<=nDurations

                title([num2str(1000*gapDurations(PSTHplotCount)) 'ms'])

            end

        end % showPSTHs

    end     % gapDurations duration

    summaryFiringRates=[summaryFiringRates allFiringRates'];

    figure(7), hold on

    semilogx(allDurations, allFiringRates/maxProbeResponse)

    ylim([0 1]), hold on

    ylim([0 inf]), xlim([min(gapDurations) max(gapDurations)])

    xlim([0.001 1])

    pause(0.1) % to allow for CTRL/C interrupts

    resultsMatrix(levelNo,:)=allFiringRates;

end          % maskerLevel

disp('delay/ rates')

disp(num2str(round( [1000*allDurations' summaryFiringRates])))

% replot with best adjustment

% figure(34), clf% use for separate plot

figure(7), clf

peakProbe=max(max(resultsMatrix));

resultsMatrix=resultsMatrix/peakProbe;

semilogx(gapDurations,HDresponse,'o'), hold on

title(['FrMsk: probe ' num2str(probeLeveldB)...

    'dB SL: peakProbe=' num2str(peakProbe,'%5.0f') ' sp/s'])

xlabel('gap duration (s)'), ylabel ('probe response')

semilogx(allDurations, resultsMatrix'), ylim([0 1])

ylim([0 inf]),

xlim([0.001 1])

legend(strvcat(num2str(maskerLevels'-thresholdAtCF)), -1)

% ------------------------------------------------- display parameters

disp(['parameter file was: ' paramsName])

fprintf('\n')

% UTIL_showStruct(inputStimulusParams, 'inputStimulusParams')

% UTIL_showStruct(outerMiddleEarParams, 'outerMiddleEarParams')

% UTIL_showStruct(DRNLParams, 'DRNLParams')

% UTIL_showStruct(IHC_VResp_VivoParams, 'IHC_VResp_VivoParams')

UTIL_showStruct(IHCpreSynapseParams, 'IHCpreSynapseParams')

UTIL_showStruct(AN_IHCsynapseParams, 'AN_IHCsynapseParams')

path(restorePath);

function testOME(paramsName, paramChanges)

savePath=path;

addpath (['..' filesep 'utilities'],['..' filesep 'MAP'])

if nargin<2

    paramChanges=[];

end

sampleRate=50000;

dt=1/sampleRate;

leveldBSPL=80; 		% dB SPL as used by Huber (may trigger AR)

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

duration=.05;

time=dt: dt: duration;

%% Comparison data (human)

% These data are taken directly from Huber 2001 (Fig. 4)

HuberFrequencies=[600	  800	 1000	 2000	 3000	4000 6000 8000];

HuberDisplacementAt80dBSPL=[1.5E-9	1.5E-09	1.5E-09	1.0E-09	7.0E-10	...

    3.0E-10	2.0E-10	1.0E-10]; % m;

% HuberVelocityAt80dBSPL= 2*pi*HuberFrequencies.*HuberDisplacementAt80dBSPL;

figure(2), clf, subplot(2,1,1)

set(2,'position',[5   349   268   327])

semilogx(HuberFrequencies, 20*log10(HuberDisplacementAt80dBSPL/1e-10),...

    'ko', 'MarkerFaceColor','k', 'Marker','o', 'markerSize',6)

hold on

%% Generate test stimulus .................................................................

% independent test using discrete frequencies

peakResponses=[];

peakTMpressure=[];

frequencies=[200 400 HuberFrequencies 10000];

for toneFrequency=frequencies

    inputSignal=amp*sin(2*pi*toneFrequency*time);

    showPlotsAndDetails=0;

    AN_spikesOrProbability='probability';

    % switch off AR & MOC (Huber's patients were deaf)

    paramChanges{1}='OMEParams.rateToAttenuationFactorProb=0;';

    paramChanges{2}='DRNLParams.rateToAttenuationFactorProb = 0;';

    global OMEoutput  OMEextEarPressure TMoutput ARattenuation

    % BF is irrelevant

    MAP1_14(inputSignal, sampleRate, -1, ...

        paramsName, AN_spikesOrProbability, paramChanges);

    peakDisplacement=max(OMEoutput(floor(end/2):end));

    peakResponses=[peakResponses peakDisplacement];

    peakTMpressure=[peakTMpressure max(OMEextEarPressure)];

    disp([' AR attenuation (dB):   ' num2str(20*log10(min(ARattenuation)))])

end

%% Report

disp('frequency displacement(m)')

% disp(num2str([frequencies' peakResponses']))

fprintf('%6.0f \t%10.3e\n',[frequencies' peakResponses']')

% stapes peak displacement

figure(2), subplot(2,1,1), hold on

semilogx(frequencies, 20*log10(peakResponses/1e-10), 'r', 'linewidth',4)

set(gca,'xScale','log')

% ylim([1e-11 1e-8])

xlim([100 10000]), ylim([0 30])

grid on

title(['stapes at ' num2str(leveldBSPL) ' (NB deaf)'])

ylabel('disp: dB re 1e-10m')

xlabel('stimulus frequency (Hz)')

legend({'Huber et al','model'},'location','southWest')

set(gcf,'name','OME')

% external ear resonance

figure(2), subplot(2,1,2),hold off

semilogx(frequencies, 20*log10(peakTMpressure/28e-6)-leveldBSPL,...

    'k', 'linewidth',2)

xlim([100 10000]) %, ylim([-10 30])

grid on

title(['External ear resonances' ])

ylabel('dB')

xlabel('frequency')

set(gcf,'name','OME: external resonances')

% ---------------------------------------------------------- display parameters

disp(['parameter file was: ' paramsName])

fprintf('\n')

path(savePath);

function testPeriphery (BMlocations, paramsName)

% testBM generates input output functions for DRNL model for any number

% of locations.

% Computations are bast on a single channel model (channelBFs=BF)

% peak displacement (peakAmp) is measured.

%  if DRNLParams.useMOC is chosen, the full model is run (slow)

%  otherwise only DRNL is computed.

% Tuning curves are generated based on a range of frequencies reletove to

% the BF of the location.

%

global    DRNLParams

savePath=path;

addpath (['..' filesep 'utilities'],['..' filesep 'MAP'])

levels=-10:10:90;   nLevels=length(levels);

% levels= 50;   nLevels=length(levels);

% relativeFrequencies=[0.25    .5   .75  1  1.25 1.5    2];

relativeFrequencies=1;

% refBMdisplacement is the displacement of the BM at threshold

% 1 nm disp at  threshold (9 kHz, Ruggero)

refBMdisplacement= 1e-8;            % adjusted for 10 nm at 1 kHz

toneDuration=.200;

rampDuration=0.01;

silenceDuration=0.01;

sampleRate=30000;

dbstop if error

figure(3), clf

% set(gcf,'position',[276    33   331   645])

set(gcf,'name','DRNL - BM')

finalSummary=[];

nBFs=length(BMlocations);

BFno=0; plotCount=0;

for BF=BMlocations

    BFno=BFno+1;

    plotCount=plotCount+nBFs;

    stimulusFrequencies=BF* relativeFrequencies;

    nFrequencies=length(stimulusFrequencies);

    peakAmpBM=zeros(nLevels,nFrequencies);

    peakAmpBMdB=NaN(nLevels,nFrequencies);

    peakEfferent=NaN(nLevels,nFrequencies);

    peakAREfferent=NaN(nLevels,nFrequencies);

    levelNo=0;

    for leveldB=levels

        disp(['level= ' num2str(leveldB)])

        levelNo=levelNo+1;

        freqNo=0;

        for frequency=stimulusFrequencies

            freqNo=freqNo+1;

            % Generate stimuli

            globalStimParams.FS=sampleRate;

            globalStimParams.overallDuration=...

                toneDuration+silenceDuration;  % s

            stim.type='tone';

            stim.phases='sin';

            stim.toneDuration=toneDuration;

            stim.frequencies=frequency;

            stim.amplitudesdB=leveldB;

            stim.beginSilence=silenceDuration;

            stim.rampOnDur=rampDuration;

            stim.rampOffDur=rampDuration;

            doPlot=0;

            inputSignal=stimulusCreate(globalStimParams, stim, doPlot);

            inputSignal=inputSignal(:,1)';

            %% run the model

            MAPparamsName=paramsName;

            AN_spikesOrProbability='probability';

            % spikes are slow but can be used to study MOC using IC units

            % AN_spikesOrProbability='spikes';

            global DRNLoutput MOCattenuation ARattenuation IHCoutput

            MAP1_14(inputSignal, sampleRate, BF, ...

                MAPparamsName, AN_spikesOrProbability);

            DRNLresponse=IHCoutput;

            peakAmp=max(max(...

                DRNLresponse(:, end-round(length(DRNLresponse)/2):end)));

            peakAmpBM(levelNo,freqNo)=peakAmp;

            if peakAmp>0

            peakAmpBMdB(levelNo,freqNo)=...

                20*log10(peakAmp/refBMdisplacement);

            else

                peakAmpBMdB(levelNo,freqNo)=peakAmp;

            end

            peakEfferent(levelNo,freqNo)=min(min(MOCattenuation));

            peakAREfferent(levelNo,freqNo)=min(min(ARattenuation));

        end  % tone frequency

    end  % level

    %% analyses results and plot

    % BM I/O plot (top panel)

    figure(3)

    subplot(3,nBFs,BFno), cla

    plot(levels,peakAmpBMdB, 'linewidth',2)

    hold on, plot(levels, repmat(refBMdisplacement,1,length(levels)))

    hold off

    title(['BF=' num2str(BF,'%5.0f') ' - ' paramsName])

    xlabel('level')

    %     set(gca,'xtick',levels),  grid on

    if length(levels)>1,xlim([min(levels) max(levels)]), end

    ylabel(['dB re:' num2str(refBMdisplacement,'%6.1e') 'm'])

    ylim([-20 50])

    set(gca,'ytick',[-10 0 10 20 40])

    %     legend({num2str(stimulusFrequencies')}, 'location', 'EastOutside')

    UTIL_printTabTable([levels' peakAmpBMdB], ...

        num2str([0 stimulusFrequencies]','%5.0f'), '%5.0f')

    finalSummary=[finalSummary peakAmpBMdB];

    % Tuning curve

    if length(relativeFrequencies)>2

        figure(3), subplot(3,nBFs, nBFs+BFno)

        %         contour(stimulusFrequencies,levels,peakAmpBM,...

        %             [refBMdisplacement refBMdisplacement],'r')

        contour(stimulusFrequencies,levels,peakAmpBM,...

            refBMdisplacement.*[1 5 10 50 100])

        title(['tuning curve at ' num2str(refBMdisplacement) 'm']);

        ylabel('level (dB) at reference')

        xlim([100 10000])

        hold on

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

    end

    % MOC contribution

    figure(3)

    subplot(3,nBFs,2*nBFs+BFno), cla

    plot(levels,20*log10(peakEfferent), 'linewidth',2)

    ylabel('MOC (dB attenuation)'), xlabel('level')

    title(['peak MOC: model= ' AN_spikesOrProbability])

    grid on

    if length(levels)>1, xlim([min(levels) max(levels)]), end

    % AR contribution

    hold on

    plot(levels,20*log10(peakAREfferent), 'r')

    hold off

end % best frequency

UTIL_showStructureSummary(DRNLParams, 'DRNLParams', 10)

    UTIL_printTabTable([levels' finalSummary], ...

        num2str([0 stimulusFrequencies]','%5.0f'), '%5.0f')

path(savePath);

function testPhaseLocking(paramsName, paramChanges)

if nargin<2

    paramChanges=[];

end

testFrequencies=[250 500 1000 2000 4000 8000];

levels=50:10:80;

figure(14), clf

set(gcf,'position', [980    36   383   321])

set(gcf,'name', 'phase locking')

allStrengths=zeros(length(testFrequencies), length(levels));

peakVectorStrength=zeros(1,length(testFrequencies));

freqCount=0;

for targetFrequency=testFrequencies;

    %single test

    freqCount=freqCount+1;

    vectorStrength=...

        testAN(targetFrequency,targetFrequency, levels,...

        paramsName, paramChanges);

    allStrengths(freqCount,:)=vectorStrength';

    peakVectorStrength(freqCount)=max(vectorStrength');

end

%% plot results

figure(14)

subplot(2,1,2)

plot(levels,allStrengths)

xlabel('levels')

ylabel('vector strength')

legend (num2str(testFrequencies'),'location','eastOutside')

subplot(2,1,1)

semilogx(testFrequencies,peakVectorStrength)

grid on

title ('peak vector strength')

xlabel('frequency')

ylabel('vector strength')

johnson=[250	0.79

500	0.82

1000	0.8

2000	0.7

4000	0.25

5500	0.05

];

hold on

plot(johnson(:,1),johnson(:,2),'o')

legend({'model','Johnson 80'},'location','eastOutside')

hold off

function testPhysiology(BF,paramsName, paramChanges)

restorePath=path;

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

if nargin<3

    paramChanges=[];

end

testOME(paramsName,paramChanges)

relativeFrequencies=[0.25    .5   .75  1  1.25 1.5    2];

testBM (BF, paramsName,relativeFrequencies,'spikes', paramChanges)

testRP(BF,paramsName,paramChanges)

testSynapse(BF,paramsName,paramChanges)

testFM(BF,paramsName,1,paramChanges)

testPhaseLocking(paramsName,paramChanges)

testAN(BF,BF, -10:10:80,paramsName,paramChanges);

figure(14)

figure(15)

path(restorePath)

function testPhysiologyProb(BF,paramsName, paramChanges)

restorePath=path;

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

if nargin<3

    paramChanges=[];

end

testOME(paramsName, paramChanges)

relativeFrequencies=[0.25    .5   .75  1  1.25 1.5    2];

testBM (BF, paramsName,relativeFrequencies,'probability', paramChanges)

testRP(BF,paramsName, paramChanges)

testSynapse(BF,paramsName, paramChanges)

testANprob(BF,BF, -10:10:80,paramsName, paramChanges);

figure(4)

path(restorePath)

function testRF

% testIHC used either for IHC I/O function or receptive field (doReceptiveFields=1)

global experiment method stimulusParameters expGUIhandles

global inputStimulusParams  IHC_ciliaParams

global IHC_VResp_VivoParams IHCpreSynapseParams  AN_IHCsynapseParams

dbstop if error

% set(expGUIhandles.pushbuttonStop, 'backgroundColor', [.941 .941 .941])

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

    ['..' filesep 'parameterStore'],  ['..' filesep 'wavFileStore'],...

    ['..' filesep 'testPrograms'])

targetFrequency=stimulusParameters.targetFrequency(1);

sampleRate=50000;

doReceptiveFields=1;

toneDuration=.05;

rampDuration=0.004;

silenceDuration=.02;

nRepeats=100;   % no. of AN fibers

plotGraphsForIHC=1;

% number of MacGregor units is set in the parameter file.

if doReceptiveFields

    % show all receptive field

    frequencies=targetFrequency*    [  0.5         0.7         0.9     1       1.1         1.3         1.6];

    levels=0:20:80; nLevels=length(levels);

    figure(14), clf

    figure(15), clf

else

    % show only I/O function at BF

    frequencies=targetFrequency;

    levels=-20:10:90;

    %         levels=10:.25:13;

    %         levels=-20:1:-15

    nLevels=length(levels);

%     figure(13), clf,

%     set (gcf, 'name', ['IHC/AN  input/output' num2str(AN_IHCsynapseParams.numFibers) ' repeats'])

%     drawnow

end

nFrequencies=length(frequencies);

IHC_RP_peak=zeros(nLevels,nFrequencies);

IHC_RP_min=zeros(nLevels,nFrequencies);

IHC_RP_dc=zeros(nLevels,nFrequencies);

AN_HSRonset=zeros(nLevels,nFrequencies);

AN_HSRsaturated=zeros(nLevels,nFrequencies);

AN_LSRonset=zeros(nLevels,nFrequencies);

AN_LSRsaturated=zeros(nLevels,nFrequencies);

CNLSRsaturated=zeros(nLevels,nFrequencies);

CNHSRsaturated=zeros(nLevels,nFrequencies);

ICHSRsaturated=zeros(nLevels,nFrequencies);

ICLSRsaturated=zeros(nLevels,nFrequencies);

levelNo=0; PSTHplotCount=0;

for leveldB=levels

    fprintf('%4.0f\t', leveldB)

    levelNo=levelNo+1;

    amp=28e-6*10^(leveldB/20);

    freqNo=0;

    for frequency=frequencies

        paramFunctionName=['method=MAPparams' experiment.name ...

            '(' num2str(targetFrequency) ');' ];

        eval(paramFunctionName);  % read parameters afresh each pass

        dt=method.dt;

        time=dt:dt:toneDuration;        

        rampTime=dt:dt:rampDuration;

        ramp=[0.5*(1+cos(2*pi*rampTime/(2*rampDuration)+pi)) ones(1,length(time)-length(rampTime))];

        ramp=ramp.*fliplr(ramp);

        silence=zeros(1,round(silenceDuration/dt));

        toneStartptr=length(silence)+1;

        toneMidptr=toneStartptr+round(toneDuration/(2*dt)) -1;

        toneEndptr=toneStartptr+round(toneDuration/dt) -1;

        % create signal (leveldB/ frequency)

        freqNo=freqNo+1;

        inputSignal=amp*sin(2*pi*frequency'*time);

        inputSignal= ramp.*inputSignal;

        inputSignal=[silence inputSignal silence];

        if doReceptiveFields  % receptive field

            method.plotGraphs=	0; % plot only PSTHs

        else

            method.plotGraphs=	plotGraphsForIHC; % show progress

        end

        targetChannelNo=1;

        % force parameters

         % the number of AN fibers at each BF

        AN_IHCsynapseParams.numFibers=	nRepeats;

        AN_IHCsynapseParams. mode= 'spikes';

        AN_IHCsynapseParams.plotSynapseContents=0;

        AN_IHCsynapseParams.PSTHbinWidth=.001;

        method.DRNLSave=1;

        method.IHC_cilia_RPSave=1;

        method.PSTHbinWidth=1e-3; % useful 1-ms default for all PSTHs

        method.AN_IHCsynapseSave=1;

        method.MacGregorMultiSave=1;

        method.MacGregorSave=1;

        method.dt=dt;

        moduleSequence=[1:8];       

        global ANdt ARAttenuation TMoutput OMEoutput ...

    DRNLoutput IHC_cilia_output IHCrestingCiliaCond IHCrestingV...

    IHCoutput ANprobRateOutput ANoutput savePavailable tauCas  ...

    CNoutput  ICoutput ICmembraneOutput ICfiberTypeRates MOCattenuation

AN_spikesOrProbability='spikes';

AN_spikesOrProbability='probability';

MAPparamsName='Normal';

MAP1_14(inputSignal, 1/dt, targetFrequency, ...

    MAPparamsName, AN_spikesOrProbability);

        % RP

        IHC_RPData=IHC_cilia_output;

        IHC_RPData=IHCoutput(targetChannelNo,:);

        IHC_RP_peak(levelNo,freqNo)=max(IHC_RPData(toneStartptr:toneEndptr));

        IHC_RP_min(levelNo,freqNo)=min(IHC_RPData(toneStartptr:toneEndptr));

        IHC_RP_dc(levelNo,freqNo)=mean(IHC_RPData(toneStartptr:toneEndptr));

        % AN next

        AN_IHCsynapseAllData=ANoutput;

        method.PSTHbinWidth=0.001;

        nTaus=length(tauCas);

        numANfibers=size(ANoutput,1);

        numLSRfibers=numANfibers/nTaus;

        %LSR (same as HSR if no LSR fibers present)

        channelPtr1=(targetChannelNo-1)*numANfibers+1;

        channelPtr2=channelPtr1+numANfibers-1;

        LSRspikes=AN_IHCsynapseAllData(channelPtr1:channelPtr2,:);

        method.dt=method.AN_IHCsynapsedt;

        PSTH=UTIL_PSTHmaker(LSRspikes, method);

        PSTH=sum(PSTH,1); % sum across fibers (HSR only)

        PSTHStartptr=round(silenceDuration/method.PSTHbinWidth)+1;

        PSTHMidptr=PSTHStartptr+round(toneDuration/(2*method.PSTHbinWidth)) -1;

        PSTHEndptr=PSTHStartptr+round(toneDuration/method.PSTHbinWidth) -1;

        AN_LSRonset(levelNo,freqNo)=max(max(PSTH))/(method.PSTHbinWidth*method.numANfibers);

        AN_LSRsaturated(levelNo,freqNo)=sum(PSTH(PSTHMidptr:PSTHEndptr))/(method.numANfibers*toneDuration/2);

        % HSR

        channelPtr1=numLSRfibers+(targetChannelNo-1)*method.numANfibers+1;

        channelPtr2=channelPtr1+method.numANfibers-1;

        HSRspikes=AN_IHCsynapseAllData(channelPtr1:channelPtr2,:);

        method.dt=method.AN_IHCsynapsedt;

        PSTH=UTIL_PSTHmaker(HSRspikes, method);

        PSTH=sum(PSTH,1); % sum across fibers (HSR only)

        PSTHStartptr=round(silenceDuration/method.PSTHbinWidth)+1;

        PSTHMidptr=PSTHStartptr+round(toneDuration/(2*method.PSTHbinWidth)) -1;

        PSTHEndptr=PSTHStartptr+round(toneDuration/method.PSTHbinWidth) -1;