MAP

function [modelResponse, MacGregorResponse]=MAPmodel( MAPplot, method)

global experiment stimulusParameters audio withinRuns

global outerMiddleEarParams DRNLParams AN_IHCsynapseParams

savePath=path;

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

modelResponse=[];

MacGregorResponse=[];

% mono only (column vector)

audio=audio(:,1)';

% if stop button pressed earlier

if experiment.stop, return, end

% -------------------------------------------------------------- run Model

MAPparamsName=experiment.name;

showPlotsAndDetails=experiment.MAPplot;

AN_spikesOrProbability='spikes';

% [response, method]=MAPsequenceSeg(audio, method, 1:8);

global ICoutput ANdt

    MAP1_14(audio, 1/method.dt, method.nonlinCF,...

        MAPparamsName, AN_spikesOrProbability);

if showPlotsAndDetails

    options.printModelParameters=0;

    options.showModelOutput=1;

    options.printFiringRates=1;

    options.showACF=0;

    options.showEfferent=1;

    UTIL_showMAP(options)

end

% No response,  probably caused by hitting 'stop' button

if isempty(ICoutput), return, end

% MacGregor response is the sum total of all final stage spiking

MacGregorResponse= sum(ICoutput,1);                 % use IC

% ---------------------------------------------------------- end model run

dt=ANdt;

time=dt:dt:dt*length(MacGregorResponse);

% group delay on unit response

MacGonsetDelay= 0.004;

MacGoffsetDelay= 0.022;

% now find the response of the MacGregor model during the target presentation + group delay

switch experiment.threshEstMethod

    case {'2I2AFC++', '2I2AFC+++'}

        idx= time>stimulusParameters.testTargetBegins+MacGonsetDelay ...

            & time<stimulusParameters.testTargetEnds+MacGoffsetDelay;

        nSpikesTrueWindow=sum(MacGregorResponse(:,idx));

        idx=find(time>stimulusParameters.testNonTargetBegins+MacGonsetDelay ...

            & time<stimulusParameters.testNonTargetEnds+MacGoffsetDelay);

        nSpikesFalseWindow=sum(MacGregorResponse(:,idx));

        % nSpikesDuringTarget is +ve when more spikes are found

        %   in the target window

        difference= nSpikesTrueWindow-nSpikesFalseWindow;

        if difference>0

            % hit

            nSpikesDuringTarget=experiment.MacGThreshold+1;

        elseif    difference<0

            % miss (wrong choice)

            nSpikesDuringTarget=experiment.MacGThreshold-1;

        else

            if rand>0.5

                % hit (random choice)

                nSpikesDuringTarget=experiment.MacGThreshold+1;

            else

                % miss (random choice)

                nSpikesDuringTarget=experiment.MacGThreshold-1;

            end

        end

        disp(['level target dummy decision: ' ...

            num2str([withinRuns.variableValue nSpikesTrueWindow ...

            nSpikesFalseWindow  nSpikesDuringTarget], '%4.0f') ] )

    otherwise

        % idx=find(time>stimulusParameters.testTargetBegins+MacGonsetDelay ...

        %         & time<stimulusParameters.testTargetEnds+MacGoffsetDelay);

        % no delay at onset

        idx=find(time>stimulusParameters.testTargetBegins +MacGonsetDelay...

            & time<stimulusParameters.testTargetEnds+MacGoffsetDelay);

        nSpikesDuringTarget=sum(MacGregorResponse(:,idx));

        % find(MacGregorResponse)*dt-stimulusParameters.stimulusDelay

        timeX=time(idx);

end

% now find the response of the MacGregor model at the end of the masker

idx2=find(time>stimulusParameters.testTargetBegins-0.02 ...

    & time<stimulusParameters.testTargetBegins);

if ~isempty(idx2)

    maskerRate=mean(mean(MacGregorResponse(idx2)));

else

    %e.g. no masker

    maskerRate=0;

end

if experiment.MAPplot

    % add vertical lines to indicate target region

    figure(99), subplot(6,1,6)

    hold on

    yL=get(gca,'YLim');

    plot([stimulusParameters.testTargetBegins + MacGonsetDelay ...

        stimulusParameters.testTargetBegins   + MacGonsetDelay],yL,'r')

    plot([stimulusParameters.testTargetEnds   + MacGoffsetDelay ...

        stimulusParameters.testTargetEnds     + MacGoffsetDelay],yL,'r')

