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