Daniel@0: function varargout = mirpeaks(orig,varargin)
Daniel@0: %   p = mirpeaks(x) detect peaks in x.
Daniel@0: %   Optional argument:
Daniel@0: %       mirpeaks(...,'Total',m): only the m highest peaks are selected.
Daniel@0: %           If m=Inf, no limitation of number of peaks.
Daniel@0: %               Default value: m = Inf
Daniel@0: %       mirpeaks(...,'Order',o): specifies the ordering of the peaks.
Daniel@0: %           Possible values for o:
Daniel@0: %               'Amplitude': orders the peaks from highest to lowest
Daniel@0: %                   (Default choice.)
Daniel@0: %               'Abscissa': orders the peaks along the abscissa axis.
Daniel@0: %       mirpeaks(...,'Contrast',cthr): a threshold value. A given local
Daniel@0: %           maximum will be considered as a peak if its distance with the
Daniel@0: %           previous and successive local minima (if any) is higher than 
Daniel@0: %           this threshold. This distance is expressed with respect to the
Daniel@0: %           total amplitude of x: a distance of 1, for instance, is
Daniel@0: %           equivalent to the distance between the maximum and the minimum
Daniel@0: %           of x.
Daniel@0: %               Default value: cthr = 0.1
Daniel@0: %           mirpeaks(...,'SelectFirst',fthr): If the 'Contrast' selection has
Daniel@0: %               been chosen, this additional option specifies that when one
Daniel@0: %               peak has to be chosen out of two candidates, and if the 
Daniel@0: %               difference of their amplitude is below the threshold fthr,
Daniel@0: %               then the most ancien one is selected.
Daniel@0: %               Option toggled off by default:
Daniel@0: %               Default value if toggled on: fthr = cthr/2
Daniel@0: %       mirpeaks(...,'Threshold',thr): a threshold value.
Daniel@0: %           A given local maximum will be considered as a peak if its
Daniel@0: %               normalized amplitude is higher than this threshold. 
Daniel@0: %           A given local minimum will be considered as a valley if its
Daniel@0: %               normalized amplitude is lower than this threshold. 
Daniel@0: %           The normalized amplitude can have value between 0 (the minimum 
Daniel@0: %               of the signal in each frame) and 1 (the maximum in each 
Daniel@0: %               frame).
Daniel@0: %           Default value: thr=0 for peaks thr = 1 for valleys
Daniel@0: %       mirpeaks(...,'Interpol',i): estimates more precisely the peak
Daniel@0: %           position and amplitude using interpolation. Performed only on
Daniel@0: %           data with numerical abscissae axis.
Daniel@0: %           Possible value for i:
Daniel@0: %               '', 'no', 'off', 0: no interpolation
Daniel@0: %               'Quadratic': quadratic interpolation. (default value).
Daniel@0: %       mirpeaks(...,'Valleys'): detect valleys instead of peaks.
Daniel@0: %       mirpeaks(...,'Reso',r): removes peaks whose distance to one or
Daniel@0: %           several higher peaks is lower than a given threshold.
Daniel@0: %           Possible value for the threshold r:
Daniel@0: %               'SemiTone': ratio between the two peak positions equal to
Daniel@0: %                   2^(1/12)
Daniel@0: %               a numerical value : difference between the two peak
Daniel@0: %                   positions equal to that value
Daniel@0: %           When two peaks are distant by an interval lower than the
Daniel@0: %               resolution, the highest of them is selected by default.
Daniel@0: %           mirpeaks(...,'Reso',r,'First') specifies on the contrary that
Daniel@0: %               the first of them is selected by default.
Daniel@0: %       mirpeaks(...,'Nearest',t,s): takes the peak nearest a given abscisse
Daniel@0: %           values t. The distance is computed either on a linear scale
Daniel@0: %           (s = 'Lin') or logarithmic scale (s = 'Log'). In this case,
Daniel@0: %           only one peak is extracted.
Daniel@0: %       mirpeaks(...,'Pref',c,std): indicates a region of preference for
Daniel@0: %           the peak picking, centered on the abscisse value c, with a
Daniel@0: %           standard deviation of std.
Daniel@0: %       mirpeaks(...,'NoBegin'): does not consider the first sample as a
Daniel@0: %           possible peak candidate.
Daniel@0: %       mirpeaks(...,'NoEnd'): does not consider the last sample as a possible
Daniel@0: %           peak candidate.
Daniel@0: %       mirpeaks(...,'Normalize',n): specifies whether frames are
Daniel@0: %           normalized globally or individually.
Daniel@0: %           Possible value for n:
Daniel@0: %               'Global': normalizes the whole frames altogether from 0 to 
Daniel@0: %                   1 (default choice).
Daniel@0: %               'Local': normalizes each frame from 0 to 1.
Daniel@0: %       mirpeaks(...,'Extract'): extracts from the initial curves all the 
Daniel@0: %           positive continuous segments (or "curve portions") where peaks
Daniel@0: %           are located.
Daniel@0: %       mirpeaks(...,'Only'): keeps from the original curve only the data
Daniel@0: %           corresponding to the peaks, and zeroes the remaining data.
Daniel@0: %       mirpeaks(...,'Track',t): tracks temporal continuities of peaks. If
Daniel@0: %           a value t is specified, the variation between successive peaks
Daniel@0: %           is tolerated up to t samples.
Daniel@0: %       mirpeaks(...,'CollapseTrack',ct): collapses tracks into one single
Daniel@0: %           track, and remove small track transitions, of length shorter
Daniel@0: %           than ct samples. Default value: ct = 7
Daniel@0: 
Daniel@0:         m.key = 'Total';
Daniel@0:         m.type = 'Integer';
Daniel@0:         m.default = Inf;
Daniel@0:     option.m = m;
Daniel@0:         
Daniel@0:         nobegin.key = 'NoBegin';
Daniel@0:         nobegin.type = 'Boolean';
Daniel@0:         nobegin.default = 0;
Daniel@0:     option.nobegin = nobegin;
Daniel@0:         
Daniel@0:         noend.key = 'NoEnd';
Daniel@0:         noend.type = 'Boolean';
Daniel@0:         noend.default = 0;
Daniel@0:     option.noend = noend;
Daniel@0:         
Daniel@0:         order.key = 'Order';
Daniel@0:         order.type = 'String';
Daniel@0:         order.choice = {'Amplitude','Abscissa'};
Daniel@0:         order.default = 'Amplitude';
Daniel@0:     option.order = order;
Daniel@0:     
Daniel@0:         chro.key = 'Chrono'; % obsolete, corresponds to 'Order','Abscissa'
Daniel@0:         chro.type = 'Boolean';
Daniel@0:         chro.default = 0;
Daniel@0:     option.chro = chro;
Daniel@0:     
Daniel@0:         ranked.key = 'Ranked'; % obsolete, corresponds to 'Order','Amplitude'
Daniel@0:         ranked.type = 'Boolean';
Daniel@0:         ranked.default = 0;
Daniel@0:     option.ranked = ranked;
Daniel@0:         
Daniel@0:         vall.key = 'Valleys';
Daniel@0:         vall.type = 'Boolean';
Daniel@0:         vall.default = 0;
Daniel@0:     option.vall = vall;
Daniel@0:     
Daniel@0:         cthr.key = 'Contrast';
Daniel@0:         cthr.type = 'Integer';
Daniel@0:         cthr.default = .1;
Daniel@0:     option.cthr = cthr;
Daniel@0:         
Daniel@0:         first.key = 'SelectFirst';
Daniel@0:         first.type = 'Integer';
Daniel@0:         first.default = 0;
Daniel@0:         first.keydefault = NaN;
Daniel@0:     option.first = first;
Daniel@0:     
Daniel@0:         thr.key = 'Threshold';
Daniel@0:         thr.type = 'Integer';
Daniel@0:         thr.default = NaN;
Daniel@0:     option.thr = thr;
Daniel@0:             
Daniel@0:         smthr.key = 'MatrixThreshold'; % to be documented in version 1.3
Daniel@0:         smthr.type = 'Integer';
Daniel@0:         smthr.default = NaN;
Daniel@0:     option.smthr = smthr;
Daniel@0:     
Daniel@0:         graph.key = 'Graph';
Daniel@0:         graph.type = 'Integer';
Daniel@0:         graph.default = 0;
Daniel@0:         graph.keydefault = .25;
Daniel@0:     option.graph = graph;
Daniel@0:         
Daniel@0:         interpol.key = 'Interpol';
Daniel@0:         interpol.type = 'String';
Daniel@0:         interpol.default = 'Quadratic';
Daniel@0:         interpol.keydefault = 'Quadratic';
Daniel@0:     option.interpol = interpol;
Daniel@0:     
Daniel@0:         reso.key = 'Reso';
Daniel@0:         %reso.type = 'String';
Daniel@0:         %reso.choice = {0,'SemiTone'};
Daniel@0:         reso.default = 0;
Daniel@0:     option.reso = reso;
Daniel@0:     
Daniel@0:         resofirst.key = 'First';
Daniel@0:         resofirst.type = 'Boolean';
Daniel@0:         resofirst.default = 0;
Daniel@0:     option.resofirst = resofirst;
Daniel@0:         
Daniel@0:         c.key = 'Pref';
Daniel@0:         c.type = 'Integer';
Daniel@0:         c.number = 2;
Daniel@0:         c.default = [0 0];
Daniel@0:     option.c = c;
Daniel@0:         
Daniel@0:         near.key = 'Nearest';
Daniel@0:         near.type = 'Integer';
Daniel@0:         near.default = NaN;
Daniel@0:     option.near = near;
Daniel@0:         
Daniel@0:         logsc.type = 'String';
Daniel@0:         logsc.choice = {'Lin','Log',0};
Daniel@0:         logsc.default = 'Lin';
Daniel@0:     option.logsc = logsc;
Daniel@0:         
Daniel@0:         normal.key = 'Normalize';
Daniel@0:         normal.type = 'String';
Daniel@0:         normal.choice = {'Local','Global'};
Daniel@0:         normal.default = 'Global';
Daniel@0:     option.normal = normal;
Daniel@0: 
Daniel@0:         extract.key = 'Extract';
Daniel@0:         extract.type = 'Boolean';
Daniel@0:         extract.default = 0;
Daniel@0:     option.extract = extract;
Daniel@0:     
Daniel@0:         only.key = 'Only';
Daniel@0:         only.type = 'Boolean';
Daniel@0:         only.default = 0;
Daniel@0:     option.only = only;
Daniel@0: 
Daniel@0:         delta.key = 'Track';
Daniel@0:         delta.type = 'Integer';
Daniel@0:         delta.default = 0;
Daniel@0:         delta.keydefault = Inf;
Daniel@0:     option.delta = delta;
Daniel@0:     
Daniel@0:         mem.key = 'TrackMem';
Daniel@0:         mem.type = 'Integer';
Daniel@0:         mem.default = Inf;
Daniel@0:     option.mem = mem;
Daniel@0: 
Daniel@0:         shorttrackthresh.key = 'CollapseTracks';
Daniel@0:         shorttrackthresh.type = 'Integer';
Daniel@0:         shorttrackthresh.default = 0;
Daniel@0:         shorttrackthresh.keydefault = 7;
Daniel@0:     option.shorttrackthresh = shorttrackthresh;
Daniel@0: 
Daniel@0:         scan.key = 'ScanForward'; % specific to mironsets(..., 'Klapuri99')
Daniel@0:         scan.default = [];
Daniel@0:     option.scan = scan;
Daniel@0:     
Daniel@0: specif.option = option;
Daniel@0: 
Daniel@0: varargout = mirfunction(@mirpeaks,orig,varargin,nargout,specif,@init,@main);
Daniel@0: 
Daniel@0: 
Daniel@0: function [x type] = init(x,option)
Daniel@0: type = mirtype(x);
Daniel@0: 
Daniel@0: 
Daniel@0: function p = main(x,option,postoption)
Daniel@0: if iscell(x)
Daniel@0:     x = x{1};
Daniel@0: end
Daniel@0: if option.chro
Daniel@0:     option.order = 'Abscissa';
Daniel@0: elseif option.ranked
Daniel@0:     option.order = 'Amplitude';
Daniel@0: end
Daniel@0: if not(isnan(option.near)) && option.m == 1
Daniel@0:     option.m = Inf;
Daniel@0: end
Daniel@0: x = purgedata(x);
Daniel@0: if option.m <= 0
Daniel@0:     p = x;
Daniel@0:     return
Daniel@0: end
Daniel@0: d = get(x,'Data');
Daniel@0: sr = get(x,'Sampling');
Daniel@0: cha = 0; % Indicates when it is possible to represent as a curve along the 
Daniel@0:          % Z-axis (channels) instead of the X-axis (initial abscissa).
Daniel@0: if isnan(option.first)
Daniel@0:     option.first = option.cthr / 2;
Daniel@0: end
Daniel@0: if isa(x,'mirscalar')
Daniel@0:     t = get(x,'FramePos');
Daniel@0:     for i = 1:length(d)
Daniel@0:         for j = 1:length(d{i})
Daniel@0:             d{i}{j} = d{i}{j}';
Daniel@0:             if size(t{i},1) == 1
Daniel@0:                 t{i}{j} = t{i}{j}(1,:,:)';
Daniel@0:             else
Daniel@0:                 t{i}{j} = (t{i}{j}(1,:,:)+t{i}{j}(2,:,:))'/2;
Daniel@0:             end
Daniel@0:         end
Daniel@0:     end
Daniel@0: elseif isa(x,'mirsimatrix')
Daniel@0:     t = get(x,'FramePos');
Daniel@0:     for i = 1:length(t)
Daniel@0:         for j = 1:length(t{i})
Daniel@0:             t{i}{j} = repmat((t{i}{j}(1,:,:)+t{i}{j}(2,:,:))'/2,...
Daniel@0:                                     [1 size(d{i}{j},2) 1]);
Daniel@0:         end
Daniel@0:     end
Daniel@0: elseif isa(x,'mirhisto')
Daniel@0:     error('ERROR IN MIRPEAKS: peaks of histogram not considered yet.');
Daniel@0: else
Daniel@0:     for i = 1:length(d)
Daniel@0:         for j = 1:length(d{i})
Daniel@0:             if iscell(d{i})
Daniel@0:                 dij = d{i}{j};
Daniel@0:                 if ~cha && j == 1 && size(dij,3) > 1 && size(dij,1) == 1
Daniel@0:                     cha = 1;
Daniel@0:                 end
Daniel@0:                 if cha && j > 1 && size(dij,1) > 1
Daniel@0:                     cha = -1;
Daniel@0:                 end
Daniel@0:             end
Daniel@0:         end
Daniel@0:         for j = 1:length(d{i})
Daniel@0:             if iscell(d{i})
Daniel@0:                 dij = d{i}{j};
Daniel@0:                 if cha == 1
Daniel@0:                     if iscell(dij)
Daniel@0:                         for k = 1:size(dij,2)
Daniel@0:                             d{i}{j}{k} = reshape(dij{k},size(dij{k},2),size(dij{k},3));
Daniel@0:                             d{i}{j}{k} = d{i}{j}{k}';
Daniel@0:                         end
Daniel@0:                     else
Daniel@0:                         d{i}{j} = reshape(dij,size(dij,2),size(dij,3));
Daniel@0:                         d{i}{j} = d{i}{j}';
Daniel@0:                     end
Daniel@0:                 end
Daniel@0:             end
Daniel@0:         end
Daniel@0:     end
Daniel@0:     if cha == 1
Daniel@0:         t = get(x,'Channels');
Daniel@0:     else
Daniel@0:         t = get(x,'Pos');
Daniel@0:     end
Daniel@0: end
Daniel@0: pp = cell(1,length(d));
Daniel@0: pv = cell(1,length(d));
Daniel@0: pm = cell(1,length(d));
Daniel@0: ppp = cell(1,length(d));
Daniel@0: ppv = cell(1,length(d));
Daniel@0: tp = cell(1,length(d));
Daniel@0: tv = cell(1,length(d));
Daniel@0: 
Daniel@0: if isnan(option.thr)
Daniel@0:     option.thr = 0;
Daniel@0: else
Daniel@0:     if option.vall
Daniel@0:         option.thr = 1-option.thr;
Daniel@0:     end
Daniel@0: end
Daniel@0: %if isnan(option.lthr)
Daniel@0: %    option.lthr = 1;
Daniel@0: %else
Daniel@0: %    if option.vall
Daniel@0: %        option.lthr = 1-option.lthr;
Daniel@0: %    end
Daniel@0: %end
Daniel@0: if isnan(option.smthr)
Daniel@0:     option.smthr = option.thr - .2;
Daniel@0: end
Daniel@0: 
Daniel@0: if not(isempty(option.scan))
Daniel@0:     pscan = get(option.scan,'PeakPos');
Daniel@0: end
Daniel@0: 
Daniel@0: interpol = get(x,'Interpolable') && not(isempty(option.interpol)) && ...
Daniel@0:                 ((isnumeric(option.interpol) && option.interpol) || ...
Daniel@0:                  (ischar(option.interpol) && not(strcmpi(option.interpol,'No')) && not(strcmpi(option.interpol,'Off'))));
Daniel@0:                 
Daniel@0: for i = 1:length(d) % For each audio file,...
Daniel@0:     di = d{i};
Daniel@0:     if cha == 1
Daniel@0:         ti = t; %sure ?
Daniel@0:     else
Daniel@0:         ti = t{i};
Daniel@0:     end
Daniel@0:     if not(iscell(di))
Daniel@0:         di = {di};
Daniel@0:         if isa(x,'mirchromagram') && not(cha)
Daniel@0:             ti = {ti};
Daniel@0:         end
Daniel@0:     end
Daniel@0:     for h = 1:length(di)    % For each segment,...
Daniel@0:         dh0 = di{h};
Daniel@0:         if option.vall
Daniel@0:             dh0 = -dh0;
Daniel@0:         end
Daniel@0:         dh2 = dh0;
Daniel@0:         [nl0 nc np nd0] = size(dh0);
Daniel@0:         if cha == 1
Daniel@0:             if iscell(ti)
Daniel@0:                 %% problem here!!!
Daniel@0:                 ti = ti{i}; %%%%%it seems that sometimes we need to use instead ti{i}(:);
Daniel@0:             end
Daniel@0:             th = repmat(ti,[1,nc,np,nd0]);
Daniel@0:         else
Daniel@0:             th = ti{h};
Daniel@0:             if iscell(th)  % Non-numerical abscissae are transformed into numerical ones. 
Daniel@0:                 th = repmat((1:size(th,1))',[1,nc,np]);
Daniel@0:             else
Daniel@0:                 if size(th,3)<np
Daniel@0:                     th = repmat(th,[1,1,np]);
Daniel@0:                 end
Daniel@0:             end
Daniel@0:         end
Daniel@0:         if option.c    % If a prefered region is specified, the data is amplified accordingly
Daniel@0:             dh0 = dh0.*exp(-(th-option.c(1)).^2/2/option.c(2)^2)...
Daniel@0:                                     /option.c(2)/sqrt(2*pi)/2;
Daniel@0:         end
Daniel@0: 
Daniel@0:         % Now the data is rescaled. the minimum is set to 0
Daniel@0:         % and the maximum to 1.
Daniel@0:         state = warning('query','MATLAB:divideByZero');
Daniel@0:         warning('off','MATLAB:divideByZero');
Daniel@0:         
Daniel@0:         % Let's first normalize all frames globally:
Daniel@0:         dh0 = (dh0-repmat(min(min(min(dh0,[],1),[],2),[],4),[nl0 nc 1 nd0]))./...
Daniel@0:             repmat(max(max(max(dh0,[],1),[],2),[],4)...
Daniel@0:                   -min(min(min(dh0,[],1),[],2),[],4),[nl0 nc 1 nd0]);
Daniel@0:         for l = 1:nd0 
Daniel@0:             [unused lowc] = find(max(dh0(:,:,:,l))<option.thr);
Daniel@0:             dh0(:,lowc,1,l) = 0;
Daniel@0:         end
Daniel@0:         
Daniel@0:         if strcmpi(option.normal,'Local')
Daniel@0:             % Normalizing each frame separately:
Daniel@0:             dh0 = (dh0-repmat(min(min(dh0,[],1),[],4),[nl0 1 1 nd0]))./... 
Daniel@0:                 repmat(max(max(dh0,[],1),[],4)...
Daniel@0:                       -min(min(dh0,[],1),[],4),[nl0 1 1 nd0]);
Daniel@0:         end
Daniel@0:         warning(state.state,'MATLAB:divideByZero');
Daniel@0: 
Daniel@0:         if option.nobegin
Daniel@0:             dh0 = [Inf(1,nc,np,nd0);dh0];   
Daniel@0:             % This infinite value at the beginning
Daniel@0:             % prevents the selection of the first sample of data
Daniel@0:             dh2 = [Inf(1,nc,np,nd0);dh2];
Daniel@0:             th = [NaN(1,nc,np,1);th];
Daniel@0:         else
Daniel@0:             dh0 = [-Inf(1,nc,np,nd0);dh0];
Daniel@0:             % This infinite negative value at the beginning
Daniel@0:             % ensures the selection of the first sample of data
Daniel@0:             dh2 = [-Inf(1,nc,np,nd0);dh2];
Daniel@0:             th = [NaN(1,nc,np,1);th];
Daniel@0:         end
Daniel@0:         if option.noend
Daniel@0:             dh0 = [dh0;Inf(1,nc,np,nd0)];
Daniel@0:             % idem for the last sample of data
Daniel@0:             dh2 = [dh2;Inf(1,nc,np,nd0)];
Daniel@0:             th = [th;NaN(1,nc,np,1)];
Daniel@0:         else
Daniel@0:             dh0 = [dh0;-Inf(1,nc,np,nd0)];
Daniel@0:             dh2 = [dh2;-Inf(1,nc,np,nd0)];
Daniel@0:             th = [th;NaN(1,nc,np,1)];
Daniel@0:         end
Daniel@0:         nl0 = nl0+2;
Daniel@0: 
Daniel@0:         % Rearrange the 4th dimension (if used) into the 1st one.
Daniel@0:         nl = nl0*nd0;
Daniel@0:         dh = zeros(nl,nc,np);
Daniel@0:         dh3 = zeros(nl,nc,np);
Daniel@0:         th2 = zeros(nl,nc,np);
Daniel@0:         for l = 1:nd0
Daniel@0:             dh((l-1)*nl0+(1:nl0)',:,:) = dh0(:,:,:,l);
Daniel@0:             dh3((l-1)*nl0+(1:nl0)',:,:) = dh2(:,:,:,l);
Daniel@0:             th2((l-1)*nl0+(1:nl0)',:,:) = th(:,:,:);
Daniel@0:         end
Daniel@0:         
Daniel@0:         th = th2; % The X-abscissa
Daniel@0: 
Daniel@0:         ddh = diff(dh);
Daniel@0:         % Let's find the local maxima
Daniel@0:         for l = 1:np
Daniel@0:             dl = dh(2:end-1,:,l);
Daniel@0:             for k = 1:nc
Daniel@0:                 dk = dl(:,k);
Daniel@0:                 mx{1,k,l} = find(and(and(dk >= option.cthr, ...
Daniel@0:                                          dk >= option.thr),...     
Daniel@0:                                          ... dk <= option.lthr)),
Daniel@0:                                      and(ddh(1:(end-1),k,l) > 0, ...
Daniel@0:                                          ddh(2:end,k,l) <= 0)))+1;
Daniel@0:             end
Daniel@0:         end
Daniel@0:         if option.cthr
Daniel@0:             for l = 1:np
Daniel@0:                 for k = 1:nc
Daniel@0:                     finalmxk = [];
Daniel@0:                     mxk = mx{1,k,l};
Daniel@0:                     if not(isempty(mxk))
Daniel@0:                         wait = 0;
Daniel@0:                         if length(mxk)>5000
Daniel@0:                             wait = waitbar(0,['Selecting peaks... (0 out of 0)']);
Daniel@0:                         end
Daniel@0:                         %if option.m < Inf
Daniel@0:                         %    [unused,idx] = sort(dh(mxk,k,l),'descend'); % The peaks are sorted in decreasing order
Daniel@0:                         %    mxk = mxk(sort(idx(1:option.m)));
Daniel@0:                         %end
Daniel@0:                         j = 1;  % Scans the peaks from begin to end.
Daniel@0:                         mxkj = mxk(j); % The current peak
Daniel@0:                         jj = j+1;
Daniel@0:                         bufmin = Inf;
Daniel@0:                         bufmax = dh(mxkj,k,l);
Daniel@0:                         oldbufmin = min(dh(1:mxk(1)-1,k,l));
Daniel@0:                         while jj <= length(mxk)
Daniel@0:                             if wait && not(mod(jj,5000))
Daniel@0:                                 waitbar(jj/length(mxk),wait,['Selecting peaks... (',num2str(length(finalmxk)),' out of ',num2str(jj),')']);
Daniel@0:                             end
Daniel@0:                             bufmin = min(bufmin, ...
Daniel@0:                                 min(dh(mxk(jj-1)+1:mxk(jj)-1,k,l)));
Daniel@0:                             if bufmax - bufmin < option.cthr
Daniel@0:                                 % There is no contrastive notch
Daniel@0:                                 if dh(mxk(jj),k,l) > bufmax && ...
Daniel@0:                                         (dh(mxk(jj),k,l) - bufmax > option.first ...
Daniel@0:                                         || (bufmax - oldbufmin < option.cthr))
Daniel@0:                                     % If the new peak is significantly
Daniel@0:                                     % higher than the previous one,
Daniel@0:                                     % The peak is transfered to the new
Daniel@0:                                     % position
Daniel@0:                                     j = jj;
Daniel@0:                                     mxkj = mxk(j); % The current peak
Daniel@0:                                     bufmax = dh(mxkj,k,l);
Daniel@0:                                     oldbufmin = min(oldbufmin,bufmin);
Daniel@0:                                     bufmin = Inf;
Daniel@0:                                 elseif dh(mxk(jj),k,l) - bufmax <= option.first
Daniel@0:                                     bufmax = max(bufmax,dh(mxk(jj),k,l));
Daniel@0:                                     oldbufmin = min(oldbufmin,bufmin);
Daniel@0:                                 end
Daniel@0:                             else
Daniel@0:                                 % There is a contrastive notch
Daniel@0:                                 if bufmax - oldbufmin < option.cthr
Daniel@0:                                     % But the previous peak candidate
Daniel@0:                                     % is too weak and therefore
Daniel@0:                                     % discarded
Daniel@0:                                     oldbufmin = min(oldbufmin,bufmin);
Daniel@0:                                 else
Daniel@0:                                     % The previous peak candidate is OK
Daniel@0:                                     % and therefore stored
Daniel@0:                                     finalmxk(end+1) = mxkj;
Daniel@0:                                     oldbufmin = bufmin;
Daniel@0:                                 end
Daniel@0:                                 bufmax = dh(mxk(jj),k,l);
Daniel@0:                                 j = jj;
Daniel@0:                                 mxkj = mxk(j); % The current peak
Daniel@0:                                 bufmin = Inf;
Daniel@0:                             end
Daniel@0:                             jj = jj+1;
Daniel@0:                         end
Daniel@0:                         if bufmax - oldbufmin >= option.cthr && ...
Daniel@0:                                 bufmax - min(dh(mxk(j)+1:end,k,l)) >= option.cthr
Daniel@0:                             % The last peak candidate is OK and stored
Daniel@0:                             finalmxk(end+1) = mxk(j);
Daniel@0:                         end
Daniel@0:                         if wait
Daniel@0:                             waitbar(1,wait);
Daniel@0:                             close(wait);
Daniel@0:                             drawnow
Daniel@0:                         end
Daniel@0:                     end
Daniel@0:                     mx{1,k,l} = finalmxk;
Daniel@0:                 end
Daniel@0:             end
Daniel@0:         end
Daniel@0:         if not(isempty(option.scan)) % 'ScanForward' option, used for 'Klapuri99' in mironsets
Daniel@0:             for l = 1:np
Daniel@0:                 for k = 1:nc   
Daniel@0:                     pscankl = pscan{i}{h}{1,k,l}; % scan seed positions
Daniel@0:                     mxkl = [];
Daniel@0:                     lp = length(pscankl); % number of peaks
Daniel@0:                     for jj = 1:lp         % for each peak
Daniel@0:                         fmx = find(mx{1,k,l}>pscankl(jj),1);
Daniel@0:                             % position of the next max following the
Daniel@0:                             % current seed
Daniel@0:                         fmx = mx{1,k,l}(fmx);
Daniel@0:                         if jj<lp && (isempty(fmx) || fmx>=pscankl(jj+1))
Daniel@0:                             [unused fmx] = max(dh(pscankl(jj):...
Daniel@0:                                                   pscankl(jj+1)-1,k,l));
Daniel@0:                             fmx = fmx+pscankl(jj)-1;
Daniel@0:                         elseif jj==lp && isempty(fmx)
Daniel@0:                             [unused fmx] = max(dh(pscankl(jj):end,k,l));
Daniel@0:                             fmx = fmx+pscankl(jj)-1;
Daniel@0:                         end
Daniel@0:                         mxkl = [mxkl fmx];
Daniel@0:                     end
Daniel@0:                     mx{1,k,l} = mxkl;
Daniel@0:                 end
Daniel@0:             end
Daniel@0:         end
Daniel@0:         if not(isequal(option.reso,0)) % Removing peaks too close to higher peaks
Daniel@0:             if ischar(option.reso) && strcmpi(option.reso,'SemiTone')
Daniel@0:                 compar = @semitone_compar;
Daniel@0:             elseif isnumeric(option.reso)
Daniel@0:                 compar = @dist_compar;
Daniel@0:             end
Daniel@0:             for l = 1:np
Daniel@0:                 for k = 1:nc
Daniel@0:                     mxlk = sort(mx{1,k,l});
Daniel@0:                     j = 1;
Daniel@0:                     while j < length(mxlk)-1
Daniel@0:                         if compar(th(mxlk(j+1),k,l),th(mxlk(j),k,l),option.reso)
Daniel@0:                             decreas = option.resofirst || ...
Daniel@0:                                 dh(mxlk(j+1),k,l)<dh(mxlk(j),k,l);
Daniel@0:                             mxlk(j + decreas) = [];
Daniel@0:                         else
Daniel@0:                             j = j+1;
Daniel@0:                         end
Daniel@0:                     end
Daniel@0:                     if length(mxlk)>1 && compar(th(mxlk(end),k,l),...
Daniel@0:                                                  th(mxlk(end-1),k,l),...
Daniel@0:                                                  option.reso)
Daniel@0:                         decreas = not(option.resofirst) &&...
Daniel@0:                             dh(mxlk(end),k,l)>dh(mxlk(end-1),k,l);
Daniel@0:                         mxlk(end-decreas) = [];
Daniel@0:                     end
Daniel@0:                     mx{1,k,l} = mxlk;
Daniel@0:                 end
Daniel@0:             end
Daniel@0:         end
Daniel@0:         if not(isnan(option.near)) % Finding a peak nearest a given prefered location
Daniel@0:             for l = 1:np
Daniel@0:                 for k = 1:nc
Daniel@0:                     mxlk = mx{1,k,l};
Daniel@0:                     if strcmp(option.logsc,'Log')
Daniel@0:                         [M I] = min(abs(log(th(mxlk,k,l)/option.near)));
Daniel@0:                     else
Daniel@0:                         [M I] = min(abs(th(mxlk,k,l)-option.near));
Daniel@0:                     end
Daniel@0:                     mx{1,k,l} = mxlk(I);
Daniel@0:                 end
Daniel@0:             end
Daniel@0:         end
Daniel@0:         if option.delta % Peak tracking
Daniel@0:             tp{i}{h} = cell(1,np);
Daniel@0:             if interpol
Daniel@0:                 tpp{i}{h} = cell(1,np);
Daniel@0:                 tpv{i}{h} = cell(1,np);
Daniel@0:             end
Daniel@0:             for l = 1:np
Daniel@0:                 
Daniel@0:                 % mxl will be the resulting track position matrix
Daniel@0:                 % and myl the related track amplitude
Daniel@0:                 % In the first frame, tracks can be identified to peaks.
Daniel@0:                 mxl = mx{1,1,l}(:)-1;        
Daniel@0:                 myl = dh(mx{1,1,l}(:),k,l); 
Daniel@0:                 
Daniel@0:                 % To each peak is associated the related track ID
Daniel@0:                 tr2 = 1:length(mx{1,1,l});
Daniel@0:                 
Daniel@0:                 grvy = []; % The graveyard...
Daniel@0:                 
Daniel@0:                 wait = 0;
Daniel@0:                 if nc-1>500
Daniel@0:                     wait = waitbar(0,['Tracking peaks...']);
Daniel@0:                 end
Daniel@0: 
Daniel@0:                 for k = 1:nc-1
Daniel@0:                     % For each successive frame...
Daniel@0:                     
Daniel@0:                     if not(isempty(grvy))
Daniel@0:                         old = find(grvy(:,2) == k-option.mem-1);
Daniel@0:                         grvy(old,:) = [];
Daniel@0:                     end
Daniel@0:                     
Daniel@0:                     if wait && not(mod(k,100))
Daniel@0:                         waitbar(k/(nc-1),wait);
Daniel@0:                     end
Daniel@0: 
Daniel@0:                     mxk1 = mx{1,k,l};   % w^k
Daniel@0:                     mxk2 = mx{1,k+1,l}; % w^{k+1}
Daniel@0:                     thk1 = th(mxk1,k,l);
Daniel@0:                     thk2 = th(mxk2,k,l);
Daniel@0:                     myk2 = dh(mx{1,k+1,l},k,l); % amplitude
Daniel@0:                     tr1 = tr2;
Daniel@0:                     tr2 = NaN(1,length(mxk2));
Daniel@0:                     
Daniel@0:                     mxl(:,k+1) = mxl(:,k);
Daniel@0:                                         
Daniel@0:                     if isempty(thk1) || isempty(thk2)
Daniel@0:                         %% IS THIS TEST NECESSARY??
Daniel@0:                         
Daniel@0:                         myl(:,k+1) = 0;
Daniel@0:                     else
Daniel@0:                         for n = 1:length(mxk1)
Daniel@0:                             % Let's check each track.
Daniel@0:                             tr = tr1(n); % Current track.
Daniel@0: 
Daniel@0:                             if not(isnan(tr))
Daniel@0:                                 % still alive...
Daniel@0: 
Daniel@0:                                 % Step 1 in Mc Aulay & Quatieri
Daniel@0:                                 [int m] = min(abs(thk2-thk1(n)));
Daniel@0:                                 if isinf(int) || int > option.delta
Daniel@0:                                     % all w^{k+1} outside matching interval:
Daniel@0:                                         % partial becomes dead
Daniel@0:                                     mxl(tr,k+1) = mxl(tr,k);
Daniel@0:                                     myl(tr,k+1) = 0;
Daniel@0:                                     grvy = [grvy; tr k]; % added to the graveyard
Daniel@0:                                 else
Daniel@0:                                     % closest w^{k+1} is tentatively selected:
Daniel@0:                                         % candidate match
Daniel@0: 
Daniel@0:                                     % Step 2 in Mc Aulay & Quatieri
Daniel@0:                                     [best mm] = min(abs(thk2(m)-th(mx{1,k,l})));
Daniel@0:                                     if mm == n
Daniel@0:                                         % no better match to remaining w^k:
Daniel@0:                                             % definite match
Daniel@0:                                         mxl(tr,k+1) = mxk2(m)-1;
Daniel@0:                                         myl(tr,k+1) = myk2(m);
Daniel@0:                                         tr2(m) = tr;
Daniel@0:                                         thk1(n) = -Inf; % selected w^k is eliminated from further consideration
Daniel@0:                                         thk2(m) = Inf;  % selected w^{k+1} is eliminated as well
Daniel@0:                                         if not(isempty(grvy))
Daniel@0:                                             zz = find ((mxl(grvy(:,1),k) >= mxl(tr,k) & ...
Daniel@0:                                                         mxl(grvy(:,1),k) <= mxl(tr,k+1)) | ...
Daniel@0:                                                        (mxl(grvy(:,1),k) <= mxl(tr,k) & ...
Daniel@0:                                                         mxl(grvy(:,1),k) >= mxl(tr,k+1)));
Daniel@0:                                             grvy(zz,:) = [];
Daniel@0:                                         end
Daniel@0:                                     else
Daniel@0:                                         % let's look at adjacent lower w^{k+1}...
Daniel@0:                                         [int mmm] = min(abs(thk2(1:m)-thk1(n)));
Daniel@0:                                         if int > best || ... % New condition added (Lartillot 16.4.2010)
Daniel@0:                                                 isinf(int) || ... % Conditions proposed in Mc Aulay & Quatieri (all w^{k+1} below matching interval)
Daniel@0:                                                 int > option.delta
Daniel@0:                                                 % partial becomes dead
Daniel@0:                                             mxl(tr,k+1) = mxl(tr,k);
Daniel@0:                                             myl(tr,k+1) = 0;
Daniel@0:                                             grvy = [grvy; tr k]; % added to the graveyard
Daniel@0:                                         else
Daniel@0:                                             % definite match
Daniel@0:                                             mxl(tr,k+1) = mxk2(mmm)-1;
Daniel@0:                                             myl(tr,k+1) = myk2(mmm);
Daniel@0:                                             tr2(mmm) = tr;
Daniel@0:                                             thk1(n) = -Inf;     % selected w^k is eliminated from further consideration
Daniel@0:                                             thk2(mmm) = Inf;    % selected w^{k+1} is eliminated as well
Daniel@0:                                             if not(isempty(grvy))
Daniel@0:                                                 zz = find ((mxl(grvy(:,1),k) >= mxl(tr,k) & ...
Daniel@0:                                                             mxl(grvy(:,1),k) <= mxl(tr,k+1)) | ...
Daniel@0:                                                            (mxl(grvy(:,1),k) <= mxl(tr,k) & ...
Daniel@0:                                                             mxl(grvy(:,1),k) >= mxl(tr,k+1)));
Daniel@0:                                                 grvy(zz,:) = [];
Daniel@0:                                             end
Daniel@0:                                         end
Daniel@0:                                     end
Daniel@0:                                 end
Daniel@0:                             end
Daniel@0:                         end
Daniel@0:                     end
Daniel@0:                     
Daniel@0:                     % Step 3 in Mc Aulay & Quatieri
Daniel@0:                     for m = 1:length(mxk2)
Daniel@0:                         if not(isinf(thk2(m)))
Daniel@0:                             % unmatched w^{k+1}
Daniel@0:                             
Daniel@0:                             if isempty(grvy)
Daniel@0:                                 int = [];
Daniel@0:                             else
Daniel@0:                                 % Let's try to reuse a zombie from the
Daniel@0:                                 % graveyard (Lartillot).
Daniel@0:                                 [int z] = min(abs(th(mxl(grvy(:,1),k+1)+1,k,l)-thk2(m)));
Daniel@0:                             end
Daniel@0:                             if isempty(int) || int > option.delta ...
Daniel@0:                                     || int > min(abs(th(mxl(:,k+1)+1,k,l)-thk2(m)))
Daniel@0:                                 % No suitable zombie.
Daniel@0:                                 % birth of a new partial (Mc Aulay &
Daniel@0:                                 % Quatieri)
Daniel@0:                                 mxl = [mxl;zeros(1,k+1)];
Daniel@0:                                 tr = size(mxl,1);
Daniel@0:                                 mxl(tr,k) = mxk2(m)-1;
Daniel@0:                             else
Daniel@0:                                 % Suitable zombie found. (Lartillot)
Daniel@0:                                 tr = grvy(z,1);
Daniel@0:                                 grvy(z,:) = [];
Daniel@0:                             end
Daniel@0:                             mxl(tr,k+1) = mxk2(m)-1;
Daniel@0:                             myl(tr,k+1) = myk2(m);
Daniel@0:                             tr2(m) = tr;
Daniel@0:                         end
Daniel@0:                     end
Daniel@0:                 end
Daniel@0:                 
Daniel@0:                 if wait
Daniel@0:                     waitbar(1,wait);
Daniel@0:                     close(wait);
Daniel@0:                     drawnow
Daniel@0:                 end
Daniel@0: 
Daniel@0:                 if size(mxl,1) > option.m
Daniel@0:                     tot = sum(myl,2);
Daniel@0:                     [tot ix] = sort(tot,'descend');
Daniel@0:                     mxl(ix(option.m+1:end),:) = [];
Daniel@0:                     myl(ix(option.m+1:end),:) = [];
Daniel@0:                 end
Daniel@0:                 
Daniel@0:                 mxl(:,not(max(myl))) = 0;
Daniel@0:                 
Daniel@0:                 if option.shorttrackthresh
Daniel@0:                     [myl bestrack] = max(myl,[],1);
Daniel@0:                     mxl = mxl(bestrack + (0:size(mxl,2)-1)*size(mxl,1));
Daniel@0:                     changes = find(not(bestrack(1:end-1) == bestrack(2:end)))+1;
Daniel@0:                     if not(isempty(changes))
Daniel@0:                         lengths = diff([1 changes nc+1]);
Daniel@0:                         shorts = find(lengths < option.shorttrackthresh);
Daniel@0:                         for k = 1:length(shorts)
Daniel@0:                             if shorts(k) == 1
Daniel@0:                                 k1 = 1;
Daniel@0:                             else
Daniel@0:                                 k1 = changes(shorts(k)-1);
Daniel@0:                             end
Daniel@0:                             k2 = k1 + lengths(shorts(k)) -1;
Daniel@0:                             myl(1,k1:k2) = 0;
Daniel@0:                             mxl(1,k1:k2) = 0;
Daniel@0:                         end
Daniel@0:                     end
Daniel@0:                 end
Daniel@0:                 
Daniel@0:                 tp{i}{h}{l} = mxl;
Daniel@0:                 tv{i}{h}{l} = myl;
Daniel@0:                 
Daniel@0:                 if interpol  
Daniel@0:                     tpv{i}{h}{l} = zeros(size(mxl));
Daniel@0:                     tpp{i}{h}{l} = zeros(size(mxl));
Daniel@0:                     for k = 1:size(mxl,2)
Daniel@0:                         for j = 1:size(mxl,1)
Daniel@0:                             mj = mxl(j,k);
Daniel@0:                             if mj>2 && mj<size(dh3,1)-1
Daniel@0:                                 % More precise peak position
Daniel@0:                                 y0 = dh3(mj,k,l);
Daniel@0:                                 ym = dh3(mj-1,k,l);
Daniel@0:                                 yp = dh3(mj+1,k,l);
Daniel@0:                                 p = (yp-ym)/(2*(2*y0-yp-ym));
Daniel@0:                                 tpv{i}{h}{l}(j,k) = y0 - 0.25*(ym-yp)*p;
Daniel@0:                                 if p >= 0
Daniel@0:                                     tpp{i}{h}{l}(j,k) = (1-p)*th(mj,k,l)+p*th(mj+1,k,l);
Daniel@0:                                 elseif p < 0
Daniel@0:                                     tpp{i}{h}{l}(j,k) = (1+p)*th(mj,k,l)-p*th(mj-1,k,l);
Daniel@0:                                 end
Daniel@0:                             elseif mj
Daniel@0:                                 tpv{i}{h}{l}(j,k) = dh3(mj,k,l);
Daniel@0:                                 tpp{i}{h}{l}(j,k) = th(mj,k,l);
Daniel@0:                             end
Daniel@0:                         end
Daniel@0:                     end
Daniel@0:                 end
Daniel@0:             end
Daniel@0:         end
Daniel@0:         if isa(x,'mirsimatrix') && option.graph
Daniel@0:             % Finding the best branch inside a graph constructed out of a
Daniel@0:             % similarity matrix
Daniel@0:             g{i}{h} = cell(1,nc,np);
Daniel@0:                 % Branch info related to each peak
Daniel@0:             br{i}{h} = {};
Daniel@0:                 % Info related to each branch
Daniel@0:             scog{i}{h} = cell(1,nc,np);
Daniel@0:                 % Score related to each peak
Daniel@0:             scob{i}{h} = [];
Daniel@0:                 % Score related to each branch
Daniel@0:             for l = 1:np
Daniel@0:                 wait = waitbar(0,['Creating peaks graph...']);
Daniel@0:                 for k = 1:nc
Daniel@0:                     g{i}{h}{1,k,l} = cell(size(mx{1,k,l}));
Daniel@0:                     scog{i}{h}{1,k,l} = zeros(size(mx{1,k,l}));
Daniel@0:                     if wait && not(mod(k,50))
Daniel@0:                         waitbar(k/(nc-1),wait);
Daniel@0:                     end
Daniel@0:                     mxk = mx{1,k,l}; % Peaks in current frame
Daniel@0:                     for j = k-1:-1:max(1,k-100) % Recent frames
Daniel@0:                         mxj = mx{1,j,l};        % Peaks in one recent frame
Daniel@0:                         for kk = 1:length(mxk)
Daniel@0:                             mxkk = mxk(kk);     % For each of current peaks
Daniel@0:                             for jj = 1:length(mxj)
Daniel@0:                                 mxjj = mxj(jj); % For each of recent peaks
Daniel@0:                                 sco = k-j - abs(mxkk-mxjj);
Daniel@0:                                     % Crossprogression from recent to
Daniel@0:                                     % current peak
Daniel@0:                                 if sco >= 0 
Daniel@0:                                         % Negative crossprogression excluded
Daniel@0:                                     dist = 0;
Daniel@0:                                     % The distance between recent and
Daniel@0:                                     % current peak is the sum of all the
Daniel@0:                                     % simatrix values when joining the two
Daniel@0:                                     % peaks with a straight line.
Daniel@0:                                     for m = j:k
Daniel@0:                                         % Each point in that straight line.
Daniel@0:                                         mxm = mxjj + (mxkk-mxjj)*(m-j)/(k-j);
Daniel@0:                                         if mxm == floor(mxm)
Daniel@0:                                             dist = dist + 1-dh(mxm,m,l);
Daniel@0:                                         else
Daniel@0:                                             dhm0 = dh(floor(mxm),m,l);
Daniel@0:                                             dhm1 = dh(ceil(mxm),m,l);
Daniel@0:                                             dist = dist + 1-...
Daniel@0:                                                 (dhm0 + ...
Daniel@0:                                                  (dhm1-dhm0)*(mxm-floor(mxm)));
Daniel@0:                                         end
Daniel@0:                                         if dist > option.graph
Daniel@0:                                             break
Daniel@0:                                         end
Daniel@0:                                     end
Daniel@0:                                     if dist < option.graph
Daniel@0:                                         % If the distance between recent
Daniel@0:                                         % and current peak is not too high,
Daniel@0:                                         % a new edge is formed between the
Daniel@0:                                         % peaks, and added to the graph.
Daniel@0:                                         gj = g{i}{h}{1,j,l}{jj};
Daniel@0:                                             % Branch information associated
Daniel@0:                                             % with recent peak
Daniel@0:                                         gk = g{i}{h}{1,k,l}{kk};
Daniel@0:                                             % Branch information associated
Daniel@0:                                             % with current peak
Daniel@0:                                         if isempty(gk) || ...
Daniel@0:                                                 sco > scog{i}{h}{1,k,l}(kk)
Daniel@0:                                             % Current peak branch to be updated
Daniel@0:                                             if isempty(gj)
Daniel@0:                                                 % New branch starting
Daniel@0:                                                 % from scratch
Daniel@0:                                                 newsco = sco;
Daniel@0:                                                 scob{i}{h}(end+1) = newsco;
Daniel@0:                                                 bid = length(scob{i}{h});
Daniel@0:                                                 g{i}{h}{1,j,l}{jj} = ...
Daniel@0:                                                     [k kk bid newsco];
Daniel@0:                                                 br{i}{h}{bid} = [j jj;k kk];
Daniel@0:                                             else
Daniel@0:                                                 newsco = scog{i}{h}{1,j,l}(jj)+sco;
Daniel@0:                                                 if length(gj) == 1
Daniel@0:                                                     % Recent peak not
Daniel@0:                                                     % associated with other
Daniel@0:                                                     % branch
Daniel@0:                                                     % -> Branch extension
Daniel@0:                                                     bid = gj;
Daniel@0:                                                     g{i}{h}{1,j,l}{jj} = ...
Daniel@0:                                                         [k kk bid newsco];
Daniel@0:                                                     br{i}{h}{bid}(end+1,:) = [k kk];
Daniel@0:                                                 else
Daniel@0:                                                     % Recent peak already
Daniel@0:                                                     % associated with other
Daniel@0:                                                     % branch
Daniel@0:                                                     % -> Branch fusion
Daniel@0:                                                     bid = length(scob{i}{h})+1;
Daniel@0:                                                     g{i}{h}{1,j,l}{jj} = ...
Daniel@0:                                                         [k kk bid newsco; gj];
Daniel@0:                                                     other = br{i}{h}{gj(1,3)};
Daniel@0:                                                         % Other branch
Daniel@0:                                                         % info
Daniel@0:                                                         % Let's copy its
Daniel@0:                                                         % prefix to the new
Daniel@0:                                                         % branch:
Daniel@0:                                                     other(other(:,1)>j,:) = [];
Daniel@0:                                                     br{i}{h}{bid} = [other;k kk];
Daniel@0:                                                 end
Daniel@0:                                                 scob{i}{h}(bid) = newsco;
Daniel@0:                                             end
Daniel@0:                                             g{i}{h}{1,k,l}{kk} = bid;
Daniel@0:                                                 % New peak associated with
Daniel@0:                                                 % branch
Daniel@0:                                             scog{i}{h}{1,k,l}(kk) = newsco;
Daniel@0:                                         end
Daniel@0:                                     end
Daniel@0:                                 end
Daniel@0:                             end
Daniel@0:                         end
Daniel@0:                     end
Daniel@0:                 end
Daniel@0:                 [scob{i}{h} IX] = sort(scob{i}{h},'descend');
Daniel@0:                     % Branch are ordered from best score to lowest
Daniel@0:                 br{i}{h} = br{i}{h}(IX);
Daniel@0:                 if wait
Daniel@0:                     waitbar(1,wait);
Daniel@0:                     close(wait);
Daniel@0:                     drawnow
Daniel@0:                 end
Daniel@0:             end
Daniel@0:         end
Daniel@0:         for l = 1:np % Orders the peaks and select the best ones
Daniel@0:             for k = 1:nc
Daniel@0:                 mxk = mx{1,k,l};
Daniel@0:                 if length(mxk) > option.m 
Daniel@0:                     [unused,idx] = sort(dh(mxk,k,l),'descend');
Daniel@0:                     idx = idx(1:option.m);
Daniel@0:                 elseif strcmpi(option.order,'Amplitude')
Daniel@0:                     [unused,idx] = sort(dh(mxk,k,l),'descend');
Daniel@0:                 else
Daniel@0:                     idx = 1:length(dh(mxk,k,l));
Daniel@0:                 end
Daniel@0:                 if strcmpi(option.order,'Abscissa')
Daniel@0:                     mx{1,k,l} = sort(mxk(idx));
Daniel@0:                 elseif strcmpi(option.order,'Amplitude')
Daniel@0:                     mx{1,k,l} = mxk(idx);
Daniel@0:                 end
Daniel@0:             end
Daniel@0:         end
Daniel@0:         if option.extract % Extracting the positive part of the curve containing the peaks
Daniel@0:             if isa(x,'mirtemporal')
Daniel@0:                 filn = floor(sr{i}/25);
Daniel@0:             else
Daniel@0:                 filn = 10;
Daniel@0:             end
Daniel@0:             if filn>1 && size(dh3,1)>5
Daniel@0:                 filn = min(filn,floor(size(dh3,1)/3));
Daniel@0:                 fild = filtfilt(ones(1,filn)/2,1,dh3(2:end-1,:,:))/filn/2;
Daniel@0:             else
Daniel@0:                 fild = dh3(2:end-1,:,:);
Daniel@0:             end
Daniel@0:             fild = [zeros(1,size(fild,2),size(fild,3));diff(fild)];
Daniel@0:             for l = 1:np
Daniel@0:                 for k = 1:nc
Daniel@0:                     idx = 1:size(dh,1);
Daniel@0:                     mxlk = sort(mx{1,k,l}-1);
Daniel@0:                     for j = 1:length(mxlk)
Daniel@0:                         
Daniel@0:                         if fild(mxlk(j),k,l) < 0
Daniel@0:                             bef0 = find(fild(1:mxlk(j)-1,k,l)>=0);
Daniel@0:                             if isempty(bef0)
Daniel@0:                                 bef0 = [];
Daniel@0:                             end
Daniel@0:                         else
Daniel@0:                             bef0 = mxlk(j)-1;
Daniel@0:                         end
Daniel@0:                         if isempty(bef0)
Daniel@0:                             bef = 0;
Daniel@0:                         else
Daniel@0:                             bef = find(fild(1:bef0(end),k,l)<=0);
Daniel@0:                             if isempty(bef)
Daniel@0:                                 bef = 0;
Daniel@0:                             end
Daniel@0:                         end
Daniel@0:                         if j>1 && bef(end)<aft(1)+2
Daniel@0:                             idx(mxlk(j-1):mxlk(j)) = 0;
Daniel@0:                             [unused btw] = min(dh3(mxlk(j-1)+1:mxlk(j)+1,k,l));
Daniel@0:                             btw = btw+mxlk(j-1);
Daniel@0:                             idx(btw-2:btw+2) = btw-2:btw+2;
Daniel@0:                             bef = btw+2;
Daniel@0:                         end
Daniel@0:                         
Daniel@0:                         if fild(mxlk(j),k,l) > 0
Daniel@0:                             aft0 = find(fild(mxlk(j)+1:end,k,l)<=0)+mxlk(j);
Daniel@0:                             if isempty(aft0)
Daniel@0:                                 aft0 = [];
Daniel@0:                             end
Daniel@0:                         else
Daniel@0:                             aft0 = mxlk(j)+1;
Daniel@0:                         end
Daniel@0:                         if isempty(aft0)
Daniel@0:                             aft = size(d{i}{h},1)+1;
Daniel@0:                         else
Daniel@0:                             aft = find(fild(aft0(1):end,k,l)>=0)+mxlk(j);
Daniel@0:                             if isempty(aft)
Daniel@0:                                 aft = size(d{i}{h},1)+1;
Daniel@0:                             end
Daniel@0:                         end
Daniel@0:                         
Daniel@0:                         idx(bef(end)+3:aft(1)-3) = 0;
Daniel@0:                     end
Daniel@0:                     idx = idx(find(idx));
Daniel@0:                     dh3(idx,k,l) = NaN;
Daniel@0:                 end
Daniel@0:             end
Daniel@0:         end
Daniel@0:         if option.vall
Daniel@0:             dh3 = -dh3;
Daniel@0:         end
Daniel@0:         mmx = cell(1,nc,np);
Daniel@0:         mmy = cell(1,nc,np);
Daniel@0:         mmv = cell(1,nc,np);
Daniel@0:         for l = 1:np
Daniel@0:             for k = 1:nc
Daniel@0:                 mmx{1,k,l} = mod(mx{1,k,l}(:,:,1),nl0)-1;
Daniel@0:                 mmy{1,k,l} = ceil(mx{1,k,l}/nl0);
Daniel@0:                 mmv{1,k,l} = dh3(mx{1,k,l}(:,:,1),k,l);
Daniel@0:             end
Daniel@0:         end
Daniel@0:         pp{i}{h} = mmx;
Daniel@0:         pm{i}{h} = mmy;
Daniel@0:         pv{i}{h} = mmv;
Daniel@0:         if not(interpol)
Daniel@0:             ppp{i}{h} = {};
Daniel@0:             ppv{i}{h} = {};
Daniel@0:         else % Interpolate to find the more exact peak positions
Daniel@0:             pih = cell(1,nc,np);
Daniel@0:             vih = cell(1,nc,np);
Daniel@0:             for l = 1:np
Daniel@0:                 for k = 1:nc
Daniel@0:                     mxlk = mx{1,k,l};
Daniel@0:                     vih{1,k,l} = zeros(length(mxlk),1);
Daniel@0:                     pih{1,k,l} = zeros(length(mxlk),1);
Daniel@0:                     for j = 1:length(mxlk)
Daniel@0:                         mj = mxlk(j); % Current values
Daniel@0:                         if strcmpi(option.interpol,'quadratic')
Daniel@0:                             if mj>2 && mj<length(dh3)-1
Daniel@0:                                 % More precise peak position
Daniel@0:                                 y0 = dh3(mj,k,l);
Daniel@0:                                 ym = dh3(mj-1,k,l);
Daniel@0:                                 yp = dh3(mj+1,k,l);
Daniel@0:                                 p = (yp-ym)/(2*(2*y0-yp-ym));
Daniel@0:                                 vih{1,k,l}(j) = y0 - 0.25*(ym-yp)*p;
Daniel@0:                                 if p >= 0
Daniel@0:                                     pih{1,k,l}(j) = (1-p)*th(mj,k,l)+p*th(mj+1,k,l);
Daniel@0:                                 elseif p < 0
Daniel@0:                                     pih{1,k,l}(j) = (1+p)*th(mj,k,l)-p*th(mj-1,k,l);
Daniel@0:                                 end
Daniel@0:                             else
Daniel@0:                                 vih{1,k,l}(j) = dh3(mj,k,l);
Daniel@0:                                 pih{1,k,l}(j) = th(mj,k,l);
Daniel@0:                             end
Daniel@0:                         end
Daniel@0:                     end
Daniel@0:                 end
Daniel@0:             end
Daniel@0:             ppp{i}{h} = pih;
Daniel@0:             ppv{i}{h} = vih;
Daniel@0:         end
Daniel@0:         if not(iscell(d{i})) % for chromagram
Daniel@0:             d{i} = dh3(2:end-1,:,:,:);
Daniel@0:         else
Daniel@0:             if cha == 1
Daniel@0:                 d{i}{h} = zeros(1,size(dh3,2),size(dh3,1)-2);
Daniel@0:                 for k = 1:size(dh3,2)
Daniel@0:                      d{i}{h}(1,k,:) = dh3(2:end-1,k);
Daniel@0:                 end
Daniel@0:             else
Daniel@0:                 d{i}{h} = dh3(2:end-1,:,:,:);
Daniel@0:             end
Daniel@0:         end
Daniel@0:         if option.only
Daniel@0:             dih = zeros(size(d{i}{h}));
Daniel@0:             for l = 1:np
Daniel@0:                 for k = 1:nc
Daniel@0:                     dih(pp{i}{h}{1,k,l},k,l) = ...
Daniel@0:                         d{i}{h}(pp{i}{h}{1,k,l},k,l);
Daniel@0:                 end
Daniel@0:             end
Daniel@0:             d{i}{h} = dih;
Daniel@0:         end
Daniel@0:     end
Daniel@0: end
Daniel@0: p = set(x,'PeakPos',pp,'PeakVal',pv,'PeakMode',pm);
Daniel@0: if interpol
Daniel@0:    p = set(p,'PeakPrecisePos',ppp,'PeakPreciseVal',ppv);
Daniel@0: end
Daniel@0: if option.extract
Daniel@0:     p = set(p,'Data',d);
Daniel@0: end
Daniel@0: empty = cell(1,length(d));
Daniel@0: if option.only
Daniel@0:     p = set(p,'Data',d,'PeakPos',empty,'PeakVal',empty,'PeakMode',empty);
Daniel@0: end
Daniel@0: if option.delta
Daniel@0:     p = set(p,'TrackPos',tp,'TrackVal',tv);
Daniel@0:     if interpol
Daniel@0:        p = set(p,'TrackPrecisePos',tpp,'TrackPreciseVal',tpv);
Daniel@0:     end
Daniel@0: end
Daniel@0: if isa(x,'mirsimatrix') && option.graph
Daniel@0:     p = set(p,'Graph',g,'Branch',br);
Daniel@0: end
Daniel@0: 
Daniel@0: 
Daniel@0: function y = semitone_compar(p1,p2,thres)
Daniel@0: y = p1/p2 < 2^(1/12);
Daniel@0: 
Daniel@0: 
Daniel@0: function y = dist_compar(p1,p2,thres)
Daniel@0: y = p1-p2 < thres;