% 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] = ssrectfr(tfsupp,TFR,freq)
% - returns reconstructed from SWFT/SWT 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 SWFT or SWT based on the
%   spacings between specified frequencies [freq] - linear (SWFT) or
%   logarithmic (SWT); ssrectfr.m is suitable only for SWFT/SWT,
%   for WFT/WT use rectfr.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 (SWFT or SWT), to which [tfsupp]
%          correspond
% freq: FNx1 vector
%        - the frequencies corresponding to the rows of [TFR]
%
% NOTE: the iamp(2,:) gives imappropriate estimates, as ridge
%       reconstruction of the amplitude is not possible in the 
%       case of syncrosqueezed TFRs.
%
%-------------------------------Examples-----------------------------------
%
% /Calculate first SWFT (see sswft function):/
% [SWFT,freq]=sswft(sig,fs,...(parameters));
% /Extract somehow the component time-frequency support from it, e.g./
% tfsupp=ssecurve(SWFT,freq);
% /Then component'a parameters can be reconstructed from its extracted
% SWFT support [tfsupp] as/
% [iamp,iphi,ifreq]=ssrectfr(tfsupp,SWFT,freq);
% /For SWT the same procedure is used./
%
%--------------------------------------------------------------------------


function [iamp,iphi,ifreq] = ssrectfr(tfsupp,TFR,freq,varargin)

[NF,L]=size(TFR); freq=freq(:);

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

%==================================SWFT====================================
if strcmpi(fres,'linear')
    for tn=tn1:tn2
        cs=TFR(:,tn); ii=tfsupp(2,tn):tfsupp(3,tn); ipeak=tfsupp(1,tn);
        
        %Direct reconstruction---------------------------------------------
        if isempty(ii(isnan(cs(ii)) | ~isfinite(cs(ii))))
            ansig=sum(cs(ii));
            iamp(1,tn)=abs(ansig); iphi(1,tn)=angle(ansig);
            cfsig=real(sum(freq(ii).*cs(ii))/ansig);
            if cfsig>=freq(ii(1)) && cfsig<=freq(ii(end)), ifreq(1,tn)=cfsig;
            else ifreq(1,tn)=sum(freq(ii).*abs(cs(ii)))/sum(abs(cs(ii))); end
        end
        
        %Ridge reconstruction----------------------------------------------
        if isfinite(cs(ipeak))
            ansig=cs(ipeak);
            iamp(2,tn)=abs(ansig); iphi(2,tn)=angle(ansig);
            ifreq(2,tn)=freq(ipeak);
        end
        
    end
end

%===================================SWT====================================
if strcmpi(fres,'log')
    for tn=tn1:tn2
        cs=TFR(:,tn); ii=tfsupp(2,tn):tfsupp(3,tn); ipeak=tfsupp(1,tn);
        
        %Direct reconstruction---------------------------------------------
        if isempty(ii(isnan(cs(ii)) | ~isfinite(cs(ii))))
            ansig=sum(cs(ii));
            iamp(1,tn)=abs(ansig); iphi(1,tn)=angle(ansig);
            cfsig=real(sum(freq(ii).*cs(ii))/ansig);
            if cfsig>=freq(ii(1)) && cfsig<=freq(ii(end)), ifreq(1,tn)=cfsig;
            else ifreq(1,tn)=sum(freq(ii).*abs(cs(ii)))/sum(abs(cs(ii))); end
        end
        
        %Ridge reconstruction----------------------------------------------
        if isfinite(cs(ipeak))
            ansig=cs(ipeak);
            iamp(2,tn)=abs(ansig); iphi(2,tn)=angle(ansig);
            ifreq(2,tn)=freq(ipeak);
        end
        
    end
end

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

end