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
= 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)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)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)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)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(zerocrossingslength(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
% 1Lmax) = 0), L =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