% Version 1.00 stable
%-------------------------------Copyright----------------------------------
%
% Author: Dmytro Iatsenko
% Information about these codes (e.g. links to the Video Instructions),
% as well as other MatLab programs and many more can be found at
% http://www.physics.lancs.ac.uk/research/nbmphysics/diats/tfr
% 
% Related articles:
% [1] D. Iatsenko, A. Stefanovska and P.V.E. McClintock,
% "Linear and synchrosqueezed time-frequency representations revisited.
%  Part I: Overview, standards of use, related issues and algorithms."
% {preprint:arXiv:1310.7215}
% [2] D. Iatsenko, A. Stefanovska and P.V.E. McClintock,
% "Linear and synchrosqueezed time-frequency representations revisited.
%  Part II: Resolution, reconstruction and concentration."
% {preprint:arXiv:1310.7274}
% [3] D. Iatsenko, A. Stefanovska and P.V.E. McClintock,
% "On the extraction of instantaneous frequencies from ridges in
%  time-frequency representations of signals."
% {preprint:arXiv:1310.7276}
%
%------------------------------Documentation-------------------------------
%
% [tfsupp,Optional:ecinfo,Skel] = ssecurve(TFR,freq,wopt,'PropertyName',PropertyValue)
% - extracts the curve (i.e. the consequence of the amplitude ridge points)
%   and its full support (the widest region of unimodal TFR amplitude
%   around them) from the given time-frequency representation [TFR] of the
%   signal. TFR can be either SWFT (calculated by sswft.m) or SWT (swt.m).
%
% OUTPUT:
% tfsupp: 3xL matrix
%    - extracted time-frequency support of the component, containing
%      frequencies of the TFR amplitude peaks (ridge points) in the
%      first row (referred as \omega_p(t)/2\pi in [3]), support lower
%      bounds (referred as \omega_-(t)/2/pi in [1]) - in the second row,
%      and the upper bounds (referred as \omega_+(t)/2/pi in [1]) - in
%      the third row.
% ecinfo: structure
%    - contains all the relevant information about the process of curve
%      extraction.
% Skel: 4x1 cell (returns empty matrix [] if 'Method' is [efreq], see below)
%    - contains the number of peaks N_p(t) in [Skel{1}], their frequency
%      indices m(t) in [Skel{2}], the corresponding frequencies
%      \nu_m(t)/2\pi in [Skel{3}], and the respective amplitudes Q_m(t)
%      in [Skel{4}] (notations according to [3]).
% 
% INPUT: ({{...}} denotes default)
% TFR: NFxL matrix (rows correspond to frequencies, columns - to time)
%    - SWFT or SWT from which to extract [tfsupp]
% freq: NFx1 vector
%    - the frequencies corresponding to the rows of [TFR]
% wopt: structure returned as third output by sswft.m or sswt.m
%    - parameters of the window/wavelet and the simulation, returned as a
%      third output by functions wft, wt, sswft and sswt; [wopt] contains
%      all the needed information, i.e. name of the TFR, sampling frequency,
%      window/wavelet characteristics etc.
%
% Properties:
% 'Method': 1|{{2}}|3|'nearest'|'max'|efreq
%    - method by which to extract the curve, as described in [3]; to use
%      frequency-based extraction, specify an L-length vector [efreq] of
%      the frequencies (in Hz), in which case the program will form the
%      ridge curve from the peaks nearest to the specified frequency
%      profile [efreq] (in fact, not nearest but lying at each time in the
%      same support, i.e. region of nonzero TFR amplitude).
% 'Param': alpha - method==1 (default 1); if alpha=0 or alpha=Inf, switches to
%                  maximum-based or frequency-based extraction, respectively.
%          [alpha,beta] - method==2 or 3 (default [1,1])
%          [] (no parameters) - method=='nearest' or 'max' or [efreq]
%    - parameters for each method, as described in [3].
% 'Display': {{'on'}}|'off'|'plot'|'plot+'
%    - to display the progress or not; if set to 'plot', additionally plots
%      all the obtained curves and parameters, e.g. within polishing stage
%      in methods 'II' and 'III'; if 'plot+', plots additionally the
%      information within each iteration, e.g. parameters convergence.
% 'Plot': 'on'|{{'off'}}
%    - if set to 'on', plots additionally the full [TFR] and shows on it
%      the extracted frequency profile and the boundaries of its support at
%      each time (don't confuse with Display='plot', which shows more
%      specific and detailed information).
% 'Skel': 4x1 cell as returned in the third output argument of ecurve.m
%    - can be specified to avoid performing search of the peaks, their
%      numbers and corresponding amplitudes each time if the procedures are
%      applied to the same TFR (e.g. one compares the performance of
%      different schemes).
%
% Warning: if 'Plot' property is set to 'on', it plots the NFxL matrix of
% TFR, which might be quite big. If it is too large, e.g. 1000x50000,
% then the computer might run out of memory (as graphical representation
% of a matrix requires an additional memory) and become "frozen", so be
% careful. The same concern 'Display' set to 'plot+' when 'Method'==3.
%
%-------------------------------Examples-----------------------------------
%
% [SWFT,freq,wopt]=sswft(sig,fs); tfsupp=ssecurve(SWFT,freq,wopt);
% - extracts the ridge curve from the synchrosqueezed windowed Fourier
%   transform [SWFT] of the signal [sig] sampled at [fs] Hz.
%
%--------------------------------------------------------------------------

function [tfsupp,varargout] = ssecurve(TFR,freq,wopt,varargin)

[NF,L]=size(TFR); tfsupp=zeros(3,L)*NaN; freq=freq(:);

