wolffd@0
|
1 function varargout = mirautocor(orig,varargin)
|
wolffd@0
|
2 % a = mirautocor(x) computes the autocorrelation function related to x.
|
wolffd@0
|
3 % Optional parameters:
|
wolffd@0
|
4 % mirautocor(...,'Min',mi) indicates the lowest delay taken into
|
wolffd@0
|
5 % consideration. The unit can be precised:
|
wolffd@0
|
6 % mirautocor(...,'Min',mi,'s') (default unit)
|
wolffd@0
|
7 % mirautocor(...,'Min',mi,'Hz')
|
wolffd@0
|
8 % Default value: 0 s.
|
wolffd@0
|
9 % mirautocor(...,'Max',ma) indicates the highest delay taken into
|
wolffd@0
|
10 % consideration. The unit can be specified as for 'Min'.
|
wolffd@0
|
11 % Default value:
|
wolffd@0
|
12 % if x is a signal, the highest delay is 0.05 s
|
wolffd@0
|
13 % (corresponding to a minimum frequency of 20 Hz).
|
wolffd@0
|
14 % if x is an envelope, the highest delay is 2 s.
|
wolffd@0
|
15 % mirautocor(...,'Resonance',r) multiplies the autocorrelation function
|
wolffd@0
|
16 % with a resonance curve:
|
wolffd@0
|
17 % Possible values:
|
wolffd@0
|
18 % 'Toiviainen' from (Toiviainen & Snyder, 2003)
|
wolffd@0
|
19 % 'vanNoorden' from (van Noorden & Moelants, 2001)
|
wolffd@0
|
20 % mirautocor(...,'Center',c) assigns the center value of the
|
wolffd@0
|
21 % resonance curve, in seconds.
|
wolffd@0
|
22 % Works mainly with 'Toiviainen' option.
|
wolffd@0
|
23 % Default value: c = 0.5
|
wolffd@0
|
24 % mirautocor(...,'Enhanced',a) reduces the effect of subharmonics.
|
wolffd@0
|
25 % The original autocorrelation function is half-wave rectified,
|
wolffd@0
|
26 % time-scaled by the factor a (which can be a factor list as
|
wolffd@0
|
27 % well), and substracted from the original clipped function.
|
wolffd@0
|
28 % (Tolonen & Karjalainen)
|
wolffd@0
|
29 % If the 'Enhanced' option is not followed by any value,
|
wolffd@0
|
30 % default value is a = 2:10
|
wolffd@0
|
31 % mirautocor(...,'Halfwave') performs a half-wave rectification on the
|
wolffd@0
|
32 % result.
|
wolffd@0
|
33 % mirautocor(...,'Freq') represents the autocorrelation function in the
|
wolffd@0
|
34 % frequency domain.
|
wolffd@0
|
35 % mirautocor(...,'NormalWindow',w): applies a window to the input
|
wolffd@0
|
36 % signal and divides the autocorrelation by the autocorrelation of
|
wolffd@0
|
37 % that window (Boersma 1993).
|
wolffd@0
|
38 % Possible values: any windowing function proposed in the Signal
|
wolffd@0
|
39 % Processing Toolbox (help window) plus 'rectangle' (no
|
wolffd@0
|
40 % windowing)
|
wolffd@0
|
41 % Default value: w = 'hanning'
|
wolffd@0
|
42 % mirautocor(...,'NormalWindow',0): toggles off this normalization
|
wolffd@0
|
43 % (which is on by default).
|
wolffd@0
|
44 % All the parameters described previously can be applied to an
|
wolffd@0
|
45 % autocorrelation function itself, in order to arrange the results
|
wolffd@0
|
46 % after the actual computation of the autocorrelation computations.
|
wolffd@0
|
47 % For instance: a = mirautocor(a,'Resonance','Enhanced')
|
wolffd@0
|
48 % Other optional parameter:
|
wolffd@0
|
49 % mirautocor(...,'Compres',k) computes the autocorrelation in the
|
wolffd@0
|
50 % frequency domain and includes a magnitude compression of the
|
wolffd@0
|
51 % spectral representation. A normal autocorrelation corresponds
|
wolffd@0
|
52 % to the value k=2, but values lower than 2 are suggested by
|
wolffd@0
|
53 % (Tolonen & Karjalainen, 2000).
|
wolffd@0
|
54 % Default value: k = 0.67
|
wolffd@0
|
55 % mirautocor(...,'Normal',n) or simply mirautocor(...,n) specifies
|
wolffd@0
|
56 % the normalization strategy. Accepted values are 'biased',
|
wolffd@0
|
57 % 'unbiased', 'coeff' (default value) and 'none'.
|
wolffd@0
|
58 % See help xcorr for an explanation.
|
wolffd@0
|
59
|
wolffd@0
|
60 min.key = 'Min';
|
wolffd@0
|
61 min.type = 'Integer';
|
wolffd@0
|
62 min.unit = {'s','Hz'};
|
wolffd@0
|
63 if isamir(orig,'mirspectrum')
|
wolffd@0
|
64 min.defaultunit = 'Hz';
|
wolffd@0
|
65 else
|
wolffd@0
|
66 min.defaultunit = 's';
|
wolffd@0
|
67 end
|
wolffd@0
|
68 min.default = 0;
|
wolffd@0
|
69 min.opposite = 'max';
|
wolffd@0
|
70 option.min = min;
|
wolffd@0
|
71
|
wolffd@0
|
72 max.key = 'Max';
|
wolffd@0
|
73 max.type = 'Integer';
|
wolffd@0
|
74 max.unit = {'s','Hz'};
|
wolffd@0
|
75 if isamir(orig,'mirspectrum')
|
wolffd@0
|
76 max.defaultunit = 'Hz';
|
wolffd@0
|
77 else
|
wolffd@0
|
78 max.defaultunit = 's';
|
wolffd@0
|
79 end
|
wolffd@0
|
80 if isamir(orig,'mirenvelope') || isamir(orig,'mirdiffenvelope')
|
wolffd@0
|
81 max.default = 2; % for envelopes, longest period: 2 seconds.
|
wolffd@0
|
82 elseif isamir(orig,'miraudio') || ischar(orig) % for audio signal,lowest frequency: 20 Hz.
|
wolffd@0
|
83 max.default = 1/20;
|
wolffd@0
|
84 else
|
wolffd@0
|
85 max.default = Inf;
|
wolffd@0
|
86 end
|
wolffd@0
|
87 max.opposite = 'min';
|
wolffd@0
|
88 option.max = max;
|
wolffd@0
|
89
|
wolffd@0
|
90 scaleoptbw.key = 'Normal'; %'Normal' keyword optional
|
wolffd@0
|
91 scaleoptbw.key = 'Boolean';
|
wolffd@0
|
92 option.scaleoptbw = scaleoptbw;
|
wolffd@0
|
93 scaleopt.type = 'String';
|
wolffd@0
|
94 scaleopt.choice = {'biased','unbiased','coeff','none'};
|
wolffd@0
|
95 scaleopt.default = 'coeff';
|
wolffd@0
|
96 option.scaleopt = scaleopt;
|
wolffd@0
|
97
|
wolffd@0
|
98 gener.key = {'Generalized','Compres'};
|
wolffd@0
|
99 gener.type = 'Integer';
|
wolffd@0
|
100 gener.default = 2;
|
wolffd@0
|
101 gener.keydefault = .67;
|
wolffd@0
|
102 option.gener = gener;
|
wolffd@0
|
103
|
wolffd@0
|
104 ni.key = 'NormalInput'; %% Normalize before frame or chunk??
|
wolffd@0
|
105 ni.type = 'Boolean';
|
wolffd@0
|
106 ni.default = 0;
|
wolffd@0
|
107 option.ni = ni;
|
wolffd@0
|
108
|
wolffd@0
|
109 reso.key = 'Resonance';
|
wolffd@0
|
110 reso.type = 'String';
|
wolffd@0
|
111 reso.choice = {'ToiviainenSnyder','Toiviainen','vanNoorden','no','off',0};
|
wolffd@0
|
112 reso.keydefault = 'Toiviainen';
|
wolffd@0
|
113 reso.when = 'After';
|
wolffd@0
|
114 reso.default = 0;
|
wolffd@0
|
115 option.reso = reso;
|
wolffd@0
|
116
|
wolffd@0
|
117 resocenter.key = {'Center','Centre'};
|
wolffd@0
|
118 resocenter.type = 'Integer';
|
wolffd@0
|
119 resocenter.when = 'After';
|
wolffd@0
|
120 option.resocenter = resocenter;
|
wolffd@0
|
121
|
wolffd@0
|
122 h.key = 'Halfwave';
|
wolffd@0
|
123 h.type = 'Boolean';
|
wolffd@0
|
124 h.when = 'After';
|
wolffd@0
|
125 h.default = 0;
|
wolffd@0
|
126 option.h = h;
|
wolffd@0
|
127
|
wolffd@0
|
128 e.key = 'Enhanced';
|
wolffd@0
|
129 e.type = 'Integers';
|
wolffd@0
|
130 e.default = [];
|
wolffd@0
|
131 e.keydefault = 2:10;
|
wolffd@0
|
132 e.when = 'After';
|
wolffd@0
|
133 option.e = e;
|
wolffd@0
|
134
|
wolffd@0
|
135 fr.key = 'Freq';
|
wolffd@0
|
136 fr.type = 'Boolean';
|
wolffd@0
|
137 fr.default = 0;
|
wolffd@0
|
138 fr.when = 'After';
|
wolffd@0
|
139 option.fr = fr;
|
wolffd@0
|
140
|
wolffd@0
|
141 nw.key = 'NormalWindow';
|
wolffd@0
|
142 nw.when = 'Both';
|
wolffd@0
|
143 if isamir(orig,'mirspectrum')
|
wolffd@0
|
144 nw.default = 0;
|
wolffd@0
|
145 elseif isamir(orig,'mirenvelope')
|
wolffd@0
|
146 nw.default = 'rectangular';
|
wolffd@0
|
147 else
|
wolffd@0
|
148 nw.default = 'hanning';
|
wolffd@0
|
149 end
|
wolffd@0
|
150 option.nw = nw;
|
wolffd@0
|
151
|
wolffd@0
|
152 win.key = 'Window';
|
wolffd@0
|
153 win.type = 'String';
|
wolffd@0
|
154 win.default = NaN;
|
wolffd@0
|
155 option.win = win;
|
wolffd@0
|
156
|
wolffd@0
|
157 specif.option = option;
|
wolffd@0
|
158
|
wolffd@0
|
159 specif.defaultframelength = 0.05;
|
wolffd@0
|
160 specif.defaultframehop = 0.5;
|
wolffd@0
|
161 specif.eachchunk = @eachchunk;
|
wolffd@0
|
162 specif.combinechunk = @combinechunk;
|
wolffd@0
|
163
|
wolffd@0
|
164 if isamir(orig,'mirscalar') || isamir(orig,'mirenvelope')
|
wolffd@0
|
165 specif.nochunk = 1;
|
wolffd@0
|
166 end
|
wolffd@0
|
167
|
wolffd@0
|
168 varargout = mirfunction(@mirautocor,orig,varargin,nargout,specif,@init,@main);
|
wolffd@0
|
169
|
wolffd@0
|
170
|
wolffd@0
|
171 function [x type] = init(x,option)
|
wolffd@0
|
172 type = 'mirautocor';
|
wolffd@0
|
173
|
wolffd@0
|
174
|
wolffd@0
|
175 function a = main(orig,option,postoption)
|
wolffd@0
|
176 if iscell(orig)
|
wolffd@0
|
177 orig = orig{1};
|
wolffd@0
|
178 end
|
wolffd@0
|
179 if isa(orig,'mirautocor')
|
wolffd@0
|
180 a = orig;
|
wolffd@0
|
181 if not(isempty(option)) && ...
|
wolffd@0
|
182 (option.min || iscell(option.max) || option.max < Inf)
|
wolffd@0
|
183 coeff = get(a,'Coeff');
|
wolffd@0
|
184 delay = get(a,'Delay');
|
wolffd@0
|
185 for h = 1:length(coeff)
|
wolffd@0
|
186 if a.freq
|
wolffd@0
|
187 mi = 1/option.max;
|
wolffd@0
|
188 ma = 1/option.min;
|
wolffd@0
|
189 else
|
wolffd@0
|
190 mi = option.min;
|
wolffd@0
|
191 ma = option.max;
|
wolffd@0
|
192 end
|
wolffd@0
|
193 for k = 1:length(coeff{h})
|
wolffd@0
|
194 range = find(and(delay{h}{k}(:,1,1) >= mi,...
|
wolffd@0
|
195 delay{h}{k}(:,1,1) <= ma));
|
wolffd@0
|
196 coeff{h}{k} = coeff{h}{k}(range,:,:);
|
wolffd@0
|
197 delay{h}{k} = delay{h}{k}(range,:,:);
|
wolffd@0
|
198 end
|
wolffd@0
|
199 end
|
wolffd@0
|
200 a = set(a,'Coeff',coeff,'Delay',delay);
|
wolffd@0
|
201 end
|
wolffd@0
|
202 if not(isempty(postoption)) && not(isequal(postoption,0))
|
wolffd@0
|
203 a = post(a,postoption);
|
wolffd@0
|
204 end
|
wolffd@0
|
205 elseif ischar(orig)
|
wolffd@0
|
206 a = mirautocor(miraudio(orig),option,postoption);
|
wolffd@0
|
207 else
|
wolffd@0
|
208 if nargin == 0
|
wolffd@0
|
209 orig = [];
|
wolffd@0
|
210 end
|
wolffd@0
|
211 a.freq = 0;
|
wolffd@0
|
212 a.ofspectrum = 0;
|
wolffd@0
|
213 a.window = {};
|
wolffd@0
|
214 a.normalwindow = 0;
|
wolffd@0
|
215 a = class(a,'mirautocor',mirdata(orig));
|
wolffd@0
|
216 a = purgedata(a);
|
wolffd@0
|
217 a = set(a,'Ord','coefficients');
|
wolffd@0
|
218 sig = get(orig,'Data');
|
wolffd@0
|
219 if isa(orig,'mirspectrum')
|
wolffd@0
|
220 a = set(a,'Title','Spectrum autocorrelation','OfSpectrum',1,...
|
wolffd@0
|
221 'Abs','frequency (Hz)');
|
wolffd@0
|
222 pos = get(orig,'Pos');
|
wolffd@0
|
223 else
|
wolffd@0
|
224 if isa(orig,'mirscalar')
|
wolffd@0
|
225 a = set(a,'Title',[get(orig,'Title') ' autocorrelation']);
|
wolffd@0
|
226 pos = get(orig,'FramePos');
|
wolffd@0
|
227 for k = 1:length(sig)
|
wolffd@0
|
228 for l = 1:length(sig{k})
|
wolffd@0
|
229 sig{k}{l} = sig{k}{l}';
|
wolffd@0
|
230 pos{k}{l} = pos{k}{l}(1,:,:)';
|
wolffd@0
|
231 end
|
wolffd@0
|
232 end
|
wolffd@0
|
233 else
|
wolffd@0
|
234 if isa(orig,'mirenvelope')
|
wolffd@0
|
235 a = set(a,'Title','Envelope autocorrelation');
|
wolffd@0
|
236 elseif not(isa(orig,'mirautocor'))
|
wolffd@0
|
237 a = set(a,'Title','Waveform autocorrelation');
|
wolffd@0
|
238 end
|
wolffd@0
|
239 pos = get(orig,'Pos');
|
wolffd@0
|
240 end
|
wolffd@0
|
241 a = set(a,'Abs','lag (s)');
|
wolffd@0
|
242 end
|
wolffd@0
|
243 f = get(orig,'Sampling');
|
wolffd@0
|
244
|
wolffd@0
|
245 if isnan(option.win)
|
wolffd@0
|
246 if isequal(option.nw,0) || ...
|
wolffd@0
|
247 strcmpi(option.nw,'Off') || strcmpi(option.nw,'No')
|
wolffd@0
|
248 option.win = 0;
|
wolffd@0
|
249 elseif isequal(option.nw,1) || strcmpi(option.nw,'On') || ...
|
wolffd@0
|
250 strcmpi(option.nw,'Yes')
|
wolffd@0
|
251 option.win = 'hanning';
|
wolffd@0
|
252 else
|
wolffd@0
|
253 option.win = postoption.nw;
|
wolffd@0
|
254 end
|
wolffd@0
|
255 end
|
wolffd@0
|
256
|
wolffd@0
|
257 coeff = cell(1,length(sig));
|
wolffd@0
|
258 lags = cell(1,length(sig));
|
wolffd@0
|
259 wind = cell(1,length(sig));
|
wolffd@0
|
260 for k = 1:length(sig)
|
wolffd@0
|
261 s = sig{k};
|
wolffd@0
|
262 p = pos{k};
|
wolffd@0
|
263 fk = f{k};
|
wolffd@0
|
264 if iscell(option.max)
|
wolffd@0
|
265 mi = option.min{k};
|
wolffd@0
|
266 ma = option.max{k};
|
wolffd@0
|
267 else
|
wolffd@0
|
268 mi = option.min;
|
wolffd@0
|
269 ma = option.max;
|
wolffd@0
|
270 end
|
wolffd@0
|
271 coeffk = cell(1,length(s));
|
wolffd@0
|
272 lagsk = cell(1,length(s));
|
wolffd@0
|
273 windk = cell(1,length(s));
|
wolffd@0
|
274 for l = 1:length(s)
|
wolffd@0
|
275 sl = s{l};
|
wolffd@0
|
276 sl(isnan(sl)) = 0;
|
wolffd@0
|
277 if option.ni
|
wolffd@0
|
278 mxsl = repmat(max(sl),[size(sl,1),1,1]);
|
wolffd@0
|
279 mnsl = repmat(min(sl),[size(sl,1),1,1]);
|
wolffd@0
|
280 sl = (sl-mnsl)./(mxsl-mnsl);
|
wolffd@0
|
281 end
|
wolffd@0
|
282 pl = p{l};
|
wolffd@0
|
283 pl = pl-repmat(pl(1,:,:),[size(pl,1),1,1]);
|
wolffd@0
|
284 ls = size(sl,1);
|
wolffd@0
|
285 if mi
|
wolffd@0
|
286 misp = find(pl(:,1,1)>=mi);
|
wolffd@0
|
287 if isempty(misp)
|
wolffd@0
|
288 warning('WARNING IN MIRAUTOCOR: The specified range of delays exceeds the temporal length of the signal.');
|
wolffd@0
|
289 disp('Minimum delay set to zero.')
|
wolffd@0
|
290 misp = 1; % misp is the lowest index of the lag range
|
wolffd@0
|
291 mi = 0;
|
wolffd@0
|
292 else
|
wolffd@0
|
293 misp = misp(1);
|
wolffd@0
|
294 end
|
wolffd@0
|
295 else
|
wolffd@0
|
296 misp = 1;
|
wolffd@0
|
297 end
|
wolffd@0
|
298 if ma
|
wolffd@0
|
299 masp = find(pl(:,1,1)>=ma);
|
wolffd@0
|
300 if isempty(masp)
|
wolffd@0
|
301 masp = Inf;
|
wolffd@0
|
302 else
|
wolffd@0
|
303 masp = masp(1);
|
wolffd@0
|
304 end
|
wolffd@0
|
305 else
|
wolffd@0
|
306 masp = Inf;
|
wolffd@0
|
307 end
|
wolffd@0
|
308 masp = min(masp,ceil(ls/2));
|
wolffd@0
|
309 if masp <= misp
|
wolffd@0
|
310 if size(sl,2) > 1
|
wolffd@0
|
311 warning('WARNING IN MIRAUTOCOR: Frame length is too small.');
|
wolffd@0
|
312 else
|
wolffd@0
|
313 warning('WARNING IN MIRAUTOCOR: The audio sequence is too small.');
|
wolffd@0
|
314 end
|
wolffd@0
|
315 display('The autocorrelation is not defined for this range of delays.');
|
wolffd@0
|
316 end
|
wolffd@0
|
317 sl = center(sl);
|
wolffd@0
|
318 if not(ischar(option.win)) || strcmpi(option.win,'Rectangular')
|
wolffd@0
|
319 kw = ones(size(sl));
|
wolffd@0
|
320 else
|
wolffd@0
|
321 N = size(sl,1);
|
wolffd@0
|
322 winf = str2func(option.win);
|
wolffd@0
|
323 try
|
wolffd@0
|
324 w = window(winf,N);
|
wolffd@0
|
325 catch
|
wolffd@0
|
326 if strcmpi(option.win,'hamming')
|
wolffd@0
|
327 disp('Signal Processing Toolbox does not seem to be installed. Recompute the hamming window manually.');
|
wolffd@0
|
328 w = 0.54 - 0.46 * cos(2*pi*(0:N-1)'/(N-1));
|
wolffd@0
|
329 else
|
wolffd@0
|
330 warning(['WARNING in MIRAUTOCOR: Unknown windowing function ',option.win,' (maybe Signal Processing Toolbox is not installed).']);
|
wolffd@0
|
331 disp('No windowing performed.')
|
wolffd@0
|
332 w = ones(size(sl,1),1);
|
wolffd@0
|
333 end
|
wolffd@0
|
334 end
|
wolffd@0
|
335 kw = repmat(w,[1,size(sl,2),size(sl,3)]);
|
wolffd@0
|
336 sl = sl.* kw;
|
wolffd@0
|
337 end
|
wolffd@0
|
338
|
wolffd@0
|
339 if strcmpi(option.scaleopt,'coeff')
|
wolffd@0
|
340 scaleopt = 'none';
|
wolffd@0
|
341 else
|
wolffd@0
|
342 scaleopt = option.scaleopt;
|
wolffd@0
|
343 end
|
wolffd@0
|
344 c = zeros(masp,size(sl,2),size(sl,3));
|
wolffd@0
|
345 for i = 1:size(sl,2)
|
wolffd@0
|
346 for j = 1:size(sl,3)
|
wolffd@0
|
347 if option.gener == 2
|
wolffd@0
|
348 cc = xcorr(sl(:,i,j),masp-1,scaleopt);
|
wolffd@0
|
349 c(:,i,j) = cc(masp:end);
|
wolffd@0
|
350 else
|
wolffd@0
|
351 ss = abs(fft(sl(:,i,j)));
|
wolffd@0
|
352 ss = ss.^option.gener;
|
wolffd@0
|
353 cc = ifft(ss);
|
wolffd@0
|
354 ll = (0:masp-1);
|
wolffd@0
|
355 c(:,i,j) = cc(ll+1);
|
wolffd@0
|
356 end
|
wolffd@0
|
357 end
|
wolffd@0
|
358 if strcmpi(option.scaleopt,'coeff') && option.gener == 2
|
wolffd@0
|
359 % to be adapted to generalized autocor
|
wolffd@0
|
360 c(:,i,:) = c(:,i,:)/xcorr(sum(sl(:,i,:),3),0);
|
wolffd@0
|
361 % This is a kind of generalization of the 'coeff'
|
wolffd@0
|
362 % normalization for multi-channels signals. In the
|
wolffd@0
|
363 % original 'coeff' option, the autocorrelation at zero
|
wolffd@0
|
364 % lag is identically 1.0. In this multi-channels
|
wolffd@0
|
365 % version, the autocorrelation at zero lag is such that
|
wolffd@0
|
366 % the sum over channels becomes identically 1.0.
|
wolffd@0
|
367 end
|
wolffd@0
|
368 end
|
wolffd@0
|
369 coeffk{l} = c(misp:end,:,:);
|
wolffd@0
|
370 pl = pl(find(pl(:,1,1) >=mi),:,:);
|
wolffd@0
|
371 lagsk{l} = pl(1:min(size(coeffk{l},1),size(pl,1)),:,:);
|
wolffd@0
|
372 windk{l} = kw;
|
wolffd@0
|
373 end
|
wolffd@0
|
374 coeff{k} = coeffk;
|
wolffd@0
|
375 lags{k} = lagsk;
|
wolffd@0
|
376 wind{k} = windk;
|
wolffd@0
|
377 end
|
wolffd@0
|
378 a = set(a,'Coeff',coeff,'Delay',lags,'Window',wind);
|
wolffd@0
|
379 if not(isempty(postoption))
|
wolffd@0
|
380 a = post(a,postoption);
|
wolffd@0
|
381 end
|
wolffd@0
|
382 end
|
wolffd@0
|
383
|
wolffd@0
|
384
|
wolffd@0
|
385 function a = post(a,option)
|
wolffd@0
|
386 debug = 0;
|
wolffd@0
|
387 coeff = get(a,'Coeff');
|
wolffd@0
|
388 lags = get(a,'Delay');
|
wolffd@0
|
389 wind = get(a,'Window');
|
wolffd@0
|
390 freq = option.fr && not(get(a,'FreqDomain'));
|
wolffd@0
|
391 if isequal(option.e,1)
|
wolffd@0
|
392 option.e = 2:10;
|
wolffd@0
|
393 end
|
wolffd@0
|
394 if max(option.e) > 1
|
wolffd@0
|
395 pa = mirpeaks(a,'NoBegin','NoEnd','Contrast',.01);
|
wolffd@0
|
396 va = mirpeaks(a,'Valleys','Contrast',.01);
|
wolffd@0
|
397 pv = get(pa,'PeakVal');
|
wolffd@0
|
398 vv = get(va,'PeakVal');
|
wolffd@0
|
399 end
|
wolffd@0
|
400 for k = 1:length(coeff)
|
wolffd@0
|
401 for l = 1:length(coeff{k})
|
wolffd@0
|
402 c = coeff{k}{l}; % Coefficients of autocorrelation
|
wolffd@0
|
403 t = lags{k}{l}; % Delays of autocorrelation
|
wolffd@0
|
404 if not(isempty(c))
|
wolffd@0
|
405 if not(isequal(option.nw,0) || strcmpi(option.nw,'No') || ...
|
wolffd@0
|
406 strcmpi(option.nw,'Off') || a.normalwindow) % 'NormalWindow' option
|
wolffd@0
|
407 xw = zeros(size(c));
|
wolffd@0
|
408 lc = size(c,1);
|
wolffd@0
|
409 for j = 1:size(c,3)
|
wolffd@0
|
410 for i = 1:size(c,2)
|
wolffd@0
|
411 xwij = xcorr(wind{k}{l}(:,i,j),lc,'coeff');
|
wolffd@0
|
412 xw(:,i,j) = xwij(lc+2:end);
|
wolffd@0
|
413 end
|
wolffd@0
|
414 end
|
wolffd@0
|
415 c = c./ xw;
|
wolffd@0
|
416 a.normalwindow = 1;
|
wolffd@0
|
417 end
|
wolffd@0
|
418 if ischar(option.reso) && ...
|
wolffd@0
|
419 (strcmpi(option.reso,'ToiviainenSnyder') || ...
|
wolffd@0
|
420 strcmpi(option.reso,'Toiviainen') || ...
|
wolffd@0
|
421 strcmpi(option.reso,'vanNoorden'))
|
wolffd@0
|
422 if isa(a,'mirautocor') && get(a,'FreqDomain')
|
wolffd@0
|
423 ll = 1./t;
|
wolffd@0
|
424 else
|
wolffd@0
|
425 ll = t;
|
wolffd@0
|
426 end
|
wolffd@0
|
427 if not(option.resocenter)
|
wolffd@0
|
428 option.resocenter = .5;
|
wolffd@0
|
429 end
|
wolffd@0
|
430 if strcmpi(option.reso,'ToiviainenSnyder') || ...
|
wolffd@0
|
431 strcmpi(option.reso,'Toiviainen')
|
wolffd@0
|
432 w = max(1 - 0.25*(log2(max(ll,1e-12)/option.resocenter)).^2, 0);
|
wolffd@0
|
433 elseif strcmpi(option.reso,'vanNoorden')
|
wolffd@0
|
434 f0=2.193; b=option.resocenter;
|
wolffd@0
|
435 f=1./ll; a1=(f0*f0-f.*f).^2+b*f.^2; a2=f0^4+f.^4;
|
wolffd@0
|
436 w=(1./sqrt(a1))-(1./sqrt(a2));
|
wolffd@0
|
437 end
|
wolffd@0
|
438 if max(w) == 0
|
wolffd@0
|
439 warning('The resonance curve, not defined for this range of delays, will not be applied.')
|
wolffd@0
|
440 else
|
wolffd@0
|
441 w = w/max(w);
|
wolffd@0
|
442 c = c.* repmat(w,[1,size(c,2),size(c,3)]);
|
wolffd@0
|
443 end
|
wolffd@0
|
444 end
|
wolffd@0
|
445 if option.h
|
wolffd@0
|
446 c = hwr(c);
|
wolffd@0
|
447 end
|
wolffd@0
|
448 if max(option.e) > 1
|
wolffd@0
|
449 if a.freq
|
wolffd@0
|
450 freq = 1;
|
wolffd@0
|
451 for i = 1:size(c,3)
|
wolffd@0
|
452 c(:,:,i) = flipud(c(:,:,i));
|
wolffd@0
|
453 end
|
wolffd@0
|
454 t = flipud(1./t);
|
wolffd@0
|
455 end
|
wolffd@0
|
456
|
wolffd@0
|
457 for g = 1:size(c,2)
|
wolffd@0
|
458 for h = 1:size(c,3)
|
wolffd@0
|
459 cgh = c(:,g,h);
|
wolffd@0
|
460 if length(cgh)>1
|
wolffd@0
|
461 pvk = pv{k}{l}{1,g,h};
|
wolffd@0
|
462 mv = [];
|
wolffd@0
|
463 if not(isempty(pvk))
|
wolffd@0
|
464 mp = min(pv{k}{l}{1,g,h}); %Lowest peak
|
wolffd@0
|
465 vvv = vv{k}{l}{1,g,h}; %Valleys
|
wolffd@0
|
466 mv = vvv(find(vvv<mp,1,'last'));
|
wolffd@0
|
467 %Highest valley below the lowest peak
|
wolffd@0
|
468
|
wolffd@0
|
469 if not(isempty(mv))
|
wolffd@0
|
470 cgh = cgh-mv;
|
wolffd@0
|
471 end
|
wolffd@0
|
472 end
|
wolffd@0
|
473 cgh2 = cgh;
|
wolffd@0
|
474 tgh2 = t(:,g,1);
|
wolffd@0
|
475 coef = cgh(2)-cgh(1); % initial slope of the autocor curve
|
wolffd@0
|
476 tcoef = tgh2(2)-tgh2(1);
|
wolffd@0
|
477 deter = 0;
|
wolffd@0
|
478 inter = 0;
|
wolffd@0
|
479
|
wolffd@0
|
480 repet = find(not(diff(tgh2))); % Avoid bug if repeated x-values
|
wolffd@0
|
481 if repet
|
wolffd@0
|
482 warning('WARNING in MIRAUTOCOR: Two successive samples have exactly same temporal position.');
|
wolffd@0
|
483 tgh2(repet+1) = tgh2(repet)+1e-12;
|
wolffd@0
|
484 end
|
wolffd@0
|
485
|
wolffd@0
|
486 if coef < 0
|
wolffd@0
|
487 % initial descending slope removed
|
wolffd@0
|
488 deter = find(diff(cgh2)>0,1)-1;
|
wolffd@0
|
489 % number of removed points
|
wolffd@0
|
490 if isempty(deter)
|
wolffd@0
|
491 deter = 0;
|
wolffd@0
|
492 end
|
wolffd@0
|
493 cgh2(1:deter) = [];
|
wolffd@0
|
494 tgh2(1:deter) = [];
|
wolffd@0
|
495 coef = cgh2(2)-cgh2(1);
|
wolffd@0
|
496 end
|
wolffd@0
|
497
|
wolffd@0
|
498 if coef > 0
|
wolffd@0
|
499 % initial ascending slope prolonged to the left
|
wolffd@0
|
500 % until it reaches the x-axis
|
wolffd@0
|
501 while cgh2(1) > 0
|
wolffd@0
|
502 coef = coef*1.1;
|
wolffd@0
|
503 % the further to the left, ...
|
wolffd@0
|
504 % the more ascending is the slope
|
wolffd@0
|
505 % (not sure it always works, though...)
|
wolffd@0
|
506 inter = inter+1;
|
wolffd@0
|
507 % number of added points
|
wolffd@0
|
508 cgh2 = [cgh2(1)-coef; cgh2];
|
wolffd@0
|
509 tgh2 = [tgh2(1)-tcoef; tgh2];
|
wolffd@0
|
510 end
|
wolffd@0
|
511 cgh2(1) = 0;
|
wolffd@0
|
512 end
|
wolffd@0
|
513
|
wolffd@0
|
514 for i = option.e % Enhancing procedure
|
wolffd@0
|
515 % option.e is the list of scaling factors
|
wolffd@0
|
516 % i is the scaling factor
|
wolffd@0
|
517 if i
|
wolffd@0
|
518 be = find(tgh2 & tgh2/i >= tgh2(1),1);
|
wolffd@0
|
519 % starting point of the substraction
|
wolffd@0
|
520 % on the X-axis
|
wolffd@0
|
521
|
wolffd@0
|
522 if not(isempty(be))
|
wolffd@0
|
523 ic = interp1(tgh2,cgh2,tgh2/i);
|
wolffd@0
|
524 % The scaled autocorrelation
|
wolffd@0
|
525 ic(1:be-1) = 0;
|
wolffd@0
|
526 ic(find(isnan(ic))) = Inf;
|
wolffd@0
|
527 % All the NaN values are changed
|
wolffd@0
|
528 % into 0 in the resulting curve
|
wolffd@0
|
529 ic = max(ic,0);
|
wolffd@0
|
530
|
wolffd@0
|
531 if debug
|
wolffd@0
|
532 hold off,plot(tgh2,cgh2)
|
wolffd@0
|
533 end
|
wolffd@0
|
534
|
wolffd@0
|
535 cgh2 = cgh2 - ic;
|
wolffd@0
|
536 % The scaled autocorrelation
|
wolffd@0
|
537 % is substracted to the initial one
|
wolffd@0
|
538
|
wolffd@0
|
539 cgh2 = max(cgh2,0);
|
wolffd@0
|
540 % Half-wave rectification
|
wolffd@0
|
541
|
wolffd@0
|
542 if debug
|
wolffd@0
|
543 hold on,plot(tgh2,ic,'r')
|
wolffd@0
|
544 hold on,plot(tgh2,cgh2,'g')
|
wolffd@0
|
545 drawnow
|
wolffd@0
|
546 figure
|
wolffd@0
|
547 end
|
wolffd@0
|
548 end
|
wolffd@0
|
549 end
|
wolffd@0
|
550 end
|
wolffd@0
|
551
|
wolffd@0
|
552 % The temporary modifications are
|
wolffd@0
|
553 % removed from the final curve
|
wolffd@0
|
554 if inter>=deter
|
wolffd@0
|
555 c(:,g,h) = cgh2(inter-deter+1:end);
|
wolffd@0
|
556 if not(isempty(mv))
|
wolffd@0
|
557 c(:,g,h) = c(:,g,h) + mv;
|
wolffd@0
|
558 end
|
wolffd@0
|
559 else
|
wolffd@0
|
560 c(:,g,h) = [zeros(deter-inter,1);cgh2];
|
wolffd@0
|
561 end
|
wolffd@0
|
562 end
|
wolffd@0
|
563 end
|
wolffd@0
|
564 end
|
wolffd@0
|
565 end
|
wolffd@0
|
566 if freq
|
wolffd@0
|
567 if t(1,1) == 0
|
wolffd@0
|
568 c = c(2:end,:,:);
|
wolffd@0
|
569 t = t(2:end,:,:);
|
wolffd@0
|
570 end
|
wolffd@0
|
571 for i = 1:size(c,3)
|
wolffd@0
|
572 c(:,:,i) = flipud(c(:,:,i));
|
wolffd@0
|
573 end
|
wolffd@0
|
574 t = flipud(1./t);
|
wolffd@0
|
575 end
|
wolffd@0
|
576 coeff{k}{l} = c;
|
wolffd@0
|
577 lags{k}{l} = t;
|
wolffd@0
|
578 end
|
wolffd@0
|
579 end
|
wolffd@0
|
580 end
|
wolffd@0
|
581 a = set(a,'Coeff',coeff,'Delay',lags,'Freq');
|
wolffd@0
|
582 if freq
|
wolffd@0
|
583 a = set(a,'FreqDomain',1,'Abs','frequency (Hz)');
|
wolffd@0
|
584 end
|
wolffd@0
|
585
|
wolffd@0
|
586
|
wolffd@0
|
587 function [y orig] = eachchunk(orig,option,missing,postchunk)
|
wolffd@0
|
588 option.scaleopt = 'none';
|
wolffd@0
|
589 y = mirautocor(orig,option);
|
wolffd@0
|
590
|
wolffd@0
|
591
|
wolffd@0
|
592 function y = combinechunk(old,new)
|
wolffd@0
|
593 do = get(old,'Data');
|
wolffd@0
|
594 do = do{1}{1};
|
wolffd@0
|
595 dn = get(new,'Data');
|
wolffd@0
|
596 dn = dn{1}{1};
|
wolffd@0
|
597 if abs(size(dn,1)-size(do,1)) <= 2 % Probleme of border fluctuation
|
wolffd@0
|
598 mi = min(size(dn,1),size(do,1));
|
wolffd@0
|
599 dn = dn(1:mi,:,:);
|
wolffd@0
|
600 do = do(1:mi,:,:);
|
wolffd@0
|
601 elseif length(dn) < length(do)
|
wolffd@0
|
602 dn(length(do),:,:) = 0; % Zero-padding
|
wolffd@0
|
603 end
|
wolffd@0
|
604 y = set(old,'ChunkData',do+dn);
|