MAP

% function MAPtwoToneDemo

%

% Demonstration of two-tone suppression in the AN using a

%   F1 tone suppressed by a 1.5*F1 kHz tone

%

global  ANprobRateOutput DRNLoutput

dbstop if error

restorePath=path;

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

    ['..' filesep 'utilities'])

figure(5), clf

disp('')

primaryToneFrequency=2000;

probedBs=[-100 20 :20: 90];

nProbeLevels=length(probedBs);

% disp(['F1 F1 level= ' num2str([primaryToneFrequency probedB])])

suppressorLevels=0:10:80;

suppressorLevels=0:10:80;

lowestBF=250; 	highestBF= 8000; 	numChannels=21;

BFlist=round(logspace(log10(lowestBF), log10(highestBF), numChannels));

suppressorFrequencies=BFlist;

BFchannel=find(BFlist==primaryToneFrequency);

sampleRate= 40000; % Hz (higher sample rate needed for BF>8000 Hz)

dt=1/sampleRate; % seconds

duration=.050;		      % seconds

tone2Duration=duration/2; % s

startSilenceDuration=0.010;

startSilence=...

    zeros(1,startSilenceDuration*sampleRate);

suppressorStartsPTR=...

    round((startSilenceDuration+duration/2)*sampleRate);

probedBCount=0;

for probedB=probedBs

    probedBCount=probedBCount+1;

    peakResponse= zeros(length(suppressorLevels),length(suppressorFrequencies));

    suppFreqCount=0;

    for suppressorFrequency=suppressorFrequencies

        suppFreqCount=suppFreqCount+1;

        suppressorLevelCount=0;

        for suppressorDB=suppressorLevels

            suppressorLevelCount=suppressorLevelCount+1;

            % primary + suppressor

            primaryLevelDB=probedB;

            suppressorLevelDB=suppressorDB;

            % primary BF tone

            time1=dt: dt: duration;

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

            inputSignal=amp*sin(2*pi*primaryToneFrequency*time1);

            rampDuration=.005; rampTime=dt:dt:rampDuration;

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

            inputSignal=inputSignal.*ramp;

            inputSignal=inputSignal.*fliplr(ramp);

            % suppressor

            time2= dt: dt: tone2Duration;

            % B: tone on

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

            inputSignal2=amp*sin(2*pi*suppressorFrequency*time2);

            rampDuration=.005; rampTime=dt:dt:rampDuration;

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

            inputSignal2=inputSignal2.*ramp;

            inputSignal2=inputSignal2.*fliplr(ramp);

            % A: initial silence (delay to suppressor)

            silence=zeros(1,length(time2));

            inputSignal2=[silence inputSignal2];

            % add tone and suppressor components

            inputSignal=inputSignal+inputSignal2;

            inputSignal=...

                [startSilence inputSignal ];

            % run MAP

            MAPparamsName='Normal';

            AN_spikesOrProbability='probability';

            paramChanges={'IHCpreSynapseParams.tauCa=80e-6;'};

            MAP1_14(inputSignal, sampleRate, BFlist, ...

                MAPparamsName, AN_spikesOrProbability, paramChanges);

            response=ANprobRateOutput;

            peakResponse(suppressorLevelCount,suppFreqCount)=...

                mean(response(BFchannel,suppressorStartsPTR:end));

            disp(['F2 level= ', num2str([suppressorFrequency suppressorDB ...

                peakResponse(suppressorLevelCount,suppFreqCount)])])

            %% make movie

            figure(6), subplot(4,1,2)

            axis tight

            set(gca,'nextplot','replacechildren');

            plot(dt:dt:dt*length(inputSignal), inputSignal, 'k')

            title('probe                 -                        suppressor')

            ylim([-.5 .5])

            figure(6), subplot(2,1,2)

            PSTHbinWidth=0.010;

            PSTH= UTIL_PSTHmakerb(...

                response, dt, PSTHbinWidth);

            [nY nX]=size(PSTH);

            time=PSTHbinWidth*(0:nX-1);

            surf(time, BFlist, PSTH)

            zlim([0 500]),

            caxis([0 500])

            shading interp

            colormap(jet)

            set(gca, 'yScale','log')

            xlim([0 max(time)+dt])

            ylim([0 max(BFlist)])

            view([0 90]) % view([-8 68])

            title('probe                 -                        suppressor')

            pause(0.05)

        end

    end

    %% plot matrix

    peakResponse=peakResponse/peakResponse(1,1);

    figure(5), subplot(2,nProbeLevels, probedBCount)

    surf(suppressorFrequencies,suppressorLevels,peakResponse)

    shading interp

    view([0 90])

    set(gca,'Xscale','log')

    xlim([min(suppressorFrequencies) max(suppressorFrequencies)])

    set(gca,'Xtick', [1000  4000],'xticklabel',{'1000', '4000'})

    ylim([min(suppressorLevels) max(suppressorLevels)])

    subplot(2,nProbeLevels,nProbeLevels+probedBCount)

    contour(suppressorFrequencies,suppressorLevels,peakResponse, ...

        [.1:.1:.9 1.1] )

    set(gca,'Xscale','log')

    set(gca,'Xtick', [1000  4000],'xticklabel',{'1000', '4000'})

    set(gcf, 'name',['primaryToneFrequency= ' num2str(primaryToneFrequency)])

    title([num2str(probedB) ' dB'])

