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