end

% specify unambiguous response

switch experiment.paradigm

    case 'gapDetection'

        gapResponse=(maskerRate-nSpikesDuringTarget)/maskerRate;

        if gapResponse>0.2

            modelResponse=2;    % gap detected

        else

            modelResponse=1;    % gap not detected

        end

        [nSpikesDuringTarget maskerRate gapResponse modelResponse]

        figure(22), plot(timeX,earObject(idx))

    otherwise

        if nSpikesDuringTarget>experiment.MacGThreshold

            modelResponse=2;    % stimulus detected

        else

            modelResponse=1;    % nothing heard (default)

        end

end

path(savePath)

function UTIL_testPhysiology(BF)

testOME('Normal')

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

testBM (BF, 'Normal',relativeFrequencies)

testRP(BF,'Normal')

testSynapse(BF,'Normal')

testFM(BF,'Normal',1)

testPhaseLocking

testAN(BF,BF, 10:10:80);

figure(14)

figure(15)

function paradigm_IFMC_8ms(handles)

global stimulusParameters experiment betweenRuns

paradigm_training(handles) % default

stimulusParameters.WRVname='maskerLevel';

stimulusParameters.WRVstartValues=50;

stimulusParameters.WRVsteps= [-10 -2];

stimulusParameters.WRVlimits=[-30 110];

experiment.psyFunSlope = -1;

withinRuns.direction='up';

betweenRuns.variableName1='maskerRelativeFrequency';

betweenRuns.variableList1=[1       0.5       1.6   .9 .7   1.1 1.3      ];

betweenRuns.variableName2='targetFrequency';

% keep old list of target frequencies

betweenRuns.variableList2=str2num(get(handles.edittargetFrequency,'string'));

betweenRuns.randomizeSequence=1; % 'random sequence'

stimulusParameters.maskerDuration=0.108;

stimulusParameters.maskerLevel=stimulusParameters.WRVstartValues(1);

stimulusParameters.maskerRelativeFrequency=betweenRuns.variableList1;

stimulusParameters.gapDuration=0.01;

stimulusParameters.targetType='tone';

stimulusParameters.targetPhase='sin';

stimulusParameters.targetFrequency=betweenRuns.variableList2(1);

stimulusParameters.targetDuration=0.016;

stimulusParameters.targetLevel=NaN;

stimulusParameters.rampDuration=0.004;

% instructions to user

% single interval up/down no cue

stimulusParameters.instructions{1}=[{'YES if you hear the added click'}, { }, { 'NO if not (or you are uncertain'}];

% single interval up/down with cue

stimulusParameters.instructions{2}=[{'count how many distinct clicks you hear'},{'ignore the tones'},{' '},...

    {'The clicks must be **clearly distinct** to count'}];

function paradigm_TMC_16ms(handles)

global stimulusParameters experiment betweenRuns

paradigm_training(handles) % default

experiment.singleIntervalMaxTrials=20;

stimulusParameters.WRVname='maskerLevel';

stimulusParameters.WRVstartValues=50;

stimulusParameters.WRVsteps= [-10 -4];

stimulusParameters.WRVlimits=[-30 110];

stimulusParameters.cueTestDifference = 10;

experiment.psyFunSlope = -1;

withinRuns.direction='up';

betweenRuns.variableName1='gapDuration';

betweenRuns.variableList1=[.01 .09 .03 .05 .07];

betweenRuns.variableName2='targetFrequency';

% retain existing targetFrequencies

betweenRuns.variableList2=str2num(get(handles.edittargetFrequency,'string'));

betweenRuns.randomizeSequence=1; % 'random sequence'

stimulusParameters.maskerType='tone';

stimulusParameters.maskerPhase='sin';

stimulusParameters.maskerDuration=0.108;

stimulusParameters.maskerLevel=stimulusParameters.WRVstartValues(1);

stimulusParameters.maskerRelativeFrequency=1;

stimulusParameters.gapDuration=betweenRuns.variableList1;

stimulusParameters.targetType='tone';

stimulusParameters.targetPhase='sin';

stimulusParameters.targetFrequency=betweenRuns.variableList2(1);

stimulusParameters.targetDuration=0.016;

stimulusParameters.targetLevel=NaN;

stimulusParameters.rampDuration=0.004;

% instructions to user

%   single interval up/down no cue