end

path=restorePath;

showMapOptions.surfAN=1;       % 3D plot of HSR response 

showMapOptions.PSTHbinwidth=0.002;      % 3D plot of HSR response 

showMapOptions.view=[-14 76];           % 3D plot of HSR response

UTIL_showMAP(showMapOptions)

numChannels=11;

paramChanges={};

showMapOptions.surfAN=0;       % 2D plot of HSR response

UTIL_showMAP(showMapOptions)

function LibermanData= LiebermanMOCdata

LibermanData=[

    2	0.2;

2.1	0.19;

2.2	0.18;

2.3	0.18;

2.4	0.16;

2.5	0.15;

2.6	0.15;

2.7	0.15;

2.8	0.12;

2.9	0.12;

3	0.1;

3.1	0.1;

3.2	0.05;

3.3	0.05;

3.4	0;

3.5	-0.1;

3.6	-0.4;

3.7	-1.2;

3.8	-1.6;

3.9	-1.8;

4	-1.85;

4.1	-2;

4.2	-2.05;

4.3	-2.05;

4.4	-2.15;

4.5	-2.2;

4.6	-2.15;

4.7	-2.1;

4.8	-2.15;

4.9	-2.2;

5	-2.1;

5.1	-2.1;

5.2	-2.25;

5.3	-2.1;

5.4	-2.15;

5.5	-2.1;

5.6	-2.15;

5.7	-2.1;

5.8	-2.2;

5.9	-2.05;

6	-2.15;

6.1	-2.05;

6.2	-2;

6.3	-2.2;

6.4	-2.1;

6.5	-2.05;

6.6	-2.05;

6.7	-2.05;

6.8	-2.2;

6.9	-2.1;

7	-2.05;

7.1	-2.05;

7.2	-0.7;

7.3	-0.1;

7.4	0;

7.5	0.1;

7.6	0.2;

7.7	0.35;

7.8	0.2;

7.9	0.15;

8	0.2;

8.1	0.15;

8.2	0.15;

8.3	0.15;

8.4	0.12;

8.5	0.1;

8.6	0.09;

8.7	0.08;

8.8	0.07;

8.9	0.06;

9	0.05;

];

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

scalar=-15;

plot(LibermanData(:,1)-2.5,scalar*LibermanData(:,2))

% function MAPtwoToneDemo

% Not for distribution but code worth keeping

% Demonstration of two-tone suppression in the AN using a

%   single channel

%

dbstop if error

% create access to all MAP 1_8 facilities

% addpath ('..\modules', '..\utilities',  '..\parameterStore',  '..\wavFileStore' , '..\testPrograms')

addpath (['..' filesep 'modules'], ['..' filesep 'utilities'],  ['..' filesep 'parameterStore'],  ['..' filesep 'wavFileStore'] , ['..' filesep 'testPrograms'])

