% 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}
%
%------------------------------Documentation-------------------------------
%
% [iamp,iphi,ifreq] = rectfr(tfsupp,TFR,freq,wopt)
% - returns reconstructed from WFT/WT amplitude [iamp], phase [iphi]
%   and frequency [ifreq], each being 2xL matrix with first row
%   reconstructed by direct method, and second one - from the amplitude
%   ridge points; determines whether TFR is WFT or WT based on the
%   spacings between specified frequencies [freq] - linear (WFT) or
%   logarithmic (WT); rectfr.m is suitable only for WFT/WT, for SWFT/SWT
%   use ssrectfr.m instead.
%
% INPUT:
% tfsupp: 3xL matrix
%        - extracted time-frequency support of the component, containing
%          frequencies of the TFR amplitude peaks (ridge points) in the
%          first row, 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.
% TFR: FNxL matrix (rows correspond to frequencies, columns - to time)
%        - time-frequency representation (WFT or WT), to which [tfsupp]
%          correspond
% freq: FNx1 vector
%        - the frequencies corresponding to the rows of [TFR]
% wopt: structure returned by function wft.m or wt.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.
%
% NOTE: if the window/wavelet does not allow direct estimation of frequency
%       from WFT/WT, i.e. [wopt.omg=Inf] for WFT or [wopt.D=Inf] for WT (as 
%       for the Morlet wavelet), corresponding to not finite \bar{\omega}_g
%       and D_\psi in [1], then uses hybrid frequency reconstruction
%       instead of the direct one.
%
%-------------------------------Examples-----------------------------------
%
% /Calculate first WFT (see wft function):/
% [WFT,freq,wopt]=wft(sig,fs,...(parameters));
% /Extract somehow the component time-frequency support from it, e.g./
% tfsupp=ecurve(WFT,freq);
% /Then component'a parameters can be reconstructed from its extracted
% WFT support [tfsupp] as/
% [iamp,iphi,ifreq]=rectfr(tfsupp,WFT,freq,wopt);
% /For WT the same procedure is used./
%
%--------------------------------------------------------------------------


function [iamp,iphi,ifreq] = rectfr(tfsupp,TFR,freq,wopt)

[FN,L]=size(TFR); freq=freq(:); fs=wopt.fs; wp=wopt.wp;

ifreq=ones(2,L)*NaN; iamp=ones(2,L)*NaN; iphi=ones(2,L)*NaN;
tn1=find(~isnan(tfsupp(1,:)),1,'first'); tn2=find(~isnan(tfsupp(1,:)),1,'last');

%Determine the frequency resolution and transform [tfsupp] from frequencies to indices
if freq(1)<=0 || std(diff(freq))<std(freq(2:end)./freq(1:end-1))
    fres='linear'; fstep=mean(diff(freq));
    tfsupp=1+floor((1/2)+(tfsupp-freq(1))/fstep);
else
    fres='log'; fstep=mean(diff(log(freq)));
    tfsupp=1+floor((1/2)+(log(tfsupp)-log(freq(1)))/fstep);
end

%Define variables to not use cells or structures
C=wp.C; ompeak=wp.ompeak; fwtmax=wp.fwtmax;
if isfield(wp,'omg'), omg=wp.omg; else D=wp.D; end
if ~iscell(wp.fwt), fwt=wp.fwt; nflag=0;
else
    nflag=1; Lf=length(wp.fwt{2});
    if strcmpi(fres,'linear'), fxi=wp.fwt{2}; fwt=wp.fwt{1};
    else
        fxi=linspace(min(wp.fwt{2}),max(wp.fwt{2}),Lf)';
        fwt=interp1(wp.fwt{2},wp.fwt{1},fxi,'spline');
    end
    wstep=mean(diff(fxi));
end