method=2; pars=[]; Display='on'; PlotMode='off'; Skel=[];
for vn=1:2:nargin-3
    if strcmpi(varargin{vn},'Method'), method=varargin{vn+1};
    elseif strcmpi(varargin{vn},'Param'), pars=varargin{vn+1};
    elseif strcmpi(varargin{vn},'Display'), Display=varargin{vn+1};
    elseif strcmpi(varargin{vn},'Plot'), PlotMode=varargin{vn+1};
    elseif strcmpi(varargin{vn},'Skel'), Skel=varargin{vn+1};
    else
        error(['There is no Property "',varargin{vn},'" (which is ',num2str(1+(vn-vst)/2,'%d'),...
            '''th out of ',num2str(ceil((nargin-1-vst)/2),'%d'),' specified)']);
    end
end
if isempty(pars) %if parameters are not specified, assign the defaults
    if ischar(method) || length(method)>1, pars=[];
    elseif method==1, pars=1;
    else pars=[1,1];
    end
else %if parameters are specified, check are they specified correctly
    if (ischar(method) || length(method)>1) && ~isempty(pars)
        error('Wrong number of parameters for the chosen method: there should be no parameters.');
    elseif method==1 && length(pars)~=1
        error('Wrong number of parameters for the chosen method: there should be one parameter.');
    elseif (method==2 || method==3) && length(pars)~=2
        error('Wrong number of parameters for the chosen method: there should be two parameters.');
    end
end
if ~ischar(method) && length(pars)==1 && pars(1)==0, method='max'; pars=[]; end %switch to maximum-based extraction if method is I(0)
if ~ischar(method) && length(pars)==1 && ~isfinite(pars(1)), method='nearest'; pars=[]; end %switch to nearest neighbour extraction if method is I(Inf)
if nargout>1, ec=struct; ec.Method=method; ec.Param=pars; end

%Determine the frequency resolution
if min(freq)<=0 || std(diff(freq))<std(diff(log(freq)))
    fres=1; fstep=mean(diff(freq));
    Dtau=wopt.wp.t2h-wopt.wp.t1h; Dxi=wopt.wp.xi2h-wopt.wp.xi1h;
    dfreq=[freq(1)-freq(end:-1:2);freq(end)-freq(end:-1:1)];
else
    fres=2; fstep=mean(diff(log(freq)));
    Dtau=wopt.wp.t2h-wopt.wp.t1h; Dxi=log(wopt.wp.xi2h/wopt.wp.xi1h);
    dfreq=[log(freq(1))-log(freq(end:-1:2));log(freq(end))-log(freq(end:-1:1))];
end

%//////////////////////////////////////////////////////////////////////////
tn1=find(~isnan(TFR(end,:)),1,'first'); tn2=find(~isnan(TFR(end,:)),1,'last'); sflag=0;
if ischar(method) || length(method)==1 %if not frequency-based extraction
    sflag=1;
    %----------------------------------------------------------------------
    %Construct matrices of peak indices, frequencies and total SWFT/SWT
    %amplitudes: [Ip],[Fp],[Qp]; [Np] - number of peaks at each time.
    %(takes much more time than similar procedure for WFT/WT in ecurve.m)
    if ~isempty(Skel)
        Np=Skel{1}; Ip=Skel{2}; Fp=Skel{3}; Qp=Skel{4}; Mp=max(Np);
    else
        if ~strcmpi(Display,'off'), fprintf('Locating supports in TFR and joining them into peaks.\n'); end
        TFR=vertcat(zeros(1,L),TFR,zeros(1,L));
        idft1=1+find(TFR(1:end-1)==0 & TFR(2:end)~=0); idft2=find(TFR(1:end-1)~=0 & TFR(2:end)==0);
        [idf1,idt1]=ind2sub(size(TFR),idft1); [idf2,idt2]=ind2sub(size(TFR),idft2);
        dind=[0;find(diff(idt1(:))>0);length(idt1)]; Mp=max(diff(dind));
        Np=zeros(1,L); for dn=1:length(dind)-1, Np(idt1(dind(dn)+1))=dind(dn+1)-dind(dn); end
        Ip=ones(Mp,L)*NaN; Fp=ones(Mp,L)*NaN; Qp=ones(Mp,L)*NaN; nid=0; tn=idt1(1);
        for kn=1:length(idft1)
            ii=(idft1(kn):idft2(kn))'; iif=-1+(idf1(kn):idf2(kn))';
            casig=sum(TFR(ii)); cfsig=real(sum(freq(iif).*TFR(ii))/casig);
            if isnan(cfsig) || cfsig<freq(iif(1)) || cfsig>freq(iif(end))
                cfsig=sum(freq(iif).*abs(TFR(ii)))/sum(abs(TFR(ii)));
            end
            if fres==1, fid=1+floor((1/2)+(cfsig-freq(1))/fstep);
            else fid=1+floor((1/2)+(log(cfsig)-log(freq(1)))/fstep); end
            fid=min([max([fid,1]),NF]);
            if idt1(kn)~=tn, nid=1; else nid=nid+1; end
            tn=idt1(kn); Ip(nid,tn)=fid; Fp(nid,tn)=freq(fid); Qp(nid,tn)=abs(casig);
        end
        clear idft1 idft2 idf1 idt1 idf2 idt2 dind;
        TFR=TFR(2:end-1,:); %recover the original TFR
    end
    if nargout>2, varargout{2}={Np,Ip,Fp,Qp}; end
    %----------------------------------------------------------------------
    
    %Find the position of the maximum from which to start
    [~,imax]=max(Qp(:)); [fimax,timax]=ind2sub([Mp,L],imax);
    if nargout>1, ec.itmax=timax; ec.ifmax=Fp(fimax,timax); ec.iqmax=Qp(fimax,timax); end
    
else %frequency-based extraction
    if length(method)~=L
        error('The specified frequency profile ("Method" property) should be of the same length as signal.');
    end
    
    if ~strcmpi(Display,'off')
        fprintf('Extracting the curve near the specified frequency profile.\n');
    end
    
    efreq=method; tn1=max([tn1,find(~isnan(efreq),1,'first')]); tn2=min([tn2,find(~isnan(efreq),1,'last')]);
    if fres==1, eind=1+floor(0.5+(efreq-freq(1))/fstep);
    else eind=1+floor(0.5+log(efreq/freq(1))/fstep); end
    eind(eind<1)=1; eind(eind>NF)=NF;
    
    %Extract the indices of the nearest peaks and time-frequency support around them
    for tn=tn1:tn2
        cind=eind(tn); cv=TFR(:,tn); cs=abs(cv);
        
        %Ridge point
        if cind>1 && cind<NF
            if cs(cind+1)==cs(cind-1)
                cpeak1=cind-1+find(cs(cind:end-1)>=cs(cind-1:end-2) & cs(cind:end-1)>cs(cind+1:end),1,'first'); cpeak1=min([cpeak1,NF]);
                cpeak2=cind+1-find(cs(cind:-1:2)>=cs(cind+1:-1:3) & cs(cind:-1:2)>cs(cind-1:-1:1),1,'first'); cpeak2=max([cpeak2,1]);
                if cs(cpeak2)==0, cpeak=cpeak1; elseif cs(cpeak1)==0, cpeak=cpeak2;
                elseif cpeak1-cind<cind-cpeak2, cpeak=cpeak1; else cpeak=cpeak2; end
            elseif cs(cind+1)>cs(cind-1)
                cpeak=cind-1+find(cs(cind:end-1)>=cs(cind-1:end-2) & cs(cind:end-1)>cs(cind+1:end),1,'first'); cpeak=min([cpeak,NF]);
            elseif cs(cind+1)<cs(cind-1)
                cpeak=cind+1-find(cs(cind:-1:2)>cs(cind-1:-1:1) & cs(cind:-1:2)>=cs(cind+1:-1:3),1,'first'); cpeak=max([cpeak,1]);
            end
        elseif cind==1
            if cs(2)<cs(1), cpeak=cind;
            else
                cpeak=1+find(cs(cind+1:end-1)>=cs(cind:end-2) & cs(cind+1:end-1)>cs(cind+2:end),1,'first'); cpeak=min([cpeak,NF]);
            end
        elseif cind==NF
            if cs(NF-1)<cs(NF), cpeak=cind;
            else
                cpeak=NF-find(cs(cind-1:-1:2)>cs(cind-2:-1:1) & cs(cind-1:-1:2)>=cs(cind:-1:3),1,'first'); cpeak=max([cpeak,1]);
            end
        end
        if cs(cpeak)==0, cpeak=cind; end
        
        %Boundaries of time-frequency support
        iup=[]; idown=[];
        if cpeak<NF-1, iup=cpeak+find(cs(cpeak+1:end)==0,1,'first'); end
        if cpeak>2, idown=cpeak-find(cs(cpeak-1:-1:1)==0,1,'first'); end
        iup=min([iup,NF]); idown=max([idown,1]);
        tfsupp(2,tn)=idown; tfsupp(3,tn)=iup;
        
        %Find central frequency
        ii=idown:iup; cfsig=real(sum(freq(ii).*cv(ii))/sum(cv(ii)));
        if cfsig>=freq(ii(1)) && cfsig<=freq(ii(end)), tfsupp(1,tn)=cfsig;
        else tfsupp(1,tn)=sum(freq(ii).*cs(ii))/sum(cs(ii)); end
    end
    %Transform to frequencies
    if fres==1, ridges=1+floor(0.5+(tfsupp(1,:)-freq(1))/fstep)';
    else ridges=1+floor(0.5+log(tfsupp(1,:)/freq(1))/fstep)'; end
    tfsupp(2:3,tn1:tn2)=freq(tfsupp(2:3,tn1:tn2));
    
    %Optional output arguments
    if nargout>1
        ec.efreq=efreq; ec.eind=eind;
        ec.pfreq=[ones(1,tn1-1)*NaN,freq(ridges(tn1:tn2))',ones(1,L-tn2)*NaN];
        ec.pind=ridges(:)'; %peak indices
        ec.pamp=zeros(1,L)*NaN; ec.pamp(tn1:tn2)=abs(TFR(sub2ind(size(TFR),ridges(tn1:tn2)',tn1:tn2)));
        varargout{1}=ec;
    end
    if nargout>2, varargout{2}=[]; end
    
    %Plotting (if needed)
    if ~isempty(strfind(Display,'plot'))
        scrsz=get(0,'ScreenSize'); figure('Position',[scrsz(3)/4,scrsz(4)/8,2*scrsz(3)/3,2*scrsz(4)/3]);
        ax=axes('Position',[0.1,0.1,0.8,0.8],'FontSize',16); hold all;
        title(ax(1),'Ridge curve \omega_p(t)/2\pi'); ylabel(ax(1),'Frequency (Hz)'); xlabel(ax(1),'Time (s)');
        plot(ax(1),(0:L-1)/wopt.fs,efreq,'--','Color',[0.5,0.5,0.5],'LineWidth',2,'DisplayName','Specified frequency profile');
        plot(ax(1),(0:L-1)/wopt.fs,tfsupp(1,:),'-k','LineWidth',2,'DisplayName','Extracted frequency profile');
        legend(ax(1),'show');
    end
    if strcmpi(PlotMode,'on'), plotfinal(tfsupp,TFR,freq,wopt); end
    
    return;
end
%//////////////////////////////////////////////////////////////////////////


%--------------------------- Global Maximum -------------------------------
if strcmpi(method,'max') || length(pars)==2
    if ~strcmpi(Display,'off')
        if sflag==0, fprintf('Extracting the curve by Global Maximum scheme.\n');
        else fprintf('Extracting positions of global maximums.\n'); end
    end
    
    ridges=zeros(L,1)*NaN; maxv=zeros(L,1)*NaN;
    if sflag==0
        for tn=tn1:tn2, [maxv(tn),ridges(tn)]=max(abs(TFR(:,tn))); end
    else
        for tn=tn1:tn2, [maxv(tn),mid]=max(Qp(1:Np(tn),tn)); ridges(tn)=Ip(mid,tn); end
    end
    idz=tn1-1+find(maxv(tn1:tn2)==0);
    if ~isempty(idz)
        idnz=tn1:tn2; idnz=idnz(~ismember(idnz,idz));
        ridges(idz)=interp1(idnz,ridges(idnz),idz,'linear','extrap');
        ridges(idz)=round(ridges(idz));
    end
    tfsupp(1,:)=ridges(:)';
    
    if nargout>1
        ec.pfreq=[ones(1,tn1-1)*NaN,freq(ridges(tn1:tn2))',ones(1,L-tn2)*NaN];
        ec.pind=ridges(:)'; %peak indices
        ec.pamp=zeros(1,L)*NaN; ec.pamp(tn1:tn2)=abs(TFR(sub2ind(size(TFR),ridges(tn1:tn2)',tn1:tn2)));
    end
    
    if ~isempty(strfind(Display,'plot')) && strcmpi(method,'max')
        scrsz=get(0,'ScreenSize'); figure('Position',[scrsz(3)/4,scrsz(4)/8,2*scrsz(3)/3,2*scrsz(4)/3]);
        ax=axes('Position',[0.1,0.1,0.8,0.8],'FontSize',16); hold all;
        title(ax(1),'Ridge curve \omega_p(t)/2\pi'); ylabel(ax(1),'Frequency (Hz)'); xlabel(ax(1),'Time (s)');
        plot(ax(1),(0:L-1)/wopt.fs,[ones(1,tn1-1)*NaN,freq(ridges(tn1:tn2))',ones(1,L-tn2)*NaN],...
            '-k','LineWidth',2,'DisplayName','Global Maximum curve');
    end
end

%------------------------- Nearest neighbour ------------------------------
if strcmpi(method,'nearest')
    if ~strcmpi(Display,'off')
        fprintf('Extracting the curve by Nearest Neighbour scheme.\n');
        fprintf('The highest peak was found at time %0.3f s and frequency %0.3f Hz (indices %d and %d, respectively).\n',...
            (timax-1)/wopt.fs,Fp(fimax,timax),timax,Ip(fimax,timax));
        fprintf('Tracing the curve forward and backward from point of maximum.\n');
    end
    
    ridges=zeros(L,1)*NaN; ridges(timax)=Ip(fimax,timax);
    for tn=timax+1:tn2
        ci=Ip(1:Np(tn),tn); cq=Qp(1:Np(tn),tn);
        if Np(tn)==1 && cq(1)==0, ridges(tn)=ridges(tn-1);
        else
            [~,mid]=min(abs(ci-ridges(tn-1)));
            ridges(tn)=ci(mid);
        end
    end
    for tn=timax-1:-1:tn1
        ci=Ip(1:Np(tn),tn); cq=Qp(1:Np(tn),tn);
        if Np(tn)==1 && cq(1)==0, ridges(tn)=ridges(tn+1);
        else
            [~,mid]=min(abs(ci-ridges(tn+1)));
            ridges(tn)=ci(mid);
        end
    end
    tfsupp(1,:)=ridges(:)';
    
    if nargout>1
        ec.pfreq=[ones(1,tn1-1)*NaN,freq(ridges(tn1:tn2))',ones(1,L-tn2)*NaN];
        ec.pind=ridges(:)'; %peak indices
        ec.pamp=zeros(1,L)*NaN; for tn=tn1:tn2, ec.pamp(tn)=Qp(Ip(:,tn)==ridges(tn),tn); end
    end
    
    if ~isempty(strfind(Display,'plot'))
        scrsz=get(0,'ScreenSize'); figure('Position',[scrsz(3)/4,scrsz(4)/8,2*scrsz(3)/3,2*scrsz(4)/3]);
        ax=axes('Position',[0.1,0.1,0.8,0.8],'FontSize',16); hold all;
        title(ax(1),'Ridge curve \omega_p(t)/2\pi'); ylabel(ax(1),'Frequency (Hz)'); xlabel(ax(1),'Time (s)');
        plot(ax(1),(0:L-1)/wopt.fs,[ones(1,tn1-1)*NaN,freq(ridges(tn1:tn2))',ones(1,L-tn2)*NaN],...
            '-k','LineWidth',2,'DisplayName','Nearest Neighbour curve');
        plot(ax(1),(timax-1)/wopt.fs,Fp(fimax,timax),'ob','LineWidth',2,'MarkerSize',8,'MarkerFaceColor','r',...
            'DisplayName','Point of maximum TFR amplitude');
    end
end

%----------------------------- Method I -----------------------------------
if length(pars)==1
    if ~strcmpi(Display,'off')
        fprintf('Extracting the curve by I scheme.\n');
    end
    
    if fres==1, logwf=-(pars(1)/2)*(2*pi*dfreq/(Dxi/2)).^2;
    else logwf=-(pars(1)/2)*(dfreq/(Dxi/2)).^2; end
    
    pridges=ecurve1(Np,Ip,Fp,Qp,logwf,timax); %just for comparison
    ridges=pathopt(Np,Ip,Fp,Qp,logwf,[],freq,wopt,Display);
    rdiff=length(find(ridges(tn1:tn2)~=pridges(tn1:tn2)))/(tn2-tn1+1);
    tfsupp(1,:)=ridges(:)';
    
    if nargout>1
        ec.pfreq=[ones(1,tn1-1)*NaN,freq(ridges(tn1:tn2))',ones(1,L-tn2)*NaN];
        ec.pind=ridges(:)'; %peak indices
        ec.pamp=zeros(1,L)*NaN; for tn=tn1:tn2, ec.pamp(tn)=Qp(Ip(:,tn)==ridges(tn),tn); end
        ec.pfreq0=[ones(1,tn1-1)*NaN,freq(pridges(:))',ones(1,L-tn2)*NaN];
        ec.pind0=pridges(:)'; %peak indices
        ec.pamp0=zeros(1,L)*NaN; for tn=tn1:tn2, ec.pamp0(tn)=Qp(Ip(:,tn)==pridges(tn),tn); end
        ec.rdiff=rdiff;
    end
    
    if ~isempty(strfind(Display,'plot'))
        scrsz=get(0,'ScreenSize'); figure('Position',[scrsz(3)/4,scrsz(4)/8,2*scrsz(3)/3,2*scrsz(4)/3]);
        ax=axes('Position',[0.1,0.1,0.8,0.8],'FontSize',16); hold all;
        title(ax(1),'Ridge curve \omega_p(t)/2\pi'); ylabel(ax(1),'Frequency (Hz)'); xlabel(ax(1),'Time (s)');
        plot(ax(1),(0:L-1)/wopt.fs,[ones(1,tn1-1)*NaN,freq(pridges(tn1:tn2))',ones(1,L-tn2)*NaN],...
            'Color',[0.5,0.5,0.5],'LineWidth',1,'DisplayName',sprintf('Without path optimization (discrepancy %0.2f)',rdiff));
        plot(ax(1),(0:L-1)/wopt.fs,[ones(1,tn1-1)*NaN,freq(ridges(tn1:tn2))',ones(1,L-tn2)*NaN],...
            '-k','LineWidth',2,'DisplayName','Optimized over path');
        plot(ax(1),(timax-1)/wopt.fs,Fp(fimax,timax),'ob','LineWidth',2,'MarkerSize',8,'MarkerFaceColor','r',...
            'DisplayName','Point of maximum TFR amplitude');
    end
end


%----------------------------- Method II ----------------------------------
if length(pars)==2 && method==2
    if ~strcmpi(Display,'off')
        fprintf('The highest peak was found at time %0.3f s and frequency %0.3f Hz (indices %d and %d, respectively).\n',...
            (timax-1)/wopt.fs,Fp(fimax,timax),timax,Ip(fimax,timax));
        fprintf('Extracting the curve by II scheme.\n');
        fprintf('Performing first run forward and backward from point of maximum.\n');
        if ~isempty(strfind(Display,'plot'))
            scrsz=get(0,'ScreenSize'); figure('Position',[scrsz(3)/4,scrsz(4)/8,2*scrsz(3)/3,2*scrsz(4)/3]);
            ax=zeros(3,1);
            ax(1)=axes('Position',[0.1,0.6,0.8,0.3],'FontSize',16); hold all;
            ax(2)=axes('Position',[0.1,0.1,0.35,0.35],'FontSize',16); hold all;
            ax(3)=axes('Position',[0.55,0.1,0.35,0.35],'FontSize',16); hold all;
            title(ax(1),'Ridge curve \omega_p(t)/2\pi'); ylabel(ax(1),'Frequency (Hz)'); xlabel(ax(1),'Time (s)');
            plot(ax(1),(0:L-1)/wopt.fs,[ones(1,tn1-1)*NaN,freq(ridges(tn1:tn2))',ones(1,L-tn2)*NaN],':','Color',[0.5,0.5,0.5],'DisplayName','Global Maximum ridges');
            mpt=plot(ax(1),(timax-1)/wopt.fs,Fp(fimax,timax),'ob','LineWidth',2,'MarkerSize',8,'MarkerFaceColor','r',...
                'DisplayName','Point of maximum TFR amplitude');
            ylabel(ax(2),'Frequency (Hz)'); xlabel(ax(2),'Run number');
            title(ax(2),'\langle\omega_p\rangle/2\pi (solid), std[\omega_p]/2\pi (dashed)');
            ylabel(ax(3),'Frequency (Hz)'); xlabel(ax(3),'Run number');
            title(ax(3),'\langle\Delta\omega_p\rangle/2\pi (solid), std[\Delta\omega_p]/2\pi (dashed)');
        end
    end
    
    %Form initial values
    iomp=2*pi*freq(ridges(tn1:tn2));
    if fres==1
        imv=[mean(iomp),mean(iomp.^2),0,(1/4)*(wopt.wp.xi2h-wopt.wp.xi1h)^2];
        nn=(wopt.fs/(Dtau/2))*ones(1,4);
    else
        imv=[mean(log(iomp)),mean(log(iomp).^2),0,(1/4)*(log(wopt.wp.xi2h/wopt.wp.xi1h))^2];
        nn=(wopt.fs/(wopt.wp.ompeak*exp(-mean(log(iomp)))*Dtau/2))*ones(1,4);
    end
    
    %First run
    [ridges,mv]=ecurve2(Np,Ip,Fp,Qp,pars,imv,nn,timax,freq,wopt,Display); nn=nn*2;
    if nargout>1
        ec.pfreq0=[ones(1,tn1-1)*NaN,freq(ridges(tn1:tn2))',ones(1,L-tn2)*NaN];
        ec.pind0=ridges(:)'; %peak indices
        ec.pamp0=zeros(1,L)*NaN; for tn=tn1:tn2, ec.pamp0(tn)=Qp(Ip(:,tn)==ridges(tn),tn); end
        ec.mv0=mv;
        ec.pfreq_iter=[]; ec.pind_iter=[]; ec.pamp_iter=[]; ec.mv_iter=[]; ec.rdiff_iter=[];
    end
    if ~isempty(strfind(Display,'plot'))
        plot(ax(1),(0:L-1)/wopt.fs,[ones(1,tn1-1)*NaN,freq(ridges(tn1:tn2))',ones(1,L-tn2)*NaN],'Color',[0.5,0.5,0.5],'LineWidth',1,'DisplayName','Run 1');
        line1=plot(ax(2),[0,1],[imv(1),mv(1)]/2/pi,'-sk','LineWidth',2,'MarkerSize',6,'MarkerFaceColor','k','DisplayName','\langle\omega_p\rangle/2\pi');
        line2=plot(ax(2),[0,1],[sqrt(imv(2)-imv(1)^2),sqrt(mv(2)-mv(1)^2)]/2/pi,'--ok','LineWidth',2,'MarkerSize',6,'MarkerFaceColor','k','DisplayName','std[\omega_p]/2\pi');
        line3=plot(ax(3),[0,1],[imv(3),mv(3)]/2/pi,'-sk','LineWidth',2,'MarkerSize',6,'MarkerFaceColor','k','DisplayName','\langle\Delta\omega_p\rangle/2\pi');
        line4=plot(ax(3),[0,1],[sqrt(imv(4)-imv(3)^2),sqrt(mv(4)-mv(3)^2)]/2/pi,'--ok','LineWidth',2,'MarkerSize',6,'MarkerFaceColor','k','DisplayName','std[\Delta\omega_p]/2\pi');
    end
    
    %Polishing
    if ~strcmpi(Display,'off'), fprintf('Polishing (run number - discrepancy with previous run): '); end
    pflag=0; itn=1;
    while pflag==0
        pridges=ridges; itn=itn+1;
        [ridges,mv]=ecurve2(Np,Ip,Fp,Qp,pars,mv,nn,timax,freq,wopt,Display); nn=nn*2;
        rdiff=length(find(ridges(tn1:tn2)~=pridges(tn1:tn2)))/(tn2-tn1+1);
        
        if nargout>1
            ec.pfreq_iter=[ec.pfreq_iter;[ones(1,tn1-1)*NaN,freq(ridges(tn1:tn2))',ones(1,L-tn2)*NaN]];
            ec.pind_iter=[ec.pind_iter;ridges(:)']; %peak indices
            cpamp=zeros(1,L)*NaN; for tn=tn1:tn2, cpamp(tn)=Qp(Ip(:,tn)==ridges(tn),tn); end
            ec.pamp_iter=[ec.pamp_iter;cpamp(:)'];
            ec.mv_iter=[ec.mv_iter;mv(:)'];
            ec.rdiff_iter=[ec.rdiff_iter;rdiff];
        end
        
        if ~strcmpi(Display,'off'), fprintf('%d - %0.2f%%; ',itn,100*rdiff); end
        if ~isempty(strfind(Display,'plot'))
            plot(ax(1),(0:L-1)/wopt.fs,[ones(1,tn1-1)*NaN,freq(ridges(tn1:tn2))',ones(1,L-tn2)*NaN],...
                'DisplayName',sprintf('Run %d (discrepancy %0.2f%%)',itn,100*rdiff));
            set(line1,'XData',0:itn,'YData',[get(line1,'YData'),mv(1)/2/pi]);
            set(line2,'XData',0:itn,'YData',[get(line2,'YData'),sqrt(mv(2)-mv(1)^2)/2/pi]);
            set(line3,'XData',0:itn,'YData',[get(line3,'YData'),mv(3)/2/pi]);
            set(line4,'XData',0:itn,'YData',[get(line4,'YData'),sqrt(mv(4)-mv(3)^2)/2/pi]);
        end
        
        if max(abs(ridges(tn1:tn2)-pridges(tn1:tn2)))==0, pflag=1; end
    end
    if ~strcmpi(Display,'off'), fprintf('\n'); end
    
    %Refine path
    if ~strcmpi(Display,'off'), fprintf('Refining the path: '); end
    pridges=ridges;
    if fres==1
        logw1=-(pars(1)/2)*((2*pi*dfreq-mv(3))/sqrt(max([mv(4)-mv(3)^2,10^(-28)]))).^2;
        logw2=-(pars(2)/2)*((2*pi*freq-mv(1))/sqrt(max([mv(2)-mv(1)^2,10^(-28)]))).^2;
    else
        logw1=-(pars(1)/2)*((dfreq-mv(3))/sqrt(max([mv(4)-mv(3)^2,10^(-28)]))).^2;
        logw2=-(pars(2)/2)*((log(2*pi*freq)-mv(1))/sqrt(max([mv(2)-mv(1)^2,10^(-28)]))).^2;
    end 
    ridges=pathopt(Np,Ip,Fp,Qp,logw1,logw2,freq,wopt,Display);
    rdiff=length(find(ridges(tn1:tn2)~=pridges(tn1:tn2)))/(tn2-tn1+1);
    if ~strcmpi(Display,'off'), fprintf('discrepancy %0.2f%%\n',100*rdiff); end
    tfsupp(1,:)=ridges(:)';
    
    if nargout>1
        ec.pfreq=[ones(1,tn1-1)*NaN,freq(ridges(tn1:tn2))',ones(1,L-tn2)*NaN];
        ec.pind=ridges(:)'; %peak indices
        ec.pamp=zeros(1,L)*NaN; for tn=tn1:tn2, ec.pamp(tn)=Qp(Ip(:,tn)==ridges(tn),tn); end
        ec.rdiff=length(find(ridges(tn1:tn2)~=pridges(tn1:tn2)))/(tn2-tn1+1);
    end
    
    if ~isempty(strfind(Display,'plot'))
        plot(ax(1),(0:L-1)/wopt.fs,[ones(1,tn1-1)*NaN,freq(ridges(tn1:tn2))',ones(1,L-tn2)*NaN],'-k','LineWidth',2,...
            'DisplayName',sprintf('Refined path (discrepancy %0.2f%%)',100*rdiff));
        uistack(mpt,'top');
        
        %Change plot if the resolution is logarithmic
        if fres==2
            set(line1,'YData',exp(2*pi*get(line1,'YData'))/2/pi,'DisplayName','exp(\langlelog\omega_p\rangle)/2\pi');
            set(line2,'YData',exp(2*pi*get(line2,'YData'))-1,'DisplayName','exp(std[log\omega_p])-1');
            set(line3,'YData',exp(2*pi*get(line3,'YData')),'DisplayName','exp(\langle\Deltalog\omega_p\rangle)');
            set(line4,'YData',exp(2*pi*get(line4,'YData'))-1,'DisplayName','exp(std[\Deltalog\omega_p])-1');
            ylabel(ax(2),'Frequency (Hz)'); xlabel(ax(2),'Run number');
            title(ax(2),'exp(\langlelog\omega_p\rangle)/2\pi (solid), exp(std[log\omega_p])-1 (dashed)');
            ylabel(ax(3),'Frequency Ratio'); xlabel(ax(3),'Run number');
            title(ax(3),'exp(\langle\Deltalog\omega_p\rangle) (solid), exp(std[\Deltalog\omega_p])-1 (dashed)');
        end
    end
end

%----------------------------- Method III ---------------------------------
if length(pars)==2 && method==3
    if ~strcmpi(Display,'off')
        fprintf('The highest peak was found at time %0.3f s and frequency %0.3f Hz (indices %d and %d, respectively).\n',...
            (timax-1)/wopt.fs,Fp(fimax,timax),timax,Ip(fimax,timax));
        fprintf('Extracting the curve by III scheme.\n');
        fprintf('Performing first run forward and backward from point of maximum.\n');
        if ~isempty(strfind(Display,'plot'))
            if fres==1, dfreq_p=dfreq; else dfreq_p=exp(dfreq); end
            scrsz=get(0,'ScreenSize'); figure('Position',[scrsz(3)/4,scrsz(4)/8,2*scrsz(3)/3,2*scrsz(4)/3]);
            ax=zeros(3,1);
            ax(1)=axes('Position',[0.1,0.6,0.8,0.3],'FontSize',16); hold all;
            ax(2)=axes('Position',[0.1,0.1,0.35,0.35],'FontSize',16,'XLim',[dfreq_p(1),dfreq_p(end)]); hold all;
            ax(3)=axes('Position',[0.55,0.1,0.35,0.35],'FontSize',16,'XLim',[freq(1),freq(end)]); hold all;
            title(ax(2),'P_1(\Delta\nu) (normalized on max)'); title(ax(3),'P_2(\nu) (normalized on max)');
            xlabel(ax(2),'\Delta\nu/2\pi (Hz)'); xlabel(ax(3),'\nu/2\pi (Hz)');
            title(ax(1),'Ridge curve \omega_p(t)/2\pi'); ylabel(ax(1),'Frequency (Hz)'); xlabel(ax(1),'Time (s)');
            plot(ax(1),(0:L-1)/wopt.fs,[ones(1,tn1-1)*NaN,freq(ridges(tn1:tn2))',ones(1,L-tn2)*NaN],':','Color',[0.5,0.5,0.5],'DisplayName','Global Maximum ridges');
            mpt=plot(ax(1),(timax-1)/wopt.fs,Fp(fimax,timax),'ob','LineWidth',2,'MarkerSize',8,'MarkerFaceColor','r',...
                'DisplayName','Point of maximum TFR amplitude');
        end
    end
    
    %Form initial values
    iomp=2*pi*freq(ridges(tn1:tn2));
    if fres==1
        imv=[mean(iomp),mean(iomp.^2)]; nn=(wopt.fs/(Dtau/2))*ones(1,2);
        iP1=exp(-(1/2)*(2*pi*dfreq/(Dxi/2)).^2);
        iP2=exp(-(1/2)*((2*pi*freq-imv(1))/sqrt(max([imv(2)-imv(1)^2,10^(-28)]))).^2);
    else
        imv=[mean(log(iomp)),mean(log(iomp).^2)]; nn=(wopt.fs/(wopt.wp.ompeak*exp(-mean(log(iomp)))*Dtau/2))*ones(1,2);
        iP1=exp(-(1/2)*(dfreq/(Dxi/2)).^2);
        iP2=exp(-(1/2)*((log(2*pi*freq)-imv(1))/sqrt(max([imv(2)-imv(1)^2,10^(-28)]))).^2);
    end
    
    %First run
    [ridges,P1,P2]=ecurve3(Np,Ip,Fp,Qp,pars,nn(1)*iP1(:)/max(iP1),nn(2)*iP2(:)/max(iP2),timax,freq,wopt,Display); nn=nn*2;
    if nargout>1
        ec.pfreq0=[ones(1,tn1-1)*NaN,freq(ridges(tn1:tn2))',ones(1,L-tn2)*NaN];
        ec.pind0=ridges(:)'; %peak indices
        ec.pamp0=zeros(1,L)*NaN; for tn=tn1:tn2, ec.pamp0(tn)=Qp(Ip(:,tn)==ridges(tn),tn); end
        ec.P10=P1(:); ec.P20=P2(:);
        ec.pfreq_iter=[]; ec.pind_iter=[]; ec.pamp_iter=[]; ec.P1_iter=[]; ec.P2_iter=[]; ec.rdiff_iter=[];
    end
    if ~isempty(strfind(Display,'plot'))
        plot(ax(1),(0:L-1)/wopt.fs,[ones(1,tn1-1)*NaN,freq(ridges(tn1:tn2))',ones(1,L-tn2)*NaN],...
            'Color',[0.5,0.5,0.5],'LineWidth',1,'DisplayName','Run 1');
        plot(ax(2),dfreq_p,iP1/max(iP1),'--k','DisplayName','Initial');
        plot(ax(3),freq,iP2/max(iP2),'--k','DisplayName','Initial');
        plot(ax(2),dfreq_p,P1/max(P1),'Color',[0.5,0.5,0.5],'DisplayName','Run 1');
        plot(ax(3),freq,P2/max(P2),'Color',[0.5,0.5,0.5],'DisplayName','Run 1');
    end
    
    %Polishing
    if ~strcmpi(Display,'off'), fprintf('Polishing (run number - discrepancy with previous run): '); end
    pflag=0; itn=1;
    while pflag==0
        pridges=ridges; itn=itn+1;
        [ridges,P1,P2]=ecurve3(Np,Ip,Fp,Qp,pars,nn(1)*P1(:)/max(P1),nn(2)*P2(:)/max(P2),timax,freq,wopt,Display); nn=nn*2;
        rdiff=length(find(ridges(tn1:tn2)~=pridges(tn1:tn2)))/(tn2-tn1+1);
        
        if ~strcmpi(Display,'off'), fprintf('%d - %0.2f%%; ',itn,100*rdiff); end
        if ~isempty(strfind(Display,'plot'))
            plot(ax(1),(0:L-1)/wopt.fs,[ones(1,tn1-1)*NaN,freq(ridges(tn1:tn2))',ones(1,L-tn2)*NaN],...
                'DisplayName',sprintf('Run %d (discrepancy %0.2f%%)',itn,100*rdiff));
            line1=plot(ax(2),dfreq_p,P1/max(P1),'DisplayName',sprintf('Run %d',itn));
            line2=plot(ax(3),freq,P2/max(P2),'DisplayName',sprintf('Run %d',itn));
        end
        
        if nargout>1
            ec.pfreq_iter=[ones(1,tn1-1)*NaN,freq(ridges(tn1:tn2))',ones(1,L-tn2)*NaN];
            ec.pind_iter=[ec.pind_iter;ridges(:)']; %peak indices
            cpamp=zeros(1,L)*NaN; for tn=tn1:tn2, cpamp(tn)=Qp(Ip(:,tn)==ridges(tn),tn); end
            ec.pamp_iter=[ec.pamp_iter;cpamp(:)'];
            ec.P1_iter=[ec.P1_iter;P1(:)']; ec.P2_iter=[ec.P2_iter;P2(:)'];
            ec.rdiff_iter=[ec.rdiff_iter;length(find(ridges(tn1:tn2)~=pridges(tn1:tn2)))/(tn2-tn1+1)];
        end
        
        if max(abs(ridges(tn1:tn2)-pridges(tn1:tn2)))==0, pflag=1; end
    end
    if ~strcmpi(Display,'off'), fprintf('\n'); end
    
    %Refine path
    if ~strcmpi(Display,'off'), fprintf('Refining the path: '); end
    pridges=ridges;
    ridges=pathopt(Np,Ip,Fp,Qp,pars(1)*log(P1(:)),pars(2)*log(P2(:)),freq,wopt,Display);
    rdiff=length(find(ridges(tn1:tn2)~=pridges(tn1:tn2)))/(tn2-tn1+1);
    if ~strcmpi(Display,'off'), fprintf('discrepancy %0.2f%%\n',100*rdiff); end
    tfsupp(1,:)=ridges(:)';
    
    if nargout>1
        ec.pfreq=[ones(1,tn1-1)*NaN,freq(ridges(tn1:tn2))',ones(1,L-tn2)*NaN];
        ec.pind=ridges(:)'; %peak indices
        ec.pamp=zeros(1,L)*NaN; for tn=tn1:tn2, ec.pamp(tn)=Qp(Ip(:,tn)==ridges(tn),tn); end
        ec.rdiff=length(find(ridges(tn1:tn2)~=pridges(tn1:tn2)))/(tn2-tn1+1);
    end
    if ~isempty(strfind(Display,'plot'))
        set([line1,line2],'Color','k','LineWidth',2);
        plot(ax(1),(0:L-1)/wopt.fs,[ones(1,tn1-1)*NaN,freq(ridges(tn1:tn2))',ones(1,L-tn2)*NaN],'-k','LineWidth',2,...
            'DisplayName',sprintf('Refined path (discrepancy %0.2f%%)',100*rdiff));
        uistack(mpt,'top');
        
        %Change plot if the resolution is logarithmic
        if fres==2
            xlabel(ax(2),'1+\Delta\nu/\nu');
            title(ax(2),'P_1(1+\Delta\nu/\nu) (normalized on max)');
            set(ax(2:3),'XScale','log');
        end
    end
    
end


%//////////////////////////////////////////////////////////////////////////
%Extract the time-frequency support around the ridge points
TFR=abs(TFR);
for tn=tn1:tn2
    cs=abs(TFR(:,tn)); cpeak=tfsupp(1,tn);
    iup=[]; idown=[];
    if cpeak<NF-1, iup=cpeak+find(cs(cpeak+1:end)==0,1,'first'); end
    if cpeak>2, idown=cpeak-find(cs(cpeak-1:-1:1)==0,1,'first'); end
    iup=min([iup,NF]); idown=max([idown,1]);
    tfsupp(2,tn)=idown; tfsupp(3,tn)=iup;
end
%//////////////////////////////////////////////////////////////////////////

%Transform to frequencies
tfsupp(:,tn1:tn2)=freq(tfsupp(:,tn1:tn2));
if nargout>1, varargout{1}=ec; end

%Plot (if needed)
if strcmpi(PlotMode,'on'), plotfinal(tfsupp,TFR,freq,wopt); end

end


%==========================================================================
%==================== Curve extraction schemes ============================
%==========================================================================


%================== Path optimization (and scheme I) ======================
function ridges=pathopt(Np,Ip,Fp,Qp,logw1,logw2,freq,wopt,Display)

[Mp,L]=size(Fp); tn1=find(Np>0,1,'first'); tn2=find(Np>0,1,'last');

%Different weighting parameters
aw1=0; aw2=0;
if ~isempty(logw1)
    aw1=1; if ~isa(logw1,'function_handle'); aw1=2; logw1=logw1(:); NF=round((length(logw1)+1)/2); end
end
if ~isempty(logw2)
    aw2=1; if ~isa(logw2,'function_handle'), aw2=2; logw2=logw2(:); NF=length(logw2); end
end

%Path optimization
q=zeros(Mp,L)*NaN; U=zeros(Mp,L)*NaN;
q(1:Np(tn1),tn1)=0; U(1:Np(tn1),tn1)=log(Qp(1:Np(tn1),tn1));
if aw2>0
    if aw2==1, nU=logw2(2*pi*Fp(1:Np(tn1),tn1)); else nU=logw2(Ip(1:Np(tn1),tn1)); end
    U(1:Np(tn1),tn1)=U(1:Np(tn1),tn1)+nU(:);
end
if aw1~=1 && aw2~=1 %JIT accelerated
    for tn=tn1+1:tn2
        for pn=1:Np(tn)
            if aw1==2, clogwf=logw1(NF+Ip(pn,tn)-Ip(1:Np(tn-1),tn-1));
            else clogwf=zeros(Np(tn-1),1); end
            if aw2==2, clogwf=clogwf+logw2(Ip(pn,tn)); end
            [U(pn,tn),q(pn,tn)]=max(log(Qp(pn,tn))+clogwf(:)+U(1:Np(tn-1),tn-1));
        end
    end
else %not accelerated by JIT
    for tn=tn1+1:tn2
        for pn=1:Np(tn)
            if aw1==1, clogwf=logw1(2*pi*Fp(pn,tn),2*pi*Fp(1:Np(tn-1),tn-1));
            elseif aw1==2, clogwf=logw1(NF+Ip(pn,tn)-Ip(1:Np(tn-1),tn-1));
            else clogwf=zeros(Np(tn-1),1);
            end
            if aw2==1, clogwf=clogwf+logw2(2*pi*Fp(pn,tn));
            elseif aw2==2, clogwf=clogwf+logw2(Ip(pn,tn));
            end
            [U(pn,tn),q(pn,tn)]=max(log(Qp(pn,tn))+clogwf(:)+U(1:Np(tn-1),tn-1));
        end
    end
end

%Recover indices
ridges=zeros(L,1)*NaN;
[~,ridges(tn2)]=max(U(:,tn2));
for tn=tn2-1:-1:tn1, ridges(tn)=q(ridges(tn+1),tn+1); end
lind=sub2ind(size(Ip),ridges(tn1:tn2),(tn1:tn2)');
ridges=[ones(tn1-1,1)*NaN;Ip(lind);ones(L-tn2,1)*NaN];

%Plot if needed
%{
if ~isempty(strfind(Display,'plot++'))
    figure; axes('FontSize',16,'Box','on'); hold all;
    rlines=zeros(Np(tn2),1);
    for pn=1:Np(tn2)
        cridges=zeros(L,1)*NaN; cridges(tn2)=pn;
        for tn=tn2-1:-1:tn1
            cridges(tn)=q(cridges(tn+1),tn+1);
        end
        lind=sub2ind(size(Ip),cridges(tn1:tn2),(tn1:tn2)');
        cridges=[ones(tn1-1,1)*NaN;Ip(lind);ones(L-tn2,1)*NaN];
        crfreq=[ones(tn1-1,1);freq(cridges(tn1:tn2));ones(L-tn2,1)];
        rlines(pn)=plot((0:L-1)/wopt.fs,crfreq,'DisplayName',['U=',num2str(U(pn,tn2))]);
    end
    [~,imax]=max(U(:,tn2)); uistack(rlines(imax),'top');
    set(rlines(imax),'Color','k','LineWidth',2);
    title('Path to each of the last peaks maximizing given functional');
    ylabel('Frequency (Hz)'); xlabel('Time (s)');
end
%}

end


%========================== Scheme I simple ===============================
function ridges=ecurve1(Np,Ip,Fp,Qp,logwf,timax)

[Mp,L]=size(Fp); tn1=find(Np>0,1,'first'); tn2=find(Np>0,1,'last');

%Different weighting parameters
awf=1; if ~isa(logwf,'function_handle'), awf=2; NF=round((length(logwf)+1)/2); end

%Simple curve following
ridges=zeros(L,1)*NaN; [~,ridges(timax)]=max(Qp(:,timax));
if awf==2 %JIT accelerated
    for tn=timax+1:tn2
        ci=Ip(1:Np(tn),tn); cq=Qp(1:Np(tn),tn);
        clogwf=logwf(NF+ci-Ip(ridges(tn-1),tn-1));
        [~,ridges(tn)]=max(log(cq(:))+clogwf(:));
    end
    for tn=timax-1:-1:tn1
        ci=Ip(1:Np(tn),tn); cq=Qp(1:Np(tn),tn);
        clogwf=logwf(NF+Ip(ridges(tn+1),tn+1)-ci);
        [~,ridges(tn)]=max(log(cq(:))+clogwf(:));
    end
else %not accelerated by JIT
    for tn=timax+1:tn2
        cxi=2*pi*Fp(1:Np(tn),tn); cq=Qp(1:Np(tn),tn);
        clogwf=logwf(cxi,2*pi*Fp(ridges(tn-1),tn-1));
        [~,ridges(tn)]=max(log(cq(:))+clogwf(:));
    end
    for tn=timax-1:-1:tn1
        cxi=2*pi*Fp(1:Np(tn),tn); cq=Qp(1:Np(tn),tn);
        clogwf=logwf(Fp(ridges(tn+1),tn+1),cxi);
        [~,ridges(tn)]=max(log(cq(:))+clogwf(:));
    end
end

%Recover indices
lind=sub2ind(size(Ip),ridges(tn1:tn2),(tn1:tn2)');
ridges=[ones(tn1-1,1)*NaN;Ip(lind);ones(L-tn2,1)*NaN];

end


%============================= Scheme II ==================================
function [ridges,mv]=ecurve2(Np,Ip,Fp,Qp,pars,imv,nn,timax,freq,wopt,Display)

[Mp,L]=size(Fp); tn1=find(Np>0,1,'first'); tn2=find(Np>0,1,'last');
dflag=0; if ~isempty(strfind(Display,'plot+')), allmv=zeros(4,L); allmv(:,timax)=imv; dflag=1; end
if min(freq)<=0 || std(diff(freq))<std(diff(log(freq))), Wp=2*pi*Fp; fres=1;
else Wp=log(2*pi*Fp); fres=2; end

%Curve extraction
ridges=zeros(L,1)*NaN;
[~,ridges(timax)]=max(log(Qp(1:Np(timax),timax))...
    -(pars(2)/2)*((Wp(1:Np(timax),timax)-imv(1))/sqrt(max([imv(2)-imv(1)^2,10^(-28)]))).^2);

mv=imv; %forward
for tn=timax+1:tn2
    cxi=Wp(1:Np(tn),tn); cq=Qp(1:Np(tn),tn); pxi=Wp(ridges(tn-1),tn-1);
    clogwf1=-(pars(1)/2)*((cxi-pxi-mv(3))/sqrt(max([mv(4)-mv(3)^2,10^(-28)]))).^2;
    clogwf2=-(pars(2)/2)*((cxi-mv(1))/sqrt(max([mv(2)-mv(1)^2,10^(-28)]))).^2;
    [~,ridges(tn)]=max(log(cq(:))+clogwf1(:)+clogwf2(:)); cxi=cxi(ridges(tn));
    %Updating
    mv(1)=((nn(1)+tn-timax)*mv(1)+cxi)/(nn(1)+tn-timax+1);
    mv(2)=((nn(2)+tn-timax)*mv(2)+cxi.^2)/(nn(2)+tn-timax+1);
    mv(3)=((nn(3)+tn-timax)*mv(3)+(cxi-pxi))/(nn(3)+tn-timax+1);
    mv(4)=((nn(4)+tn-timax)*mv(4)+(cxi-pxi)^2)/(nn(4)+tn-timax+1);
    if dflag==1, allmv(:,tn)=mv(:); end
end
mvF=mv;

mv=imv; %backward
for tn=timax-1:-1:tn1
    cxi=Wp(1:Np(tn),tn); cq=Qp(1:Np(tn),tn); pxi=Wp(ridges(tn+1),tn+1);
    clogwf1=-(pars(1)/2)*((pxi-cxi-mv(3))/sqrt(max([mv(4)-mv(3)^2,10^(-28)]))).^2;
    clogwf2=-(pars(2)/2)*((cxi-mv(1))/sqrt(max([mv(2)-mv(1)^2,10^(-28)]))).^2;
    [~,ridges(tn)]=max(log(cq(:))+clogwf1(:)+clogwf2(:)); cxi=cxi(ridges(tn));
    %Updating
    mv(1)=((nn(1)+timax-tn)*mv(1)+cxi)/(nn(1)+timax-tn+1);
    mv(2)=((nn(2)+timax-tn)*mv(2)+cxi.^2)/(nn(2)+timax-tn+1);
    mv(3)=((nn(3)+timax-tn)*mv(3)+(pxi-cxi))/(nn(3)+timax-tn+1);
    mv(4)=((nn(4)+timax-tn)*mv(4)+(pxi-cxi)^2)/(nn(4)+timax-tn+1);
    if dflag==1, allmv(:,tn)=mv(:); end
end
mvB=mv;

%Recover indices
lind=sub2ind(size(Ip),ridges(tn1:tn2),(tn1:tn2)');
ridges=[ones(tn1-1,1)*NaN;Ip(lind);ones(L-tn2,1)*NaN];

%Update the time-averages
mv=((timax-tn1)*mvB+(tn2-timax)*mvF)/(tn2-tn1);

if dflag==1
    scrsz=get(0,'ScreenSize');
    figure('Position',[scrsz(3)/4,scrsz(4)/8,2*scrsz(3)/3,3*scrsz(4)/4]);
    ax=zeros(3,1);
    ax(1)=axes('Position',[0.1,0.7,0.8,0.2],'FontSize',16); hold all;
    ax(2)=axes('Position',[0.1,0.4,0.8,0.2],'FontSize',16); hold all;
    ax(3)=axes('Position',[0.1,0.1,0.8,0.2],'FontSize',16); hold all;
    tt=(0:L-1)/wopt.fs; rfreq=[ones(tn1-1,1)*NaN;freq(ridges(tn1:tn2));ones(L-tn2,1)*NaN];
    axes(ax(1)); plot(tt,rfreq,'-k');
    plot(tt(timax),rfreq(timax),'ok','LineWidth',1,'MarkerSize',6,'MarkerFaceColor','r',...
                'DisplayName','Point of maximum TFR amplitude');
    title('Extracted ridge curve \omega_p(t)/2\pi'); ylabel('Frequency (Hz)');
    
    if fres==1 %for WFT
        axes(ax(2)); plot(tt,allmv(1,:)/2/pi,'-k','DisplayName','\langle\omega_p\rangle/2\pi');
        plot(tt,sqrt(allmv(2,:)-allmv(1,:).^2)/2/pi,'--k','DisplayName','std[\omega_p]/2\pi');
        plot(tt(timax),allmv(1,timax)/2/pi,'ok','LineWidth',1,'MarkerSize',6,'MarkerFaceColor','r',...
            'DisplayName','Point of maximum TFR amplitude');
        plot(tt(timax),sqrt(allmv(2,timax)-allmv(1,timax)^2)/2/pi,'ok','LineWidth',1,'MarkerSize',6,'MarkerFaceColor','r',...
            'DisplayName','Point of maximum TFR amplitude');
        title('\langle\omega_p\rangle/2\pi (solid), std[\omega_p]/2\pi (dashed)'); ylabel('Frequency (Hz)');
        axes(ax(3)); plot(tt,allmv(3,:)/2/pi,'-k','DisplayName','\langle\Delta\omega_p\rangle/2\pi');
        plot(tt,sqrt(allmv(4,:)-allmv(3,:).^2)/2/pi,'--k','DisplayName','std[\Delta\omega_p]/2\pi');
        plot(tt(timax),allmv(3,timax)/2/pi,'ok','LineWidth',1,'MarkerSize',6,'MarkerFaceColor','r',...
            'DisplayName','Point of maximum TFR amplitude');
        plot(tt(timax),sqrt(allmv(4,timax)-allmv(3,timax)^2)/2/pi,'ok','LineWidth',1,'MarkerSize',6,'MarkerFaceColor','r',...
            'DisplayName','Point of maximum TFR amplitude');
        title('\langle\Delta\omega_p\rangle/2\pi (solid), std[\Delta\omega_p]/2\pi (dashed)'); ylabel('Frequency (Hz)');
        xlabel('Time (s)');
    else %for WT
        axes(ax(2)); plot(tt,exp(allmv(1,:))/2/pi,'-k','DisplayName','exp(\langlelog\omega_p\rangle)/2\pi');
        plot(tt,exp(sqrt(allmv(2,:)-allmv(1,:).^2))-1,'--k','DisplayName','exp(std[log\omega_p])-1');
        plot(tt(timax),exp(allmv(1,timax))/2/pi,'ok','LineWidth',1,'MarkerSize',6,'MarkerFaceColor','r',...
            'DisplayName','Point of maximum TFR amplitude');
        plot(tt(timax),exp(sqrt(allmv(2,timax)-allmv(1,timax)^2))-1,'ok','LineWidth',1,'MarkerSize',6,'MarkerFaceColor','r',...
            'DisplayName','Point of maximum TFR amplitude');
        title('exp(\langlelog\omega_p\rangle)/2\pi (solid), exp(std[log\omega_p])-1 (dashed)'); ylabel('Frequency (Hz)');
        axes(ax(3)); plot(tt,exp(allmv(3,:)),'-k','DisplayName','exp(\langle\Deltalog\omega_p\rangle)');
        plot(tt,exp(sqrt(allmv(4,:)-allmv(3,:).^2))-1,'--k','DisplayName','exp(std[\Deltalog\omega_p])-1');
        plot(tt(timax),exp(allmv(3,timax)),'ok','LineWidth',1,'MarkerSize',6,'MarkerFaceColor','r',...
            'DisplayName','Point of maximum TFR amplitude');
        plot(tt(timax),exp(sqrt(allmv(4,timax)-allmv(3,timax)^2))-1,'ok','LineWidth',1,'MarkerSize',6,'MarkerFaceColor','r',...
            'DisplayName','Point of maximum TFR amplitude');
        title('exp(\langle\Deltalog\omega_p\rangle) (solid), exp(std[\Deltalog\omega_p])-1 (dashed)'); ylabel('Frequency Ratio');
        xlabel('Time (s)');
    end
end
    
end


%============================= Scheme III =================================
function [ridges,P1,P2]=ecurve3(Np,Ip,Fp,Qp,pars,iP1,iP2,timax,freq,wopt,Display)

[Mp,L]=size(Fp); NF=length(freq);
tn1=find(Np>0,1,'first'); tn2=find(Np>0,1,'last');
if min(freq)<=0 || std(diff(freq))<std(diff(log(freq)))
    fres=1; Dxi=wopt.wp.xi2h-wopt.wp.xi1h;
    Wp=2*pi*Fp; omg=2*pi*freq;
else
    fres=2; Dxi=log(wopt.wp.xi2h/wopt.wp.xi1h);
    Wp=log(2*pi*Fp); omg=log(2*pi*freq);
end
domg=[omg(1)-omg(end:-1:2);omg(end)-omg(end:-1:1)];

dflag=0;
if ~isempty(strfind(Display,'plot+'))
    dflag=1;
    allP1=zeros(2*NF-1,L); allP1(:,timax)=iP1(:)/max(iP1);
    allP2=zeros(NF,L); allP2(:,timax)=iP2(:)/max(iP2);
end

%Curve extraction
ridges=zeros(L,1)*NaN;
[~,ridges(timax)]=max(log(Qp(1:Np(timax),timax))...
    +pars(2)*log(iP2(Ip(1:Np(timax),timax))));

P1=iP1; P2=iP2; %forward
for tn=timax+1:tn2
    cxi=Wp(1:Np(tn),tn); ci=Ip(1:Np(tn),tn); cq=Qp(1:Np(tn),tn);
    pid=Ip(ridges(tn-1),tn-1); pxi=Wp(ridges(tn-1),tn-1);
    logw1=pars(1)*log(P1(NF+ci-pid)); logw2=pars(2)*log(P2(ci));
    [~,ridges(tn)]=max(log(cq(:))+logw1(:)+logw2(:)); cxi=cxi(ridges(tn));
    %Updating
    P1=P1+exp(-(1/2)*((domg-(cxi-pxi))/(Dxi/16)).^2);
    P2=P2+exp(-(1/2)*((omg-cxi)/(Dxi/16)).^2);
    if dflag==1, allP1(:,tn)=P1(:)/max(P1); allP2(:,tn)=P2(:)/max(P2); end
end
P1F=P1; P2F=P2;

P1=iP1; P2=iP2; %backward
for tn=timax-1:-1:tn1
    cxi=Wp(1:Np(tn),tn); ci=Ip(1:Np(tn),tn); cq=Qp(1:Np(tn),tn);
    pid=Ip(ridges(tn+1),tn+1); pxi=Wp(ridges(tn+1),tn+1);
    logw1=pars(1)*log(P1(NF+pid-ci)); logw2=pars(2)*log(P2(ci));
    [~,ridges(tn)]=max(log(cq(:))+logw1(:)+logw2(:)); cxi=cxi(ridges(tn));
    %Updating
    P1=P1+exp(-(1/2)*((domg-(pxi-cxi))/(Dxi/16)).^2);
    P2=P2+exp(-(1/2)*((omg-cxi)/(Dxi/16)).^2);
    if dflag==1, allP1(:,tn)=P1(:)/max(P1); allP2(:,tn)=P2(:)/max(P2); end
end
P1B=P1; P2B=P2;

%Recover indices
lind=sub2ind(size(Ip),ridges(tn1:tn2),(tn1:tn2)');
ridges=[ones(tn1-1,1)*NaN;Ip(lind);ones(L-tn2,1)*NaN];

%Update the time-averages
P1=((timax-tn1)*P1B+(tn2-timax)*P1F)/(tn2-tn1);
P2=((timax-tn1)*P2B+(tn2-timax)*P2F)/(tn2-tn1);

if dflag==1
    scrsz=get(0,'ScreenSize');
    figure('Position',[scrsz(3)/4,scrsz(4)/8,2*scrsz(3)/3,3*scrsz(4)/4]);
    ax=zeros(3,1);
    ax(1)=axes('Position',[0.1,0.7,0.8,0.2],'FontSize',16); hold all;
    ax(2)=axes('Position',[0.1,0.4,0.8,0.2],'FontSize',16); hold all;
    ax(3)=axes('Position',[0.1,0.1,0.8,0.2],'FontSize',16); hold all;
    tt=(0:L-1)/wopt.fs; rfreq=[ones(tn1-1,1)*NaN;freq(ridges(tn1:tn2));ones(L-tn2,1)*NaN];
    axes(ax(1)); plot(tt,rfreq,'-k');
    plot(tt(timax),rfreq(timax),'ok','LineWidth',1,'MarkerSize',6,'MarkerFaceColor','r',...
                'DisplayName','Point of maximum TFR amplitude');
    title('Extracted ridge curve \omega_p(t)/2\pi'); ylabel('Frequency (Hz)');
    
    if fres==1 %for WFT
        axes(ax(2)); surf(tt,domg/2/pi,allP1,'LineStyle','none'); view(2);
        line(tt(timax)*ones(2*NF-1,1),domg/2/pi,allP1(:,timax),'Color','r','LineWidth',1);
        line(tt(timax)*ones(2*NF-1,1),domg/2/pi,ones(2*NF-1,1),'Color','r','LineWidth',1);
        line([tt(timax),tt(timax)],[domg(1),domg(1)]/2/pi,[allP1(1,timax),1],'Color','r','LineWidth',1);
        line([tt(timax),tt(timax)],[domg(end),domg(end)]/2/pi,[allP1(end,timax),1],'Color','r','LineWidth',1);
        title('Evolution of P_1(\Delta\nu) (normalized on max)'); ylabel('\Delta\nu/2\pi (Hz)');
        xlim([tt(1),tt(end)]); ylim([domg(1)/2/pi,domg(end)/2/pi]); zlim([0,1]);
        axes(ax(3)); surf(tt,omg/2/pi,allP2,'LineStyle','none'); view(2);
        line(tt(timax)*ones(NF,1),omg/2/pi,allP2(:,timax),'Color','r','LineWidth',1);
        line(tt(timax)*ones(NF,1),omg/2/pi,ones(NF,1),'Color','r','LineWidth',1);
        line([tt(timax),tt(timax)],[omg(1),omg(1)]/2/pi,[allP2(1,timax),1],'Color','r','LineWidth',1);
        line([tt(timax),tt(timax)],[omg(end),omg(end)]/2/pi,[allP2(end,timax),1],'Color','r','LineWidth',1);
        title('Evolution of P_2(\nu) (normalized on max)'); ylabel('\nu/2\pi (Hz)');
        xlim([tt(1),tt(end)]); ylim([omg(1)/2/pi,omg(end)/2/pi]); zlim([0,1]);
        xlabel('Time (s)');
    else %for WT
        axes(ax(2)); surf(tt,exp(domg),allP1,'LineStyle','none'); view(2);
        line(tt(timax)*ones(2*NF-1,1),exp(domg),allP1(:,timax),'Color','r','LineWidth',1);
        line(tt(timax)*ones(2*NF-1,1),exp(domg),ones(2*NF-1,1),'Color','r','LineWidth',1);
        line([tt(timax),tt(timax)],exp([domg(1),domg(1)]),[allP1(1,timax),1],'Color','r','LineWidth',1);
        line([tt(timax),tt(timax)],exp([domg(end),domg(end)]),[allP1(end,timax),1],'Color','r','LineWidth',1);
        title('Evolution of P_1(1+\Delta\nu/\nu) (normalized on max)'); ylabel('1+\Delta\nu/\nu');
        xlim([tt(1),tt(end)]); ylim([exp(domg(1)),exp(domg(end))]); zlim([0,1]);
        axes(ax(3)); surf(tt,exp(omg)/2/pi,allP2,'LineStyle','none'); view(2);
        line(tt(timax)*ones(NF,1),exp(omg)/2/pi,allP2(:,timax),'Color','r','LineWidth',1);
        line(tt(timax)*ones(NF,1),exp(omg)/2/pi,ones(NF,timax),'Color','r','LineWidth',1);
        line([tt(timax),tt(timax)],exp([omg(1),omg(1)])/2/pi,[allP2(1,timax),1],'Color','r','LineWidth',1);
        line([tt(timax),tt(timax)],exp([omg(end),omg(end)])/2/pi,[allP2(end,timax),1],'Color','r','LineWidth',1);
        title('Evolution of P_2(\nu) (normalized on max)'); ylabel('\nu/2\pi (Hz)');
        xlim([tt(1),tt(end)]); ylim([exp(omg(1))/2/pi,exp(omg(end))/2/pi]); zlim([0,1]);
        xlabel('Time (s)');
        set(ax(2:3),'YScale','log');
    end
        
end

end


%======================== Function for plotting ===========================
function plotfinal(tfsupp,TFR,freq,wopt)

[NF,L]=size(TFR); scrsz=get(0,'ScreenSize');
vmax=max(abs(TFR(:))); tt=(0:(L-1))/wopt.fs;
if min(freq)<=0 || std(diff(freq))<std(diff(log(freq)))
    fres=1; fstep=mean(diff(freq)); eind=1+floor(0.5+(tfsupp-freq(1))/fstep);
else
    fres=2; fstep=mean(diff(log(freq))); eind=1+floor(0.5+log(tfsupp/freq(1))/fstep);
end
eind(eind<1)=1; eind(eind>NF)=NF;

figure('Position',[scrsz(3)/4,scrsz(4)/4,scrsz(3)/2,2*scrsz(4)/3]);
axes('Position',[0.1,0.1,0.8,0.75],'NextPlot','Add','FontSize',18,...
    'XLim',[tt(1),tt(end)],'YLim',[freq(1),freq(end)],'ZLim',[0,vmax]);
xlabel('Time (s)'); ylabel('Frequency (Hz)');
title({'TFR amplitude', '(black: extracted peaks; gray: their supports)'});
if fres==2, set(gca,'YScale','log'); end
surf(tt,freq,abs(TFR),'LineStyle','none'); view(2);

cc={'k',[0.5,0.5,0.5],[0.5,0.5,0.5]};
for kn=1:3
    line(tt,tfsupp(kn,:),vmax*ones(1,L),'Color',cc{kn},'LineWidth',1);
    line([tt(1),tt(1)],[tfsupp(kn,1),tfsupp(kn,1)],[abs(TFR(eind(kn,1),1)),vmax],'Color',cc{kn},'LineWidth',1);
    line([tt(end),tt(end)],[tfsupp(kn,end),tfsupp(kn,end)],[abs(TFR(eind(kn,end),end)),vmax],'Color',cc{kn},'LineWidth',1);
    line(tt,tfsupp(kn,:),abs(TFR(sub2ind(size(TFR),eind(kn,:),1:L))),'Color',cc{kn},'LineWidth',1);
end

end