stimulusParameters.instructions{1}=[{'YES if you hear the added click'}, { }, { 'NO if not (or you are uncertain'}];

%   single interval up/down with cue

stimulusParameters.instructions{2}=[{'count how many distinct clicks you hear'},{'ignore the tones'},{' '},...

    {'The clicks must be **clearly distinct** to count'}];

function paradigm_TMCmodel(handles)

global stimulusParameters experiment betweenRuns

paradigm_training(handles) % default

stimulusParameters.WRVname='maskerLevel';

stimulusParameters.WRVstartValues=30;

stimulusParameters.WRVsteps= [-10 -4];

stimulusParameters.WRVlimits=[-30 110];

stimulusParameters.cueTestDifference = 10;

experiment.psyFunSlope = -1;

withinRuns.direction='up';

betweenRuns.variableName1='gapDuration';

betweenRuns.variableList1=[ 0.09 0.01:0.02:0.07 0.02:0.02:.08 0.005];

betweenRuns.variableName2='targetFrequency';

% retain existing targetFrequencies

betweenRuns.variableList2=str2num(get(handles.edittargetFrequency,'string'));

betweenRuns.randomizeSequence=1; % 'random sequence'

stimulusParameters.maskerType='tone';

stimulusParameters.maskerPhase='sin';

stimulusParameters.maskerDuration=0.108;

stimulusParameters.maskerLevel=stimulusParameters.WRVstartValues(1);

stimulusParameters.maskerRelativeFrequency=1;

experiment.singleIntervalMaxTrials=10;

stimulusParameters.gapDuration=betweenRuns.variableList1;

stimulusParameters.targetType='tone';

stimulusParameters.targetPhase='sin';

stimulusParameters.targetFrequency=betweenRuns.variableList2(1);

stimulusParameters.targetDuration=0.016;

stimulusParameters.targetLevel=NaN;

stimulusParameters.rampDuration=0.004;

% instructions to user

%   single interval up/down no cue

stimulusParameters.instructions{1}=[{'YES if you hear the added click'}, { }, { 'NO if not (or you are uncertain'}];

%   single interval up/down with cue

stimulusParameters.instructions{2}=[{'count how many distinct clicks you hear'},{'ignore the tones'},{' '},...

    {'The clicks must be **clearly distinct** to count'}];

function paradigm_absThreshold_8(handles)

global stimulusParameters experiment betweenRuns

paradigm_training(handles) % default

betweenRuns.variableName1='targetFrequency';

betweenRuns.variableList1=1000;

betweenRuns.variableList1=str2num(get(handles.edittargetFrequency,'string'));

betweenRuns.variableName2='targetDuration';

betweenRuns.variableList2=0.008 ;

betweenRuns.randomizeSequence=1; % 'random sequence'

% delay > masker > gap > target

stimulusParameters.targetType='tone';

stimulusParameters.targetPhase='sin';

stimulusParameters.targetFrequency=1000;

stimulusParameters.targetDuration=0.008;

stimulusParameters.targetLevel=stimulusParameters.WRVstartValues(1);

stimulusParameters.rampDuration=0.004;

% forced choice window interval

stimulusParameters.AFCsilenceDuration=0.5;

% instructions to user

%   single interval up/down no cue

stimulusParameters.instructions{1}= [{'YES if you hear the tone clearly'}, { }, { 'NO if not (or you are uncertain'}];

%   single interval up/down with cue

stimulusParameters.instructions{2}= [{'count the tones you hear clearly'}, { }, { 'ignore indistinct tones'}];

function paradigm_thr_IFMC(handles)

global stimulusParameters experiment betweenRuns

paradigm_training(handles) % default

betweenRuns.variableName1='targetFrequency';

betweenRuns.variableList1=1000;

betweenRuns.variableList1=str2num(get(handles.edittargetFrequency,'string'));

betweenRuns.variableName2='targetDuration';

betweenRuns.variableList2=0.016;

betweenRuns.randomizeSequence=1; % 'random sequence'

% delay > masker > gap > target

stimulusParameters.targetType='tone';

stimulusParameters.targetPhase='sin';

stimulusParameters.targetFrequency=1000;

stimulusParameters.targetDuration=0.016;

stimulusParameters.targetLevel=stimulusParameters.WRVstartValues(1);

stimulusParameters.rampDuration=0.004;

experiment.stopCriteria2IFC=[75 3 5];

