
function [EMDdenoised,EMDstd,IMF,noise]=EMD_filtering(signal,swh,ia)
%%%%%%%% Output and input data 
%EMDdenoised  : denoised signal
%EMDstd       : uncertainty on the denoised signal 
%IMF          : IMF series for the input signal
%noise        : high-frequency noise from wavelet processing of the first IMF
%signal       : noisy signal (SLA or SWH), should contain at least 8 continuous measurements
%swh          : significant wave height, for possible tuning of wavelet processing for high sea state 
%ia           : identification of altimeter: 1=AltiKa; 2=Sentinel-3; 3=Jason-3 
%%%%%%%% Parameters to be set
%method       : EMDdenoising method 
%               'IT' for Interval Thresholding (see [2]),
%               'IIT' for Iterative Interval Thresholding (see [2]),
%               'CIIT' for Clear first Iterative Interval Thresholding [2],
%iterations   : Number of averaging iterations for IIT and CIIT methods
%altermethod  : Noise altering method
%               'circ' for random circulations
%               'perm' for random permutations
%nofsifts     : Number of sifting iterations 
%threstype    : Thresholding method
%               'hard' for hard Thresholding,
%               'soft' for soft Thresholding,
%               'softSCAD' for smoothly clipped absolute deviation
%                          (SCAD) penalty Thresholding.
%T_mult       : Multiplication factor of the universal threshold. 

% TO REFERENCE AS:
% [1] Quilfen Y., Chapron B.: On denoising satellite altimeter measurements for high-resolution geophysical signal analysis. Advances in Space Research, 68. https://doi.org/10.1016/j.asr.2020.01.005, 2021.
% OTHER REFERENCES
% [2] Y. Kopsinis, S. McLaughlin, Development of EMD-based Denoising Methods Inspired by Wavelet Thresholding.
% IEEE Trans. on Signal Processing, VOL. 57, NO. 4, APRIL 2009.
% [3] Quilfen, Y., Piolle, J.-F., and Chapron, B.: Towards improved analysis of short mesoscale sea level signals from
% satellite altimetry, Earth Syst. Sci. Data Discuss. [preprint], https://doi.org/10.5194/essd-2021-352,2022.


%%%%%%%%%%% Parameters setting
method='CIIT';
iterations=20; 
if strcmp(method,'IT')==1
    iterations=1;
end
altermethod='perm';
nofsifts=8;
threstype='hard';
if (ia==1|ia==2)
T_mult=1.925;
elseif ia==3
2.400;
end
M1=1;%Rank of the first IMF to be filtered, IMFs of lower rank are discarded in the signal reconstruction
M2=M1+4;%Rank of the last IMF to be filtered, IMFs of higher rank are fully kept in the signal reconstruction
NW=2;%Parameter to set the Maximum level to retain in the reconstruction in wavelet analysis of IMF1
%sw=prctile(swh,80);
%if sw>4.5
%  NW=3;%More filtering in IMF1 wavelet processing
%end
%%%%%%%%%%%

n=length(signal);
if n >= 8
t=1:n;


%%%%%%%%%%%%% Outliers detection
%%% Initial EMD processing for outlier detection 
[IMF,localmean] = emdfull(signal,t,nofsifts);

slan=signal;
kj=0;
if (size(IMF,1)>0 & n>16)
st=std(IMF(1,:));

ik=1;
while ik>0
clear mea dd db ki
dd=ones(length(slan),1)*0;
mea=ones(length(slan),1)*NaN;
mea(1)=slan(1);
mea(end)=slan(end);
mea(2)=slan(2);
mea(end-1)=slan(end-1);
for i=3:length(slan)-2
me=(slan(i-1)+slan(i+1))/2;
dd(i)= abs(slan(i)-me);
mea(i)=(slan(i-1)+slan(i+1)+slan(i-2)+slan(i+2))/4;
end
[db,ki]=sort(dd,'descend');
if db(1) > 4.5*st & kj~=ki(1)
ik=1;
slan(ki(1))=mea(ki(1));
kj=ki(1);
else
ik=0;
end
end
signal=slan;
end
%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%
% Second EMD processing for the EMD filtering process
[IMF,localmean] = emdfull(signal,t,nofsifts);
%%%%%%%%%%%%%%%%%%%

nofIMFs=size(IMF,1);
clear localmean
if(nofIMFs<=M2+1)
M2=nofIMFs-1;
end


if strcmp(method,'CIIT')==1
    clearfirst=1;
else
    clearfirst=0;
end

if(nofIMFs>=M1&&M2>=M1)% If not (too short segments), the processing is not performed
%%%%%%%%%%%%%%%%%%%
%%Wavelet analysis of IMF1 and generation of an iterations number of new noisy data segments
[triggered,noise,sigIMF1]=alteringnoise(IMF,iterations,nofsifts,altermethod,clearfirst,NW);
%%%%%%%%%%%%%%%%%%

