rmeddis@0: function testRF
rmeddis@0: % testIHC used either for IHC I/O function or receptive field (doReceptiveFields=1)
rmeddis@0: 
rmeddis@0: global experiment method stimulusParameters expGUIhandles
rmeddis@0: global inputStimulusParams  IHC_ciliaParams
rmeddis@0: global IHC_VResp_VivoParams IHCpreSynapseParams  AN_IHCsynapseParams
rmeddis@0: dbstop if error
rmeddis@0: % set(expGUIhandles.pushbuttonStop, 'backgroundColor', [.941 .941 .941])
rmeddis@0: 
rmeddis@0: addpath (['..' filesep 'MAP'], ['..' filesep 'utilities'], ...
rmeddis@0:     ['..' filesep 'parameterStore'],  ['..' filesep 'wavFileStore'],...
rmeddis@0:     ['..' filesep 'testPrograms'])
rmeddis@0: 
rmeddis@0: targetFrequency=stimulusParameters.targetFrequency(1);
rmeddis@0: 
rmeddis@0: sampleRate=50000;
rmeddis@0: doReceptiveFields=1;
rmeddis@0: 
rmeddis@0: toneDuration=.05;
rmeddis@0: rampDuration=0.004;
rmeddis@0: silenceDuration=.02;
rmeddis@0: 
rmeddis@0: nRepeats=100;   % no. of AN fibers
rmeddis@0: 
rmeddis@0: plotGraphsForIHC=1;
rmeddis@0: % number of MacGregor units is set in the parameter file.
rmeddis@0: 
rmeddis@0: if doReceptiveFields
rmeddis@0:     % show all receptive field
rmeddis@0:     frequencies=targetFrequency*    [  0.5         0.7         0.9     1       1.1         1.3         1.6];
rmeddis@0:     levels=0:20:80; nLevels=length(levels);
rmeddis@0:     figure(14), clf
rmeddis@0:     figure(15), clf
rmeddis@0: else
rmeddis@0:     % show only I/O function at BF
rmeddis@0:     frequencies=targetFrequency;
rmeddis@0:     levels=-20:10:90;
rmeddis@0:     %         levels=10:.25:13;
rmeddis@0:     %         levels=-20:1:-15
rmeddis@0:     nLevels=length(levels);
rmeddis@0: %     figure(13), clf,
rmeddis@0: %     set (gcf, 'name', ['IHC/AN  input/output' num2str(AN_IHCsynapseParams.numFibers) ' repeats'])
rmeddis@0: %     drawnow
rmeddis@0: end
rmeddis@0: nFrequencies=length(frequencies);
rmeddis@0: 
rmeddis@0: IHC_RP_peak=zeros(nLevels,nFrequencies);
rmeddis@0: IHC_RP_min=zeros(nLevels,nFrequencies);
rmeddis@0: IHC_RP_dc=zeros(nLevels,nFrequencies);
rmeddis@0: AN_HSRonset=zeros(nLevels,nFrequencies);
rmeddis@0: AN_HSRsaturated=zeros(nLevels,nFrequencies);
rmeddis@0: AN_LSRonset=zeros(nLevels,nFrequencies);
rmeddis@0: AN_LSRsaturated=zeros(nLevels,nFrequencies);
rmeddis@0: CNLSRsaturated=zeros(nLevels,nFrequencies);
rmeddis@0: CNHSRsaturated=zeros(nLevels,nFrequencies);
rmeddis@0: ICHSRsaturated=zeros(nLevels,nFrequencies);
rmeddis@0: ICLSRsaturated=zeros(nLevels,nFrequencies);
rmeddis@0: 
rmeddis@0: 
rmeddis@0: levelNo=0; PSTHplotCount=0;
rmeddis@0: for leveldB=levels
rmeddis@0:     fprintf('%4.0f\t', leveldB)
rmeddis@0:     levelNo=levelNo+1;
rmeddis@0:     amp=28e-6*10^(leveldB/20);
rmeddis@0:     
rmeddis@0:     freqNo=0;
rmeddis@0:     for frequency=frequencies
rmeddis@0: 
rmeddis@0:         paramFunctionName=['method=MAPparams' experiment.name ...
rmeddis@0:             '(' num2str(targetFrequency) ');' ];
rmeddis@0:         eval(paramFunctionName);  % read parameters afresh each pass
rmeddis@0: 
rmeddis@0:         dt=method.dt;
rmeddis@0:         time=dt:dt:toneDuration;        
rmeddis@0:         rampTime=dt:dt:rampDuration;
rmeddis@0:         ramp=[0.5*(1+cos(2*pi*rampTime/(2*rampDuration)+pi)) ones(1,length(time)-length(rampTime))];
rmeddis@0:         ramp=ramp.*fliplr(ramp);
rmeddis@0:         
rmeddis@0:         silence=zeros(1,round(silenceDuration/dt));
rmeddis@0:         
rmeddis@0:         toneStartptr=length(silence)+1;
rmeddis@0:         toneMidptr=toneStartptr+round(toneDuration/(2*dt)) -1;
rmeddis@0:         toneEndptr=toneStartptr+round(toneDuration/dt) -1;
rmeddis@0:         
rmeddis@0:         % create signal (leveldB/ frequency)
rmeddis@0:         freqNo=freqNo+1;
rmeddis@0:         inputSignal=amp*sin(2*pi*frequency'*time);
rmeddis@0:         inputSignal= ramp.*inputSignal;
rmeddis@0:         inputSignal=[silence inputSignal silence];
rmeddis@0:         
rmeddis@0:         if doReceptiveFields  % receptive field
rmeddis@0:             method.plotGraphs=	0; % plot only PSTHs
rmeddis@0:         else
rmeddis@0:             method.plotGraphs=	plotGraphsForIHC; % show progress
rmeddis@0:         end
rmeddis@0:         
rmeddis@0:         targetChannelNo=1;
rmeddis@0:         
rmeddis@0:         % force parameters
rmeddis@0:          % the number of AN fibers at each BF
rmeddis@0:         AN_IHCsynapseParams.numFibers=	nRepeats;
rmeddis@0:         AN_IHCsynapseParams. mode= 'spikes';
rmeddis@0:         AN_IHCsynapseParams.plotSynapseContents=0;
rmeddis@0:         AN_IHCsynapseParams.PSTHbinWidth=.001;
rmeddis@0:         
rmeddis@0:         method.DRNLSave=1;
rmeddis@0:         method.IHC_cilia_RPSave=1;
rmeddis@0:         method.PSTHbinWidth=1e-3; % useful 1-ms default for all PSTHs
rmeddis@0:         method.AN_IHCsynapseSave=1;
rmeddis@0:         method.MacGregorMultiSave=1;
rmeddis@0:         method.MacGregorSave=1;
rmeddis@0:         method.dt=dt;
rmeddis@0:               
rmeddis@0:         moduleSequence=[1:8];       
rmeddis@0: 
rmeddis@0:         global ANdt ARAttenuation TMoutput OMEoutput ...
rmeddis@0:     DRNLoutput IHC_cilia_output IHCrestingCiliaCond IHCrestingV...
rmeddis@0:     IHCoutput ANprobRateOutput ANoutput savePavailable tauCas  ...
rmeddis@0:     CNoutput  ICoutput ICmembraneOutput ICfiberTypeRates MOCattenuation
rmeddis@0: 
rmeddis@0: AN_spikesOrProbability='spikes';
rmeddis@0: AN_spikesOrProbability='probability';
rmeddis@0: MAPparamsName='Normal';
rmeddis@0: 
rmeddis@0: MAP1_14(inputSignal, 1/dt, targetFrequency, ...
rmeddis@0:     MAPparamsName, AN_spikesOrProbability);
rmeddis@0:         
rmeddis@0:         % RP
rmeddis@0:         IHC_RPData=IHC_cilia_output;
rmeddis@0:         IHC_RPData=IHCoutput(targetChannelNo,:);
rmeddis@0:         IHC_RP_peak(levelNo,freqNo)=max(IHC_RPData(toneStartptr:toneEndptr));
rmeddis@0:         IHC_RP_min(levelNo,freqNo)=min(IHC_RPData(toneStartptr:toneEndptr));
rmeddis@0:         IHC_RP_dc(levelNo,freqNo)=mean(IHC_RPData(toneStartptr:toneEndptr));
rmeddis@0:         
rmeddis@0:         % AN next
rmeddis@0:         AN_IHCsynapseAllData=ANoutput;
rmeddis@0:         method.PSTHbinWidth=0.001;
rmeddis@0:         
rmeddis@0:         nTaus=length(tauCas);
rmeddis@0:         numANfibers=size(ANoutput,1);
rmeddis@0:         numLSRfibers=numANfibers/nTaus;
rmeddis@0:         
rmeddis@0:         %LSR (same as HSR if no LSR fibers present)
rmeddis@0:         channelPtr1=(targetChannelNo-1)*numANfibers+1;
rmeddis@0:         channelPtr2=channelPtr1+numANfibers-1;
rmeddis@0:         LSRspikes=AN_IHCsynapseAllData(channelPtr1:channelPtr2,:);
rmeddis@0:         method.dt=method.AN_IHCsynapsedt;
rmeddis@0:         PSTH=UTIL_PSTHmaker(LSRspikes, method);
rmeddis@0:         PSTH=sum(PSTH,1); % sum across fibers (HSR only)
rmeddis@0:         PSTHStartptr=round(silenceDuration/method.PSTHbinWidth)+1;
rmeddis@0:         PSTHMidptr=PSTHStartptr+round(toneDuration/(2*method.PSTHbinWidth)) -1;
rmeddis@0:         PSTHEndptr=PSTHStartptr+round(toneDuration/method.PSTHbinWidth) -1;
rmeddis@0:         AN_LSRonset(levelNo,freqNo)=max(max(PSTH))/(method.PSTHbinWidth*method.numANfibers);
rmeddis@0:         AN_LSRsaturated(levelNo,freqNo)=sum(PSTH(PSTHMidptr:PSTHEndptr))/(method.numANfibers*toneDuration/2);
rmeddis@0:         
rmeddis@0:         % HSR
rmeddis@0:         channelPtr1=numLSRfibers+(targetChannelNo-1)*method.numANfibers+1;
rmeddis@0:         channelPtr2=channelPtr1+method.numANfibers-1;
rmeddis@0:         HSRspikes=AN_IHCsynapseAllData(channelPtr1:channelPtr2,:);
rmeddis@0:         method.dt=method.AN_IHCsynapsedt;
rmeddis@0:         PSTH=UTIL_PSTHmaker(HSRspikes, method);
rmeddis@0:         PSTH=sum(PSTH,1); % sum across fibers (HSR only)
rmeddis@0:         PSTHStartptr=round(silenceDuration/method.PSTHbinWidth)+1;
rmeddis@0:         PSTHMidptr=PSTHStartptr+round(toneDuration/(2*method.PSTHbinWidth)) -1;
rmeddis@0:         PSTHEndptr=PSTHStartptr+round(toneDuration/method.PSTHbinWidth) -1;
rmeddis@0:         AN_HSRonset(levelNo,freqNo)=max(max(PSTH))/(method.PSTHbinWidth*method.numANfibers);
rmeddis@0:         AN_HSRsaturated(levelNo,freqNo)=sum(PSTH(PSTHMidptr:PSTHEndptr))/(method.numANfibers*toneDuration/2);
rmeddis@0:         [cvANHSR, cvTimes, allTimeStamps, allISIs]=  UTIL_CV(HSRspikes, method.AN_IHCsynapsedt);
rmeddis@0:         
rmeddis@0:         PSTHplotCount=PSTHplotCount+1;
rmeddis@0:         if doReceptiveFields  % receptive field for HSR only
rmeddis@0:             figure(14), set(gcf,'name','AN')
rmeddis@0:             plotReceptiveFields(method, PSTH, PSTHplotCount, levels, frequencies)
rmeddis@0:             ylim([0 method.numANfibers])
rmeddis@0:             xlabel(['CV= ' num2str(max(cvANHSR),'%4.2f')],'fontsize',8)
rmeddis@0:         end % doReceptiveFields
rmeddis@0:         
rmeddis@0:         % CN
rmeddis@0:         MacGregorMultiAllData=method.MacGregorMultiData;
rmeddis@0:         numLSRfibers=method.McGMultinNeuronsPerBF*length(method.nonlinCF)* (nTaus-1);
rmeddis@0:         
rmeddis@0:         %LSR (same as HSR if no LSR fibers present)
rmeddis@0:         channelPtr1=(targetChannelNo-1)*method.McGMultinNeuronsPerBF+1;
rmeddis@0:         channelPtr2=channelPtr1+method.McGMultinNeuronsPerBF-1;
rmeddis@0:         MacGregorMultiLSRspikes=MacGregorMultiAllData(channelPtr1:channelPtr2,:);
rmeddis@0:         method.dt=method.MacGregorMultidt;
rmeddis@0:         PSTH=UTIL_PSTHmaker(MacGregorMultiLSRspikes, method);
rmeddis@0:         PSTH=sum(PSTH,1); % sum across fibers (HSR only)
rmeddis@0:         PSTHStartptr=round(silenceDuration/method.PSTHbinWidth)+1;
rmeddis@0:         PSTHMidptr=PSTHStartptr+round(toneDuration/(2*method.PSTHbinWidth)) -1;
rmeddis@0:         PSTHEndptr=PSTHStartptr+round(toneDuration/method.PSTHbinWidth) -1;
rmeddis@0:         CNLSRsaturated(levelNo,freqNo)=sum(PSTH(PSTHMidptr:PSTHEndptr));
rmeddis@0:         CNLSRsaturated(levelNo,freqNo)=CNLSRsaturated(levelNo,freqNo)...
rmeddis@0:             /((toneDuration/2)*method.McGMultinNeuronsPerBF);
rmeddis@0:         
rmeddis@0:         %HSR
rmeddis@0:         channelPtr1=numLSRfibers+(targetChannelNo-1)*method.McGMultinNeuronsPerBF+1;
rmeddis@0:         channelPtr2=channelPtr1+method.McGMultinNeuronsPerBF-1;
rmeddis@0:         MacGregorMultiHSRspikes=MacGregorMultiAllData(channelPtr1:channelPtr2,:);
rmeddis@0:         method.dt=method.MacGregorMultidt;
rmeddis@0:         PSTH=UTIL_PSTHmaker(MacGregorMultiHSRspikes, method);
rmeddis@0:         PSTH=sum(PSTH,1); % sum across fibers (HSR only)
rmeddis@0:         PSTHStartptr=round(silenceDuration/method.PSTHbinWidth)+1;
rmeddis@0:         PSTHMidptr=PSTHStartptr+round(toneDuration/(2*method.PSTHbinWidth)) -1;
rmeddis@0:         PSTHEndptr=PSTHStartptr+round(toneDuration/method.PSTHbinWidth) -1;
rmeddis@0:         CNHSRsaturated(levelNo,freqNo)=sum(PSTH(PSTHMidptr:PSTHEndptr));
rmeddis@0:         CNHSRsaturated(levelNo,freqNo)=CNHSRsaturated(levelNo,freqNo)...
rmeddis@0:             /((toneDuration/2)*method.McGMultinNeuronsPerBF);
rmeddis@0:         [cvMMHSR, cvTimes, allTimeStamps, allISIs]=  UTIL_CV(MacGregorMultiHSRspikes, method.MacGregorMultidt);
rmeddis@0:         
rmeddis@0:         if doReceptiveFields  % receptive field
rmeddis@0:             figure(15), set(gcf,'name','CN HSR input')
rmeddis@0:             plotReceptiveFields(method, PSTH, PSTHplotCount, levels, frequencies)
rmeddis@0:             ylim([0 method.McGMultinNeuronsPerBF])
rmeddis@0:             xlabel(['CV= ' num2str(max(cvMMHSR),'%4.2f')],'fontsize',8)
rmeddis@0:         end
rmeddis@0:         
rmeddis@0:         MacGregorAllData=method.MacGregorData;
rmeddis@0:         numLSRfibers=length(method.nonlinCF)* (nTaus-1);
rmeddis@0:         
rmeddis@0:         %LSR (same as HSR if no LSR fibers present)
rmeddis@0:         channelPtr1=targetChannelNo;
rmeddis@0:         MacGregorLSR=MacGregorAllData(channelPtr1,:);
rmeddis@0:         method.dt=method.MacGregordt;
rmeddis@0:         PSTH=UTIL_PSTHmaker(MacGregorLSR, method);
rmeddis@0:         PSTHStartptr=round(silenceDuration/method.PSTHbinWidth)+1;
rmeddis@0:         PSTHMidptr=PSTHStartptr+round(toneDuration/(2*method.PSTHbinWidth)) -1;
rmeddis@0:         PSTHEndptr=PSTHStartptr+round(toneDuration/method.PSTHbinWidth) -1;
rmeddis@0:         ICLSRsaturated(levelNo,freqNo)=sum(PSTH(PSTHMidptr:PSTHEndptr));
rmeddis@0:         ICLSRsaturated(levelNo,freqNo)=ICLSRsaturated(levelNo,freqNo)/(toneDuration/2);
rmeddis@0:         
rmeddis@0:         %LSR (same as HSR if no LSR fibers present)
rmeddis@0:         channelPtr1=numLSRfibers+targetChannelNo;
rmeddis@0:         MacGregorHSR=MacGregorAllData(channelPtr1,:);
rmeddis@0:         method.dt=method.MacGregordt;
rmeddis@0:         PSTH=UTIL_PSTHmaker(MacGregorHSR, method);
rmeddis@0:         PSTHStartptr=round(silenceDuration/method.PSTHbinWidth)+1;
rmeddis@0:         PSTHMidptr=PSTHStartptr+round(toneDuration/(2*method.PSTHbinWidth)) -1;
rmeddis@0:         PSTHEndptr=PSTHStartptr+round(toneDuration/method.PSTHbinWidth) -1;
rmeddis@0:         ICHSRsaturated(levelNo,freqNo)=sum(PSTH(PSTHMidptr:PSTHEndptr));
rmeddis@0:         ICHSRsaturated(levelNo,freqNo)=ICHSRsaturated(levelNo,freqNo)/(toneDuration/2);
rmeddis@0:         [cvICHSR, cvTimes, allTimeStamps, allISIs]=  UTIL_CV(MacGregorHSR, method.MacGregordt);
rmeddis@0:         
rmeddis@0: %         if doReceptiveFields  % receptive field
rmeddis@0: %             figure(16), set(gcf,'name','IC HSR input')
rmeddis@0: %             plotReceptiveFields(method, PSTH, PSTHplotCount, levels, frequencies)
rmeddis@0: %             ylim([0 method.McGMultinNeuronsPerBF])
rmeddis@0: %             xlabel(['CV= ' num2str(max(cvICHSR),'%4.2f')],'fontsize',8)
rmeddis@0: %         end
rmeddis@0:     end % frequency
rmeddis@0: end % level
rmeddis@0: fprintf('\n')
rmeddis@0: toneDuration=2;
rmeddis@0: rampDuration=0.004;
rmeddis@0: silenceDuration=.02;
rmeddis@0: nRepeats=200;   % no. of AN fibers
rmeddis@0: fprintf('toneDuration  %6.3f\n', toneDuration)
rmeddis@0: fprintf(' %6.0f  AN fibers (repeats)\n', nRepeats)
rmeddis@0: fprintf('levels')
rmeddis@0: fprintf('%6.2f\t', levels)
rmeddis@0: fprintf('\n')
rmeddis@0: 
rmeddis@0: 
rmeddis@0: % ---------------------------------------------------------- display parameters
rmeddis@0: disp(['parameter file was: ' experiment.name])
rmeddis@0: fprintf('\n')
rmeddis@0: UTIL_showStruct(IHC_VResp_VivoParams, 'IHC_cilia_RPParams')
rmeddis@0: UTIL_showStruct(IHCpreSynapseParams, 'IHCpreSynapseParams')
rmeddis@0: UTIL_showStruct(AN_IHCsynapseParams, 'AN_IHCsynapseParams')
rmeddis@0: 
rmeddis@0: 
rmeddis@0: 
rmeddis@0: function plotReceptiveFields(method, PSTH, PSTHplotCount,  levels, frequencies)
rmeddis@0: 
rmeddis@0: % show PSTH for each level/frequency combination
rmeddis@0: nLevels=length(levels);
rmeddis@0: nFrequencies=length(frequencies);
rmeddis@0: 
rmeddis@0: PSTHtime=method.PSTHbinWidth:method.PSTHbinWidth:method.PSTHbinWidth*length(PSTH);
rmeddis@0: subplot(nLevels,nFrequencies,PSTHplotCount)
rmeddis@0: bar(PSTHtime, PSTH)
rmeddis@0: xlim([0 max(PSTHtime)])
rmeddis@0: % write axis labels only at left and bottom
rmeddis@0: if PSTHplotCount< (nLevels-1) * nFrequencies+1
rmeddis@0:     set(gca,'xticklabel',[])
rmeddis@0: end
rmeddis@0: if ~isequal(mod(PSTHplotCount,nFrequencies),1)
rmeddis@0:     set(gca,'yticklabel',[])
rmeddis@0: else
rmeddis@0:     ylabel([num2str(levels(round(PSTHplotCount/nFrequencies) +1)) ' dB'])
rmeddis@0: end
rmeddis@0: % add titles only on top row
rmeddis@0: if PSTHplotCount<=nFrequencies
rmeddis@0:     title([num2str(frequencies(PSTHplotCount)) ' Hz'])
rmeddis@0: end