annotate testPrograms/temp.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 3ea506487b3b
children
rev   line source
rmeddis@38 1 function vectorStrength=testAN(probeFrequency,BFlist, levels, ...
rmeddis@38 2 paramsName,paramChanges)
rmeddis@38 3 % testAN generates rate/level functions for AN and brainstem units.
rmeddis@38 4 % also other information like PSTHs, MOC efferent activity levels.
rmeddis@38 5 % Vector strength calculations require the computation of period
rmeddis@38 6 % histograms. for this reason the sample rate must always be an integer
rmeddis@38 7 % multiple of the tone frequency. This applies to both the sampleRate and
rmeddis@38 8 % the spikes sampling rate.
rmeddis@38 9 % A full 'spikes' model is used.
rmeddis@38 10 % paramChanges is a cell array of strings containing MATLAB statements that
rmeddis@38 11 % change model parameters. See MAPparamsNormal for examples.
rmeddis@38 12 % e.g.
rmeddis@38 13 % testAN(1000,1000, -10:10:80,'Normal',{});
rmeddis@35 14
rmeddis@38 15
rmeddis@38 16 global IHC_VResp_VivoParams IHC_cilia_RPParams IHCpreSynapseParams
rmeddis@38 17 global AN_IHCsynapseParams
rmeddis@38 18 global ANoutput dtSpikes CNoutput ICoutput ICmembraneOutput ANtauCas
rmeddis@38 19 global ARattenuation MOCattenuation
rmeddis@38 20 tic
rmeddis@36 21 dbstop if error
rmeddis@36 22 restorePath=path;
rmeddis@38 23 addpath (['..' filesep 'MAP'], ['..' filesep 'utilities'], ...
rmeddis@38 24 ['..' filesep 'parameterStore'], ['..' filesep 'wavFileStore'],...
rmeddis@38 25 ['..' filesep 'testPrograms'])
rmeddis@35 26
rmeddis@38 27 if nargin<5, paramChanges={'IHCpreSynapseParams.tauCa= [30e-6 90e-6];'}; end
rmeddis@38 28 if nargin<4, paramsName='PL'; end
rmeddis@38 29 if nargin<3, levels=-10:10:80; end
rmeddis@38 30 if nargin==0, probeFrequency=1000; BFlist=1000; end
rmeddis@38 31 nLevels=length(levels);
rmeddis@36 32
rmeddis@38 33 toneDuration=.2; rampDuration=0.005; silenceDuration=.02;
rmeddis@38 34 localPSTHbinwidth=0.001;
rmeddis@36 35
rmeddis@38 36 %% guarantee that the sample rate is an interger multiple
rmeddis@38 37 % of the probe frequency
rmeddis@38 38 % we want 5 bins per period for spikes
rmeddis@38 39 spikesSampleRate=5*probeFrequency;
rmeddis@38 40 % model sample rate must be an integer multiple of this and in the region
rmeddis@38 41 % of 50000
rmeddis@38 42 sampleRate=spikesSampleRate*round(50000/spikesSampleRate);
rmeddis@38 43 dt=1/sampleRate;
rmeddis@38 44 % avoid very slow spikes sampling rate
rmeddis@38 45 spikesSampleRate=spikesSampleRate*ceil(10000/spikesSampleRate);
rmeddis@36 46
rmeddis@38 47 %% pre-allocate storage
rmeddis@38 48 AN_HSRonset=zeros(nLevels,1);
rmeddis@38 49 AN_HSRsaturated=zeros(nLevels,1);
rmeddis@38 50 AN_LSRonset=zeros(nLevels,1);
rmeddis@38 51 AN_LSRsaturated=zeros(nLevels,1);
rmeddis@38 52 CNLSRrate=zeros(nLevels,1);
rmeddis@38 53 CNHSRsaturated=zeros(nLevels,1);
rmeddis@38 54 ICHSRsaturated=zeros(nLevels,1);
rmeddis@38 55 ICLSRsaturated=zeros(nLevels,1);
rmeddis@38 56 vectorStrength=zeros(nLevels,1);
rmeddis@36 57
rmeddis@38 58 AR=zeros(nLevels,1);
rmeddis@38 59 MOC=zeros(nLevels,1);
rmeddis@36 60
rmeddis@38 61 figure(15), clf
rmeddis@38 62 set(gcf,'position',[980 356 401 321])
rmeddis@38 63 figure(5), clf
rmeddis@38 64 set(gcf,'position', [980 34 400 295])
rmeddis@38 65 drawnow
rmeddis@36 66
rmeddis@36 67
rmeddis@38 68 %% main computational loop (vary level)
rmeddis@38 69 levelNo=0;
rmeddis@38 70 for leveldB=levels
rmeddis@38 71 levelNo=levelNo+1;
rmeddis@38 72 amp=28e-6*10^(leveldB/20);
rmeddis@38 73 fprintf('%4.0f\t', leveldB)
rmeddis@36 74
rmeddis@38 75 %% generate tone and silences
rmeddis@38 76 time=dt:dt:toneDuration;
rmeddis@38 77 rampTime=dt:dt:rampDuration;
rmeddis@38 78 ramp=[0.5*(1+cos(2*pi*rampTime/(2*rampDuration)+pi)) ...
rmeddis@38 79 ones(1,length(time)-length(rampTime))];
rmeddis@38 80 ramp=ramp.*fliplr(ramp);
rmeddis@38 81
rmeddis@38 82 silence=zeros(1,round(silenceDuration/dt));
rmeddis@38 83
rmeddis@38 84 inputSignal=amp*sin(2*pi*probeFrequency'*time);
rmeddis@38 85 inputSignal= ramp.*inputSignal;
rmeddis@38 86 inputSignal=[silence inputSignal];
rmeddis@38 87
rmeddis@38 88 %% run the model
rmeddis@38 89 AN_spikesOrProbability='spikes';
rmeddis@38 90 nExistingParamChanges=length(paramChanges);
rmeddis@38 91 paramChanges{nExistingParamChanges+1}=...
rmeddis@38 92 ['AN_IHCsynapseParams.spikesTargetSampleRate=' ...
rmeddis@38 93 num2str(spikesSampleRate) ';'];
rmeddis@36 94
rmeddis@38 95 MAP1_14(inputSignal, 1/dt, BFlist, ...
rmeddis@38 96 paramsName, AN_spikesOrProbability, paramChanges);
rmeddis@38 97
rmeddis@38 98 nTaus=length(ANtauCas);
rmeddis@38 99
rmeddis@38 100 %% Auditory nerve evaluate and display (Fig. 5)
rmeddis@38 101 %LSR (same as HSR if no LSR fibers present)
rmeddis@38 102 [nANFibers nTimePoints]=size(ANoutput);
rmeddis@38 103 numLSRfibers=nANFibers/nTaus;
rmeddis@38 104 numHSRfibers=numLSRfibers;
rmeddis@38 105
rmeddis@38 106 LSRspikes=ANoutput(1:numLSRfibers,:);
rmeddis@38 107 PSTH=UTIL_PSTHmaker(LSRspikes, dtSpikes, localPSTHbinwidth);
rmeddis@38 108 PSTHLSR=mean(PSTH,1)/localPSTHbinwidth; % across fibers rates
rmeddis@38 109 PSTHtime=localPSTHbinwidth:localPSTHbinwidth:...
rmeddis@38 110 localPSTHbinwidth*length(PSTH);
rmeddis@38 111 AN_LSRonset(levelNo)= max(PSTHLSR); % peak in 5 ms window
rmeddis@38 112 AN_LSRsaturated(levelNo)= mean(PSTHLSR(round(length(PSTH)/2):end));
rmeddis@38 113
rmeddis@38 114 % AN HSR
rmeddis@38 115 HSRspikes= ANoutput(end- numHSRfibers+1:end, :);
rmeddis@38 116 PSTH=UTIL_PSTHmaker(HSRspikes, dtSpikes, localPSTHbinwidth);
rmeddis@38 117 PSTH=mean(PSTH,1)/localPSTHbinwidth; % sum across fibers (HSR only)
rmeddis@38 118 AN_HSRonset(levelNo)= max(PSTH);
rmeddis@38 119 AN_HSRsaturated(levelNo)= mean(PSTH(round(length(PSTH)/2): end));
rmeddis@38 120
rmeddis@38 121 figure(5), subplot(2,2,2)
rmeddis@38 122 hold off, bar(PSTHtime,PSTH, 'b')
rmeddis@38 123 hold on, bar(PSTHtime,PSTHLSR,'r')
rmeddis@38 124 ylim([0 1000])
rmeddis@38 125 xlim([0 length(PSTH)*localPSTHbinwidth])
rmeddis@38 126 grid on
rmeddis@38 127 set(gcf,'name',['PSTH: ' num2str(BFlist), ' Hz: ' num2str(leveldB) ' dB']);
rmeddis@38 128
rmeddis@38 129 % AN - CV
rmeddis@38 130 % CV is computed 5 times. Use the middle one (3) as most typical
rmeddis@38 131 cvANHSR= UTIL_CV(HSRspikes, dtSpikes);
rmeddis@38 132
rmeddis@38 133 % AN - vector strength
rmeddis@38 134 PSTH=sum(CNoutput (11:20,:));
rmeddis@38 135 [PH, binTimes]=UTIL_periodHistogram...
rmeddis@38 136 (PSTH, dtSpikes, probeFrequency);
rmeddis@38 137 VS=UTIL_vectorStrength(PH);
rmeddis@38 138 vectorStrength(levelNo)=VS;
rmeddis@38 139 disp(['sat rate= ' num2str(AN_HSRsaturated(levelNo)) ...
rmeddis@38 140 '; phase-locking VS = ' num2str(VS)])
rmeddis@38 141 title(['AN HSR: CV=' num2str(cvANHSR(3),'%5.2f') ...
rmeddis@38 142 'VS=' num2str(VS,'%5.2f')])
rmeddis@38 143
rmeddis@38 144 %% CN - first-order neurons
rmeddis@38 145
rmeddis@38 146 % CN driven by LSR fibers
rmeddis@38 147 [nCNneurons c]=size(CNoutput);
rmeddis@38 148 nLSRneurons=round(nCNneurons/nTaus);
rmeddis@38 149 CNLSRspikes=CNoutput(1:nLSRneurons,:);
rmeddis@38 150 PSTH=UTIL_PSTHmaker(CNLSRspikes, dtSpikes, localPSTHbinwidth);
rmeddis@38 151 PSTH=sum(PSTH)/nLSRneurons;
rmeddis@38 152 % CNLSRrate(levelNo)=mean(PSTH(round(length(PSTH)/2):end))/localPSTHbinwidth;
rmeddis@38 153 CNLSRrate(levelNo)=mean(PSTH)/localPSTHbinwidth;
rmeddis@38 154
rmeddis@38 155 %CN HSR
rmeddis@38 156 MacGregorMultiHSRspikes=...
rmeddis@38 157 CNoutput(end-nLSRneurons+1:end,:);
rmeddis@38 158 PSTH=UTIL_PSTHmaker(MacGregorMultiHSRspikes, dtSpikes, localPSTHbinwidth);
rmeddis@38 159 PSTH=sum(PSTH)/nLSRneurons;
rmeddis@38 160 PSTH=mean(PSTH,1)/localPSTHbinwidth; % sum across fibers (HSR only)
rmeddis@38 161
rmeddis@38 162 % CNHSRsaturated(levelNo)=mean(PSTH(length(PSTH)/2:end));
rmeddis@38 163 CNHSRsaturated(levelNo)=mean(PSTH);
rmeddis@38 164
rmeddis@38 165 figure(5), subplot(2,2,3)
rmeddis@38 166 bar(PSTHtime,PSTH)
rmeddis@38 167 ylim([0 1000])
rmeddis@38 168 xlim([0 length(PSTH)*localPSTHbinwidth])
rmeddis@38 169 cvMMHSR= UTIL_CV(MacGregorMultiHSRspikes, dtSpikes);
rmeddis@38 170 title(['CN CV= ' num2str(cvMMHSR(3),'%5.2f')])
rmeddis@38 171
rmeddis@38 172 %% IC LSR
rmeddis@38 173 [nICneurons c]=size(ICoutput);
rmeddis@38 174 nLSRneurons=round(nICneurons/nTaus);
rmeddis@38 175 ICLSRspikes=ICoutput(1:nLSRneurons,:);
rmeddis@38 176 PSTH=UTIL_PSTHmaker(ICLSRspikes, dtSpikes, localPSTHbinwidth);
rmeddis@38 177 % ICLSRsaturated(levelNo)=mean(PSTH(round(length(PSTH)/2):end))/localPSTHbinwidth;
rmeddis@38 178 ICLSRsaturated(levelNo)=mean(PSTH)/localPSTHbinwidth;
rmeddis@38 179
rmeddis@38 180 %IC HSR
rmeddis@38 181 MacGregorMultiHSRspikes=...
rmeddis@38 182 ICoutput(end-nLSRneurons+1:end,:);
rmeddis@38 183 PSTH=UTIL_PSTHmaker(MacGregorMultiHSRspikes, dtSpikes, localPSTHbinwidth);
rmeddis@38 184 ICHSRsaturated(levelNo)= (sum(PSTH)/nLSRneurons)/toneDuration;
rmeddis@36 185
rmeddis@38 186 AR(levelNo)=min(ARattenuation);
rmeddis@38 187 MOC(levelNo)=min(MOCattenuation(length(MOCattenuation)/2:end));
rmeddis@38 188
rmeddis@38 189 time=dt:dt:dt*size(ICmembraneOutput,2);
rmeddis@38 190 figure(5), subplot(2,2,4)
rmeddis@38 191 % plot HSR (last of two)
rmeddis@38 192 plot(time,ICmembraneOutput(end-nLSRneurons+1, 1:end),'k')
rmeddis@38 193 ylim([-0.07 0])
rmeddis@38 194 xlim([0 max(time)])
rmeddis@38 195 title(['IC ' num2str(leveldB,'%4.0f') 'dB'])
rmeddis@38 196 drawnow
rmeddis@38 197
rmeddis@38 198 figure(5), subplot(2,2,1)
rmeddis@38 199 plot(20*log10(MOC), 'k'), hold on
rmeddis@38 200 plot(20*log10(AR), 'r'), hold off
rmeddis@38 201 title(' MOC'), ylabel('dB attenuation')
rmeddis@38 202 ylim([-30 0])
rmeddis@38 203
rmeddis@38 204 % x=input('prompt')
rmeddis@38 205 end % level
rmeddis@36 206
rmeddis@38 207 %% plot with levels on x-axis
rmeddis@38 208 figure(5), subplot(2,2,1)
rmeddis@38 209 plot(levels,20*log10(MOC), 'k'),hold on
rmeddis@38 210 plot(levels,20*log10(AR), 'r'), hold off
rmeddis@38 211 title(' MOC'), ylabel('dB attenuation')
rmeddis@38 212 ylim([-30 0])
rmeddis@38 213 xlim([0 max(levels)])
rmeddis@36 214
rmeddis@38 215 fprintf('\n')
rmeddis@38 216 toneDuration=2;
rmeddis@38 217 rampDuration=0.004;
rmeddis@38 218 silenceDuration=.02;
rmeddis@38 219 nRepeats=200; % no. of AN fibers
rmeddis@38 220 fprintf('toneDuration %6.3f\n', toneDuration)
rmeddis@38 221 fprintf(' %6.0f AN fibers (repeats)\n', nRepeats)
rmeddis@38 222 fprintf('levels')
rmeddis@38 223 fprintf('%6.2f\t', levels)
rmeddis@38 224 fprintf('\n')
rmeddis@36 225
rmeddis@38 226 %% ---------------------- display summary results (Fig 15)
rmeddis@38 227 figure(15), clf
rmeddis@38 228 nRows=2; nCols=2;
rmeddis@36 229
rmeddis@38 230 % AN rate - level ONSET functions
rmeddis@38 231 subplot(nRows,nCols,1)
rmeddis@38 232 plot(levels,AN_LSRonset,'ro'), hold on
rmeddis@38 233 plot(levels,AN_HSRonset,'ko'), hold off
rmeddis@38 234 ylim([0 1000])
rmeddis@38 235 if length(levels)>1
rmeddis@38 236 xlim([min(levels) max(levels)])
rmeddis@36 237 end
rmeddis@36 238
rmeddis@38 239 ttl=['tauCa= ' num2str(IHCpreSynapseParams.tauCa)];
rmeddis@38 240 title( ttl)
rmeddis@38 241 xlabel('level dB SPL'), ylabel('peak rate (sp/s)'), grid on
rmeddis@38 242 text(0, 800, 'AN onset', 'fontsize', 14)
rmeddis@38 243
rmeddis@38 244 % AN rate - level ADAPTED function
rmeddis@38 245 subplot(nRows,nCols,2)
rmeddis@38 246 plot(levels,AN_LSRsaturated, 'ro'), hold on
rmeddis@38 247 plot(levels,AN_HSRsaturated, 'ko'), hold off
rmeddis@38 248 ylim([0 400])
rmeddis@38 249 set(gca,'ytick',0:50:300)
rmeddis@38 250 if length(levels)>1
rmeddis@38 251 xlim([min(levels) max(levels)])
rmeddis@36 252 end
rmeddis@38 253 set(gca,'xtick',[levels(1):20:levels(end)])
rmeddis@38 254 % grid on
rmeddis@38 255 ttl=[ 'spont=' num2str(mean(AN_HSRsaturated(1,:)),'%4.0f')...
rmeddis@38 256 ' sat=' num2str(mean(AN_HSRsaturated(end,1)),'%4.0f')];
rmeddis@38 257 title( ttl)
rmeddis@38 258 xlabel('level dB SPL'), ylabel ('adapted rate (sp/s)')
rmeddis@38 259 text(0, 340, 'AN adapted', 'fontsize', 14), grid on
rmeddis@36 260
rmeddis@38 261 % CN rate - level function
rmeddis@38 262 subplot(nRows,nCols,3)
rmeddis@38 263 plot(levels,CNLSRrate, 'ro'), hold on
rmeddis@38 264 plot(levels,CNHSRsaturated, 'ko'), hold off
rmeddis@38 265 ylim([0 400])
rmeddis@38 266 set(gca,'ytick',0:50:300)
rmeddis@38 267 if length(levels)>1
rmeddis@38 268 xlim([min(levels) max(levels)])
rmeddis@38 269 end
rmeddis@38 270 set(gca,'xtick',[levels(1):20:levels(end)])
rmeddis@38 271 % grid on
rmeddis@38 272 ttl=[ 'spont=' num2str(mean(CNHSRsaturated(1,:)),'%4.0f') ' sat=' ...
rmeddis@38 273 num2str(mean(CNHSRsaturated(end,1)),'%4.0f')];
rmeddis@38 274 title( ttl)
rmeddis@38 275 xlabel('level dB SPL'), ylabel ('adapted rate (sp/s)')
rmeddis@38 276 text(0, 350, 'CN', 'fontsize', 14), grid on
rmeddis@36 277
rmeddis@38 278 % IC rate - level function
rmeddis@38 279 subplot(nRows,nCols,4)
rmeddis@38 280 plot(levels,ICLSRsaturated, 'ro'), hold on
rmeddis@38 281 plot(levels,ICHSRsaturated, 'ko'), hold off
rmeddis@38 282 ylim([0 400])
rmeddis@38 283 set(gca,'ytick',0:50:300)
rmeddis@38 284 if length(levels)>1
rmeddis@38 285 xlim([min(levels) max(levels)])
rmeddis@38 286 end
rmeddis@38 287 set(gca,'xtick',[levels(1):20:levels(end)]), grid on
rmeddis@38 288 ttl=['spont=' num2str(mean(ICHSRsaturated(1,:)),'%4.0f') ...
rmeddis@38 289 ' sat=' num2str(mean(ICHSRsaturated(end,1)),'%4.0f')];
rmeddis@38 290 title( ttl)
rmeddis@38 291 xlabel('level dB SPL'), ylabel ('adapted rate (sp/s)')
rmeddis@38 292 text(0, 350, 'IC', 'fontsize', 14)
rmeddis@38 293 set(gcf,'name',' AN CN IC rate/level')
rmeddis@36 294
rmeddis@38 295 fprintf('\n')
rmeddis@38 296 disp('levels vectorStrength')
rmeddis@38 297 fprintf('%3.0f \t %6.4f \n', [levels; vectorStrength'])
rmeddis@38 298 fprintf('\n')
rmeddis@38 299 fprintf('Phase locking, max vector strength=\t %6.4f\n\n',...
rmeddis@38 300 max(vectorStrength))
rmeddis@36 301
rmeddis@38 302 allData=[ levels' AN_HSRonset AN_HSRsaturated...
rmeddis@38 303 AN_LSRonset AN_LSRsaturated ...
rmeddis@38 304 CNHSRsaturated CNLSRrate...
rmeddis@38 305 ICHSRsaturated ICLSRsaturated];
rmeddis@38 306 fprintf('\n levels \tANHSR Onset \tANHSR adapted\tANLSR Onset \tANLSR adapted\tCNHSR\tCNLSR\tICHSR \tICLSR \n');
rmeddis@38 307 UTIL_printTabTable(round(allData))
rmeddis@38 308 fprintf('VS (phase locking)= \t%6.4f\n\n',...
rmeddis@38 309 max(vectorStrength))
rmeddis@36 310
rmeddis@38 311 UTIL_showStruct(IHC_cilia_RPParams, 'IHC_cilia_RPParams')
rmeddis@38 312 UTIL_showStruct(IHCpreSynapseParams, 'IHCpreSynapseParams')
rmeddis@38 313 UTIL_showStruct(AN_IHCsynapseParams, 'AN_IHCsynapseParams')
rmeddis@38 314
rmeddis@35 315 fprintf('\n')
rmeddis@38 316 disp('levels vectorStrength')
rmeddis@38 317 fprintf('%3.0f \t %6.4f \n', [levels; vectorStrength'])
rmeddis@38 318 fprintf('\n')
rmeddis@38 319 fprintf('Phase locking, max vector strength= \t%6.4f\n\n',...
rmeddis@38 320 max(vectorStrength))
rmeddis@35 321
rmeddis@38 322 allData=[ levels' AN_HSRonset AN_HSRsaturated...
rmeddis@38 323 AN_LSRonset AN_LSRsaturated ...
rmeddis@38 324 CNHSRsaturated CNLSRrate...
rmeddis@38 325 ICHSRsaturated ICLSRsaturated];
rmeddis@38 326 fprintf('\n levels \tANHSR Onset \tANHSR adapted\tANLSR Onset \tANLSR adapted\tCNHSR\tCNLSR\tICHSR \tICLSR \n');
rmeddis@38 327 UTIL_printTabTable(round(allData))
rmeddis@38 328 fprintf('VS (phase locking)= \t%6.4f\n\n',...
rmeddis@38 329 max(vectorStrength))
rmeddis@36 330
rmeddis@38 331 path(restorePath)
rmeddis@36 332 toc