if(clearfirst==1)
%%%%%%%%%%%%%%%%%%
%%Setting of noise energy variance
if(nofIMFs>2&n>iterations&NW==2)%The general case
     if(var(noise(5:end-4))/var(IMF(1,5:end-4))>1/2)%The optimal case
      estimenergy_F(1)=(median(abs(noise))/0.6745)^2; 
     else
      estimenergy_F(1)=(2*(median(abs(noise))/0.6745)^2+(median(abs(IMF(1,:)))/0.6745)^2)/3;
     end
 else
      estimenergy_F(1)=(median(abs(IMF(1,:)))/0.6745)^2; 
 end
else
  estimenergy_F(1)=(median(abs(IMF(1,:)))/0.6745)^2; 
end

%Parameters for a coloured noise
h=0.8;
beta=1.025;
h=1;
beta=1.2;
rhoh=2.01+0.2*(h-1/2)+0.12*(h-1/2)^2;

%Setting of noise variance in IMFs of rank > 1
for k=2:size(IMF,1)+3
    estimenergy_F(k)=estimenergy_F(1)/0.719*2.01^(-k);

%   estimenergy_F(k)=estimenergy_F(1)/beta*rhoh^(-2*(1-h)*k);%For a coloured noise
end
%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%
%%Calculation of denoised data and associated uncertainty
[EMDdenoised,EMDstd]=EMDdenoise_averaging(IMF,triggered.IMFpros,triggered.zcpos_pr,triggered.extrema_pr,estimenergy_F,threstype,T_mult,M1,M2);
%%%%%%%%%%%%%%%%%%%%%%%

else
EMDdenoised=signal;
EMDstd=signal*NaN;
noise=signal*NaN;
end
end
end


function [out_aver,out_std]=EMDdenoise_averaging(IMF,IMFpros,zcpos,extrema,IMFenergy,threstype,mult,M1,M2)
% IMF: Matrix containing the IMFs
% IMFpros: Cell array containing matrices of IMFs of each different random realization
% IMFenergy: Vector containing the estimated variances of each IMF. If
% isempty [], then the energies of the n_a processed IMFpros are estimated
% separately. If it is equal to 0 then the average of all the IMFpros{} is
% adopted as the final variance
%
% mult: The multiplication factor for the Threshold computation. T=std*mult
% Lcount: The number up to which the first IMFs are excluded from the sum
% Lcount_end: In the output summations, the last Lcount_end IMFs are taken
% not from the thresholded IMFs but from the original ones
%
% out_aver: Cell array {Lcount,noiseaverages} containing the denoised signal that
% corresponds to the specific i-th Lcount and averaging up to n_a, i.e.  1:n_a




noiseaverages=length(IMFpros);

for n_a=1:noiseaverages

    out_all=EMDthres_mine(IMFpros{n_a},zcpos{n_a},extrema{n_a},IMFenergy,threstype,mult);

        out(n_a,:)=IMFsubsums(out_all,M1,IMFpros{n_a},M2);
end
out_aver=mean(out,1);
out_std=std(out,1);
end

%% IMFsubsums
function out=IMFsubsums(IMFthresholded,M1,IMF,M2)
% IMFthresholded: Matrix with the thresholded IMFs.
% Lcount: The number up to which the first IMFs are excluded from the sum
% IMF: the original IMFs
% Lcount_end, the IMF number up to the end that the IMFs of the orginal
% ones are participating in the final sum.
%
% Produces a matrix, out{1,:), out(2,:), ..., out(Lcount,:), with the
% corresponding subsums of IMFthresholded. If nargin=4 then the last
% Lcount_end summed IMFs are the original and not the thresholded ones.
is=size(IMFthresholded,1);
if (is<=M2)
    out=sum(IMFthresholded(M1:end,:),1); 
else 
    out=sum(IMFthresholded(M1:M2,:),1); 
    out=out+sum(IMF(M2+1:end,:),1); 
end 

end
%% EMDthres_mine
function IMF_thresholded=EMDthres_mine(IMF,zcpos,extrema_Pr,IMFenergy,threstype,mult)
% IMF: Matrix containing the IMFs
% IMFenergy: Vector containing the estimated variances of each IMF
% mult: The multiplication factor for the Threshold computation. T=std*mult
extrema=extrema_Pr.extrema;
IMFflag=extrema_Pr.IMFflag;

nofIMFs=size(IMF,1);
for i=1:nofIMFs
    if i<=length(IMFenergy)
        T = sqrt(IMFenergy(i))*mult;
    else
        T=0;
    end

    md=min(diff(zcpos{i}));
    if isempty(md)
        md=0;
    end
    if md<=2 | IMFflag{i}==0  %i==nofIMFs
        IMF_thresholded(i,:)=thresholdIMF_no_iterp(IMF(i,:),T,threstype,zcpos{i});
    else
        IMF_thresholded(i,:)=thresholdIMF_no_iterp2(IMF(i,:),T,threstype,zcpos{i},extrema{i});

    end

end

end

%%
function [triggered,noise,sigIMF1]=alteringnoise(IMF,noiseaverages,maxiter,altermethod,clearfirst,NW)

