MAP

%

% 1. Measurements using real listeners

global checkForPreviousGUI  % holds screen positioning across repeated calls

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

%Paradigm: read in all relevant parameters

set(handles.popupmenuRandomize,'value', betweenRuns.randomizeSequence)

global   stimulusParameters betweenRuns variableNames

    end

        showModelPushButtons(handles, 0)

        showModelPushButtons(handles, 1)

        showModelPushButtons(handles, 0)

global checkForPreviousGUI % holds screen positioning across repeated calls

        if strcmp(errormsg,'manually stopped')

        if strcmp(errormsg,'manually stopped')

        pause(.1)

        if strcmp(errormsg,'manually stopped')

        pause(.1)

        IFMCs=resultsTable(2:end,2:end);

        offBFs=resultsTable(2:end,1);

        % reset original paradigm

        optionNo=strmatch('profile',paradigmNames);

        set(handles.popupmenuParadigm,'value',optionNo);

        aParadigmSelection(handles)

        save profile longTone shortTone gaps BFs TMC offBFs IFMCs

        plotProfile(longTone,shortTone,gaps,BFs,TMC,offBFs,IFMCs)

        if strcmp(errormsg,'manually stopped')

        end

% set(handles.editparamChanges,'visible','off')

global experiment

global variableNames  LevittControl paramChanges

            'MAPmodelSingleCh', 'MAPmodelListen'}

        % MAPmodel writes forced parameter settings to the screen

%  this is used to *reduce* the output signal from what it otherwise

% calibrationdB adjusts these values to compensate for equipment

%  this will change the range. e.g. a 7 dB calibration will yield

        LevittControl.steadyLevittStep= stimulusParameters.WRVsteps(2);

if length(stimulusParameters.WRVstartValues)==1

        end

% identify model parameter changes if any

paramChanges=get(handles.editparamChanges,'string');

eval(paramChanges);

        set(handles.editcalibrationdB,'string','0')

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

global experiment paramChanges

testOME(experiment.name, paramChanges);

global  stimulusParameters experiment paramChanges

AN_spikesOrProbability='probability';

    experiment.name,relativeFrequencies, AN_spikesOrProbability, ...

    paramChanges);

global stimulusParameters experiment paramChanges

testAN(targetFrequency,BFlist,-10:10:90,experiment.name, paramChanges);

global experiment

testPhaseLocking(experiment.name)

global stimulusParameters experiment paramChanges

% now carry out tests

showPSTHs=0;

targetFrequency=stimulusParameters.targetFrequency(1);

BFlist=targetFrequency;

testSynapse(BFlist,experiment.name, paramChanges)

global stimulusParameters experiment paramChanges

testFM(stimulusParameters.targetFrequency(1),experiment.name, ...

    showPSTHs, paramChanges)

eval(paramChanges);

global experiment stimulusParameters paramChanges

testRP(stimulusParameters.targetFrequency,experiment.name, paramChanges)

% if experiment.stop

%     disp('******** experiment manually stopped  *****************')

%     experiment.status= 'waitingForStart';

%     errormsg='manually stopped';

%     addToMsg(errormsg,1)

%     return

% end

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

    options.printModelParameters=0;

    UTIL_showMAP(options)

stimulusParameters.WRVstartValues=50;

stimulusParameters.WRVstartValues=50;

stimulusParameters.targetFrequency=[250 500 1000 2000 4000 8000 ];

% comparison data (e.g. participants)

% rows are BFs

% -------------------------------------------JSan

compareBFs=[250	500	1000	2000	3000

];

compareLongTone=  [67.46485	56.95655	65.01985	61.46655	73.33265

];

compareShortTone=[	72.3185	63.2818	69.0373	65.2853	76

];

compareGaps=[0.02 0.05 0.08];