moduleSequence= 1:7;  	% up to the AN

figure(3), clf

primaryToneFrequency=2000;

suppressorFrequency=primaryToneFrequency*1.5;

primaryDB=30;

suppressors=20:10:80;

frameCount=0;

for suppressorDB=suppressors

    frameCount=frameCount+1;

    lowestBF=1000; 	highestBF= 5000; 	numChannels=30;

    BFlist=round(logspace(log10(lowestBF), log10(highestBF), numChannels));

    BFlist=2000;

    duration=.020;		      % seconds

    sampleRate= 40000; % Hz (higher sample rate needed for BF>8000 Hz)

    dt=1/sampleRate; % seconds

    % for conditionNo=1:3  % probe alone/ suppressor alone/ combined

    for conditionNo=1:3  % probe alone/ suppressor alone/ combined

        switch conditionNo

            case 1

                primaryLevelDB=primaryDB;

                suppressorLevelDB= -100;

                plotGuide.subPlotNo=	1;     % initialize subplot count

                figureTitle=['probe alone:  ' num2str(primaryToneFrequency) ' Hz;   ' num2str(primaryLevelDB) ' dB SPL'];

            case 2

                primaryLevelDB=-100;

                suppressorLevelDB=suppressorDB;

                plotGuide.subPlotNo=	2;     % initialize subplot count

                figureTitle=['suppressor alone:  ' num2str(suppressorFrequency) ' Hz;   ' num2str(suppressorLevelDB) ' dB SPL'];

            case 3

                primaryLevelDB=primaryDB;

                suppressorLevelDB=suppressorDB;

                plotGuide.subPlotNo=	3;     % initialize subplot count

                figureTitle=['probe + suppressor'];

        end

        % primary BF tone

        time1=dt: dt: duration;

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

        inputSignal=amp*sin(2*pi*primaryToneFrequency*time1);

        rampDuration=.005; rampTime=dt:dt:rampDuration;

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

        inputSignal=inputSignal.*ramp;

        inputSignal=inputSignal.*fliplr(ramp);

        % suppressor

        tone2Duration=duration/2; % s

        time2= dt: dt: tone2Duration;

        % B: tone on

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

        inputSignal2=amp*sin(2*pi*suppressorFrequency*time2);

        rampDuration=.005; rampTime=dt:dt:rampDuration;

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

        inputSignal2=inputSignal2.*ramp;

        % A: initial silence (delay to suppressor)

        silence=zeros(1,length(time2));

        inputSignal2=[silence inputSignal2];

        % add tone and suppressor components

        inputSignal=inputSignal+inputSignal2;

        % specify model parameters

        method=MAPparamsDEMO(BFlist, sampleRate);

        % parameter change (must be global to take effect)

        global   AN_IHCsynapseParams

        AN_IHCsynapseParams.mode=	'probability';

%         method.useEfferent=0;

        method.plotGraphs=	1;	   % please plot

        figure(4), subplot(3,1,conditionNo)

        plot(time1, inputSignal, 'k')% *************

        [ANresponse, method, A]=MAPsequenceSeg(inputSignal, method, moduleSequence);

        response{conditionNo}=ANresponse;

% 	F(frameCount) = getframe(gcf);

    end

%     peakResponse=max(max([response{1} response{2} response{3}]));

%     figure(3), subplot(4,1,1)

%     mesh((response{1}), [0 peakResponse])

%     title(['primary: ' num2str(primaryDB) ' dB SPL'])

%     zlim([0 5000]), view([-28 36])

%     

%     figure(3), subplot(4,1,2)

%     mesh((response{2}), [0 peakResponse])

%     title(['suppressor: ' num2str(suppressorDB) ' dB SPL'])

%     zlim([0 5000]),  view([-28 36])

%     

%     figure(3), subplot(4,1,3)

%     mesh((response{3}), [0 peakResponse])

%     title([' primary + suppressor '])

%     zlim([0 5000]), view([-28 36])

%     

%     figure(3), subplot(4,1,4)

