Mercurial > hg > camir-aes2014

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