compareTMC=	[

    84	69	77	75	93

88	73	81	79	97

95	79	85	83	98

compareMaskerFreqs=[0.7  0.9 1 1.1  1.3 ];

compareIFMCs=[83.1698	77.3165	79.8474	82.9074	82.3294

80.9667	73.6653	80.9446	80.7005	79.0022

82.0135	71.2284	78.7345	74.3342	84.126

79.3348	70.3347	78.5769	79.7274	91.9849

79.2308	76.3164	83.6881	86.2105	NaN

  ];

% % -------------------------------------------JE

% BFs=[250 500 1000 2000 4000 8000];

% compareLongTone=  [32 30	31	40	54 NaN];

% compareShortTone=[	49	50	47	56	63	NaN];

% 

% compareGaps=0.01:0.01:0.09;

% compareTMC=	[

%     69	83	82	NaN	NaN	NaN	NaN	NaN	NaN

%     61	68	79	88	93	NaN	NaN	NaN	NaN

%     63	69	79	84	92	NaN	NaN	NaN	NaN

%     67	71	75	80	82	84	88	93	NaN

%     82	82	86	86	NaN	83	88	90	75

%     NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN

%     ];

% compareTMC=compareTMC';

% 

% compareMaskerFreqs=[0.5	0.7	0.9	1	1.1	1.3	1.6];

% compareIFMCs=[

%     64	60	58	58	57	60	64

%     65	63	58	55	54	59	69

%     68	64	60	59	62	73	79

%     76	75	71	67	68	71	77

%     79	71	68	69	73	75	77

%     76	73	75	75	76	80	NaN

%     ];

% compareIFMCs=compareIFMCs';

% -------------------------------------------CMR

% CMR 

% BFs=[250 500 1000 2000 4000 8000];

% compareLongTone=[	11.4	1.55	-13.5	-6.35	-6.4	7.45];

% compareShortTone=[	23.85	18.9	9.85	10.6	9.55	21.9];

% 

% compareGaps=0.01:0.01:0.09;

% compareTMC=	[

%     28.5	35.0	49.3	70.1	80.5	85.5    NaN      NaN    NaN;

%     31.1	44.3	48.4	59.5	56.4	76.7	70.2	82.4	76.3;

%     33.4	38.4	48.8	55.8	64.5	78.7	84.2	88.3	90.3;

%     25.4	37.0	49.2	49.7	58.2	69.6	87.7	95.8	93.0;

%     18.2	23.5	27.4	41.5	64.3	82.1	86.7	91.2	NaN;

%     32.5	35.8	43.5	52.1	69.1	78.6	86.6	86.0	NaN;

%     ];

% compareTMC=compareTMC';

% 

% compareMaskerFreqs=[0.5	0.7	0.9	1	1.1	1.3	1.6];

% compareIFMCs=[

%     50	42	34	35	34	33	37;

%     58	51	38	33	28	41	49;

%     57	41	27	20	28	37	66;

%     61	49	27	20	34	68	79;

%     67	45	27	22	46	74	87;

%     62	62	43	22	47	56	83;

%     ];

% compareIFMCs=compareIFMCs';

semilogx(compareBFs,compareLongTone,'ko:'); hold on

semilogx(compareBFs,compareShortTone,'bo:'); hold on

    plot(gaps,TMC(:,BFno),'b','lineWidth',3), hold on

    ylim([-10 110])

    grid on

%         text(0.02,80,' TMC','backgroundColor','w')

    title([num2str(BFs(BFno)) ' Hz'])

    set(gca,'XTick',[ 0.1],'xTickLabel', { '0.1'})

for BFno=1:length(compareBFs)

    subplot(2,6,BFno)

    plot(compareGaps,compareTMC(:,BFno)-longTone(BFno),'k:')

    plot(compareGaps,compareTMC(:,BFno),'k:')

    ylim([-10 110])

    xlim([0.01 0.1])

    grid on

    if BFno==1

        ylabel('masker dB SL')

        xlabel('gap')

%         text(0.02,80,' TMC','backgroundColor','w')

    end

    title([num2str(BFs(BFno)) ' Hz'])

    set(gca,'XTick',[ 0.1],'xTickLabel', { '0.1'})

end

% IFMCs

for BFno=1:length(compareBFs)

    compareFreq=compareMaskerFreqs'*BFs(BFno);

    subplot(2,1,2)

    semilogx(compareFreq,compareIFMCs(:,BFno),'k:')

    ylim([0 100])

    xlim([100 12000])

    grid on

end

global paramChanges

        saveFileName=...

            ['savedData/' experiment.name '_' ...

            experiment.date '_' experiment.paradigm];

disp(paramChanges)

disp(paramChanges)

    % show sample Rate on GUI; it must be set in MAPparams ##??

    if experiment.singleShot % ##??

% global outerMiddleEarParams DRNLParams AN_IHCsynapseParams

global ICoutput ANdt dt savedBFlist ANprobRateOutput expGUIhandles

global paramChanges

% important buried constant ##??

% convert from string to a cell array

eval(paramChanges);

function method=MAPparamsJE ...

    (BFlist, sampleRate, showParams, paramChanges)

% MAPparams<> establishes a complete set of MAP parameters

% Parameter file names must be of the form <MAPparams> <name>

%

% input arguments

%  BFlist     (optional) specifies the desired list of channel BFs

%    otherwise defaults set below

%  sampleRate (optional), default is 50000.

%  showParams (optional) =1 prints out the complete set of parameters

% output argument

%  method passes a miscelleny of values

global inputStimulusParams OMEParams DRNLParams IHC_cilia_RPParams

global IHC_VResp_VivoParams IHCpreSynapseParams  AN_IHCsynapseParams

global MacGregorParams MacGregorMultiParams  filteredSACFParams

global experiment % used by calls from multiThreshold only

currentFile=mfilename;                      % i.e. the name of this mfile

method.parameterSource=currentFile(10:end); % for the record

efferentDelay=0.010;

method.segmentDuration=efferentDelay;

if nargin<3, showParams=0; end

if nargin<2, sampleRate=50000; end

if nargin<1 || BFlist(1)<0 % if BFlist= -1, set BFlist to default

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

    % 21 chs (250-8k)includes BFs at 250 500 1000 2000 4000 8000

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

end

% BFlist=1000;

% preserve for backward campatibility

method.nonlinCF=BFlist; 

method.dt=1/sampleRate; 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% set  model parameters

%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%  #1 inputStimulus

inputStimulusParams=[];

inputStimulusParams.sampleRate= sampleRate; 

%%  #2 outerMiddleEar

OMEParams=[];  % clear the structure first

% outer ear resonances band pass filter  [gain lp order  hp]

OMEParams.externalResonanceFilters=      [ 10 1 1000 4000];

% highpass stapes filter  

%  Huber gives 2e-9 m at 80 dB and 1 kHz (2e-13 at 0 dB SPL)

OMEParams.OMEstapesLPcutoff= 1000;

OMEParams.stapesScalar=	     45e-9;

% Acoustic reflex: maximum attenuation should be around 25 dB Price (1966)

% i.e. a minimum ratio of 0.056.

% 'spikes' model: AR based on brainstem spiking activity (LSR)

OMEParams.rateToAttenuationFactor=0.006;   % * N(all ICspikes)

% OMEParams.rateToAttenuationFactor=0;   % * N(all ICspikes)

% 'probability model': Ar based on AN firing probabilities (LSR)

OMEParams.rateToAttenuationFactorProb=0.01;% * N(all ANrates)

% OMEParams.rateToAttenuationFactorProb=0;% * N(all ANrates)

% asymptote should be around 100-200 ms

OMEParams.ARtau=.05; % AR smoothing function

% delay must be longer than the segment length

OMEParams.ARdelay=efferentDelay;  %Moss gives 8.5 ms latency

OMEParams.ARrateThreshold=0;

%%  #3 DRNL

DRNLParams=[];  % clear the structure first

DRNLParams.BFlist=BFlist;

% DRNL nonlinear path

DRNLParams.a=5e4;     % DRNL.a=0 means no OHCs (no nonlinear path)

DRNLParams.a=2e4;     % DRNL.a=0 means no OHCs (no nonlinear path)

DRNLParams.b=8e-6;    % *compression threshold raised compression

% DRNLParams.b=1;    % b=1 means no compression

DRNLParams.c=0.2;     % compression exponent

% nonlinear filters

DRNLParams.nonlinCFs=BFlist;

DRNLParams.nonlinOrder=	3;           % order of nonlinear gammatone filters

p=0.2895;   q=170;  % human  (% p=0.14;   q=366;  % cat)

DRNLParams.nlBWs=  p * BFlist + q;

DRNLParams.p=p;   DRNLParams.q=q;   % save p and q for printing only

% DRNL linear path:

DRNLParams.g=1000;     % linear path gain factor

% linCF is not necessarily the same as nonlinCF

minLinCF=153.13; coeffLinCF=0.7341;   % linCF>nonlinBF for BF < 1 kHz

DRNLParams.linCFs=minLinCF+coeffLinCF*BFlist;

DRNLParams.linOrder=	3;           % order of linear gammatone filters

minLinBW=100; coeffLinBW=0.6531;

DRNLParams.linBWs=minLinBW + coeffLinBW*BFlist; % bandwidths of linear  filters

% DRNL MOC efferents

DRNLParams.MOCdelay = efferentDelay;            % must be < segment length!

% 'spikes' model: MOC based on brainstem spiking activity (HSR)

DRNLParams.rateToAttenuationFactor = .01;  % strength of MOC

%      DRNLParams.rateToAttenuationFactor = 0;  % strength of MOC

% 'probability' model: MOC based on AN spiking activity (HSR)

DRNLParams.rateToAttenuationFactorProb = .0055;  % strength of MOC

% DRNLParams.rateToAttenuationFactorProb = .0;  % strength of MOC

DRNLParams.MOCrateThresholdProb =70;                % spikes/s probability only

DRNLParams.MOCtau =.1;                         % smoothing for MOC

%% #4 IHC_cilia_RPParams

IHC_cilia_RPParams.tc=	0.0003;   % 0.0003 filter time simulates viscocity

% IHC_cilia_RPParams.tc=	0.0005;   % 0.0003 filter time simulates viscocity

IHC_cilia_RPParams.C=	0.03;      % 0.1 scalar (C_cilia ) 

IHC_cilia_RPParams.u0=	5e-9;       

IHC_cilia_RPParams.s0=	30e-9;

IHC_cilia_RPParams.u1=	1e-9;

IHC_cilia_RPParams.s1=	1e-9;

IHC_cilia_RPParams.Gmax= 6e-9;    % 2.5e-9 maximum conductance (Siemens)

IHC_cilia_RPParams.Ga=	1e-9;  % 4.3e-9 fixed apical membrane conductance

IHC_cilia_RPParams.Ga=	.8e-9;  % 4.3e-9 fixed apical membrane conductance

%  #5 IHC_RP

IHC_cilia_RPParams.Cab=	4e-012;         % IHC capacitance (F)

% IHC_cilia_RPParams.Cab=	1e-012;         % IHC capacitance (F)

IHC_cilia_RPParams.Et=	0.05;          % endocochlear potential (V)

IHC_cilia_RPParams.Gk=	2e-008;         % 1e-8 potassium conductance (S)

IHC_cilia_RPParams.Ek=	-0.08;          % -0.084 K equilibrium potential

IHC_cilia_RPParams.Rpc=	0.04;           % combined resistances

%%  #5 IHCpreSynapse

IHCpreSynapseParams=[];

IHCpreSynapseParams.GmaxCa=	14e-9;% maximum calcium conductance

IHCpreSynapseParams.GmaxCa=	12e-9;% maximum calcium conductance

IHCpreSynapseParams.ECa=	0.066;  % calcium equilibrium potential

IHCpreSynapseParams.beta=	400;	% determine Ca channel opening

IHCpreSynapseParams.gamma=	100;	% determine Ca channel opening

IHCpreSynapseParams.tauM=	0.00005;	% membrane time constant ?0.1ms

IHCpreSynapseParams.power=	3;

% reminder: changing z has a strong effect on HF thresholds (like Et)

IHCpreSynapseParams.z=	    2e42;   % scalar Ca -> vesicle release rate

LSRtauCa=35e-6;            HSRtauCa=85e-6;            % seconds

% LSRtauCa=35e-6;            HSRtauCa=70e-6;            % seconds

IHCpreSynapseParams.tauCa= [LSRtauCa HSRtauCa]; %LSR and HSR fiber

%%  #6 AN_IHCsynapse

% c=kym/(y(l+r)+kl)	(spontaneous rate)

% c=(approx)  ym/l  (saturated rate)

AN_IHCsynapseParams=[];             % clear the structure first

AN_IHCsynapseParams.M=	12;         % maximum vesicles at synapse

AN_IHCsynapseParams.y=	4;          % depleted vesicle replacement rate

AN_IHCsynapseParams.y=	6;          % depleted vesicle replacement rate

AN_IHCsynapseParams.x=	30;         % replenishment from re-uptake store

AN_IHCsynapseParams.x=	60;         % replenishment from re-uptake store

% reduce l to increase saturated rate

AN_IHCsynapseParams.l=	100; % *loss rate of vesicles from the cleft

AN_IHCsynapseParams.l=	250; % *loss rate of vesicles from the cleft

AN_IHCsynapseParams.r=	500; % *reuptake rate from cleft into cell

% AN_IHCsynapseParams.r=	300; % *reuptake rate from cleft into cell

AN_IHCsynapseParams.refractory_period=	0.00075;

% number of AN fibers at each BF (used only for spike generation)

AN_IHCsynapseParams.numFibers=	100; 

AN_IHCsynapseParams.TWdelay=0.004;  % ?delay before stimulus first spike

AN_IHCsynapseParams.ANspeedUpFactor=5; % longer epochs for computing spikes.

%%  #7 MacGregorMulti (first order brainstem neurons)

MacGregorMultiParams=[];

MacGregorMultiType='chopper'; % MacGregorMultiType='primary-like'; %choose

switch MacGregorMultiType

    case 'primary-like'

        MacGregorMultiParams.nNeuronsPerBF=	10;   % N neurons per BF

        MacGregorMultiParams.type = 'primary-like cell';

        MacGregorMultiParams.fibersPerNeuron=4;   % N input fibers

        MacGregorMultiParams.dendriteLPfreq=200;  % dendritic filter

        MacGregorMultiParams.currentPerSpike=0.11e-6; % (A) per spike

        MacGregorMultiParams.Cap=4.55e-9;   % cell capacitance (Siemens)

        MacGregorMultiParams.tauM=5e-4;     % membrane time constant (s)

        MacGregorMultiParams.Ek=-0.01;      % K+ eq. potential (V)

        MacGregorMultiParams.dGkSpike=3.64e-5; % K+ cond.shift on spike,S

        MacGregorMultiParams.tauGk=	0.0012; % K+ conductance tau (s)

        MacGregorMultiParams.Th0=	0.01;   % equilibrium threshold (V)

        MacGregorMultiParams.c=	0.01;       % threshold shift on spike, (V)

        MacGregorMultiParams.tauTh=	0.015;  % variable threshold tau

        MacGregorMultiParams.Er=-0.06;      % resting potential (V)

        MacGregorMultiParams.Eb=0.06;       % spike height (V)

    case 'chopper'

        MacGregorMultiParams.nNeuronsPerBF=	10;   % N neurons per BF

        MacGregorMultiParams.type = 'chopper cell';

        MacGregorMultiParams.fibersPerNeuron=10;  % N input fibers

%         MacGregorMultiParams.fibersPerNeuron=6;  % N input fibers

        MacGregorMultiParams.dendriteLPfreq=50;   % dendritic filter

        MacGregorMultiParams.currentPerSpike=35e-9; % *per spike

%         MacGregorMultiParams.currentPerSpike=30e-9; % *per spike

        MacGregorMultiParams.Cap=1.67e-8; % ??cell capacitance (Siemens)

        MacGregorMultiParams.tauM=0.002;  % membrane time constant (s)

        MacGregorMultiParams.Ek=-0.01;    % K+ eq. potential (V)

        MacGregorMultiParams.dGkSpike=1.33e-4; % K+ cond.shift on spike,S

        MacGregorMultiParams.tauGk=	0.0005;% K+ conductance tau (s)

        MacGregorMultiParams.Th0=	0.01; % equilibrium threshold (V)

        MacGregorMultiParams.c=	0;        % threshold shift on spike, (V)

        MacGregorMultiParams.tauTh=	0.02; % variable threshold tau

        MacGregorMultiParams.Er=-0.06;    % resting potential (V)

        MacGregorMultiParams.Eb=0.06;     % spike height (V)

        MacGregorMultiParams.PSTHbinWidth=	1e-4;

end

%%  #8 MacGregor (second-order neuron). Only one per channel

MacGregorParams=[];                 % clear the structure first

MacGregorParams.type = 'chopper cell';

MacGregorParams.fibersPerNeuron=10; % N input fibers

MacGregorParams.dendriteLPfreq=100; % dendritic filter

MacGregorParams.currentPerSpike=120e-9;% *(A) per spike

MacGregorParams.currentPerSpike=40e-9;% *(A) per spike

MacGregorParams.Cap=16.7e-9;        % cell capacitance (Siemens)

MacGregorParams.tauM=0.002;         % membrane time constant (s)

MacGregorParams.Ek=-0.01;           % K+ eq. potential (V)

MacGregorParams.dGkSpike=1.33e-4;   % K+ cond.shift on spike,S

MacGregorParams.tauGk=	0.0005;     % K+ conductance tau (s)

MacGregorParams.Th0=	0.01;       % equilibrium threshold (V)

MacGregorParams.c=	0;              % threshold shift on spike, (V)

MacGregorParams.tauTh=	0.02;       % variable threshold tau

MacGregorParams.Er=-0.06;           % resting potential (V)

MacGregorParams.Eb=0.06;            % spike height (V)

MacGregorParams.debugging=0;        % (special)

% wideband accepts input from all channels (of same fiber type)

% use wideband to create inhibitory units

MacGregorParams.wideband=0;         % special for wideband units

% MacGregorParams.saveAllData=0;

%%  #9 filteredSACF

minPitch=	300; maxPitch=	3000; numPitches=60;    % specify lags

pitches=100*log10(logspace(minPitch/100, maxPitch/100, numPitches));

filteredSACFParams.lags=1./pitches;     % autocorrelation lags vector

filteredSACFParams.acfTau=	.003;       % time constant of running ACF

filteredSACFParams.lambda=	0.12;       % slower filter to smooth ACF

filteredSACFParams.plotFilteredSACF=1;  % 0 plots unfiltered ACFs

filteredSACFParams.plotACFs=0;          % special plot (see code)

%  filteredSACFParams.usePressnitzer=0; % attenuates ACF at  long lags

filteredSACFParams.lagsProcedure=  'useAllLags';

% filteredSACFParams.lagsProcedure=  'useBernsteinLagWeights';

% filteredSACFParams.lagsProcedure=  'omitShortLags';

filteredSACFParams.criterionForOmittingLags=3;

% checks

if AN_IHCsynapseParams.numFibers<MacGregorMultiParams.fibersPerNeuron

    error('MacGregorMulti: too few input fibers for input to MacG unit')

end

%% now accept last minute parameter changes required by the calling program

% paramChanges

if nargin>3 && ~isempty(paramChanges)

    nChanges=length(paramChanges);

    for idx=1:nChanges

        eval(paramChanges{idx})

    end

end

%% write all parameters to the command window

% showParams is currently set at the top of htis function

if showParams

    fprintf('\n %%%%%%%%\n')

    fprintf('\n%s\n', method.parameterSource)

    fprintf('\n')

    nm=UTIL_paramsList(whos);

    for i=1:length(nm)

        %         eval(['UTIL_showStruct(' nm{i} ', ''' nm{i} ''')'])

        if ~strcmp(nm(i), 'method')

            eval(['UTIL_showStructureSummary(' nm{i} ', ''' nm{i} ''', 10)'])

        end

    end

    % highlight parameter changes made locally

    if nargin>3 && ~isempty(paramChanges)

        fprintf('\n Local parameter changes:\n')

        for i=1:length(paramChanges)

            disp(paramChanges{i})

        end

    end

end

% for backward compatibility

experiment.comparisonData=[];

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