%     sum= response{3}-response{1}-response{2};

%     mesh(sum)

%     zlim([-2000 2000])

%     title(' difference matrix (combined - primary - suppressor)')

%     view([-28 20])

%     disp( [num2str([ suppressorDB primaryDB max(max(sum))])] )

%     

%     colormap(jet)

    pause(1)

end

figure(4)

xlabel('time (s)')

figure(3)

xlabel('time (s)')

%  figure(99), movie(F, 1, 1)

% UTIL_showAllMAPStructures

global OMEParams DRNLParams IHC_cilia_RPParams IHCpreSynapseParams

global AN_IHCsynapseParams MacGregorParams MacGregorMultiParams

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

global  savedParamChanges savedBFlist saveAN_spikesOrProbability ...

    saveMAPparamsName savedInputSignal dt dtSpikes ...

    OMEextEarPressure TMoutput OMEoutput DRNLoutput...

    IHC_cilia_output IHCrestingCiliaCond IHCrestingV...

    IHCoutput ANprobRateOutput ANoutput savePavailable saveNavailable ...

    ANtauCas CNtauGk CNoutput ICoutput ICmembraneOutput ICfiberTypeRates...

    MOCattenuation ARattenuation 

function runMAP1_14

% runMAP1_14 is a general purpose test routine that can be adjusted to

% explore different applications of MAP1_14

%

% It is also designed as 'starter' code for building new applications.

%

% A range of options are supplied in the early part of the program

%

% #1

% Identify the file (in 'MAPparamsName') containing the model parameters

%

% #2

% Identify the kind of model required (in 'AN_spikesOrProbability').

%  A full brainstem model ('spikes') can be computed or a shorter model

%  ('probability') that computes only so far as the auditory nerve

%

% #3

% Choose between a tone signal or file input (in 'signalType')

%

% #4

% Set the signal rms level (in leveldBSPL)

%

% #5

% Identify the channels in terms of their best frequencies in the vector

%  BFlist.

%

% #6

% Last minute changes to the model parameters can be made using

%  a cell array of strings, 'paramChanges'.

dbstop if error

restorePath=path;

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

    ['..' filesep 'utilities'])

%%  #1 parameter file name

MAPparamsName='Normal';

%% #2 probability (fast) or spikes (slow) representation: select one

% AN_spikesOrProbability='spikes';

%   or

AN_spikesOrProbability='probability';

%% #3 A. pure tone, B. harmonic sequence or C. speech file input

% comment out unwanted code

% A. tone

sampleRate= 95000;

signalType= 'tones';

toneFrequency= 1000;            % or a pure tone (Hz)

duration=0.500;                 % seconds

beginSilence=0.010;

endSilence=0.020;

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

%   or

% B. harmonic tone (Hz) - useful to demonstrate a broadband sound

% sampleRate= 44100;

% signalType= 'tones';

% toneFrequency= F0:F0:8000;

% duration=0.500;                 % seconds

% beginSilence=0.250;

% endSilence=0.250;

% F0=210;

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

%   or

% C. signalType= 'file';

% fileName='twister_44kHz';

%% #4 rms level

% signal details

leveldBSPL= 80;                  % dB SPL (80 for Lieberman)

%% #5 number of channels in the model

%   21-channel model (log spacing)

numChannels=21;

lowestBF=100; 	highestBF= 6000;

BFlist=round(logspace(log10(lowestBF), log10(highestBF), numChannels));

%   or specify your own channel BFs

% numChannels=1;

% BFlist=toneFrequency;

%% #6 change model parameters

paramChanges={};    % no changes

% Parameter changes can be used to change one or more model parameters

%  *after* the MAPparams file has been read

%  Each string must have the same format as the corresponding line in the

%   file identified in 'MAPparamsName'

% This example declares only one fiber type with a calcium clearance time

% constant of 80e-6 s (HSR fiber) when the probability option is selected.

%   paramChanges={'AN_IHCsynapseParams.ANspeedUpFactor=5;', ...

%     'IHCpreSynapseParams.tauCa=86e-6; '};