n=length(IMF);
nofIMFs=size(IMF,1);
t=1:n;

    if clearfirst==0
        noise=IMF(1,:);
        signal=sum(IMF(2:nofIMFs,:),1);
    else
        IMF1=IMF(1,:);
        qmf=MakeONFilter('Symmlet',8);
        siz=length(IMF1);
        iexp=ceil(log(siz)/log(2));
        idif=2^iexp-siz;
        if(idif>0)
            ss=wextend(1,'symw',IMF1,ceil(idif/2),'l');
            sss=wextend(1,'symw',ss,floor(idif/2),'r');
        else
            sss=IMF1;
        end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        N=NW;%To estimate the signal, level J-1 is considered as full noise

        IMF1cleaned = recdecompsh(sss,qmf,N);
        sshff=IMF1cleaned(ceil(idif/2)+1:ceil(idif/2)+siz);
        signal=sum(IMF(2:nofIMFs,:),1)+sshff;
        sigIMF1=sshff;

        N=1;%To estimate the noise, levels up to J-1 are processed to remove signal
        IMF1cleaned = recdecompsh(sss,qmf,N);
        sshff=IMF1cleaned(ceil(idif/2)+1:ceil(idif/2)+siz);
        noise=IMF1-sshff; 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    end

for i=1:size(IMF,1)
    [maxim,minim,zc]=extremes(IMF(i,:),t);
    zcpos{i}=zc;
    extrema{i}=union(maxim,minim);
    if any(IMF(i,maxim)<0) | any(IMF(i,minim)>0)
        IMFflag{i}=0;
    else
        IMFflag{i}=1;
    end
    zerocrossings(i)=length(zc);
    energy(i)=cov(IMF(i,:));
    energyrobust(i)=(median(abs(IMF(i,:)))/0.6745)^2;
end
zcpos_pr{1}=zcpos;
IMFpros{1}=IMF;
extrema_pr{1}.extrema=extrema;
extrema_pr{1}.IMFflag=IMFflag;
zerocrossings_pr{1}=zerocrossings;