%==================================WFT=====================================
if strcmpi(fres,'linear')
    for tn=tn1:tn2
        cs=TFR(:,tn); ii=tfsupp(2,tn):tfsupp(3,tn); ipeak=tfsupp(1,tn);
        if ~isfinite(omg)
            if tn>tn1 && tn<tn2, cw=angle(TFR(ii,tn+1))-angle(TFR(ii,tn-1)); cw(cw<0)=cw(cw<0)+2*pi; cw=cw*fs/2;
            elseif tn==1, cw=angle(TFR(ii,tn+1))-angle(cs(ii)); cw(cw<0)=cw(cw<0)+2*pi; cw=cw*fs;
            else cw=angle(cs(ii))-angle(TFR(ii,tn-1)); cw(cw<0)=cw(cw<0)+2*pi; cw=cw*fs; end
            cw=cw/(2*pi);
        end
        
        %Direct reconstruction---------------------------------------------
        if isempty(ii(isnan(cs(ii)) | ~isfinite(cs(ii))))
            ansig=(1/C)*sum(cs(ii)*(2*pi*fstep));
            iamp(1,tn)=abs(ansig); iphi(1,tn)=angle(ansig);
            if isfinite(omg)
                ifreq(1,tn)=real(-omg+(1/C)*sum(freq(ii).*cs(ii)*(2*pi*fstep))/ansig);
            else
                ifreq(1,tn)=real((1/C)*sum(cw.*cs(ii)*(2*pi*fstep))/ansig);
            end
        end
        
        %Ridge reconstruction----------------------------------------------
        if isfinite(cs(ipeak))
            ifreq(2,tn)=freq(ipeak)-ompeak/2/pi;
            if ipeak>1 && ipeak<FN %quadratic interpolation
                a1=abs(cs(ipeak-1)); a2=abs(cs(ipeak)); a3=abs(cs(ipeak+1));
                p=(1/2)*(a1-a3)/(a1-2*a2+a3);
                if abs(p)<=1, ifreq(2,tn)=ifreq(2,tn)+p*fstep; end
            end
            ximax=2*pi*(freq(ipeak)-ifreq(2,tn));
            if nflag==0 %if window FT is known in analytic form
                cmax=fwt(ximax); if isnan(cmax), cmax=fwt(ximax+10^(-14)); end
                if isnan(cmax), cmax=fwtmax; end
            else %if window FT is numerically estimated
                cid1=1+floor((ximax-fxi(1))/wstep); cid2=cid1+1;
                cid1=min([max([cid1,1]),Lf]); cid2=min([max([cid2,1]),Lf]);
                if cid1==cid2, cmax=fwt(cid1);
                else
                    cmax=fwt(cid1)+(fwt(cid2)-fwt(cid1))*(ximax-fxi(cid1))/(fxi(cid2)-fxi(cid1));
                end
            end
            ansig=2*cs(ipeak)/cmax;
            iamp(2,tn)=abs(ansig); iphi(2,tn)=angle(ansig);
        end
        
    end
end

%===================================WT=====================================
if strcmpi(fres,'log')
    for tn=tn1:tn2
        cs=TFR(:,tn); ii=tfsupp(2,tn):tfsupp(3,tn); ipeak=tfsupp(1,tn);
        if ~isfinite(D)
            if tn>tn1 && tn<tn2, cw=angle(TFR(ii,tn+1))-angle(TFR(ii,tn-1)); cw(cw<0)=cw(cw<0)+2*pi; cw=cw*fs/2;
            elseif tn==1, cw=angle(TFR(ii,tn+1))-angle(cs(ii)); cw(cw<0)=cw(cw<0)+2*pi; cw=cw*fs;
            else cw=angle(cs(ii))-angle(TFR(ii,tn-1)); cw(cw<0)=cw(cw<0)+2*pi; cw=cw*fs; end
            cw=cw/(2*pi);
        end
        
        %Direct reconstruction---------------------------------------------
        if isempty(ii(isnan(cs(ii)) | ~isfinite(cs(ii))))
            ansig=(1/C)*sum(cs(ii)*fstep);
            iamp(1,tn)=abs(ansig); iphi(1,tn)=angle(ansig);
            if isfinite(D)
                ifreq(1,tn)=real((1/D)*sum(freq(ii).*cs(ii)*fstep)/ansig);
            else
                ifreq(1,tn)=real((1/C)*sum(cw.*cs(ii)*fstep)/ansig);
            end
        end
        
        %Ridge reconstruction----------------------------------------------
        if isfinite(cs(ipeak))
            ifreq(2,tn)=freq(ipeak);
            if ipeak>1 && ipeak<FN %quadratic interpolation
                a1=abs(cs(ipeak-1)); a2=abs(cs(ipeak)); a3=abs(cs(ipeak+1));
                p=(1/2)*(a1-a3)/(a1-2*a2+a3);
                if abs(p)<=1, ifreq(2,tn)=exp(log(ifreq(2,tn))+p*fstep); end
            end
            ximax=ompeak*ifreq(2,tn)/freq(ipeak);
            if nflag==0 %if wavelet FT is known in analytic form
                cmax=fwt(ximax); if isnan(cmax), cmax=fwt(ximax*(1+10^(-14))); end
                if isnan(cmax), cmax=fwtmax; end
            else %if wavelet FT is numerically estimated
                cid1=1+floor((ximax-fxi(1))/wstep); cid2=cid1+1;
                cid1=min([max([cid1,1]),Lf]); cid2=min([max([cid2,1]),Lf]);
                if cid1==cid2, cmax=fwt(cid1);
                else
                    cmax=fwt(cid1)+(fwt(cid2)-fwt(cid1))*(ximax-fxi(cid1))/(fxi(cid2)-fxi(cid1));
                end
            end
            ansig=2*cs(ipeak)/cmax;
            iamp(2,tn)=abs(ansig); iphi(2,tn)=angle(ansig);
        end
        
    end
end

%Unwrap phases at the end
iphi(1,:)=unwrap(iphi(1,:)); iphi(2,:)=unwrap(iphi(2,:));

end