Mercurial > hg > map

comparison userProgramsRM/testACF.m @ 38:c2204b18f4a2 tip

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
End nov big change
author Ray Meddis <rmeddis@essex.ac.uk>
date Mon, 28 Nov 2011 13:34:28 +0000
parents
children
comparison
equal deleted inserted replaced
37:771a643d5c29 38:c2204b18f4a2
1 % function [LP_SACF dt lags SACF]= testACF
2 % testACF is a *script* to demonstrate the smoothed ACF of
3 % Balaguer-Ballestera, E. Denham, S.L. and Meddis, R. (2008).
4 %
5 % Convert this to a *function* by uncommenting the first line
6 % The function returns the LP_SACF matrix plotted in Figure 96.
7 % If a function is used, the following outputs are returned:
8 % LP_SACF: smoothed SACF (lags x time matrix)
9 % dt: time interval between successive columns of LP_SACF
10 % lags: lags used in computing LP_SACF
11 % SACF: unsmoothed SACFs
12 %
13 % A range of options are supplied in the early part of the program
14 %
15 % #1
16 % Identify the model parameter file (in 'MAPparamsName')
17 %
18 % #2
19 % Identify the kind of model required (in 'AN_spikesOrProbability')
20 % 'probability' is recommended for ACF work
21 %
22 % #3
23 % Choose between a harmonic complex or file input
24 % by commenting out unwanted code
25 %
26 % #4
27 % Set the signal rms level (in leveldBSPL)
28 %
29 % #5
30 % Identify the model channel BFs in the vector 'BFlist'.
31 %
32 % #6
33 % Last minute changes to the model parameters can be made using
34 % the cell array of strings 'paramChanges'.
35 % This is used here to control the details of the ACF computations
36 % Read the notes in this section for more information
37 %
38 % displays:
39 % Figure 97 shows the AN response to the stimulus. this is a channel x time
40 % display. The z-axis (and colour) is the AN fiber firing rate
41 %
42 % Figure 96 shows the LP_SACF-matrix, the smoothed SACF.
43 %
44 % Figure 89 shows a summary of the evolution of the unsmoothed SACF
45 % over time. If you wish to take a snapshot of the LP_SACF-matrix at a
46 % particular time, this figure can help identify when to take it.
47 % The index on the y-axis, identifies the required row numbers
48 % of the LP_SACF or SACF matrix, e.g. LP_SACF(:,2000)
49 %
50 % On request, (filteredSACFParams.plotACFs=1) Figure 89 shows the channel
51 % by channel ACFs at intervals during the computation as a movie.
52 % The number of ACF displays is controlled by 'plotACFsInterval'
53 % and the movie can be slowed or speeded up using 'plotMoviePauses'
54 % (see paramChanges section below).
55
56 % - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
57 % This global will find results from MAP1_14
58 global savedInputSignal ANprobRateOutput ANoutput dt dtSpikes savedBFlist
59 % This global,from model parameter file
60 global filteredSACFParams
61
62 % User sets up requirements
63 %% #1 parameter file name
64 MAPparamsName='Normal'; % recommended
65
66
67 %% #2 probability (fast) or spikes (slow) representation: select one
68 % AN_spikesOrProbability='spikes';
69 % or
70 AN_spikesOrProbability='probability'; % recommended
71
72 %% #3 A. harmonic sequence or B. speech file input
73 % Comment out unwanted code
74 % A. harmonic tone (Hz) - useful to demonstrate a broadband sound
75 sampleRate= 44100; % recommended 44100
76 signalType= 'tones';
77 duration=0.100; % seconds
78 beginSilence=0.020;
79 endSilence=0.020;
80 rampDuration=.005; % raised cosine ramp (seconds)
81
82 % toneFrequency is a vector of component frequencies
83 F0=120;
84 toneFrequency= [3*F0 4*F0 5*F0];
85
86 % or
87 % B. file input
88 % signalType= 'file';
89 % fileName='Oh No';
90 % fileName='twister_44kHz';
91
92 %% #4 rms level
93 leveldBSPL= 100; % dB SPL (80 for Lieberman)
94
95 %% #5 number of channels in the model
96 % 21-channel model (log spacing of BFs)
97 numChannels=21;
98 lowestBF=250; highestBF= 5000;
99 BFlist=round(logspace(log10(lowestBF), log10(highestBF), numChannels));
100
101 %% #6 change model parameters
102 % Parameter changes can be used to change one or more model parameters
103 % *after* the MAPparams file has been read (see manual)
104
105 % Take control of ACF parameters
106 % The filteredACF parameters are set in the MAPparamsNormal file
107 % However, it is convenient to change them here leving the file intacta
108 minPitch= 400; maxPitch= 3000; numPitches=200;
109 maxLag=1/minPitch; minLag=1/maxPitch;
110 lags= linspace(minLag, maxLag, numPitches);
111
112 paramChanges={...
113 'filteredSACFParams.lags=lags; % autocorrelation lags vector;',...
114 'filteredSACFParams.acfTau= 2; % (Wiegrebe) time constant ACF;',...
115 'filteredSACFParams.lambda= 0.12; % slower filter to smooth ACF;',...
116 'filteredSACFParams.plotACFs=1; % plot ACFs while computing;',...
117 'filteredSACFParams.plotACFsInterval=0.01;',...
118 'filteredSACFParams.plotMoviePauses=.1; ',...
119 'filteredSACFParams.usePressnitzer=0; % attenuates ACF at long lags;',...
120 'filteredSACFParams.lagsProcedure= ''useAllLags'';',...
121 };
122
123 % Notes:
124 % acfTau: time constant of unsmoothed ACF
125 % lambda: time constant of smoothed ACFS
126 % plotACFs: plot ACFs during computation (0 to switch off, for speed)
127 % plotACFsInterval: sampling interval for plots
128 % plotMoviePauses: pause duration between frames to allow viewing
129 % usePressnitzer: gives low weights to long lags
130 % lagsProcedure: used to fiddle with output (ignore)
131
132 %% delare 'showMap' options to control graphical output
133 % see UTIL_showMAP for more options
134 showMapOptions=[];
135 % showMapOptions.showModelOutput=0; % plot of all stages
136 showMapOptions.surfAN=1; % surface plot of HSR response
137 showMapOptions.PSTHbinwidth=0.001; % smoothing for PSTH
138
139 if exist('fileName','var')
140 % needed for labeling plot
141 showMapOptions.fileName=fileName;
142 end
143
144 %% Generate stimuli
145 switch signalType
146 case 'tones'
147 % Create tone stimulus
148 dt=1/sampleRate; % seconds
149 time=dt: dt: duration;
150 inputSignal=sum(sin(2*pi*toneFrequency'*time), 1);
151 amp=10^(leveldBSPL/20)*28e-6; % converts to Pascals (peak)
152 inputSignal=amp*inputSignal;
153 % apply ramps
154 % catch rampTime error
155 if rampDuration>0.5*duration, rampDuration=duration/2; end
156 rampTime=dt:dt:rampDuration;
157 ramp=[0.5*(1+cos(2*pi*rampTime/(2*rampDuration)+pi)) ...
158 ones(1,length(time)-length(rampTime))];
159 inputSignal=inputSignal.*ramp;
160 ramp=fliplr(ramp);
161 inputSignal=inputSignal.*ramp;
162 % add silence
163 intialSilence= zeros(1,round(beginSilence/dt));
164 finalSilence= zeros(1,round(endSilence/dt));
165 inputSignal= [intialSilence inputSignal finalSilence];
166
167 case 'file'
168 %% file input simple or mixed
169 [inputSignal sampleRate]=wavread(fileName);
170 dt=1/sampleRate;
171 inputSignal=inputSignal(:,1);
172 targetRMS=20e-6*10^(leveldBSPL/20);
173 rms=(mean(inputSignal.^2))^0.5;
174 amp=targetRMS/rms;
175 inputSignal=inputSignal*amp;
176 end
177
178 wavplay(inputSignal, sampleRate)
179
180 %% run the model
181 dbstop if error
182 restorePath=path;
183 addpath (['..' filesep 'MAP'], ['..' filesep 'wavFileStore'], ...
184 ['..' filesep 'utilities'])
185
186 fprintf('\n')
187 disp(['Signal duration= ' num2str(length(inputSignal)/sampleRate)])
188 disp([num2str(numChannels) ' channel model: ' AN_spikesOrProbability])
189 disp('Computing MAP ...')
190
191 MAP1_14(inputSignal, sampleRate, BFlist, ...
192 MAPparamsName, AN_spikesOrProbability, paramChanges);
193
194
195 %% The model run is now complete. Now display the results
196 % display the AN response
197 UTIL_showMAP(showMapOptions)
198
199 % compute ACF
200 switch AN_spikesOrProbability
201 case 'probability'
202 % use only HSR fibers
203 inputToACF=ANprobRateOutput(end-length(savedBFlist)+1:end,:);
204 otherwise
205 inputToACF=ANoutput;
206 dt=dtSpikes;
207 end
208
209 disp ('computing ACF...')
210
211 % read paramChanges to get new filteredSACFParams
212 for i=1:length(paramChanges)
213 eval(paramChanges{i});
214 end
215
216 [LP_SACF BFlist SACF]= filteredSACF(inputToACF, dt, savedBFlist, ...
217 filteredSACFParams);
218 disp(' ACF done.')
219
220 %% plot original waveform on summary/smoothed ACF plot
221 figure(96), clf
222 subplot(3,1,3)
223 t=dt*(1:length(savedInputSignal));
224 plot(t,savedInputSignal, 'k')
225 xlim([0 t(end)])
226 title(['stimulus: ' num2str(leveldBSPL, '%4.0f') ' dB SPL']);
227
228 % plot SACF
229 figure(96)
230 subplot(2,1,1)
231 imagesc(LP_SACF)
232 colormap bone
233 ylabel('periodicities (Hz)'), xlabel('time (s)')
234 title(['smoothed SACF. (periodicity x time)'])
235 % y-axis specifies pitches (1/lags)
236 % Force MATLAB to show the lowest pitch
237 postedYvalues=[1 get(gca,'ytick')]; set(gca,'ytick',postedYvalues)
238 pitches=1./filteredSACFParams.lags;
239 set(gca,'ytickLabel', round(pitches(postedYvalues)))
240 % x-axis is time at which LP_SACF is samples
241 [nCH nTimes]=size(LP_SACF);
242 t=dt:dt:dt*nTimes;
243 tt=get(gca,'xtick');
244 set(gca,'xtickLabel', round(100*t(tt))/100)
245
246 %% On a new figure show a cascade of SACFs
247 figure(89), clf
248 % select 100 samples;
249 [r c]=size(SACF);
250 step=round(c/100);
251 idx=step:step:c;
252
253 UTIL_cascadePlot(SACF(:,idx)', 1./pitches)
254
255 xlabel('lag (s)'), ylabel('time pointer -->')
256 title(' SACF summary over time')
257 yValues=get(gca,'yTick');
258 set(gca,'yTickLabel', num2str(yValues'*100))
259
260 path(restorePath)
261