rmeddis@32: function testFM(BFlist,paramsName, AN_spikesOrProbability,...
rmeddis@32:     paramChanges)
rmeddis@32: % testFM(1000,'Normal','probability', [])
rmeddis@32: % testFM(1000,'Normal','spikes', [])
rmeddis@32: 
rmeddis@32: % Demonstration is based on Harris and Dallos
rmeddis@29: %   specify whether you want AN 'probability' or 'spikes'
rmeddis@29: %       spikes is more realistic but takes longer
rmeddis@29: %       refractory effect is included only for spikes
rmeddis@29: %
rmeddis@29: 
rmeddis@29: 
rmeddis@32: global inputStimulusParams outerMiddleEarParams DRNLParams
rmeddis@29: global IHC_VResp_VivoParams IHCpreSynapseParams  AN_IHCsynapseParams
rmeddis@38: global  ANprobRateOutput  ANoutput ANtauCas  dtSpikes
rmeddis@29: dbstop if error
rmeddis@29: restorePath=path;
rmeddis@29: addpath (['..' filesep 'MAP'], ['..' filesep 'utilities'], ...
rmeddis@29:     ['..' filesep 'parameterStore'],  ['..' filesep 'wavFileStore'],...
rmeddis@29:     ['..' filesep 'testPrograms'])
rmeddis@29: 
rmeddis@35: if nargin==0
rmeddis@35:     BFlist=1000;
rmeddis@35:     paramsName=('Normal');
rmeddis@35:     AN_spikesOrProbability='spikes';
rmeddis@29:     paramChanges=[];
rmeddis@35: else
rmeddis@35:     if nargin<3
rmeddis@38:         paramChanges=[];
rmeddis@35:     end
rmeddis@29: end
rmeddis@29: 
rmeddis@29: % masker and probe levels are relative to this threshold
rmeddis@29: thresholdAtCF=10; % dB SPL
rmeddis@33: maskerLevels=[-80   10 20 30 40 60 ] + thresholdAtCF;
rmeddis@29: 
rmeddis@32: showPSTHs=1;
rmeddis@29: 
rmeddis@29: sampleRate=50000;
rmeddis@32: dt=1/sampleRate;
rmeddis@29: 
rmeddis@29: % fetch BF from GUI: use only the first target frequency
rmeddis@29: maskerFrequency=BFlist;
rmeddis@29: maskerDuration=.1;
rmeddis@29: 
rmeddis@29: targetFrequency=maskerFrequency;
rmeddis@29: probeLeveldB=20+thresholdAtCF;	% H&D use 20 dB SL/ TMC uses 10 dB SL
rmeddis@29: probeDuration=0.015; % HD use 15 ms probe (fig 3).
rmeddis@29: 
rmeddis@29: rampDuration=.001;  % HD use 1 ms linear ramp
rmeddis@29: initialSilenceDuration=0.02;
rmeddis@29: finalSilenceDuration=0.05;      % useful for showing recovery
rmeddis@29: 
rmeddis@29: figure(7), clf
rmeddis@29: set(gcf,'position',[613    36   360   247])
rmeddis@29: set(gcf,'name', ['forward masking: thresholdAtCF=' num2str(thresholdAtCF)])
rmeddis@29: 
rmeddis@29: if showPSTHs
rmeddis@29:     figure(8), clf
rmeddis@29:     set(gcf,'name', 'Harris and Dallos simulation')
rmeddis@29:     set(gcf,'position',[980    36   380   249])
rmeddis@29: end
rmeddis@29: 
rmeddis@29: % Plot Harris and Dallos result for comparison
rmeddis@32: figure(7)
rmeddis@29: gapDurations=[0.001	0.002	0.005	0.01	0.02	0.05	0.1	0.3];
rmeddis@29: HDmaskerLevels=[+10	+20	+30	+40	+60];
rmeddis@29: HDresponse=[
rmeddis@29:     1 1 1 1 1 1 1 1;
rmeddis@29:     0.8  	0.82	0.81	0.83	0.87	0.95	1	    1;
rmeddis@29:     0.48	0.5	    0.54	0.58	0.7	    0.85	0.95	1;
rmeddis@29:     0.3	    0.31	0.35	0.4	    0.5	    0.68	0.82	0.94;
rmeddis@29:     0.2 	0.27	0.27	0.29	0.42	0.64	0.75	0.92;
rmeddis@29:     0.15	0.17	0.18	0.23	0.3	     0.5	0.6	    0.82];
rmeddis@29: semilogx(gapDurations,HDresponse,'o'), hold on
rmeddis@29: legend(strvcat(num2str(maskerLevels')),'location','southeast')
rmeddis@29: title([ 'masker dB: ' num2str(HDmaskerLevels)])
rmeddis@29: 
rmeddis@29: %% Run the trials
rmeddis@29: maxProbeResponse=0;
rmeddis@29: levelNo=0;
rmeddis@29: resultsMatrix=zeros(length(maskerLevels), length(gapDurations));
rmeddis@29: summaryFiringRates=[];
rmeddis@29: 
rmeddis@29: % initial silence
rmeddis@29: time=dt: dt: initialSilenceDuration;
rmeddis@29: initialSilence=zeros(1,length(time));
rmeddis@29: 
rmeddis@29: % probe
rmeddis@29: time=dt: dt: probeDuration;
rmeddis@29: amp=28e-6*10^(probeLeveldB/20);
rmeddis@29: probe=amp*sin(2*pi.*targetFrequency*time);
rmeddis@29: % ramps
rmeddis@29: % catch rampTime error
rmeddis@29: if rampDuration>0.5*probeDuration, rampDuration=probeDuration/2; end
rmeddis@29: rampTime=dt:dt:rampDuration;
rmeddis@29: % raised cosine ramp
rmeddis@29: ramp=[0.5*(1+cos(2*pi*rampTime/(2*rampDuration)+pi)) ...
rmeddis@29:     ones(1,length(time)-length(rampTime))];
rmeddis@29: %  onset ramp
rmeddis@29: probe=probe.*ramp;
rmeddis@29: %  offset ramp makes complete ramp for probe
rmeddis@29: ramp=fliplr(ramp);
rmeddis@29: % apply ramp mask to probe. Probe does not change below
rmeddis@29: probe=probe.*ramp;
rmeddis@29: 
rmeddis@29: % final silence
rmeddis@29: time=dt: dt: finalSilenceDuration;
rmeddis@29: finalSilence=zeros(1,length(time));
rmeddis@29: 
rmeddis@29: PSTHplotCount=0;
rmeddis@29: longestSignalDuration=initialSilenceDuration + maskerDuration +...
rmeddis@29:     max(gapDurations) + probeDuration + finalSilenceDuration ;
rmeddis@29: for maskerLeveldB=maskerLevels
rmeddis@29:     levelNo=levelNo+1;
rmeddis@29:     allDurations=[];
rmeddis@29:     allFiringRates=[];
rmeddis@29: 
rmeddis@29:     %masker
rmeddis@29:     time=dt: dt: maskerDuration;
rmeddis@29:     masker=sin(2*pi.*maskerFrequency*time);
rmeddis@29:     % masker ramps
rmeddis@29:     if rampDuration>0.5*maskerDuration
rmeddis@29:         % catch ramp duration error
rmeddis@29:         rampDuration=maskerDuration/2;
rmeddis@29:     end
rmeddis@29:     % NB masker ramp (not probe ramp)
rmeddis@29:     rampTime=dt:dt:rampDuration;
rmeddis@29:     % raised cosine ramp
rmeddis@29:     ramp=[0.5*(1+cos(2*pi*rampTime/(2*rampDuration)+pi))...
rmeddis@29:         ones(1,length(time)-length(rampTime))];
rmeddis@29:     %  onset ramp
rmeddis@29:     masker=masker.*ramp;
rmeddis@29:     %  offset ramp
rmeddis@29:     ramp=fliplr(ramp);
rmeddis@29:     % apply ramp
rmeddis@29:     masker=masker.*ramp;
rmeddis@29:     amp=28e-6*10^(maskerLeveldB/20);
rmeddis@29:     maskerPa=amp*masker;
rmeddis@29: 
rmeddis@29:     for gapDuration=gapDurations
rmeddis@29:         time=dt: dt: gapDuration;
rmeddis@29:         gap=zeros(1,length(time));
rmeddis@29: 
rmeddis@29:         inputSignal=...
rmeddis@29:             [initialSilence maskerPa gap probe finalSilence];
rmeddis@38:         %         time=dt:dt:length(inputSignal)*dt;
rmeddis@38:         %         figure(99), plot(time,inputSignal)
rmeddis@29: 
rmeddis@29:         % **********************************  run MAP model
rmeddis@38:         %         showPlotsAndDetails=0;
rmeddis@38:         nChanges=length(paramChanges);
rmeddis@38:         paramChanges{nChanges+1}='AN_IHCsynapseParams.numFibers=	500;';
rmeddis@32:         MAP1_14(inputSignal, 1/dt, targetFrequency, ...
rmeddis@32:             paramsName, AN_spikesOrProbability, paramChanges);
rmeddis@29: 
rmeddis@32:         if strcmp(AN_spikesOrProbability,'probability')
rmeddis@32:             [nFibers c]=size(ANprobRateOutput);
rmeddis@32:             nLSRfibers=nFibers/length(ANtauCas);
rmeddis@32:             ANresponse=ANprobRateOutput(end-nLSRfibers:end,:);
rmeddis@38:             dtSpikes=dt; % no adjustment for spikes speedup
rmeddis@32:         else
rmeddis@32:             [nFibers c]=size(ANoutput);
rmeddis@32:             nLSRfibers=nFibers/length(ANtauCas);
rmeddis@32:             ANresponse=ANoutput(end-nLSRfibers:end,:);
rmeddis@32:         end
rmeddis@29: 
rmeddis@34:         ANresponse=sum(ANresponse);
rmeddis@38:         %         ANresponseTimes=dtSpikes:dtSpikes:length(ANresponse)*dtSpikes;
rmeddis@38:         %         figure(99), plot(ANresponseTimes,ANresponse)
rmeddis@32: 
rmeddis@29:         % analyse results
rmeddis@29:         probeStart=initialSilenceDuration+maskerDuration+gapDuration;
rmeddis@38:         PSTHbinWidth=dtSpikes;
rmeddis@29:         responseDelay=0.005;
rmeddis@29:         probeTimes=probeStart+responseDelay:...
rmeddis@29:             PSTHbinWidth:probeStart+probeDuration+responseDelay;
rmeddis@29:         probeIDX=round(probeTimes/PSTHbinWidth);
rmeddis@29:         probePSTH=ANresponse(probeIDX);
rmeddis@29:         firingRate=mean(probePSTH);
rmeddis@29:         % NB this only works if you start with the lowest level masker
rmeddis@29:         maxProbeResponse=max([maxProbeResponse firingRate]);
rmeddis@29:         allDurations=[allDurations gapDuration];
rmeddis@29:         allFiringRates=[allFiringRates firingRate];
rmeddis@32: 
rmeddis@29:         %% show PSTHs
rmeddis@29:         if showPSTHs
rmeddis@29:             nLevels=length(maskerLevels);
rmeddis@29:             nDurations=length(gapDurations);
rmeddis@29:             figure(8)
rmeddis@29:             PSTHbinWidth=1e-3;
rmeddis@38:             PSTH=UTIL_PSTHmaker(ANresponse, dtSpikes, PSTHbinWidth);
rmeddis@38:             PSTH=PSTH*dtSpikes/PSTHbinWidth;
rmeddis@29:             PSTHplotCount=PSTHplotCount+1;
rmeddis@29:             subplot(nLevels,nDurations,PSTHplotCount)
rmeddis@33:             PSTHtime=PSTHbinWidth:PSTHbinWidth:...
rmeddis@29:                 PSTHbinWidth*length(PSTH);
rmeddis@34:             if strcmp(AN_spikesOrProbability, 'spikes')
rmeddis@38:                 bar(PSTHtime, PSTH/PSTHbinWidth/nFibers)
rmeddis@38:                 %                 ylim([0 500])
rmeddis@34:             else
rmeddis@38:                 bar(PSTHtime, PSTH)
rmeddis@29:                 ylim([0 500])
rmeddis@29:             end
rmeddis@38:             xlim([0 longestSignalDuration])
rmeddis@29:             grid on
rmeddis@29: 
rmeddis@29:             if PSTHplotCount< (nLevels-1) * nDurations+1
rmeddis@29:                 set(gca,'xticklabel',[])
rmeddis@29:             end
rmeddis@29: 
rmeddis@29:             if ~isequal(mod(PSTHplotCount,nDurations),1)
rmeddis@29:                 set(gca,'yticklabel',[])
rmeddis@29:             else
rmeddis@29:                 ylabel([num2str(maskerLevels...
rmeddis@29:                     (round(PSTHplotCount/nDurations) +1))])
rmeddis@29:             end
rmeddis@29: 
rmeddis@29:             if PSTHplotCount<=nDurations
rmeddis@29:                 title([num2str(1000*gapDurations(PSTHplotCount)) 'ms'])
rmeddis@29:             end
rmeddis@33: 
rmeddis@38:             %         figure(99),            bar(PSTHtime, PSTH)
rmeddis@33: 
rmeddis@29:         end % showPSTHs
rmeddis@29: 
rmeddis@29:     end     % gapDurations duration
rmeddis@29:     summaryFiringRates=[summaryFiringRates allFiringRates'];
rmeddis@29: 
rmeddis@29:     figure(7), hold on
rmeddis@29:     semilogx(allDurations, allFiringRates/maxProbeResponse)
rmeddis@29:     ylim([0 1]), hold on
rmeddis@29:     ylim([0 inf]), xlim([min(gapDurations) max(gapDurations)])
rmeddis@29:     xlim([0.001 1])
rmeddis@29:     pause(0.1) % to allow for CTRL/C interrupts
rmeddis@29:     resultsMatrix(levelNo,:)=allFiringRates;
rmeddis@29: end          % maskerLevel
rmeddis@29: 
rmeddis@29: disp('delay/ rates')
rmeddis@29: disp(num2str(round( [1000*allDurations' summaryFiringRates])))
rmeddis@29: 
rmeddis@29: % replot with best adjustment
rmeddis@29: % figure(34), clf% use for separate plot
rmeddis@29: figure(7), clf
rmeddis@29: peakProbe=max(max(resultsMatrix));
rmeddis@29: resultsMatrix=resultsMatrix/peakProbe;
rmeddis@29: semilogx(gapDurations,HDresponse,'o'), hold on
rmeddis@29: title(['FrMsk: probe ' num2str(probeLeveldB)...
rmeddis@29:     'dB SL: peakProbe=' num2str(peakProbe,'%5.0f') ' sp/s'])
rmeddis@29: xlabel('gap duration (s)'), ylabel ('probe response')
rmeddis@29: semilogx(allDurations, resultsMatrix'), ylim([0 1])
rmeddis@29: ylim([0 inf]),
rmeddis@29: xlim([0.001 1])
rmeddis@29: legend(strvcat(num2str(maskerLevels'-thresholdAtCF)), -1)
rmeddis@29: 
rmeddis@29: % ------------------------------------------------- display parameters
rmeddis@29: disp(['parameter file was: ' paramsName])
rmeddis@29: fprintf('\n')
rmeddis@29: % UTIL_showStruct(inputStimulusParams, 'inputStimulusParams')
rmeddis@29: % UTIL_showStruct(outerMiddleEarParams, 'outerMiddleEarParams')
rmeddis@29: % UTIL_showStruct(DRNLParams, 'DRNLParams')
rmeddis@29: % UTIL_showStruct(IHC_VResp_VivoParams, 'IHC_VResp_VivoParams')
rmeddis@29: UTIL_showStruct(IHCpreSynapseParams, 'IHCpreSynapseParams')
rmeddis@29: UTIL_showStruct(AN_IHCsynapseParams, 'AN_IHCsynapseParams')
rmeddis@29: 
rmeddis@29: 
rmeddis@32: path(restorePath);