nit=noiseaverages;
for n_a=2:noiseaverages
    if strcmp(altermethod,'circ')==1
    noisetrig=circshift(noise',round(rand*length(noise)))';
    elseif strcmp(altermethod,'perm')==1
	nl=length(noise);
	il=floor(nl/nit);
	in=il*nit;
	if il==0
	   noisetrig=noise(randperm(length(noise)));
	elseif in==nl
	   for ii=1:il
	      nois=noise((ii-1)*nit+1:ii*nit);
	      noisetrig((ii-1)*nit+1:ii*nit)=nois(randperm(nit));
	   end
	else
	   for ii=1:il-1
	      nois=noise((ii-1)*nit+1:ii*nit);
	      noisetrig((ii-1)*nit+1:ii*nit)=nois(randperm(length(nois)));
	   end 
	   nois=noise(nl-nit+1:nl);
	   noisetrig(nl-nit+1:nl)=nois(randperm(length(nois)));
	end

    end
    
    signaltriggered=noisetrig+signal;
    [IMFtrig,localmean] = emdfull(signaltriggered,t,maxiter);
    nofIMFs_trig=size(IMFtrig,1);

    clear localmean
    for i=1:size(IMFtrig,1)
        [maxim,minim,zc]=extremes(IMFtrig(i,:),t);
        zcpos_tr{i}=zc;
        extrema_tr{i}=union(maxim,minim);
        if any(IMFtrig(i,maxim)<0) | any(IMFtrig(i,minim)>0)
            IMFflag_tr{i}=0;
        else
            IMFflag_tr{i}=1;
        end
        zerocrossings_tr(i)=length(zc);
    end

    IMFpros{n_a}=IMFtrig;
    clear IMFtrig
    zcpos_pr{n_a}=zcpos_tr;
    extrema_pr{n_a}.extrema=extrema_tr;
    extrema_pr{n_a}.IMFflag=IMFflag_tr;
    zerocrossings_pr{n_a}=zerocrossings_tr;

end
triggered.IMFpros=IMFpros;
triggered.zcpos_pr=zcpos_pr;
triggered.extrema_pr=extrema_pr;

end

function [IMF,localmean] = emdfull(signal,t,maxiter);

if size(signal,2)==1, signal=signal.'; end
if size(t,2)==1, t=t'; end
nofextr = 10000;
imfnum = 1;
if imfnum==1
    axx=[min(signal),max(signal)];
end

while nofextr>2
    iter = 1;
    options.iter=iter;
    [m,siftdetails]=sifting(signal,t);
    if isempty(siftdetails)
        ico=0;
        localmean(imfnum,:) = m;
        break 
    else
        ico=1;

    end

    localmean(imfnum,:) = m;

    h = signal - m;
    if maxiter>=iter
        crettotal = 1;
    else
        crettotal = 0;
    end

    iter = iter+1;
    h1=h;

    while crettotal
        [m,siftdetails]=sifting(h,t);
        if isempty(siftdetails)
            break
        end
        if maxiter>=iter
            crettotal = 1;
        else
            crettotal = 0;
        end

        if crettotal==1
            h = h - m;
            localmean(imfnum,:) = localmean(imfnum,:) + m;
        end

        iter = iter+1;
    end
    if(imfnum==1)
      ik=find(h>0&h1>0&h>h1);
      h(ik)=h1(ik);
      ik=find(h<0&h1<0&h<h1);
      h(ik)=h1(ik);
    end
    IMF(imfnum,:) = h;
    signal = signal - h;
    imfnum = imfnum + 1;
end
IMF(imfnum,:)=signal;


end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [m,siftdetails]=sifting(signal,t);
interporder=3;
[maxi, mini, zerocross] = extremes(signal,t);

nofextr = length(maxi)+length(mini);

if length(maxi) >= 2 & length(mini) >= 2
    [maxp_ex, maxv_ex, minp_ex, minv_ex] = edge_extrapolation(maxi, mini, signal, t, 3, 'mirror');
    if length(maxp_ex) <= interporder | length(minp_ex) <= interporder
        m=zeros(1,length(signal));
        siftdetails=[];
        return
    end
else
    m=zeros(1,length(signal));
    siftdetails=[];
    return
end


envmax = interp1(maxp_ex,maxv_ex,t,'spline'); % envelop of the max
envmin = interp1(minp_ex,minv_ex,t,'spline'); % envelop of the min



m = 1/2*(envmax + envmin);


siftdetails.envmax=envmax;
siftdetails.envmin=envmin;
siftdetails.maxi=maxi;
siftdetails.mini=mini;
siftdetails.zerocross=zerocross;
siftdetails.nofextr=nofextr;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%





function [maxima , minima, zerocross] = extremes(signal,t)

dif = diff(signal); %approximate derivative
dif = sign(dif);

% It considers the appearence of single zeros only.
if any(dif==0)
    zdif=find(dif==0);      %in case of zero the sign of the previous 
    if zdif(end)>1
    if zdif(1)==1, zdif(1)=[]; end 
    if zdif(end)==length(dif), zdif(end)=[]; end 
    else
    if zdif(1)==1, zdif(1)=[]; end
    end
    dif(zdif)=dif(zdif-1);  %point is adopted. 
end

N = length(dif);
sdif = dif(1:N-1).*dif(2:N);
extremespos = find(sdif==-1);
extremestype = dif(extremespos);

maxima = extremespos(find(extremestype>0))+1;
minima = extremespos(find(extremestype<0))+1;

zerocross = find(signal(1:N-1).*signal(2:N)<0);

end

function [maxp_ex, maxv_ex, minp_ex, minv_ex] = edge_extrapolation(maxi, mini, signal, t, N_extr, method);

% This function contains some parts (possibly modified) from the EMD.m function made by
%
% G. Rilling, July 2002
%
% that computes EMD (Empirical Mode Decomposition) according to:
%
% N. E. Huang et al., "The empirical mode decomposition and the
% Hilbert spectrum for non-linear and non stationary time series analysis,"
% Proc. Royal Soc. London A, Vol. 454, pp. 903-995, 1998
%
% with variations reported in:
%
% G. Rilling, P. Flandrin and P. Gon?alv?s
% "On Empirical Mode Decomposition and its algorithms"
% IEEE-EURASIP Workshop on Nonlinear Signal and Image Processing
% NSIP-03, Grado (I), June 2003




% Extrapolates a number of N_extr extemes around the edges of the signal
%
% method: 'mirror' referes to the method of Rilling et al., "on empirical
%          mode decomposition and its algorithms",
%          sortly, new extrapolation methods will be added.
%
% maxi (mini): the index positions (1,2,...,N) of the maxima (minima)
% N_extr: The number of extrapolating points

% maxp_ex (minp_ex): the time positions of the maxima (minima) including
%                    the extrapolated points at the edges
% maxv_ex (minv_ex): the values of the maxima (minima) including
%                    the values of the extrapolated points at the edges

maxp = t(maxi); % maxp (minp): the time positions of the maxima (minima)
minp = t(mini);
maxv = signal(maxi); % maxv (minv): the values of the maximuma (minimuma)
minv = signal(mini);


if strcmp(method,'mirror')==1
    indmax = maxi;
    indmin = mini;
    m = signal;
    NBSYM = N_extr;
    lx = length(signal);
    %% copy from the Rilling - Flandrin m-file.
    if indmax(1) < indmin(1)
        if m(1) > m(indmin(1))
            lmax = fliplr(indmax(2:min(end,NBSYM+1)));
            lmin = fliplr(indmin(1:min(end,NBSYM)));
            lsym = indmax(1);
        else
            lmax = fliplr(indmax(1:min(end,NBSYM)));
            lmin = [fliplr(indmin(1:min(end,NBSYM-1))),1];
            lsym = 1;
        end
    else
        if m(1) < m(indmax(1))
            lmax = fliplr(indmax(1:min(end,NBSYM)));
            lmin = fliplr(indmin(2:min(end,NBSYM+1)));
            lsym = indmin(1);
        else
            lmax = [fliplr(indmax(1:min(end,NBSYM-1))),1];
            lmin = fliplr(indmin(1:min(end,NBSYM)));
            lsym = 1;
        end
    end

    if indmax(end) < indmin(end)
        if m(end) < m(indmax(end))
            rmax = fliplr(indmax(max(end-NBSYM+1,1):end));
            rmin = fliplr(indmin(max(end-NBSYM,1):end-1));
            rsym = indmin(end);
        else
            rmax = [lx,fliplr(indmax(max(end-NBSYM+2,1):end))];
            rmin = fliplr(indmin(max(end-NBSYM+1,1):end));
            rsym = lx;
        end
    else
        if m(end) > m(indmin(end))
            rmax = fliplr(indmax(max(end-NBSYM,1):end-1));
            rmin = fliplr(indmin(max(end-NBSYM+1,1):end));
            rsym = indmax(end);
        else
            rmax = fliplr(indmax(max(end-NBSYM+1,1):end));
            rmin = [lx,fliplr(indmin(max(end-NBSYM+2,1):end))];
            rsym = lx;
        end
    end

    tlmin = 2*t(lsym)-t(lmin);
    tlmax = 2*t(lsym)-t(lmax);
    trmin = 2*t(rsym)-t(rmin);
    trmax = 2*t(rsym)-t(rmax);

    % in case symmetrized parts do not extend enough
    if tlmin(1) > t(1) | tlmax(1) > t(1)
        if lsym == indmax(1)
            lmax = fliplr(indmax(1:min(end,NBSYM)));
        else
            lmin = fliplr(indmin(1:min(end,NBSYM)));
        end
        if lsym == 1
            error('bug')
        end
        lsym = 1;
        tlmin = 2*t(lsym)-t(lmin);
        tlmax = 2*t(lsym)-t(lmax);
    end

    if trmin(end) < t(lx) | trmax(end) < t(lx)
        if rsym == indmax(end)
            rmax = fliplr(indmax(max(end-NBSYM+1,1):end));
        else
            rmin = fliplr(indmin(max(end-NBSYM+1,1):end));
        end
        if rsym == lx
            error('bug')
        end
        rsym = lx;
        trmin = 2*t(rsym)-t(rmin);
        trmax = 2*t(rsym)-t(rmax);
    end

    mlmax =m(lmax);
    mlmin =m(lmin);
    mrmax =m(rmax);
    mrmin =m(rmin);


    if length(mlmax)<NBSYM
        KK=NBSYM-length(mlmax);
        if length(mlmax)>1
            dml=mlmax(1)-mlmax(2);
            dtl=abs(tlmax(1)-tlmax(2));
        else
            dml=mlmax-maxv(1);
            dtl=abs(tlmax-maxp(1));
        end
        mlmax=[fliplr(mlmax(1)+dml*(1:KK)) mlmax];
        tlmax=[fliplr(tlmax(1)-dtl*(1:KK)) tlmax];
    end


    if length(mrmax)<NBSYM
        KK=NBSYM-length(mrmax);
        if length(mrmax)>1
            dmr=mrmax(end)-mrmax(end-1);
            dtr=abs(trmax(end)-trmax(end-1));
        else
            dmr=mrmax-maxv(end);
            dtr=abs(trmax-maxp(end));
        end

        mrmax=[mrmax mrmax(end)+dmr*(1:KK)];
        trmax=[trmax trmax(end)+dtr*(1:KK)];
    end


    if length(mlmin)<NBSYM
        KK=NBSYM-length(mlmin);
        if length(mlmin)>1
            dml=mlmin(1)-mlmin(2);
            dtl=abs(tlmin(1)-tlmin(2));
        else
            dml=mlmin-minv(1);
            dtl=abs(tlmin-minp(1));
        end
        mlmin=[fliplr(mlmin(1)+dml*(1:KK)) mlmin];
        tlmin=[fliplr(tlmin(1)-dtl*(1:KK)) tlmin];
    end


    if length(mrmin)<NBSYM
        KK=NBSYM-length(mrmin);
        if length(mrmin)>1
            dmr=mrmin(end)-mrmin(end-1);
            dtr=abs(trmin(end)-trmin(end-1));
        else
            dmr=mrmin-minv(end);
            dtr=abs(trmin-minp(end));
        end
        mrmin=[mrmin mrmin(end)+dmr*(1:KK)];
        trmin=[trmin trmin(end)+dtr*(1:KK)];
    end
    maxv_ex = [mlmax maxv mrmax];
    maxp_ex = [tlmax maxp trmax];
    minv_ex = [mlmin minv mrmin];
    minp_ex = [tlmin minp trmin];

else
    fprintf('Unknown extrapolation method')
end

end


%% thresholdIMF_no_iterp
function thresholded=thresholdIMF_no_iterp(IMF,T,threstype,zcpos)
% It does not uses the final interpolation.

n=length(IMF);
temp_out=zeros(1,n);

temp_out(find(IMF > T))=1;
DD=diff([0 temp_out]);
signalpos_plus=find(DD==1);
temp_out=zeros(1,n);
temp_out(find(IMF < -T))=1;
DD=diff([0 temp_out]);
signalpos_minus=find(DD==1);
signalpos=sort([signalpos_plus signalpos_minus]);

temp_out=zeros(1,n);
zerocrossings=[0 zcpos n]+0.5;
nodesexclude=[];
for fsp=1:length(signalpos)
    ff=find(zerocrossings<signalpos(fsp));
    ff=length(ff);
    if zerocrossings(ff)~=n
        if strcmp(threstype,'hard')==1
            temp_out(ceil(zerocrossings(ff)):floor(zerocrossings(ff+1)))=1;
        elseif strcmp(threstype,'soft')==1
            pick=max(abs(IMF((ceil(zerocrossings(ff)):floor(zerocrossings(ff+1))))));
            temp_out(ceil(zerocrossings(ff)):floor(zerocrossings(ff+1)))=(pick-T)/pick;
        elseif strcmp(threstype,'softSCAD')==1
            pick=max(abs(IMF((ceil(zerocrossings(ff)):floor(zerocrossings(ff+1))))));
            alpha=3.7;
            if pick<=2*T
                temp_out(ceil(zerocrossings(ff)):floor(zerocrossings(ff+1)))=max(0,(pick-T))/pick;
            elseif pick<=alpha*T
                temp_out(ceil(zerocrossings(ff)):floor(zerocrossings(ff+1)))=(((alpha-1)*pick-alpha*T)/(alpha-2))/pick;
            else
                temp_out(ceil(zerocrossings(ff)):floor(zerocrossings(ff+1)))=1;
            end
        end
    end
end
temp_out1=IMF.*temp_out;
thresholded=temp_out1;
end

%% thresholdIMF_no_iterp2
function thresholded=thresholdIMF_no_iterp2(IMF,T,threstype,zcpos,extremapos)
% It does not uses the final interpolation.

if length(extremapos)>length(zcpos)+1 | length(zcpos)<=1 | length(extremapos)<2 %The trend IMF
    thresholded=IMF;
else

    n=length(IMF);

    not_thresholded=find(abs(IMF(extremapos)) >= T);

    if extremapos(1)<=zcpos(1)
        if extremapos(2)<=zcpos(1)
            zcpos=[1 round((extremapos(2)+extremapos(1))/2) zcpos];
        else
            zcpos=[1 zcpos];
        end
    end
    if extremapos(end)>=zcpos(end)
        if extremapos(end-1)>=zcpos(end)
            zcpos=[zcpos round((extremapos(end)+extremapos(end-1))/2) n];
        else
            zcpos=[zcpos n];
        end
    end

    temp_out=zeros(1,n);
    for fsp=1:length(not_thresholded)
        if strcmp(threstype,'hard')==1
            if(length(zcpos)>=not_thresholded(fsp)+1)
            temp_out(zcpos(not_thresholded(fsp))+1:zcpos(not_thresholded(fsp)+1))=1;
            end
        elseif strcmp(threstype,'soft')==1
            pick=max(abs(IMF(zcpos(not_thresholded(fsp))+1:zcpos(not_thresholded(fsp)+1))));
            temp_out(zcpos(not_thresholded(fsp))+1:zcpos(not_thresholded(fsp)+1))=(pick-T)/pick;

        elseif strcmp(threstype,'softSCAD')==1
            pick=max(abs(IMF(zcpos(not_thresholded(fsp))+1:zcpos(not_thresholded(fsp)+1))));
            alpha=3.7;
            if pick<=2*T
                temp_out(zcpos(not_thresholded(fsp))+1:zcpos(not_thresholded(fsp)+1))=max(0,(pick-T))/pick;
            elseif pick<=alpha*T
                temp_out(zcpos(not_thresholded(fsp))+1:zcpos(not_thresholded(fsp)+1))=(((alpha-1)*pick-alpha*T)/(alpha-2))/pick;
            else
                temp_out(zcpos(not_thresholded(fsp))+1:zcpos(not_thresholded(fsp)+1))=1;
            end
        end
    end

    % In order to let the fist and last parts of IMF if they are over the
    % threshold
    if abs(IMF(1)) >= T;
        temp_out(1:zcpos(1))=1;
    end
    if abs(IMF(end)) >= T;
        temp_out(zcpos(end)+1:n)=1;
    end

    thresholded=IMF.*temp_out;
end
end

function [f,ns] = recdecompsh(signal,h,N)

% recdecompsh:  Smoothing method using a deterministic/stochastic 
%               approach. It calls the DECOMPSH procedure.
% Usage
%               f = recdecompsh(signal,h)
% Inputs
%   signal	    1-d Noisy signal, length(signal)= 2^J
%   h 		    Quadrature Mirror Filter for Wavelet Transform
%		          Optional, Default = Symmlet 8
%   N           Parameter to set the Maximum level to retain in the reconstruction: Lmax = J - N 
% Outputs
%   f		       Estimate, obtained by applying thresholding on the 
%               wavelet coefficients.
% Reference
%               Huang, H.-C. & Cressie, N. (2000). Deterministic/Stochastic 
%               Wavelet Decomposition for Recovery of Signal from Noisy Data. 
%               Technometrics, 42, 262-276.

if nargin < 2,
	h = MakeONFilter('Symmlet',8);
end

%initialisations
n=length(signal);
lev=floor(log2(log(n)))+1;

%Extraction of the wavelet coefficients
J = log2(n);
[reconstruct,yd,ys,ns, beta] = decompsh(signal,lev,h,J-N,0);

f = reconstruct;
end
% Copyright (c) 2001
%
% Anestis Antoniadis, Jeremie Bigot
% Laboratoire IMAG-LMC
% University Joseph Fourier
% BP 53, 38041 Grenoble Cedex 9
% France.
%
% mailto: Anestis.Antoniadis@imag.fr
% mailto: Jeremie.Bigot@imag.fr
%
% and
%
% Theofanis Sapatinas
% Department of Mathematics and Statistics
% University of Cyprus
% P.O. Box 20537  
% CY 1678 Nicosia
% Cyprus.
%
% mailto: T.Sapatinas@ucy.ac.cy 
  
function [y,yd,ys,ns, beta]  = decompsh(x,L,qmf,Lmax,show)

% decompsh:  Signal is decomposed by an orthogonal wavelet tranform:  w = Tx
%            The tranformed signal is modeled as:  w = beta + epsilon, where 
%            epsilon is a Gaussian noise. The expected denoised component b is 
%            decomposed as: beta = mu + eta, where mu is the deterministic 
%            component and eta the stochastic component. This procedure is called 
%            by the RECDECOMPSH procedure.
% Usage 
%            [y,yd,ys,ns, beta] = decompsh(x,L,qmf,Lmax,show)
% Inputs
%    x       1-d signal. length(x)= 2^J
%    L       Level of coarsest scale L
%            1<L<<J
%    qmf     Quadrature Mirror Filter for Wavelet Transform
%            Optional, Default = Symmlet 8
%    Lmax    Maximum level to retain in the reconstruction
%            (wcoef(L>Lmax) = 0),  L  <Lmax <J
%            Lmax(Default)=J-2:    finest details are supposed to be only noise
%    show    if 1 then show figures showing the process, if 0 don't
% Outputs 
%    y       Final (stochastic + deterministic) denoised estimate
%    yd      Deterministic structure estimate
%    ys      Stochastic component estimate
%    ns      Noise component
%    beta    Wavelet transform of estimate
% Notes
%            Variance estimation in finest scale is not done with the 
%            semivariogram method, but with traditional median absolute
%            deviation.
% REFERENCES
% Huang, H.-C. and Cressie, N. (2000). 
% Deterministic/stochastic wavelet decomposition for recovery of signal from noisy data. 
% Technometrics, 42, 262-276.
%
% ACKNOWLEDGEMENT
% This function was contributed by Camilo La Rota (Lab TIMC, Grenoble). 

x  = ShapeAsRow(x);
[n,J] = dyadlength(x) ;

  if (nargin < 2)|(isempty(L)),
      L=3;
  end
  if ((nargin < 3)|(isempty(qmf))),
      qmf = makeonfilter('Symmlet',8);
  end
  if ((nargin < 4) | (isempty(Lmax)))
      Lmax =J-2;  %J-1 level is supposed to be only noise for biological signals
  end
  if (nargin < 5)
      show=1;  %show decomposition process 
  end
%

%
%-------------------Wavelet Tranform ------------------------
  wcoef         = FWT_PO(x,L,qmf) ;
%
% noise variance estimation with high resolution wavelet coefficients
  wfine    = wcoef(dyad(J-1));
  %noisevar = varvar(wfine);  % variance estimation with Huang method is not as robust as they say ...
  noisevar=(mad2(wfine))^2;  % estimate with donoho method, 

  if show 
       fig1=figure; 
       set(fig1,'Name','Deterministic Component Estimation','Units','normalized','Position',[0.05,0.55,0.4,0.4]); 
  end          

% ------------------------mean structure computation (mu)------------------------------------------------


  mu            = wcoef;
  
  %cancel high resolution coefficients dues to noise
  for j=J-1:-1:Lmax+1,
        mu(dyad(j))       = zeros(size(dyad(j)));
       if show subplot(221); plot(1:n,x,'b',1:n,IWT_PO(mu,L,qmf),'r');title(['Aproximation up to level ' num2str(j)]);axis tight; pause;end         
  end

  for j=Lmax:-1:L,             %loop through scales
   wj= wcoef(dyad(j));

   if show   %show quantile plot
     subplot(2,2,2);qqplot(wj);title('qnqplot and estimated threshold'); axis tight;
   end          

   nj= length(wj);
   [wjs,l]=sort(wj);
   qj= norminv((1:nj)/(nj+1),0,1);
   Tj                = median(abs(wj))/0.6745;   %slope of quantiles-normalquantiles plot
   qwj(l)            = qj;                 
   resj              = wj-(Tj*qwj);              %residual of wavelet coefficient to corresponding quantile
   for i=1:nj
      if (abs(wj(i))>=abs(Tj*qwj(i)))
           qwj(i)=0;
      end
   end
   lambdaj         = Tj*max(abs(qwj));          %estimation of threshold for scale j

   if show   %show quantile plot
     t=axis;
     t=t(1):(t(2)-t(1))/100:t(2);
     thres=lambdaj*ones(size(t));
     hold;plot(t,thres,'r:',t,thres*(-1),'r:'); hold;
   end          
   
   muj       = ((resj.^2)./(Tj^2 + resj.^2)).*wj;      %deterministic shrinkage, determinstic component estimation


    for k = 1:2^j,
     if abs(wj(k)) <= lambdaj
        muj(k)       = 0;         %adaptive hardthreshold of deterministic component
     end
    end

    mu(dyad(j))       = muj;

    if show 
       subplot(2,1,2);
       plot(1:nj,wj,'b*',1:nj,muj,'ko');title(['Wavelet coefficients shrinkage, Level:' num2str(j)]);axis tight;hold;  
       plot(1:nj,muj,'k');    %show estimated deterministic component for scale j
       t=axis;
       t=t(1):(t(2)-t(1))/100:t(2);
       thres=lambdaj*ones(size(t));
       plot(t,thres,'r:',t,thres*(-1),'r:'); hold;
    end

   clear qwj resj muj; 
   if show subplot(221); plot(1:n,x,'b',1:n,IWT_PO(mu,L,qmf),'r');title(['Aproximation up to level ' num2str(j)]);axis tight; pause; end         

 end       %end loop trough scales

  yd    = IWT_PO(mu,L,qmf);          % intermediate result, deterministic component
  yd   = ShapeLike(yd,x);
%---------------------------end deterministic estimation-----------------------------
  if show 
       fig2=figure; 
       set(fig2,'Name','Stochastic Component Estimation','Units','normalized','Position',[0.55,0.55,0.4,0.4]); 
       subplot(2,1,1); plot(1:n,x,'b',1:n,yd,'r');title('Estimated Deterministic Component'); axis tight; 
  
       fig3=figure; 
       set(fig3,'Name','Estimated components','Units','normalized','Position',[0.55,0.05,0.4,0.4]); 
       subplot(2,1,1); plot(1:n,x,'b',1:n,yd,'k');title('Estimated Deterministic Component'); axis tight; pause;
  end          

  

%  --------------Stochastic component estimation (eta),(scale independence hypothesis, decompI)----------------------
  if show 
  end          
  beta = mu;

  for j=Lmax:-1:L,        %loop trough scales
     muj = mu(dyad(j));
     nj= length(muj);
     %stochastic process covariance estimation
     Dj       = wcoef(dyad(j)) - muj;
     Dj2      = Dj*Dj';
     Dvar     =  (Dj2)/(2^j) - noisevar;
     sigma2j = max(Dvar,0);

    %  ----------Bayesian shrinkage, estimation with stochastic and mean components-----------
     dbj          = (sigma2j/(sigma2j+noisevar)) * Dj;
     beta(dyad(j)) = muj + dbj;
    
     if show 
          figure(fig2); 
          subplot(2,1,1); plot(1:n,x,'b',1:n,yd,'r',1:n,IWT_PO(beta,L,qmf),'k');title(['Deterministic+Stochastic Aproximation up to level ' num2str(j)]); axis tight; 
          

          subplot(2,1,2);plot(1:nj,wcoef(dyad(j)),'b:' ,1:nj,muj,'r:'); axis tight; hold; 
          plot(1:nj,wcoef(dyad(j)),'b*' ,1:nj,muj,'ro'); axis tight;  
          plot(1:nj,beta(dyad(j)),'k+',1:nj,beta(dyad(j)),'k' ); hold; 
          title(['Stochastic wavelet coefficient shrinkage at level:' num2str(j)]);
          pause;
     end
     
     clear muj;
     clear Dj;
  end                 %end loop trough scales


 %------------------- end stochastic component estimation ------------------------ 
   

% Inverse transformation
  y    = IWT_PO(beta,L,qmf);      
  y    = ShapeLike(y,x);
  ys   = y-yd;
  ns   = y-x;
end  
 
% Copyright (c) 2001
%
% Anestis Antoniadis, Jeremie Bigot
% Laboratoire IMAG-LMC
% University Joseph Fourier
% BP 53, 38041 Grenoble Cedex 9
% France.
%
% mailto: Anestis.Antoniadis@imag.fr
% mailto: Jeremie.Bigot@imag.fr
%
% and
%
% Theofanis Sapatinas
% Department of Mathematics and Statistics
% University of Cyprus
% P.O. Box 20537  
% CY 1678 Nicosia
% Cyprus.
%
% mailto: T.Sapatinas@ucy.ac.cy 