paramChanges={ 'IHCpreSynapseParams.tauCa=86e-6; ',...

    'DRNLParams.rateToAttenuationFactorProb = 0;'};

%% delare 'showMap' options to control graphical output

% see UTIL_showMAP for more options

showMapOptions.printModelParameters=1;   % prints all parameters

showMapOptions.showModelOutput=0;       % plot of all stages

showMapOptions.printFiringRates=1;      % prints stage activity levels

showMapOptions.showEfferent=0;          % tracks of AR and MOC

showMapOptions.surfAN=1;       % 2D plot of HSR response

if strcmp(signalType, 'file')

    % needed for labeling plot

    showMapOptions.fileName=fileName;

else

    showMapOptions.fileName=[];

end

%% Generate stimuli

switch signalType

    case 'tones'

        % Create pure tone stimulus

        dt=1/sampleRate; % seconds

        time=dt: dt: duration;

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

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

        inputSignal=amp*inputSignal;

        % apply ramps

        % catch rampTime error

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

        rampTime=dt:dt:rampDuration;

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

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

        inputSignal=inputSignal.*ramp;

        ramp=fliplr(ramp);

        inputSignal=inputSignal.*ramp;

        % add silence

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

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

        inputSignal= [intialSilence inputSignal finalSilence];

    case 'file'

        %% file input simple or mixed

        [inputSignal sampleRate]=wavread(fileName);

        dt=1/sampleRate;

        inputSignal=inputSignal(:,1);

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

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

        amp=targetRMS/rms;

        inputSignal=inputSignal*amp;

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

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

        inputSignal= [intialSilence inputSignal' finalSilence];

end

%% run the model

tic

fprintf('\n')

disp(['Signal duration= ' num2str(length(inputSignal)/sampleRate)])

disp([num2str(numChannels) ' channel model: ' AN_spikesOrProbability])

disp('Computing ...')

MAP1_14(inputSignal, sampleRate, BFlist, ...

    MAPparamsName, AN_spikesOrProbability, paramChanges);

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

UTIL_showMAP(showMapOptions)

if strcmp(signalType,'tones')

    disp(['duration=' num2str(duration)])

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

    disp(['toneFrequency=' num2str(toneFrequency)])

    global DRNLParams

    disp(['attenuation factor =' ...

        num2str(DRNLParams.rateToAttenuationFactor, '%5.3f') ])

    disp(['attenuation factor (probability)=' ...

        num2str(DRNLParams.rateToAttenuationFactorProb, '%5.3f') ])

    disp(AN_spikesOrProbability)

end

disp('paramChanges')

for i=1:length(paramChanges)

    disp(paramChanges{i})

end

toc

path(restorePath)

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

    paramsName,paramChanges)

% testAN generates rate/level functions for AN and brainstem units.

%  also other information like PSTHs, MOC efferent activity levels.

% Vector strength calculations require the computation of period

% histograms. for this reason the sample rate must always be an integer

% multiple of the tone frequency. This applies to both the sampleRate and

% the spikes sampling rate.

% A full 'spikes' model is used.

% paramChanges is a cell array of strings containing MATLAB statements that

% change model parameters. See MAPparamsNormal for examples.

% e.g.

% testAN(1000,1000, -10:10:80,'Normal',{});

global IHC_VResp_VivoParams  IHC_cilia_RPParams IHCpreSynapseParams

global AN_IHCsynapseParams

global ANoutput dtSpikes CNoutput ICoutput ICmembraneOutput ANtauCas

global ARattenuation MOCattenuation

tic

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

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

    ['..' filesep 'testPrograms'])

if nargin<5, paramChanges={'IHCpreSynapseParams.tauCa= [30e-6 90e-6];'}; end

if nargin<4, paramsName='PL'; end

if nargin<3, levels=-10:10:80; end

if nargin==0, probeFrequency=1000; BFlist=1000; end

nLevels=length(levels);

toneDuration=.2;   rampDuration=0.005;   silenceDuration=.02;

localPSTHbinwidth=0.001;

%% guarantee that the sample rate is an interger multiple 

%   of the probe frequency

% we want 5 bins per period for spikes