experiment.singleIntervalMaxTrials=[20];

% forced choice window interval

stimulusParameters.AFCsilenceDuration=0.5;

% instructions to user

%   single interval up/down no cue

stimulusParameters.instructions{1}= [{'YES if you hear the tone clearly'}, { }, { 'NO if not (or you are uncertain'}];

%   single interval up/down with cue

stimulusParameters.instructions{2}= [{'count the tones you hear clearly'}, { }, { 'ignore indistinct tones'}];

function paradigm_thr_TMC(handles)

global stimulusParameters experiment betweenRuns

paradigm_training(handles) % default

betweenRuns.variableName1='targetFrequency';

betweenRuns.variableList1=1000;

betweenRuns.variableList1=str2num(get(handles.edittargetFrequency,'string'));

betweenRuns.variableName2='targetDuration';

betweenRuns.variableList2=0.016;

betweenRuns.randomizeSequence=1; % 'random sequence'

% delay > masker > gap > target

stimulusParameters.targetType='tone';

stimulusParameters.targetPhase='sin';

stimulusParameters.targetFrequency=1000;

stimulusParameters.targetDuration=0.016;

stimulusParameters.targetLevel=stimulusParameters.WRVstartValues(1);

stimulusParameters.rampDuration=0.004;

experiment.stopCriteria2IFC=[75 3 5];

experiment.singleIntervalMaxTrials=[20];

% forced choice window interval

stimulusParameters.AFCsilenceDuration=0.5;

% instructions to user

%   single interval up/down no cue

stimulusParameters.instructions{1}= [{'YES if you hear the tone clearly'}, { }, { 'NO if not (or you are uncertain'}];

%   single interval up/down with cue

stimulusParameters.instructions{2}= [{'count the tones you hear clearly'}, { }, { 'ignore indistinct tones'}];

function experiment=OHIOthresholds(experiment)

experiment.OHIOfrequencies=[494, 663, 870, 1125, 1442, 1838, 2338, 2957, 3725, 4680, 5866,  7334]; %Hz.

% User must specify abs thresholds (dB SPL) of each tone frequency

experiment.OHIOthresholds= [

    9.1	8.6	8.1	7.9	9.8	10.5	13.5	15.0	17.4	19.4	22.6	25.2

];

function paradigm_OHIOabs(handles)

global stimulusParameters experiment betweenRuns

paradigm_training(handles) % default

% find the  threshold for a tonecomplex  consisting of a sequence of 10-ms tones

%   whose frequencies are chosen at random from a list (OHIOfrequencies)

% All tones are presented at the same level (SL) computed using absolute

%   threshols specified in OHIOthresholds;

% The duration of the complex is increased across runs and the number of tones is

