annotate testPrograms/testRF.m @ 38:c2204b18f4a2 tip

End nov big change
author Ray Meddis <rmeddis@essex.ac.uk>
date Mon, 28 Nov 2011 13:34:28 +0000
parents 82fb37eb430e
children
rev   line source
rmeddis@29 1 function testRF
rmeddis@29 2 % testIHC used either for IHC I/O function or receptive field (doReceptiveFields=1)
rmeddis@29 3
rmeddis@29 4 global experiment method stimulusParameters expGUIhandles
rmeddis@29 5 global inputStimulusParams IHC_ciliaParams
rmeddis@29 6 global IHC_VResp_VivoParams IHCpreSynapseParams AN_IHCsynapseParams
rmeddis@29 7 dbstop if error
rmeddis@29 8 % set(expGUIhandles.pushbuttonStop, 'backgroundColor', [.941 .941 .941])
rmeddis@29 9
rmeddis@29 10 addpath (['..' filesep 'MAP'], ['..' filesep 'utilities'], ...
rmeddis@29 11 ['..' filesep 'parameterStore'], ['..' filesep 'wavFileStore'],...
rmeddis@29 12 ['..' filesep 'testPrograms'])
rmeddis@29 13
rmeddis@29 14 targetFrequency=stimulusParameters.targetFrequency(1);
rmeddis@29 15
rmeddis@29 16 sampleRate=50000;
rmeddis@29 17 doReceptiveFields=1;
rmeddis@29 18
rmeddis@29 19 toneDuration=.05;
rmeddis@29 20 rampDuration=0.004;
rmeddis@29 21 silenceDuration=.02;
rmeddis@29 22
rmeddis@29 23 nRepeats=100; % no. of AN fibers
rmeddis@29 24
rmeddis@29 25 plotGraphsForIHC=1;
rmeddis@29 26 % number of MacGregor units is set in the parameter file.
rmeddis@29 27
rmeddis@29 28 if doReceptiveFields
rmeddis@29 29 % show all receptive field
rmeddis@29 30 frequencies=targetFrequency* [ 0.5 0.7 0.9 1 1.1 1.3 1.6];
rmeddis@29 31 levels=0:20:80; nLevels=length(levels);
rmeddis@29 32 figure(14), clf
rmeddis@29 33 figure(15), clf
rmeddis@29 34 else
rmeddis@29 35 % show only I/O function at BF
rmeddis@29 36 frequencies=targetFrequency;
rmeddis@29 37 levels=-20:10:90;
rmeddis@29 38 % levels=10:.25:13;
rmeddis@29 39 % levels=-20:1:-15
rmeddis@29 40 nLevels=length(levels);
rmeddis@29 41 % figure(13), clf,
rmeddis@29 42 % set (gcf, 'name', ['IHC/AN input/output' num2str(AN_IHCsynapseParams.numFibers) ' repeats'])
rmeddis@29 43 % drawnow
rmeddis@29 44 end
rmeddis@29 45 nFrequencies=length(frequencies);
rmeddis@29 46
rmeddis@29 47 IHC_RP_peak=zeros(nLevels,nFrequencies);
rmeddis@29 48 IHC_RP_min=zeros(nLevels,nFrequencies);
rmeddis@29 49 IHC_RP_dc=zeros(nLevels,nFrequencies);
rmeddis@29 50 AN_HSRonset=zeros(nLevels,nFrequencies);
rmeddis@29 51 AN_HSRsaturated=zeros(nLevels,nFrequencies);
rmeddis@29 52 AN_LSRonset=zeros(nLevels,nFrequencies);
rmeddis@29 53 AN_LSRsaturated=zeros(nLevels,nFrequencies);
rmeddis@29 54 CNLSRsaturated=zeros(nLevels,nFrequencies);
rmeddis@29 55 CNHSRsaturated=zeros(nLevels,nFrequencies);
rmeddis@29 56 ICHSRsaturated=zeros(nLevels,nFrequencies);
rmeddis@29 57 ICLSRsaturated=zeros(nLevels,nFrequencies);
rmeddis@29 58
rmeddis@29 59
rmeddis@29 60 levelNo=0; PSTHplotCount=0;
rmeddis@29 61 for leveldB=levels
rmeddis@29 62 fprintf('%4.0f\t', leveldB)
rmeddis@29 63 levelNo=levelNo+1;
rmeddis@29 64 amp=28e-6*10^(leveldB/20);
rmeddis@29 65
rmeddis@29 66 freqNo=0;
rmeddis@29 67 for frequency=frequencies
rmeddis@29 68
rmeddis@29 69 paramFunctionName=['method=MAPparams' experiment.name ...
rmeddis@29 70 '(' num2str(targetFrequency) ');' ];
rmeddis@29 71 eval(paramFunctionName); % read parameters afresh each pass
rmeddis@29 72
rmeddis@29 73 dt=method.dt;
rmeddis@29 74 time=dt:dt:toneDuration;
rmeddis@29 75 rampTime=dt:dt:rampDuration;
rmeddis@29 76 ramp=[0.5*(1+cos(2*pi*rampTime/(2*rampDuration)+pi)) ones(1,length(time)-length(rampTime))];
rmeddis@29 77 ramp=ramp.*fliplr(ramp);
rmeddis@29 78
rmeddis@29 79 silence=zeros(1,round(silenceDuration/dt));
rmeddis@29 80
rmeddis@29 81 toneStartptr=length(silence)+1;
rmeddis@29 82 toneMidptr=toneStartptr+round(toneDuration/(2*dt)) -1;
rmeddis@29 83 toneEndptr=toneStartptr+round(toneDuration/dt) -1;
rmeddis@29 84
rmeddis@29 85 % create signal (leveldB/ frequency)
rmeddis@29 86 freqNo=freqNo+1;
rmeddis@29 87 inputSignal=amp*sin(2*pi*frequency'*time);
rmeddis@29 88 inputSignal= ramp.*inputSignal;
rmeddis@29 89 inputSignal=[silence inputSignal silence];
rmeddis@29 90
rmeddis@29 91 if doReceptiveFields % receptive field
rmeddis@29 92 method.plotGraphs= 0; % plot only PSTHs
rmeddis@29 93 else
rmeddis@29 94 method.plotGraphs= plotGraphsForIHC; % show progress
rmeddis@29 95 end
rmeddis@29 96
rmeddis@29 97 targetChannelNo=1;
rmeddis@29 98
rmeddis@29 99 % force parameters
rmeddis@29 100 % the number of AN fibers at each BF
rmeddis@29 101 AN_IHCsynapseParams.numFibers= nRepeats;
rmeddis@29 102 AN_IHCsynapseParams. mode= 'spikes';
rmeddis@29 103 AN_IHCsynapseParams.plotSynapseContents=0;
rmeddis@29 104 AN_IHCsynapseParams.PSTHbinWidth=.001;
rmeddis@29 105
rmeddis@29 106 method.DRNLSave=1;
rmeddis@29 107 method.IHC_cilia_RPSave=1;
rmeddis@29 108 method.PSTHbinWidth=1e-3; % useful 1-ms default for all PSTHs
rmeddis@29 109 method.AN_IHCsynapseSave=1;
rmeddis@29 110 method.MacGregorMultiSave=1;
rmeddis@29 111 method.MacGregorSave=1;
rmeddis@29 112 method.dt=dt;
rmeddis@29 113
rmeddis@29 114 moduleSequence=[1:8];
rmeddis@29 115
rmeddis@38 116 global dtSpikes ARAttenuation TMoutput OMEoutput ...
rmeddis@29 117 DRNLoutput IHC_cilia_output IHCrestingCiliaCond IHCrestingV...
rmeddis@32 118 IHCoutput ANprobRateOutput ANoutput savePavailable ANtauCas ...
rmeddis@29 119 CNoutput ICoutput ICmembraneOutput ICfiberTypeRates MOCattenuation
rmeddis@29 120
rmeddis@29 121 AN_spikesOrProbability='spikes';
rmeddis@29 122 AN_spikesOrProbability='probability';
rmeddis@29 123 MAPparamsName='Normal';
rmeddis@29 124
rmeddis@29 125 MAP1_14(inputSignal, 1/dt, targetFrequency, ...
rmeddis@29 126 MAPparamsName, AN_spikesOrProbability);
rmeddis@29 127
rmeddis@29 128 % RP
rmeddis@29 129 IHC_RPData=IHC_cilia_output;
rmeddis@29 130 IHC_RPData=IHCoutput(targetChannelNo,:);
rmeddis@29 131 IHC_RP_peak(levelNo,freqNo)=max(IHC_RPData(toneStartptr:toneEndptr));
rmeddis@29 132 IHC_RP_min(levelNo,freqNo)=min(IHC_RPData(toneStartptr:toneEndptr));
rmeddis@29 133 IHC_RP_dc(levelNo,freqNo)=mean(IHC_RPData(toneStartptr:toneEndptr));
rmeddis@29 134
rmeddis@29 135 % AN next
rmeddis@29 136 AN_IHCsynapseAllData=ANoutput;
rmeddis@29 137 method.PSTHbinWidth=0.001;
rmeddis@29 138
rmeddis@32 139 nTaus=length(ANtauCas);
rmeddis@29 140 numANfibers=size(ANoutput,1);
rmeddis@29 141 numLSRfibers=numANfibers/nTaus;
rmeddis@29 142
rmeddis@29 143 %LSR (same as HSR if no LSR fibers present)
rmeddis@29 144 channelPtr1=(targetChannelNo-1)*numANfibers+1;
rmeddis@29 145 channelPtr2=channelPtr1+numANfibers-1;
rmeddis@29 146 LSRspikes=AN_IHCsynapseAllData(channelPtr1:channelPtr2,:);
rmeddis@29 147 method.dt=method.AN_IHCsynapsedt;
rmeddis@29 148 PSTH=UTIL_PSTHmaker(LSRspikes, method);
rmeddis@29 149 PSTH=sum(PSTH,1); % sum across fibers (HSR only)
rmeddis@29 150 PSTHStartptr=round(silenceDuration/method.PSTHbinWidth)+1;
rmeddis@29 151 PSTHMidptr=PSTHStartptr+round(toneDuration/(2*method.PSTHbinWidth)) -1;
rmeddis@29 152 PSTHEndptr=PSTHStartptr+round(toneDuration/method.PSTHbinWidth) -1;
rmeddis@29 153 AN_LSRonset(levelNo,freqNo)=max(max(PSTH))/(method.PSTHbinWidth*method.numANfibers);
rmeddis@29 154 AN_LSRsaturated(levelNo,freqNo)=sum(PSTH(PSTHMidptr:PSTHEndptr))/(method.numANfibers*toneDuration/2);
rmeddis@29 155
rmeddis@29 156 % HSR
rmeddis@29 157 channelPtr1=numLSRfibers+(targetChannelNo-1)*method.numANfibers+1;
rmeddis@29 158 channelPtr2=channelPtr1+method.numANfibers-1;
rmeddis@29 159 HSRspikes=AN_IHCsynapseAllData(channelPtr1:channelPtr2,:);
rmeddis@29 160 method.dt=method.AN_IHCsynapsedt;
rmeddis@29 161 PSTH=UTIL_PSTHmaker(HSRspikes, method);
rmeddis@29 162 PSTH=sum(PSTH,1); % sum across fibers (HSR only)
rmeddis@29 163 PSTHStartptr=round(silenceDuration/method.PSTHbinWidth)+1;
rmeddis@29 164 PSTHMidptr=PSTHStartptr+round(toneDuration/(2*method.PSTHbinWidth)) -1;
rmeddis@29 165 PSTHEndptr=PSTHStartptr+round(toneDuration/method.PSTHbinWidth) -1;
rmeddis@29 166 AN_HSRonset(levelNo,freqNo)=max(max(PSTH))/(method.PSTHbinWidth*method.numANfibers);
rmeddis@29 167 AN_HSRsaturated(levelNo,freqNo)=sum(PSTH(PSTHMidptr:PSTHEndptr))/(method.numANfibers*toneDuration/2);
rmeddis@29 168 [cvANHSR, cvTimes, allTimeStamps, allISIs]= UTIL_CV(HSRspikes, method.AN_IHCsynapsedt);
rmeddis@29 169
rmeddis@29 170 PSTHplotCount=PSTHplotCount+1;
rmeddis@29 171 if doReceptiveFields % receptive field for HSR only
rmeddis@29 172 figure(14), set(gcf,'name','AN')
rmeddis@29 173 plotReceptiveFields(method, PSTH, PSTHplotCount, levels, frequencies)
rmeddis@29 174 ylim([0 method.numANfibers])
rmeddis@29 175 xlabel(['CV= ' num2str(max(cvANHSR),'%4.2f')],'fontsize',8)
rmeddis@29 176 end % doReceptiveFields
rmeddis@29 177
rmeddis@29 178 % CN
rmeddis@29 179 MacGregorMultiAllData=method.MacGregorMultiData;
rmeddis@29 180 numLSRfibers=method.McGMultinNeuronsPerBF*length(method.nonlinCF)* (nTaus-1);
rmeddis@29 181
rmeddis@29 182 %LSR (same as HSR if no LSR fibers present)
rmeddis@29 183 channelPtr1=(targetChannelNo-1)*method.McGMultinNeuronsPerBF+1;
rmeddis@29 184 channelPtr2=channelPtr1+method.McGMultinNeuronsPerBF-1;
rmeddis@29 185 MacGregorMultiLSRspikes=MacGregorMultiAllData(channelPtr1:channelPtr2,:);
rmeddis@29 186 method.dt=method.MacGregorMultidt;
rmeddis@29 187 PSTH=UTIL_PSTHmaker(MacGregorMultiLSRspikes, method);
rmeddis@29 188 PSTH=sum(PSTH,1); % sum across fibers (HSR only)
rmeddis@29 189 PSTHStartptr=round(silenceDuration/method.PSTHbinWidth)+1;
rmeddis@29 190 PSTHMidptr=PSTHStartptr+round(toneDuration/(2*method.PSTHbinWidth)) -1;
rmeddis@29 191 PSTHEndptr=PSTHStartptr+round(toneDuration/method.PSTHbinWidth) -1;
rmeddis@29 192 CNLSRsaturated(levelNo,freqNo)=sum(PSTH(PSTHMidptr:PSTHEndptr));
rmeddis@29 193 CNLSRsaturated(levelNo,freqNo)=CNLSRsaturated(levelNo,freqNo)...
rmeddis@29 194 /((toneDuration/2)*method.McGMultinNeuronsPerBF);
rmeddis@29 195
rmeddis@29 196 %HSR
rmeddis@29 197 channelPtr1=numLSRfibers+(targetChannelNo-1)*method.McGMultinNeuronsPerBF+1;
rmeddis@29 198 channelPtr2=channelPtr1+method.McGMultinNeuronsPerBF-1;
rmeddis@29 199 MacGregorMultiHSRspikes=MacGregorMultiAllData(channelPtr1:channelPtr2,:);
rmeddis@29 200 method.dt=method.MacGregorMultidt;
rmeddis@29 201 PSTH=UTIL_PSTHmaker(MacGregorMultiHSRspikes, method);
rmeddis@29 202 PSTH=sum(PSTH,1); % sum across fibers (HSR only)
rmeddis@29 203 PSTHStartptr=round(silenceDuration/method.PSTHbinWidth)+1;
rmeddis@29 204 PSTHMidptr=PSTHStartptr+round(toneDuration/(2*method.PSTHbinWidth)) -1;
rmeddis@29 205 PSTHEndptr=PSTHStartptr+round(toneDuration/method.PSTHbinWidth) -1;
rmeddis@29 206 CNHSRsaturated(levelNo,freqNo)=sum(PSTH(PSTHMidptr:PSTHEndptr));
rmeddis@29 207 CNHSRsaturated(levelNo,freqNo)=CNHSRsaturated(levelNo,freqNo)...
rmeddis@29 208 /((toneDuration/2)*method.McGMultinNeuronsPerBF);
rmeddis@29 209 [cvMMHSR, cvTimes, allTimeStamps, allISIs]= UTIL_CV(MacGregorMultiHSRspikes, method.MacGregorMultidt);
rmeddis@29 210
rmeddis@29 211 if doReceptiveFields % receptive field
rmeddis@29 212 figure(15), set(gcf,'name','CN HSR input')
rmeddis@29 213 plotReceptiveFields(method, PSTH, PSTHplotCount, levels, frequencies)
rmeddis@29 214 ylim([0 method.McGMultinNeuronsPerBF])
rmeddis@29 215 xlabel(['CV= ' num2str(max(cvMMHSR),'%4.2f')],'fontsize',8)
rmeddis@29 216 end
rmeddis@29 217
rmeddis@29 218 MacGregorAllData=method.MacGregorData;
rmeddis@29 219 numLSRfibers=length(method.nonlinCF)* (nTaus-1);
rmeddis@29 220
rmeddis@29 221 %LSR (same as HSR if no LSR fibers present)
rmeddis@29 222 channelPtr1=targetChannelNo;
rmeddis@29 223 MacGregorLSR=MacGregorAllData(channelPtr1,:);
rmeddis@29 224 method.dt=method.MacGregordt;
rmeddis@29 225 PSTH=UTIL_PSTHmaker(MacGregorLSR, method);
rmeddis@29 226 PSTHStartptr=round(silenceDuration/method.PSTHbinWidth)+1;
rmeddis@29 227 PSTHMidptr=PSTHStartptr+round(toneDuration/(2*method.PSTHbinWidth)) -1;
rmeddis@29 228 PSTHEndptr=PSTHStartptr+round(toneDuration/method.PSTHbinWidth) -1;
rmeddis@29 229 ICLSRsaturated(levelNo,freqNo)=sum(PSTH(PSTHMidptr:PSTHEndptr));
rmeddis@29 230 ICLSRsaturated(levelNo,freqNo)=ICLSRsaturated(levelNo,freqNo)/(toneDuration/2);
rmeddis@29 231
rmeddis@29 232 %LSR (same as HSR if no LSR fibers present)
rmeddis@29 233 channelPtr1=numLSRfibers+targetChannelNo;
rmeddis@29 234 MacGregorHSR=MacGregorAllData(channelPtr1,:);
rmeddis@29 235 method.dt=method.MacGregordt;
rmeddis@29 236 PSTH=UTIL_PSTHmaker(MacGregorHSR, method);
rmeddis@29 237 PSTHStartptr=round(silenceDuration/method.PSTHbinWidth)+1;
rmeddis@29 238 PSTHMidptr=PSTHStartptr+round(toneDuration/(2*method.PSTHbinWidth)) -1;
rmeddis@29 239 PSTHEndptr=PSTHStartptr+round(toneDuration/method.PSTHbinWidth) -1;
rmeddis@29 240 ICHSRsaturated(levelNo,freqNo)=sum(PSTH(PSTHMidptr:PSTHEndptr));
rmeddis@29 241 ICHSRsaturated(levelNo,freqNo)=ICHSRsaturated(levelNo,freqNo)/(toneDuration/2);
rmeddis@29 242 [cvICHSR, cvTimes, allTimeStamps, allISIs]= UTIL_CV(MacGregorHSR, method.MacGregordt);
rmeddis@29 243
rmeddis@29 244 % if doReceptiveFields % receptive field
rmeddis@29 245 % figure(16), set(gcf,'name','IC HSR input')
rmeddis@29 246 % plotReceptiveFields(method, PSTH, PSTHplotCount, levels, frequencies)
rmeddis@29 247 % ylim([0 method.McGMultinNeuronsPerBF])
rmeddis@29 248 % xlabel(['CV= ' num2str(max(cvICHSR),'%4.2f')],'fontsize',8)
rmeddis@29 249 % end
rmeddis@29 250 end % frequency
rmeddis@29 251 end % level
rmeddis@29 252 fprintf('\n')
rmeddis@29 253 toneDuration=2;
rmeddis@29 254 rampDuration=0.004;
rmeddis@29 255 silenceDuration=.02;
rmeddis@29 256 nRepeats=200; % no. of AN fibers
rmeddis@29 257 fprintf('toneDuration %6.3f\n', toneDuration)
rmeddis@29 258 fprintf(' %6.0f AN fibers (repeats)\n', nRepeats)
rmeddis@29 259 fprintf('levels')
rmeddis@29 260 fprintf('%6.2f\t', levels)
rmeddis@29 261 fprintf('\n')
rmeddis@29 262
rmeddis@29 263
rmeddis@29 264 % ---------------------------------------------------------- display parameters
rmeddis@29 265 disp(['parameter file was: ' experiment.name])
rmeddis@29 266 fprintf('\n')
rmeddis@29 267 UTIL_showStruct(IHC_VResp_VivoParams, 'IHC_cilia_RPParams')
rmeddis@29 268 UTIL_showStruct(IHCpreSynapseParams, 'IHCpreSynapseParams')
rmeddis@29 269 UTIL_showStruct(AN_IHCsynapseParams, 'AN_IHCsynapseParams')
rmeddis@29 270
rmeddis@29 271
rmeddis@29 272
rmeddis@29 273 function plotReceptiveFields(method, PSTH, PSTHplotCount, levels, frequencies)
rmeddis@29 274
rmeddis@29 275 % show PSTH for each level/frequency combination
rmeddis@29 276 nLevels=length(levels);
rmeddis@29 277 nFrequencies=length(frequencies);
rmeddis@29 278
rmeddis@29 279 PSTHtime=method.PSTHbinWidth:method.PSTHbinWidth:method.PSTHbinWidth*length(PSTH);
rmeddis@29 280 subplot(nLevels,nFrequencies,PSTHplotCount)
rmeddis@29 281 bar(PSTHtime, PSTH)
rmeddis@29 282 xlim([0 max(PSTHtime)])
rmeddis@29 283 % write axis labels only at left and bottom
rmeddis@29 284 if PSTHplotCount< (nLevels-1) * nFrequencies+1
rmeddis@29 285 set(gca,'xticklabel',[])
rmeddis@29 286 end
rmeddis@29 287 if ~isequal(mod(PSTHplotCount,nFrequencies),1)
rmeddis@29 288 set(gca,'yticklabel',[])
rmeddis@29 289 else
rmeddis@29 290 ylabel([num2str(levels(round(PSTHplotCount/nFrequencies) +1)) ' dB'])
rmeddis@29 291 end
rmeddis@29 292 % add titles only on top row
rmeddis@29 293 if PSTHplotCount<=nFrequencies
rmeddis@29 294 title([num2str(frequencies(PSTHplotCount)) ' Hz'])
rmeddis@29 295 end