spikesSampleRate=5*probeFrequency;

% model sample rate must be an integer multiple of this and in the region

% of 50000

sampleRate=spikesSampleRate*round(50000/spikesSampleRate);

dt=1/sampleRate;

% avoid very slow spikes sampling rate

spikesSampleRate=spikesSampleRate*ceil(10000/spikesSampleRate);

%% pre-allocate storage

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

figure(15), clf

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

figure(5), clf

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

drawnow

%% main computational loop (vary level)

levelNo=0;

for leveldB=levels

    levelNo=levelNo+1;

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

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

    %% generate tone and silences

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

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

    inputSignal= ramp.*inputSignal;

    inputSignal=[silence inputSignal];

    %% run the model

    AN_spikesOrProbability='spikes';  

    nExistingParamChanges=length(paramChanges);

    paramChanges{nExistingParamChanges+1}=...

        ['AN_IHCsynapseParams.spikesTargetSampleRate=' ...

        num2str(spikesSampleRate) ';'];

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

        paramsName, AN_spikesOrProbability, paramChanges);

    nTaus=length(ANtauCas);

    %% Auditory nerve evaluate and display (Fig. 5)

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

    % AN HSR

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

    PSTH=UTIL_PSTHmaker(HSRspikes, dtSpikes, 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])

    grid on

    set(gcf,'name',['PSTH: ' 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, dtSpikes);

    % AN - vector strength

    PSTH=sum(CNoutput (11:20,:));

    [PH, binTimes]=UTIL_periodHistogram...

        (PSTH, dtSpikes, probeFrequency);

    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 driven by LSR fibers

    [nCNneurons c]=size(CNoutput);

    nLSRneurons=round(nCNneurons/nTaus);

    CNLSRspikes=CNoutput(1:nLSRneurons,:);

    PSTH=UTIL_PSTHmaker(CNLSRspikes, dtSpikes, localPSTHbinwidth);

    PSTH=sum(PSTH)/nLSRneurons;

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

    CNLSRrate(levelNo)=mean(PSTH)/localPSTHbinwidth;

    %CN HSR

    MacGregorMultiHSRspikes=...

        CNoutput(end-nLSRneurons+1:end,:);

    PSTH=UTIL_PSTHmaker(MacGregorMultiHSRspikes, dtSpikes, localPSTHbinwidth);

    PSTH=sum(PSTH)/nLSRneurons;

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

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

    CNHSRsaturated(levelNo)=mean(PSTH);

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

    bar(PSTHtime,PSTH)

    ylim([0 1000])

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

    cvMMHSR= UTIL_CV(MacGregorMultiHSRspikes, dtSpikes);

    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, dtSpikes, localPSTHbinwidth);

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

    ICLSRsaturated(levelNo)=mean(PSTH)/localPSTHbinwidth;

    %IC HSR

    MacGregorMultiHSRspikes=...

        ICoutput(end-nLSRneurons+1:end,:);

    PSTH=UTIL_PSTHmaker(MacGregorMultiHSRspikes, dtSpikes, localPSTHbinwidth);

    ICHSRsaturated(levelNo)= (sum(PSTH)/nLSRneurons)/toneDuration;

    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 HSR (last of two)

    plot(time,ICmembraneOutput(end-nLSRneurons+1, 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'), hold on

    plot(20*log10(AR), 'r'), hold off

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

    ylim([-30 0])

%     x=input('prompt')

end % level

%% plot with levels on x-axis

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

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

plot(levels,20*log10(AR), 'r'), hold off

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 summary results (Fig 15)

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

if length(levels)>1

    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', 14)

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

if length(levels)>1

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', 14), grid on

% CN rate - level 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)

if length(levels)>1

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

end

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', 14), grid on

% IC rate - level 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)

if length(levels)>1

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

end

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', 14)

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

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

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 testPhaseLocking(paramsName, paramChanges)

if nargin<2

    paramChanges=[];

end

if nargin<1

    paramsName='PL';

end

testFrequencies=[250 500 1000 2000];

levels=0:10:100;

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=...

        temp(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