Mercurial > hg > camir-aes2014

comparison toolboxes/MIRtoolbox1.3.2/MIRToolbox/mirpeaks.m @ 0:e9a9cd732c1e tip

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
first hg version after svn
author wolffd
date Tue, 10 Feb 2015 15:05:51 +0000
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:e9a9cd732c1e
1 function varargout = mirpeaks(orig,varargin)
2 % p = mirpeaks(x) detect peaks in x.
3 % Optional argument:
4 % mirpeaks(...,'Total',m): only the m highest peaks are selected.
5 % If m=Inf, no limitation of number of peaks.
6 % Default value: m = Inf
7 % mirpeaks(...,'Order',o): specifies the ordering of the peaks.
8 % Possible values for o:
9 % 'Amplitude': orders the peaks from highest to lowest
10 % (Default choice.)
11 % 'Abscissa': orders the peaks along the abscissa axis.
12 % mirpeaks(...,'Contrast',cthr): a threshold value. A given local
13 % maximum will be considered as a peak if its distance with the
14 % previous and successive local minima (if any) is higher than
15 % this threshold. This distance is expressed with respect to the
16 % total amplitude of x: a distance of 1, for instance, is
17 % equivalent to the distance between the maximum and the minimum
18 % of x.
19 % Default value: cthr = 0.1
20 % mirpeaks(...,'SelectFirst',fthr): If the 'Contrast' selection has
21 % been chosen, this additional option specifies that when one
22 % peak has to be chosen out of two candidates, and if the
23 % difference of their amplitude is below the threshold fthr,
24 % then the most ancien one is selected.
25 % Option toggled off by default:
26 % Default value if toggled on: fthr = cthr/2
27 % mirpeaks(...,'Threshold',thr): a threshold value.
28 % A given local maximum will be considered as a peak if its
29 % normalized amplitude is higher than this threshold.
30 % A given local minimum will be considered as a valley if its
31 % normalized amplitude is lower than this threshold.
32 % The normalized amplitude can have value between 0 (the minimum
33 % of the signal in each frame) and 1 (the maximum in each
34 % frame).
35 % Default value: thr=0 for peaks thr = 1 for valleys
36 % mirpeaks(...,'Interpol',i): estimates more precisely the peak
37 % position and amplitude using interpolation. Performed only on
38 % data with numerical abscissae axis.
39 % Possible value for i:
40 % '', 'no', 'off', 0: no interpolation
41 % 'Quadratic': quadratic interpolation. (default value).
42 % mirpeaks(...,'Valleys'): detect valleys instead of peaks.
43 % mirpeaks(...,'Reso',r): removes peaks whose distance to one or
44 % several higher peaks is lower than a given threshold.
45 % Possible value for the threshold r:
46 % 'SemiTone': ratio between the two peak positions equal to
47 % 2^(1/12)
48 % a numerical value : difference between the two peak
49 % positions equal to that value
50 % When two peaks are distant by an interval lower than the
51 % resolution, the highest of them is selected by default.
52 % mirpeaks(...,'Reso',r,'First') specifies on the contrary that
53 % the first of them is selected by default.
54 % mirpeaks(...,'Nearest',t,s): takes the peak nearest a given abscisse
55 % values t. The distance is computed either on a linear scale
56 % (s = 'Lin') or logarithmic scale (s = 'Log'). In this case,
57 % only one peak is extracted.
58 % mirpeaks(...,'Pref',c,std): indicates a region of preference for
59 % the peak picking, centered on the abscisse value c, with a
60 % standard deviation of std.
61 % mirpeaks(...,'NoBegin'): does not consider the first sample as a
62 % possible peak candidate.
63 % mirpeaks(...,'NoEnd'): does not consider the last sample as a possible
64 % peak candidate.
65 % mirpeaks(...,'Normalize',n): specifies whether frames are
66 % normalized globally or individually.
67 % Possible value for n:
68 % 'Global': normalizes the whole frames altogether from 0 to
69 % 1 (default choice).
70 % 'Local': normalizes each frame from 0 to 1.
71 % mirpeaks(...,'Extract'): extracts from the initial curves all the
72 % positive continuous segments (or "curve portions") where peaks
73 % are located.
74 % mirpeaks(...,'Only'): keeps from the original curve only the data
75 % corresponding to the peaks, and zeroes the remaining data.
76 % mirpeaks(...,'Track',t): tracks temporal continuities of peaks. If
77 % a value t is specified, the variation between successive peaks
78 % is tolerated up to t samples.
79 % mirpeaks(...,'CollapseTrack',ct): collapses tracks into one single
80 % track, and remove small track transitions, of length shorter
81 % than ct samples. Default value: ct = 7
82
83 m.key = 'Total';
84 m.type = 'Integer';
85 m.default = Inf;
86 option.m = m;
87
88 nobegin.key = 'NoBegin';
89 nobegin.type = 'Boolean';
90 nobegin.default = 0;
91 option.nobegin = nobegin;
92
93 noend.key = 'NoEnd';
94 noend.type = 'Boolean';
95 noend.default = 0;
96 option.noend = noend;
97
98 order.key = 'Order';
99 order.type = 'String';
100 order.choice = {'Amplitude','Abscissa'};
101 order.default = 'Amplitude';
102 option.order = order;
103
104 chro.key = 'Chrono'; % obsolete, corresponds to 'Order','Abscissa'
105 chro.type = 'Boolean';
106 chro.default = 0;
107 option.chro = chro;
108
109 ranked.key = 'Ranked'; % obsolete, corresponds to 'Order','Amplitude'
110 ranked.type = 'Boolean';
111 ranked.default = 0;
112 option.ranked = ranked;
113
114 vall.key = 'Valleys';
115 vall.type = 'Boolean';
116 vall.default = 0;
117 option.vall = vall;
118
119 cthr.key = 'Contrast';
120 cthr.type = 'Integer';
121 cthr.default = .1;
122 option.cthr = cthr;
123
124 first.key = 'SelectFirst';
125 first.type = 'Integer';
126 first.default = 0;
127 first.keydefault = NaN;
128 option.first = first;
129
130 thr.key = 'Threshold';
131 thr.type = 'Integer';
132 thr.default = NaN;
133 option.thr = thr;
134
135 smthr.key = 'MatrixThreshold'; % to be documented in version 1.3
136 smthr.type = 'Integer';
137 smthr.default = NaN;
138 option.smthr = smthr;
139
140 graph.key = 'Graph';
141 graph.type = 'Integer';
142 graph.default = 0;
143 graph.keydefault = .25;
144 option.graph = graph;
145
146 interpol.key = 'Interpol';
147 interpol.type = 'String';
148 interpol.default = 'Quadratic';
149 interpol.keydefault = 'Quadratic';
150 option.interpol = interpol;
151
152 reso.key = 'Reso';
153 %reso.type = 'String';
154 %reso.choice = {0,'SemiTone'};
155 reso.default = 0;
156 option.reso = reso;
157
158 resofirst.key = 'First';
159 resofirst.type = 'Boolean';
160 resofirst.default = 0;
161 option.resofirst = resofirst;
162
163 c.key = 'Pref';
164 c.type = 'Integer';
165 c.number = 2;
166 c.default = [0 0];
167 option.c = c;
168
169 near.key = 'Nearest';
170 near.type = 'Integer';
171 near.default = NaN;
172 option.near = near;
173
174 logsc.type = 'String';
175 logsc.choice = {'Lin','Log',0};
176 logsc.default = 'Lin';
177 option.logsc = logsc;
178
179 normal.key = 'Normalize';
180 normal.type = 'String';
181 normal.choice = {'Local','Global'};
182 normal.default = 'Global';
183 option.normal = normal;
184
185 extract.key = 'Extract';
186 extract.type = 'Boolean';
187 extract.default = 0;
188 option.extract = extract;
189
190 only.key = 'Only';
191 only.type = 'Boolean';
192 only.default = 0;
193 option.only = only;
194
195 delta.key = 'Track';
196 delta.type = 'Integer';
197 delta.default = 0;
198 delta.keydefault = Inf;
199 option.delta = delta;
200
201 mem.key = 'TrackMem';
202 mem.type = 'Integer';
203 mem.default = Inf;
204 option.mem = mem;
205
206 shorttrackthresh.key = 'CollapseTracks';
207 shorttrackthresh.type = 'Integer';
208 shorttrackthresh.default = 0;
209 shorttrackthresh.keydefault = 7;
210 option.shorttrackthresh = shorttrackthresh;
211
212 scan.key = 'ScanForward'; % specific to mironsets(..., 'Klapuri99')
213 scan.default = [];
214 option.scan = scan;
215
216 specif.option = option;
217
218 varargout = mirfunction(@mirpeaks,orig,varargin,nargout,specif,@init,@main);
219
220
221 function [x type] = init(x,option)
222 type = mirtype(x);
223
224
225 function p = main(x,option,postoption)
226 if iscell(x)
227 x = x{1};
228 end
229 if option.chro
230 option.order = 'Abscissa';
231 elseif option.ranked
232 option.order = 'Amplitude';
233 end
234 if not(isnan(option.near)) && option.m == 1
235 option.m = Inf;
236 end
237 x = purgedata(x);
238 if option.m <= 0
239 p = x;
240 return
241 end
242 d = get(x,'Data');
243 sr = get(x,'Sampling');
244 cha = 0; % Indicates when it is possible to represent as a curve along the
245 % Z-axis (channels) instead of the X-axis (initial abscissa).
246 if isnan(option.first)
247 option.first = option.cthr / 2;
248 end
249 if isa(x,'mirscalar')
250 t = get(x,'FramePos');
251 for i = 1:length(d)
252 for j = 1:length(d{i})
253 d{i}{j} = d{i}{j}';
254 if size(t{i},1) == 1
255 t{i}{j} = t{i}{j}(1,:,:)';
256 else
257 t{i}{j} = (t{i}{j}(1,:,:)+t{i}{j}(2,:,:))'/2;
258 end
259 end
260 end
261 elseif isa(x,'mirsimatrix')
262 t = get(x,'FramePos');
263 for i = 1:length(t)
264 for j = 1:length(t{i})
265 t{i}{j} = repmat((t{i}{j}(1,:,:)+t{i}{j}(2,:,:))'/2,...
266 [1 size(d{i}{j},2) 1]);
267 end
268 end
269 elseif isa(x,'mirhisto')
270 error('ERROR IN MIRPEAKS: peaks of histogram not considered yet.');
271 else
272 for i = 1:length(d)
273 for j = 1:length(d{i})
274 if iscell(d{i})
275 dij = d{i}{j};
276 if ~cha && j == 1 && size(dij,3) > 1 && size(dij,1) == 1
277 cha = 1;
278 end
279 if cha && j > 1 && size(dij,1) > 1
280 cha = -1;
281 end
282 end
283 end
284 for j = 1:length(d{i})
285 if iscell(d{i})
286 dij = d{i}{j};
287 if cha == 1
288 if iscell(dij)
289 for k = 1:size(dij,2)
290 d{i}{j}{k} = reshape(dij{k},size(dij{k},2),size(dij{k},3));
291 d{i}{j}{k} = d{i}{j}{k}';
292 end
293 else
294 d{i}{j} = reshape(dij,size(dij,2),size(dij,3));
295 d{i}{j} = d{i}{j}';
296 end
297 end
298 end
299 end
300 end
301 if cha == 1
302 t = get(x,'Channels');
303 else
304 t = get(x,'Pos');
305 end
306 end
307 pp = cell(1,length(d));
308 pv = cell(1,length(d));
309 pm = cell(1,length(d));
310 ppp = cell(1,length(d));
311 ppv = cell(1,length(d));
312 tp = cell(1,length(d));
313 tv = cell(1,length(d));
314
315 if isnan(option.thr)
316 option.thr = 0;
317 else
318 if option.vall
319 option.thr = 1-option.thr;
320 end
321 end
322 %if isnan(option.lthr)
323 % option.lthr = 1;
324 %else
325 % if option.vall
326 % option.lthr = 1-option.lthr;
327 % end
328 %end
329 if isnan(option.smthr)
330 option.smthr = option.thr - .2;
331 end
332
333 if not(isempty(option.scan))
334 pscan = get(option.scan,'PeakPos');
335 end
336
337 interpol = get(x,'Interpolable') && not(isempty(option.interpol)) && ...
338 ((isnumeric(option.interpol) && option.interpol) || ...
339 (ischar(option.interpol) && not(strcmpi(option.interpol,'No')) && not(strcmpi(option.interpol,'Off'))));
340
341 for i = 1:length(d) % For each audio file,...
342 di = d{i};
343 if cha == 1
344 ti = t; %sure ?
345 else
346 ti = t{i};
347 end
348 if not(iscell(di))
349 di = {di};
350 if isa(x,'mirchromagram') && not(cha)
351 ti = {ti};
352 end
353 end
354 for h = 1:length(di) % For each segment,...
355 dh0 = di{h};
356 if option.vall
357 dh0 = -dh0;
358 end
359 dh2 = dh0;
360 [nl0 nc np nd0] = size(dh0);
361 if cha == 1
362 if iscell(ti)
363 %% problem here!!!
364 ti = ti{i}; %%%%%it seems that sometimes we need to use instead ti{i}(:);
365 end
366 th = repmat(ti,[1,nc,np,nd0]);
367 else
368 th = ti{h};
369 if iscell(th) % Non-numerical abscissae are transformed into numerical ones.
370 th = repmat((1:size(th,1))',[1,nc,np]);
371 else
372 if size(th,3)<np
373 th = repmat(th,[1,1,np]);
374 end
375 end
376 end
377 if option.c % If a prefered region is specified, the data is amplified accordingly
378 dh0 = dh0.*exp(-(th-option.c(1)).^2/2/option.c(2)^2)...
379 /option.c(2)/sqrt(2*pi)/2;
380 end
381
382 % Now the data is rescaled. the minimum is set to 0
383 % and the maximum to 1.
384 state = warning('query','MATLAB:divideByZero');
385 warning('off','MATLAB:divideByZero');
386
387 % Let's first normalize all frames globally:
388 dh0 = (dh0-repmat(min(min(min(dh0,[],1),[],2),[],4),[nl0 nc 1 nd0]))./...
389 repmat(max(max(max(dh0,[],1),[],2),[],4)...
390 -min(min(min(dh0,[],1),[],2),[],4),[nl0 nc 1 nd0]);
391 for l = 1:nd0
392 [unused lowc] = find(max(dh0(:,:,:,l))<option.thr);
393 dh0(:,lowc,1,l) = 0;
394 end
395
396 if strcmpi(option.normal,'Local')
397 % Normalizing each frame separately:
398 dh0 = (dh0-repmat(min(min(dh0,[],1),[],4),[nl0 1 1 nd0]))./...
399 repmat(max(max(dh0,[],1),[],4)...
400 -min(min(dh0,[],1),[],4),[nl0 1 1 nd0]);
401 end
402 warning(state.state,'MATLAB:divideByZero');
403
404 if option.nobegin
405 dh0 = [Inf(1,nc,np,nd0);dh0];
406 % This infinite value at the beginning
407 % prevents the selection of the first sample of data
408 dh2 = [Inf(1,nc,np,nd0);dh2];
409 th = [NaN(1,nc,np,1);th];
410 else
411 dh0 = [-Inf(1,nc,np,nd0);dh0];
412 % This infinite negative value at the beginning
413 % ensures the selection of the first sample of data
414 dh2 = [-Inf(1,nc,np,nd0);dh2];
415 th = [NaN(1,nc,np,1);th];
416 end
417 if option.noend
418 dh0 = [dh0;Inf(1,nc,np,nd0)];
419 % idem for the last sample of data
420 dh2 = [dh2;Inf(1,nc,np,nd0)];
421 th = [th;NaN(1,nc,np,1)];
422 else
423 dh0 = [dh0;-Inf(1,nc,np,nd0)];
424 dh2 = [dh2;-Inf(1,nc,np,nd0)];
425 th = [th;NaN(1,nc,np,1)];
426 end
427 nl0 = nl0+2;
428
429 % Rearrange the 4th dimension (if used) into the 1st one.
430 nl = nl0*nd0;
431 dh = zeros(nl,nc,np);
432 dh3 = zeros(nl,nc,np);
433 th2 = zeros(nl,nc,np);
434 for l = 1:nd0
435 dh((l-1)*nl0+(1:nl0)',:,:) = dh0(:,:,:,l);
436 dh3((l-1)*nl0+(1:nl0)',:,:) = dh2(:,:,:,l);
437 th2((l-1)*nl0+(1:nl0)',:,:) = th(:,:,:);
438 end
439
440 th = th2; % The X-abscissa
441
442 ddh = diff(dh);
443 % Let's find the local maxima
444 for l = 1:np
445 dl = dh(2:end-1,:,l);
446 for k = 1:nc
447 dk = dl(:,k);
448 mx{1,k,l} = find(and(and(dk >= option.cthr, ...
449 dk >= option.thr),...
450 ... dk <= option.lthr)),
451 and(ddh(1:(end-1),k,l) > 0, ...
452 ddh(2:end,k,l) <= 0)))+1;
453 end
454 end
455 if option.cthr
456 for l = 1:np
457 for k = 1:nc
458 finalmxk = [];
459 mxk = mx{1,k,l};
460 if not(isempty(mxk))
461 wait = 0;
462 if length(mxk)>5000
463 wait = waitbar(0,['Selecting peaks... (0 out of 0)']);
464 end
465 %if option.m < Inf
466 % [unused,idx] = sort(dh(mxk,k,l),'descend'); % The peaks are sorted in decreasing order
467 % mxk = mxk(sort(idx(1:option.m)));
468 %end
469 j = 1; % Scans the peaks from begin to end.
470 mxkj = mxk(j); % The current peak
471 jj = j+1;
472 bufmin = Inf;
473 bufmax = dh(mxkj,k,l);
474 oldbufmin = min(dh(1:mxk(1)-1,k,l));
475 while jj <= length(mxk)
476 if wait && not(mod(jj,5000))
477 waitbar(jj/length(mxk),wait,['Selecting peaks... (',num2str(length(finalmxk)),' out of ',num2str(jj),')']);
478 end
479 bufmin = min(bufmin, ...
480 min(dh(mxk(jj-1)+1:mxk(jj)-1,k,l)));
481 if bufmax - bufmin < option.cthr
482 % There is no contrastive notch
483 if dh(mxk(jj),k,l) > bufmax && ...
484 (dh(mxk(jj),k,l) - bufmax > option.first ...
485 || (bufmax - oldbufmin < option.cthr))
486 % If the new peak is significantly
487 % higher than the previous one,
488 % The peak is transfered to the new
489 % position
490 j = jj;
491 mxkj = mxk(j); % The current peak
492 bufmax = dh(mxkj,k,l);
493 oldbufmin = min(oldbufmin,bufmin);
494 bufmin = Inf;
495 elseif dh(mxk(jj),k,l) - bufmax <= option.first
496 bufmax = max(bufmax,dh(mxk(jj),k,l));
497 oldbufmin = min(oldbufmin,bufmin);
498 end
499 else
500 % There is a contrastive notch
501 if bufmax - oldbufmin < option.cthr
502 % But the previous peak candidate
503 % is too weak and therefore
504 % discarded
505 oldbufmin = min(oldbufmin,bufmin);
506 else
507 % The previous peak candidate is OK
508 % and therefore stored
509 finalmxk(end+1) = mxkj;
510 oldbufmin = bufmin;
511 end
512 bufmax = dh(mxk(jj),k,l);
513 j = jj;
514 mxkj = mxk(j); % The current peak
515 bufmin = Inf;
516 end
517 jj = jj+1;
518 end
519 if bufmax - oldbufmin >= option.cthr && ...
520 bufmax - min(dh(mxk(j)+1:end,k,l)) >= option.cthr
521 % The last peak candidate is OK and stored
522 finalmxk(end+1) = mxk(j);
523 end
524 if wait
525 waitbar(1,wait);
526 close(wait);
527 drawnow
528 end
529 end
530 mx{1,k,l} = finalmxk;
531 end
532 end
533 end
534 if not(isempty(option.scan)) % 'ScanForward' option, used for 'Klapuri99' in mironsets
535 for l = 1:np
536 for k = 1:nc
537 pscankl = pscan{i}{h}{1,k,l}; % scan seed positions
538 mxkl = [];
539 lp = length(pscankl); % number of peaks
540 for jj = 1:lp % for each peak
541 fmx = find(mx{1,k,l}>pscankl(jj),1);
542 % position of the next max following the
543 % current seed
544 fmx = mx{1,k,l}(fmx);
545 if jj<lp && (isempty(fmx) || fmx>=pscankl(jj+1))
546 [unused fmx] = max(dh(pscankl(jj):...
547 pscankl(jj+1)-1,k,l));
548 fmx = fmx+pscankl(jj)-1;
549 elseif jj==lp && isempty(fmx)
550 [unused fmx] = max(dh(pscankl(jj):end,k,l));
551 fmx = fmx+pscankl(jj)-1;
552 end
553 mxkl = [mxkl fmx];
554 end
555 mx{1,k,l} = mxkl;
556 end
557 end
558 end
559 if not(isequal(option.reso,0)) % Removing peaks too close to higher peaks
560 if ischar(option.reso) && strcmpi(option.reso,'SemiTone')
561 compar = @semitone_compar;
562 elseif isnumeric(option.reso)
563 compar = @dist_compar;
564 end
565 for l = 1:np
566 for k = 1:nc
567 mxlk = sort(mx{1,k,l});
568 j = 1;
569 while j < length(mxlk)-1
570 if compar(th(mxlk(j+1),k,l),th(mxlk(j),k,l),option.reso)
571 decreas = option.resofirst || ...
572 dh(mxlk(j+1),k,l)<dh(mxlk(j),k,l);
573 mxlk(j + decreas) = [];
574 else
575 j = j+1;
576 end
577 end
578 if length(mxlk)>1 && compar(th(mxlk(end),k,l),...
579 th(mxlk(end-1),k,l),...
580 option.reso)
581 decreas = not(option.resofirst) &&...
582 dh(mxlk(end),k,l)>dh(mxlk(end-1),k,l);
583 mxlk(end-decreas) = [];
584 end
585 mx{1,k,l} = mxlk;
586 end
587 end
588 end
589 if not(isnan(option.near)) % Finding a peak nearest a given prefered location
590 for l = 1:np
591 for k = 1:nc
592 mxlk = mx{1,k,l};
593 if strcmp(option.logsc,'Log')
594 [M I] = min(abs(log(th(mxlk,k,l)/option.near)));
595 else
596 [M I] = min(abs(th(mxlk,k,l)-option.near));
597 end
598 mx{1,k,l} = mxlk(I);
599 end
600 end
601 end
602 if option.delta % Peak tracking
603 tp{i}{h} = cell(1,np);
604 if interpol
605 tpp{i}{h} = cell(1,np);
606 tpv{i}{h} = cell(1,np);
607 end
608 for l = 1:np
609
610 % mxl will be the resulting track position matrix
611 % and myl the related track amplitude
612 % In the first frame, tracks can be identified to peaks.
613 mxl = mx{1,1,l}(:)-1;
614 myl = dh(mx{1,1,l}(:),k,l);
615
616 % To each peak is associated the related track ID
617 tr2 = 1:length(mx{1,1,l});
618
619 grvy = []; % The graveyard...
620
621 wait = 0;
622 if nc-1>500
623 wait = waitbar(0,['Tracking peaks...']);
624 end
625
626 for k = 1:nc-1
627 % For each successive frame...
628
629 if not(isempty(grvy))
630 old = find(grvy(:,2) == k-option.mem-1);
631 grvy(old,:) = [];
632 end
633
634 if wait && not(mod(k,100))
635 waitbar(k/(nc-1),wait);
636 end
637
638 mxk1 = mx{1,k,l}; % w^k
639 mxk2 = mx{1,k+1,l}; % w^{k+1}
640 thk1 = th(mxk1,k,l);
641 thk2 = th(mxk2,k,l);
642 myk2 = dh(mx{1,k+1,l},k,l); % amplitude
643 tr1 = tr2;
644 tr2 = NaN(1,length(mxk2));
645
646 mxl(:,k+1) = mxl(:,k);
647
648 if isempty(thk1) || isempty(thk2)
649 %% IS THIS TEST NECESSARY??
650
651 myl(:,k+1) = 0;
652 else
653 for n = 1:length(mxk1)
654 % Let's check each track.
655 tr = tr1(n); % Current track.
656
657 if not(isnan(tr))
658 % still alive...
659
660 % Step 1 in Mc Aulay & Quatieri
661 [int m] = min(abs(thk2-thk1(n)));
662 if isinf(int) || int > option.delta
663 % all w^{k+1} outside matching interval:
664 % partial becomes dead
665 mxl(tr,k+1) = mxl(tr,k);
666 myl(tr,k+1) = 0;
667 grvy = [grvy; tr k]; % added to the graveyard
668 else
669 % closest w^{k+1} is tentatively selected:
670 % candidate match
671
672 % Step 2 in Mc Aulay & Quatieri
673 [best mm] = min(abs(thk2(m)-th(mx{1,k,l})));
674 if mm == n
675 % no better match to remaining w^k:
676 % definite match
677 mxl(tr,k+1) = mxk2(m)-1;
678 myl(tr,k+1) = myk2(m);
679 tr2(m) = tr;
680 thk1(n) = -Inf; % selected w^k is eliminated from further consideration
681 thk2(m) = Inf; % selected w^{k+1} is eliminated as well
682 if not(isempty(grvy))
683 zz = find ((mxl(grvy(:,1),k) >= mxl(tr,k) & ...
684 mxl(grvy(:,1),k) <= mxl(tr,k+1)) | ...
685 (mxl(grvy(:,1),k) <= mxl(tr,k) & ...
686 mxl(grvy(:,1),k) >= mxl(tr,k+1)));
687 grvy(zz,:) = [];
688 end
689 else
690 % let's look at adjacent lower w^{k+1}...
691 [int mmm] = min(abs(thk2(1:m)-thk1(n)));
692 if int > best || ... % New condition added (Lartillot 16.4.2010)
693 isinf(int) || ... % Conditions proposed in Mc Aulay & Quatieri (all w^{k+1} below matching interval)
694 int > option.delta
695 % partial becomes dead
696 mxl(tr,k+1) = mxl(tr,k);
697 myl(tr,k+1) = 0;
698 grvy = [grvy; tr k]; % added to the graveyard
699 else
700 % definite match
701 mxl(tr,k+1) = mxk2(mmm)-1;
702 myl(tr,k+1) = myk2(mmm);
703 tr2(mmm) = tr;
704 thk1(n) = -Inf; % selected w^k is eliminated from further consideration
705 thk2(mmm) = Inf; % selected w^{k+1} is eliminated as well
706 if not(isempty(grvy))
707 zz = find ((mxl(grvy(:,1),k) >= mxl(tr,k) & ...
708 mxl(grvy(:,1),k) <= mxl(tr,k+1)) | ...
709 (mxl(grvy(:,1),k) <= mxl(tr,k) & ...
710 mxl(grvy(:,1),k) >= mxl(tr,k+1)));
711 grvy(zz,:) = [];
712 end
713 end
714 end
715 end
716 end
717 end
718 end
719
720 % Step 3 in Mc Aulay & Quatieri
721 for m = 1:length(mxk2)
722 if not(isinf(thk2(m)))
723 % unmatched w^{k+1}
724
725 if isempty(grvy)
726 int = [];
727 else
728 % Let's try to reuse a zombie from the
729 % graveyard (Lartillot).
730 [int z] = min(abs(th(mxl(grvy(:,1),k+1)+1,k,l)-thk2(m)));
731 end
732 if isempty(int) || int > option.delta ...
733 || int > min(abs(th(mxl(:,k+1)+1,k,l)-thk2(m)))
734 % No suitable zombie.
735 % birth of a new partial (Mc Aulay &
736 % Quatieri)
737 mxl = [mxl;zeros(1,k+1)];
738 tr = size(mxl,1);
739 mxl(tr,k) = mxk2(m)-1;
740 else
741 % Suitable zombie found. (Lartillot)
742 tr = grvy(z,1);
743 grvy(z,:) = [];
744 end
745 mxl(tr,k+1) = mxk2(m)-1;
746 myl(tr,k+1) = myk2(m);
747 tr2(m) = tr;
748 end
749 end
750 end
751
752 if wait
753 waitbar(1,wait);
754 close(wait);
755 drawnow
756 end
757
758 if size(mxl,1) > option.m
759 tot = sum(myl,2);
760 [tot ix] = sort(tot,'descend');
761 mxl(ix(option.m+1:end),:) = [];
762 myl(ix(option.m+1:end),:) = [];
763 end
764
765 mxl(:,not(max(myl))) = 0;
766
767 if option.shorttrackthresh
768 [myl bestrack] = max(myl,[],1);
769 mxl = mxl(bestrack + (0:size(mxl,2)-1)*size(mxl,1));
770 changes = find(not(bestrack(1:end-1) == bestrack(2:end)))+1;
771 if not(isempty(changes))
772 lengths = diff([1 changes nc+1]);
773 shorts = find(lengths < option.shorttrackthresh);
774 for k = 1:length(shorts)
775 if shorts(k) == 1
776 k1 = 1;
777 else
778 k1 = changes(shorts(k)-1);
779 end
780 k2 = k1 + lengths(shorts(k)) -1;
781 myl(1,k1:k2) = 0;
782 mxl(1,k1:k2) = 0;
783 end
784 end
785 end
786
787 tp{i}{h}{l} = mxl;
788 tv{i}{h}{l} = myl;
789
790 if interpol
791 tpv{i}{h}{l} = zeros(size(mxl));
792 tpp{i}{h}{l} = zeros(size(mxl));
793 for k = 1:size(mxl,2)
794 for j = 1:size(mxl,1)
795 mj = mxl(j,k);
796 if mj>2 && mj<size(dh3,1)-1
797 % More precise peak position
798 y0 = dh3(mj,k,l);
799 ym = dh3(mj-1,k,l);
800 yp = dh3(mj+1,k,l);
801 p = (yp-ym)/(2*(2*y0-yp-ym));
802 tpv{i}{h}{l}(j,k) = y0 - 0.25*(ym-yp)*p;
803 if p >= 0
804 tpp{i}{h}{l}(j,k) = (1-p)*th(mj,k,l)+p*th(mj+1,k,l);
805 elseif p < 0
806 tpp{i}{h}{l}(j,k) = (1+p)*th(mj,k,l)-p*th(mj-1,k,l);
807 end
808 elseif mj
809 tpv{i}{h}{l}(j,k) = dh3(mj,k,l);
810 tpp{i}{h}{l}(j,k) = th(mj,k,l);
811 end
812 end
813 end
814 end
815 end
816 end
817 if isa(x,'mirsimatrix') && option.graph
818 % Finding the best branch inside a graph constructed out of a
819 % similarity matrix
820 g{i}{h} = cell(1,nc,np);
821 % Branch info related to each peak
822 br{i}{h} = {};
823 % Info related to each branch
824 scog{i}{h} = cell(1,nc,np);
825 % Score related to each peak
826 scob{i}{h} = [];
827 % Score related to each branch
828 for l = 1:np
829 wait = waitbar(0,['Creating peaks graph...']);
830 for k = 1:nc
831 g{i}{h}{1,k,l} = cell(size(mx{1,k,l}));
832 scog{i}{h}{1,k,l} = zeros(size(mx{1,k,l}));
833 if wait && not(mod(k,50))
834 waitbar(k/(nc-1),wait);
835 end
836 mxk = mx{1,k,l}; % Peaks in current frame
837 for j = k-1:-1:max(1,k-100) % Recent frames
838 mxj = mx{1,j,l}; % Peaks in one recent frame
839 for kk = 1:length(mxk)
840 mxkk = mxk(kk); % For each of current peaks
841 for jj = 1:length(mxj)
842 mxjj = mxj(jj); % For each of recent peaks
843 sco = k-j - abs(mxkk-mxjj);
844 % Crossprogression from recent to
845 % current peak
846 if sco >= 0
847 % Negative crossprogression excluded
848 dist = 0;
849 % The distance between recent and
850 % current peak is the sum of all the
851 % simatrix values when joining the two
852 % peaks with a straight line.
853 for m = j:k
854 % Each point in that straight line.
855 mxm = mxjj + (mxkk-mxjj)*(m-j)/(k-j);
856 if mxm == floor(mxm)
857 dist = dist + 1-dh(mxm,m,l);
858 else
859 dhm0 = dh(floor(mxm),m,l);
860 dhm1 = dh(ceil(mxm),m,l);
861 dist = dist + 1-...
862 (dhm0 + ...
863 (dhm1-dhm0)*(mxm-floor(mxm)));
864 end
865 if dist > option.graph
866 break
867 end
868 end
869 if dist < option.graph
870 % If the distance between recent
871 % and current peak is not too high,
872 % a new edge is formed between the
873 % peaks, and added to the graph.
874 gj = g{i}{h}{1,j,l}{jj};
875 % Branch information associated
876 % with recent peak
877 gk = g{i}{h}{1,k,l}{kk};
878 % Branch information associated
879 % with current peak
880 if isempty(gk) || ...
881 sco > scog{i}{h}{1,k,l}(kk)
882 % Current peak branch to be updated
883 if isempty(gj)
884 % New branch starting
885 % from scratch
886 newsco = sco;
887 scob{i}{h}(end+1) = newsco;
888 bid = length(scob{i}{h});
889 g{i}{h}{1,j,l}{jj} = ...
890 [k kk bid newsco];
891 br{i}{h}{bid} = [j jj;k kk];
892 else
893 newsco = scog{i}{h}{1,j,l}(jj)+sco;
894 if length(gj) == 1
895 % Recent peak not
896 % associated with other
897 % branch
898 % -> Branch extension
899 bid = gj;
900 g{i}{h}{1,j,l}{jj} = ...
901 [k kk bid newsco];
902 br{i}{h}{bid}(end+1,:) = [k kk];
903 else
904 % Recent peak already
905 % associated with other
906 % branch
907 % -> Branch fusion
908 bid = length(scob{i}{h})+1;
909 g{i}{h}{1,j,l}{jj} = ...
910 [k kk bid newsco; gj];
911 other = br{i}{h}{gj(1,3)};
912 % Other branch
913 % info
914 % Let's copy its
915 % prefix to the new
916 % branch:
917 other(other(:,1)>j,:) = [];
918 br{i}{h}{bid} = [other;k kk];
919 end
920 scob{i}{h}(bid) = newsco;
921 end
922 g{i}{h}{1,k,l}{kk} = bid;
923 % New peak associated with
924 % branch
925 scog{i}{h}{1,k,l}(kk) = newsco;
926 end
927 end
928 end
929 end
930 end
931 end
932 end
933 [scob{i}{h} IX] = sort(scob{i}{h},'descend');
934 % Branch are ordered from best score to lowest
935 br{i}{h} = br{i}{h}(IX);
936 if wait
937 waitbar(1,wait);
938 close(wait);
939 drawnow
940 end
941 end
942 end
943 for l = 1:np % Orders the peaks and select the best ones
944 for k = 1:nc
945 mxk = mx{1,k,l};
946 if length(mxk) > option.m
947 [unused,idx] = sort(dh(mxk,k,l),'descend');
948 idx = idx(1:option.m);
949 elseif strcmpi(option.order,'Amplitude')
950 [unused,idx] = sort(dh(mxk,k,l),'descend');
951 else
952 idx = 1:length(dh(mxk,k,l));
953 end
954 if strcmpi(option.order,'Abscissa')
955 mx{1,k,l} = sort(mxk(idx));
956 elseif strcmpi(option.order,'Amplitude')
957 mx{1,k,l} = mxk(idx);
958 end
959 end
960 end
961 if option.extract % Extracting the positive part of the curve containing the peaks
962 if isa(x,'mirtemporal')
963 filn = floor(sr{i}/25);
964 else
965 filn = 10;
966 end
967 if filn>1 && size(dh3,1)>5
968 filn = min(filn,floor(size(dh3,1)/3));
969 fild = filtfilt(ones(1,filn)/2,1,dh3(2:end-1,:,:))/filn/2;
970 else
971 fild = dh3(2:end-1,:,:);
972 end
973 fild = [zeros(1,size(fild,2),size(fild,3));diff(fild)];
974 for l = 1:np
975 for k = 1:nc
976 idx = 1:size(dh,1);
977 mxlk = sort(mx{1,k,l}-1);
978 for j = 1:length(mxlk)
979
980 if fild(mxlk(j),k,l) < 0
981 bef0 = find(fild(1:mxlk(j)-1,k,l)>=0);
982 if isempty(bef0)
983 bef0 = [];
984 end
985 else
986 bef0 = mxlk(j)-1;
987 end
988 if isempty(bef0)
989 bef = 0;
990 else
991 bef = find(fild(1:bef0(end),k,l)<=0);
992 if isempty(bef)
993 bef = 0;
994 end
995 end
996 if j>1 && bef(end)<aft(1)+2
997 idx(mxlk(j-1):mxlk(j)) = 0;
998 [unused btw] = min(dh3(mxlk(j-1)+1:mxlk(j)+1,k,l));
999 btw = btw+mxlk(j-1);
1000 idx(btw-2:btw+2) = btw-2:btw+2;
1001 bef = btw+2;
1002 end
1003
1004 if fild(mxlk(j),k,l) > 0
1005 aft0 = find(fild(mxlk(j)+1:end,k,l)<=0)+mxlk(j);
1006 if isempty(aft0)
1007 aft0 = [];
1008 end
1009 else
1010 aft0 = mxlk(j)+1;
1011 end
1012 if isempty(aft0)
1013 aft = size(d{i}{h},1)+1;
1014 else
1015 aft = find(fild(aft0(1):end,k,l)>=0)+mxlk(j);
1016 if isempty(aft)
1017 aft = size(d{i}{h},1)+1;
1018 end
1019 end
1020
1021 idx(bef(end)+3:aft(1)-3) = 0;
1022 end
1023 idx = idx(find(idx));
1024 dh3(idx,k,l) = NaN;
1025 end
1026 end
1027 end
1028 if option.vall
1029 dh3 = -dh3;
1030 end
1031 mmx = cell(1,nc,np);
1032 mmy = cell(1,nc,np);
1033 mmv = cell(1,nc,np);
1034 for l = 1:np
1035 for k = 1:nc
1036 mmx{1,k,l} = mod(mx{1,k,l}(:,:,1),nl0)-1;
1037 mmy{1,k,l} = ceil(mx{1,k,l}/nl0);
1038 mmv{1,k,l} = dh3(mx{1,k,l}(:,:,1),k,l);
1039 end
1040 end
1041 pp{i}{h} = mmx;
1042 pm{i}{h} = mmy;
1043 pv{i}{h} = mmv;
1044 if not(interpol)
1045 ppp{i}{h} = {};
1046 ppv{i}{h} = {};
1047 else % Interpolate to find the more exact peak positions
1048 pih = cell(1,nc,np);
1049 vih = cell(1,nc,np);
1050 for l = 1:np
1051 for k = 1:nc
1052 mxlk = mx{1,k,l};
1053 vih{1,k,l} = zeros(length(mxlk),1);
1054 pih{1,k,l} = zeros(length(mxlk),1);
1055 for j = 1:length(mxlk)
1056 mj = mxlk(j); % Current values
1057 if strcmpi(option.interpol,'quadratic')
1058 if mj>2 && mj<length(dh3)-1
1059 % More precise peak position
1060 y0 = dh3(mj,k,l);
1061 ym = dh3(mj-1,k,l);
1062 yp = dh3(mj+1,k,l);
1063 p = (yp-ym)/(2*(2*y0-yp-ym));
1064 vih{1,k,l}(j) = y0 - 0.25*(ym-yp)*p;
1065 if p >= 0
1066 pih{1,k,l}(j) = (1-p)*th(mj,k,l)+p*th(mj+1,k,l);
1067 elseif p < 0
1068 pih{1,k,l}(j) = (1+p)*th(mj,k,l)-p*th(mj-1,k,l);
1069 end
1070 else
1071 vih{1,k,l}(j) = dh3(mj,k,l);
1072 pih{1,k,l}(j) = th(mj,k,l);
1073 end
1074 end
1075 end
1076 end
1077 end
1078 ppp{i}{h} = pih;
1079 ppv{i}{h} = vih;
1080 end
1081 if not(iscell(d{i})) % for chromagram
1082 d{i} = dh3(2:end-1,:,:,:);
1083 else
1084 if cha == 1
1085 d{i}{h} = zeros(1,size(dh3,2),size(dh3,1)-2);
1086 for k = 1:size(dh3,2)
1087 d{i}{h}(1,k,:) = dh3(2:end-1,k);
1088 end
1089 else
1090 d{i}{h} = dh3(2:end-1,:,:,:);
1091 end
1092 end
1093 if option.only
1094 dih = zeros(size(d{i}{h}));
1095 for l = 1:np
1096 for k = 1:nc
1097 dih(pp{i}{h}{1,k,l},k,l) = ...
1098 d{i}{h}(pp{i}{h}{1,k,l},k,l);
1099 end
1100 end
1101 d{i}{h} = dih;
1102 end
1103 end
1104 end
1105 p = set(x,'PeakPos',pp,'PeakVal',pv,'PeakMode',pm);
1106 if interpol
1107 p = set(p,'PeakPrecisePos',ppp,'PeakPreciseVal',ppv);
1108 end
1109 if option.extract
1110 p = set(p,'Data',d);
1111 end
1112 empty = cell(1,length(d));
1113 if option.only
1114 p = set(p,'Data',d,'PeakPos',empty,'PeakVal',empty,'PeakMode',empty);
1115 end
1116 if option.delta
1117 p = set(p,'TrackPos',tp,'TrackVal',tv);
1118 if interpol
1119 p = set(p,'TrackPrecisePos',tpp,'TrackPreciseVal',tpv);
1120 end
1121 end
1122 if isa(x,'mirsimatrix') && option.graph
1123 p = set(p,'Graph',g,'Branch',br);
1124 end
1125
1126
1127 function y = semitone_compar(p1,p2,thres)
1128 y = p1/p2 < 2^(1/12);
1129
1130
1131 function y = dist_compar(p1,p2,thres)
1132 y = p1-p2 < thres;