%   controlled by OHIOdurations, (for each 20 ms a further tone is added.

% The frequency of the tones is changed on each trial

experiment.OHIOfrequencies=[494, 663, 870, 1125, 1442, 1838, 2338, 2957, 3725, 4680, 5866,  7334]; %Hz.

% User must specify abs thresholds (dB SPL) of each tone frequency

% experiment.OHIOthresholds= [18	16	16	19	20	22	24	26	27	30	32	35];

%  assessment method

% {'oneIntervalUpDown', 'MaxLikelihood', '2I2AFC++', '2I2AFC+++'}

experiment.threshEstMethod='oneIntervalUpDown';

% {'cued', 'noCue'};

stimulusParameters.WRVname='targetLevel';

stimulusParameters.WRVstartValues=20 ;

stimulusParameters.WRVsteps=[10 2];

stimulusParameters.WRVlimits=[-30 110];

betweenRuns.variableName1='numOHIOtones';

betweenRuns.variableList1= 1:12; % i.e. the frequency to be used

betweenRuns.variableName2='stimulusDelay';

betweenRuns.variableList2=0.05;

betweenRuns.randomizeSequence=2; % not random sequence

stimulusParameters.targetType='OHIO';

stimulusParameters.targetPhase='sin';

stimulusParameters.targetFrequency=experiment.OHIOfrequencies;

stimulusParameters.targetDuration=0.01; % overruled by OHIO program

stimulusParameters.targetLevel=stimulusParameters.WRVstartValues(1);

stimulusParameters.rampDuration=0.005;

% instructions to user

%   single interval up/down no cue

stimulusParameters.instructions{1}= [{'YES if you hear the tone clearly'}, { }, { 'NO if not (or you are uncertain'}];

%   single interval up/down with cue

stimulusParameters.instructions{2}= [{'count the tones you hear clearly'}, { }, { 'ignore indistinct tones'}];

function paradigm_OHIOrand(handles)

global stimulusParameters experiment betweenRuns

paradigm_training(handles) % default

% find the  threshold for a tonecomplex  consisting of a sequence of 10-ms tones

%   whose frequencies are chosen at random from a list (OHIOfrequencies)

% All tones are presented at the same level (SL) computed using absolute

%   threshols specified in OHIOthresholds;

% The duration of the complex is increased across runs and the number of tones is

%   controlled by OHIOdurations, (for each 20 ms a further tone is added.

% The frequency of the tones is changed on each trial

% fetch thresholds and frequencies

experiment=OHIOthresholds(experiment);

stimulusParameters.WRVname='targetLevel';

stimulusParameters.WRVstartValues=0 ;

stimulusParameters.WRVsteps=[10 2];

stimulusParameters.WRVlimits=[-30 110];

betweenRuns.variableName1='numOHIOtones';

betweenRuns.variableList1= [1 2 4 8 12];

betweenRuns.variableName2='stimulusDelay';

betweenRuns.variableList2=0.1;

betweenRuns.randomizeSequence=2; % not random sequence

stimulusParameters.targetType='OHIO';

stimulusParameters.targetPhase='sin';

stimulusParameters.targetFrequency=experiment.OHIOfrequencies;

stimulusParameters.targetDuration=betweenRuns.variableList2;

stimulusParameters.targetLevel=stimulusParameters.WRVstartValues(1);

stimulusParameters.rampDuration=0.005;

% instructions to user

%   single interval up/down no cue

stimulusParameters.instructions{1}= [{'YES if you hear the tone clearly'}, { }, { 'NO if not (or you are uncertain'}];

%   single interval up/down with cue

stimulusParameters.instructions{2}= [{'count the tones you hear clearly'}, { }, { 'ignore indistinct tones'}];

function paradigm_OHIOspect(handles)

global stimulusParameters experiment betweenRuns

paradigm_training(handles) % default

% find the  threshold for a tonecomplex  consisting of a sequence of 10-ms tones

%   whose frequencies are chosen at random from a list (OHIOfrequencies)

% All tones are presented at the same level (SL) computed using absolute

%   threshols specified in OHIOthresholds;

% The duration of the complex is increased across runs and the number of tones is

%   controlled by OHIOdurations, (for each 20 ms a further tone is added.

% The frequency of the tones is changed on each trial

% fetch thresholds and frequencies

experiment=OHIOthresholds(experiment);

stimulusParameters.WRVname='targetLevel';

stimulusParameters.WRVstartValues=0 ;

stimulusParameters.WRVsteps=[10 2];

stimulusParameters.WRVlimits=[-30 110];

betweenRuns.variableName1='numOHIOtones';

betweenRuns.variableList1= [1 2 4 8 12];

betweenRuns.variableName2='stimulusDelay';

betweenRuns.variableList2=0.1;

betweenRuns.randomizeSequence=2; % not random sequence

stimulusParameters.targetType='OHIO';

stimulusParameters.targetPhase='sin';

stimulusParameters.targetFrequency=experiment.OHIOfrequencies;

stimulusParameters.targetDuration=betweenRuns.variableList2;

stimulusParameters.targetLevel=stimulusParameters.WRVstartValues(1);

stimulusParameters.rampDuration=0.005;

% instructions to user

%   single interval up/down no cue

stimulusParameters.instructions{1}= [{'YES if you hear the tone clearly'}, { }, { 'NO if not (or you are uncertain'}];

%   single interval up/down with cue

stimulusParameters.instructions{2}= [{'count the tones you hear clearly'}, { }, { 'ignore indistinct tones'}];

function paradigm_OHIOspectemp(handles)

global stimulusParameters experiment betweenRuns

paradigm_training(handles) % default

% find the  threshold for a tonecomplex  consisting of a sequence of 10-ms tones

%   whose frequencies are chosen at random from a list (OHIOfrequencies)

% All tones are presented at the same level (SL) computed using absolute

%   threshols specified in OHIOthresholds;

% The duration of the complex is increased across runs and the number of tones is

%   controlled by OHIOdurations, (for each 20 ms a further tone is added.

% The frequency of the tones is changed on each trial

% fetch thresholds and frequencies

experiment=OHIOthresholds(experiment);

stimulusParameters.WRVname='targetLevel';

stimulusParameters.WRVstartValues=0 ;

stimulusParameters.WRVsteps=[10 2];

stimulusParameters.WRVlimits=[-30 110];

betweenRuns.variableName1='numOHIOtones';

betweenRuns.variableList1= [1 2 4 8 12];

betweenRuns.variableName2='stimulusDelay';

betweenRuns.variableList2=0.1;

betweenRuns.randomizeSequence=2; % not random sequence

stimulusParameters.targetType='OHIO';

stimulusParameters.targetPhase='sin';

stimulusParameters.targetFrequency=experiment.OHIOfrequencies;

stimulusParameters.targetDuration=betweenRuns.variableList2;

stimulusParameters.targetLevel=stimulusParameters.WRVstartValues(1);

stimulusParameters.rampDuration=0.005;

% instructions to user

%   single interval up/down no cue

stimulusParameters.instructions{1}= [{'YES if you hear the tone clearly'}, { }, { 'NO if not (or you are uncertain'}];

%   single interval up/down with cue

stimulusParameters.instructions{2}= [{'count the tones you hear clearly'}, { }, { 'ignore indistinct tones'}];

function paradigm_OHIOtemp(handles)

global stimulusParameters experiment betweenRuns

paradigm_training(handles) % default

% find the  threshold for a tonecomplex  consisting of a sequence of 10-ms tones

%   whose frequencies are chosen at random from a list (OHIOfrequencies)

% All tones are presented at the same level (SL) computed using absolute

%   threshols specified in OHIOthresholds;

% The duration of the complex is increased across runs and the number of tones is

%   controlled by OHIOdurations, (for each 20 ms a further tone is added.

% The frequency of the tones is changed on each trial

% fetch thresholds and frequencies

experiment=OHIOthresholds(experiment);

stimulusParameters.WRVname='targetLevel';

stimulusParameters.WRVstartValues=0 ;

stimulusParameters.WRVsteps=[10 2];

stimulusParameters.WRVlimits=[-30 110];

% target variable: slope=1, start going down.

stimulusParameters.cueTestDifference=10;

experiment.psyFunSlope= 1;

withinRuns.direction='down';

betweenRuns.variableName1='numOHIOtones';

betweenRuns.variableList1= [1 2 4 8 12];

betweenRuns.variableName2='stimulusDelay';

betweenRuns.variableList2=0.1;

betweenRuns.randomizeSequence=2; % not random sequence

stimulusParameters.targetType='OHIO';

stimulusParameters.targetPhase='sin';

stimulusParameters.targetFrequency=experiment.OHIOfrequencies;

stimulusParameters.targetDuration=betweenRuns.variableList2;

stimulusParameters.targetLevel=stimulusParameters.WRVstartValues(1);

stimulusParameters.rampDuration=0.005;

% instructions to user

%   single interval up/down no cue

stimulusParameters.instructions{1}= [{'YES if you hear the tone clearly'}, { }, { 'NO if not (or you are uncertain'}];

%   single interval up/down with cue

stimulusParameters.instructions{2}= [{'count the tones you hear clearly'}, { }, { 'ignore indistinct tones'}];

function vectorStrength=testAN(targetFrequency,BFlist, levels)

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

%  or receptive field (doReceptiveFields=1)

global experiment method stimulusParameters

global IHC_VResp_VivoParams  IHC_cilia_RPParams IHCpreSynapseParams

global AN_IHCsynapseParams

    global ANoutput ANdt CNoutput ICoutput ICmembraneOutput tauCas

    global ARattenuation MOCattenuation

dbstop if error

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

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

    ['..' filesep 'testPrograms'])

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';

    MAPparamsName=experiment.name;

    showPlotsAndDetails=0;

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

        MAPparamsName, AN_spikesOrProbability);

    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))

function testBM (BMlocations, paramsName,relativeFrequencies)

% 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);

% 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',[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;

            AN_spikesOrProbability='probability';

%             AN_spikesOrProbability='spikes';

            global DRNLoutput MOCattenuation ARattenuation

            MAP1_14(inputSignal, sampleRate, BF, ...

                MAPparamsName, AN_spikesOrProbability);

            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)

%   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

% 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);

    [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')

function testOME(paramsName)

savePath=path;

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

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