309 lines
13 KiB
Matlab
309 lines
13 KiB
Matlab
% This function was extracted from https://nl.mathworks.com/matlabcentral/fileexchange/21409-empirical-mode-decomposition
|
|
% Authors: Raul Rato (rtr@uninova.DOT.pt) and Manuel Ortigueira (mdortigueira@uninova.pt or mdo@fct.unl.pt)
|
|
|
|
function rParabEmd = rParabEmd__L (x, qResol, qResid, qAlfa)
|
|
dbstop if warning
|
|
if(nargin~=4), error('rParabEmd__L: Use with 4 inputs.'), end
|
|
if(nargout>1), error('rParabEmd__L: Use with just one output.'), end
|
|
ArgCheck_s(x, qResol, qResid, qAlfa)
|
|
% Actual computation -------------------------------------
|
|
kc = x(:); % ket copy of the input signal
|
|
Wx= kc'*kc; % Original signal energy
|
|
quntN = length(kc); % Signal length
|
|
% loop to decompose the input signal into successive IMFs
|
|
rParabEmd= []; % Matrix which will contain the successive IMFs, and the residue
|
|
rParabEmdCnt= 0;
|
|
qDbResid= 0; %Equal energies at start
|
|
quntOscCnt= quntNOsc_s(kc);
|
|
while ((qDbResid<qResid) && (quntOscCnt>2) ) % c has some energy and oscilates
|
|
kImf = kc; % at the beginning of the sifting process, kImf is the signal
|
|
rPMOri= rGetPMaxs_s(kImf); % rPM= [xM(M), yM(M)];
|
|
rPmOri= rGetPMins_s(kImf); % rPm= [xm(m), ym(m)];
|
|
rPM= rPMaxExtrapol_s(rPMOri, rPmOri, quntN);
|
|
rPm= rPMinExtrapol_s(rPMOri, rPmOri, quntN);
|
|
quntLM= length(rPM); quntLm= length(rPm);
|
|
% if (abs(quntLM-quntLm)>2), disp('Debug: Max-Min count mismatch.'),keyboard,end;
|
|
if (abs(quntLM-quntLm)>2), disp('Debug: Max-Min count mismatch.'),end;
|
|
if(sum(abs(diff(sign(rPM(1:min(quntLM,quntLm),1)- rPm(1:min(quntLM,quntLm),1)))))>0)
|
|
% disp('Debug: Max-Min sequence mismatch.'),keyboard;
|
|
disp('Debug: Max-Min sequence mismatch.');
|
|
end
|
|
if(sum(abs(diff(sign(rPm(1:min(quntLM,quntLm),1)- rPM(1:min(quntLM,quntLm),1)))))>0)
|
|
% disp('Debug: Max-Min reverse sequence mismatch.'),keyboard;
|
|
disp('Debug: Max-Min reverse sequence mismatch.');
|
|
end
|
|
bTenv= spline(rPM(:,1), rPM(:,2), 1:quntN); % Top envelop: bTenv[n];
|
|
bDenv= spline(rPm(:,1), rPm(:,2), 1:quntN); % Down envelop: bDenv[n];
|
|
bBias= (bTenv+bDenv)/2; % first bias estimate
|
|
while true(1) % inner loop to find each IMF
|
|
WImf= kImf'*kImf; %current IMF energy
|
|
WBias= bBias*bBias'; %bias energy
|
|
if WBias*WImf<0 , warning('rParabEmd__L: Ooops, negative energy detected.'), end
|
|
if WBias> 0, DbqResol= 10*log10(WImf/WBias); else DbqResol= Inf; end
|
|
if (DbqResol>qResol), break, end %Resolution reached
|
|
%Resolution not reached. More work is needed
|
|
kImf = kImf- qAlfa*bBias'; % subtract qAlfa bias from kImf
|
|
rPMOri= rGetPMaxs_s(kImf); % rPM= [xM(M), yM(M)];
|
|
rPmOri= rGetPMins_s(kImf); % rPm= [xm(m), ym(m)];
|
|
rPM= rPMaxExtrapol_s(rPMOri, rPmOri, quntN);
|
|
rPm= rPMinExtrapol_s(rPMOri, rPmOri, quntN);
|
|
bTenv= spline(rPM(:,1), rPM(:,2), 1:quntN); % Top envelop: bTenv[n];
|
|
bDenv= spline(rPm(:,1), rPm(:,2), 1:quntN); % Down envelop: bDenv[n];
|
|
bBias= (bTenv+bDenv)/2; % new bias estimate
|
|
end % Wend true
|
|
%
|
|
rParabEmd = [rParabEmd; kImf']; % store the extracted rParabEmd in the matrix rParabEmd
|
|
kc = kc - kImf; % subtract the extracted rParabEmd from the signal
|
|
quntOscCnt= quntNOsc_s(kc);
|
|
rParabEmdCnt=rParabEmdCnt+1;
|
|
if (kc'*kc)>0
|
|
qDbResid= 10*log10(Wx/(kc'*kc));
|
|
else
|
|
qDbResid = Inf
|
|
end
|
|
%
|
|
end % Wend ((DbR... ))
|
|
if ((kc'*kc)/Wx)>(10^-12)
|
|
rParabEmd=[rParabEmd; kc']; %The residual is the last IMF
|
|
rParabEmdCnt=rParabEmdCnt+1;
|
|
NumOscqResiduais= quntNOsc_s(kc);
|
|
end
|
|
rParabEmd= rParabEmd';
|
|
|
|
end %main function
|
|
%SubFunctions ------------------------------------------------------------
|
|
%-------------------------------------------------------------------------
|
|
function ArgCheck_s(x, qResol, qResid, qAlfa)
|
|
[qL, qC] = size(x);
|
|
if ((qL*qC)~= max(qL,qC)), error('rParabEmd__L: Input signal must be a one dim vector.'), end
|
|
if ((qL*qC)<= 1), error('rParabEmd__L: Input signal must be a vector.'), end
|
|
[qL,qC] = size(qResol);
|
|
if ( ~((qL==1)&(qC==1)) ), error('rParabEmd__L: Input resolution must be a scalar.'), end
|
|
if ( qResol<=0 ), error('rParabEmd__L: Input resolution must strictly positive.'), end
|
|
[qL,qC] = size(qResid);
|
|
if ( ~((qL==1)&(qC==1)) ), error('rParabEmd__L: Input residual must be a scalar.'), end
|
|
if ( qResid<=0 ), error('rParabEmd__L: Input residual must strictly positive.'), end
|
|
[qL,qC] = size(qAlfa);
|
|
if ( ~((qL==1)&(qC==1)) ), error('rParabEmd__L: qAlfa step must be a scalar.'), end
|
|
if ( qAlfa<=0 ), error('rParabEmd__L: qAlfa step must be strictly positive.'), end
|
|
end
|
|
%--------------------------------------------------------------------------
|
|
%---------- make at 17-Jul-07 10:16:59.44
|
|
% quntNOsc_s v1.01
|
|
% build 20070409001
|
|
% Returns the oscilation count, no steps
|
|
function quntNOsc = quntNOsc_s (x)
|
|
y=0; qisTop= false; qisDown= false;
|
|
for i=2:(length(x)-1)
|
|
if( ((x(i-1)) < (x(i))) && ((x(i+1))< (x(i))) ) %Max /-\
|
|
y=y+1;
|
|
end
|
|
if( ((x(i-1)) > (x(i))) && ((x(i+1))> (x(i))) ) %min \_/
|
|
y=y+1;
|
|
end
|
|
|
|
%Top
|
|
if( ((x(i-1)) < (x(i))) && ((x(i+1))== (x(i))) ) %StepL /-
|
|
qisTop= true; qisDown= false;
|
|
end
|
|
if( ((x(i-1)) == (x(i))) && ((x(i+1))< (x(i))) ) %stepR -\
|
|
if qisTop; y=y+1; end;
|
|
qisTop= false;
|
|
end
|
|
|
|
%Downs
|
|
if( ((x(i-1)) > (x(i))) && ((x(i+1))== (x(i))) ) %stepL \_
|
|
qisTop= false; qisDown= true;
|
|
end
|
|
if( ((x(i-1)) == (x(i))) && ((x(i+1))> (x(i))) ) %StepR _/
|
|
if qisDown; y=y+1; end
|
|
qisDown=false;
|
|
end
|
|
end % for i=2:(length(x)-1)
|
|
quntNOsc= y;
|
|
end % function y = quntNOsc_s (x)
|
|
%---------- make at 17-Jul-07 10:16:59.44
|
|
function rPMaxExtrapol= rPMaxExtrapol_s(rPM, rPm, quntL)
|
|
%rPMaxExtrapol_s V1.00
|
|
% build 2007407001
|
|
% Time-mirrored top extrema (Parabolic Maxs) extrapolation
|
|
%Init ------------------------------------
|
|
rPM= sortrows(rPM); %assumes nothing on rPM sort order
|
|
rPm= sortrows(rPm); %assumes nothing on rPm sort order
|
|
kTopTim1= rPM(:,1); kTopVal= rPM(:,2);
|
|
kDwnTim1= rPm(:,1); kDwnVal= rPm(:,2);
|
|
%Start extrapolation ---------------------
|
|
if ( (kTopTim1(1)== 1) && (kDwnTim1(1)== 1) )
|
|
disp ('rPMaxExtrapol_s: Poliextrema at signal''s start');
|
|
elseif ( (kTopTim1(1)<1) || (kDwnTim1(1)< 1) )
|
|
disp ('rPMaxExtrapol_s: Invalid extrema at signal''s start');
|
|
else
|
|
kTopTim1=[2-kDwnTim1(1); kTopTim1]; % New first Top at the (one based) specular Min
|
|
kTopVal=[kTopVal(1); kTopVal]; % Same Val as old first Top
|
|
end
|
|
% End extrapolation -----------------------
|
|
if ( (kTopTim1(end)== quntL) && (kDwnTim1(end)== quntL) )
|
|
disp ('rPMaxExtrapol_s: Poliextrema at signal''s end');
|
|
elseif ( (kTopTim1(end)> quntL) || (kDwnTim1(end)> quntL) )
|
|
disp ('rPMaxExtrapol_s: Invalid extrema at signal''s end');
|
|
else
|
|
kTopTim1=[kTopTim1; (2*quntL - kDwnTim1(end))]; % New last Top at the specular Min
|
|
kTopVal=[ kTopVal; kTopVal(end)]; % Same Val as old last Top
|
|
end
|
|
% return value ------------------------
|
|
rPMaxExtrapol= sortrows([kTopTim1, kTopVal]);
|
|
end
|
|
%-------------------------------------------------------------------------
|
|
%---------- make at 17-Jul-07 10:16:59.44
|
|
function rPMinExtrapol= rPMinExtrapol_s(rPM, rPm, quntL)
|
|
%rPMinExtrapol_s V1.00
|
|
% build 2007407001
|
|
% Time-mirrored down extrema (Parabolic Mins) extrapolation
|
|
%Init ------------------------------------
|
|
rPM= sortrows(rPM); %assumes nothing on rPM sort order
|
|
rPm= sortrows(rPm); %assumes nothing on rPm sort order
|
|
kTopTim1= rPM(:,1); kTopVal= rPM(:,2);
|
|
kDwnTim1= rPm(:,1); kDwnVal= rPm(:,2);
|
|
%Start extrapolation ---------------------
|
|
if ( (kTopTim1(1)== 1) && (kDwnTim1(1)== 1) )
|
|
disp ('rPMinExtrapol_s: Poliextrema at signal''s start');
|
|
elseif ( (kTopTim1(1)<1) || (kDwnTim1(1)< 1) )
|
|
disp ('rPMinExtrapol_s: Invalid extrema at signal''s start');
|
|
else
|
|
kDwnTim1=[2-kTopTim1(1); kDwnTim1]; % New first Dwn at the (one based) specular Max
|
|
kDwnVal=[kDwnVal(1); kDwnVal]; % Same Val as old first Dwn
|
|
end
|
|
% End extrapolation -----------------------
|
|
if ( (kTopTim1(end)== quntL) && (kDwnTim1(end)== quntL) )
|
|
disp ('rPMinExtrapol_s: Poliextrema at signal''s end');
|
|
elseif ( (kTopTim1(end)> quntL) || (kDwnTim1(end)> quntL) )
|
|
disp ('rPMinExtrapol_s: Invalid extrema at signal''s end');
|
|
else
|
|
kDwnTim1=[kDwnTim1; (2*quntL - kTopTim1(end))]; % New last Dwn at the specular Max
|
|
kDwnVal=[ kDwnVal; kDwnVal(end)]; % Same Val as old last Dwn
|
|
end
|
|
% return value ------------------------
|
|
rPMinExtrapol= sortrows([kDwnTim1, kDwnVal]);
|
|
end
|
|
%-------------------------------------------------------------------------
|
|
%---------- make at 17-Jul-07 10:16:59.44
|
|
function rPMax= rGetPMaxs_s(aS) %Get Parabolic Maxs, plateaus out
|
|
% build 20070612001
|
|
kS= aS(:);
|
|
quntLenS=length(kS);
|
|
quntMaxCnt=0;
|
|
kSMNdx1= []; kSMVal=[]; %signal S Maxima indices and values
|
|
kSPMTim1= []; kSPMVal=[]; %signal S Parabolic Maxima times and values
|
|
if (quntLenS>2) %if signal has enough length
|
|
for Cnt=2:(quntLenS-1) %search the Maxs
|
|
if ( ((kS(Cnt) > kS(Cnt+1))) && ((kS(Cnt) >= kS(Cnt-1))) || ((kS(Cnt) >= kS(Cnt+1))) && ((kS(Cnt) > kS(Cnt-1))) )
|
|
quntMaxCnt=quntMaxCnt+1;
|
|
kSMNdx1= [kSMNdx1; Cnt]; kSMVal=[kSMVal; kS(Cnt)];
|
|
end
|
|
end
|
|
end
|
|
% Now we have the Maxs, lets get the Parabolic Maxs
|
|
oldxv= -Inf; oldyv= -Inf;
|
|
intGapMax= max(kS)-min(kS);
|
|
for jj=1:quntMaxCnt %for all Maxs
|
|
%xa= -1; xb= 0; xc= 1;
|
|
ya= kS(kSMNdx1(jj)-1); % Sample point before
|
|
yb= kS(kSMNdx1(jj)); % Sample point, == kSMVal(jj)
|
|
yc= kS(kSMNdx1(jj)+1); % Sample point after
|
|
D= (-4*yb+2*ya+2*yc);
|
|
if (D==0), xv= kSMNdx1(jj);
|
|
else xv= kSMNdx1(jj)+(ya-yc)/D; end; % Vertix abscissa
|
|
D= (-16*yb+ 8*ya+ 8*yc);
|
|
if (D==0), yv= yb;
|
|
else yv= yb+ (2*yc*ya- ya*ya- yc*yc)/D; end;
|
|
% Lets check for double maxima
|
|
if ( (xv==oldxv)||(abs(yv-oldyv)/abs(xv-oldxv))> (2*intGapMax) )
|
|
xv= (xv+ oldxv)/2; yv= max(yv,oldyv); %Double found
|
|
kSPMTim1(length(kSPMTim1))= xv; kSPMVal(length(kSPMVal))= yv;
|
|
else
|
|
kSPMTim1= [kSPMTim1; xv]; kSPMVal=[kSPMVal; yv];
|
|
end
|
|
oldxv= xv; oldyv= yv;
|
|
end % for jj=1:quntMaxCnt
|
|
if quntMaxCnt>0
|
|
if ( kS(1) >= kSPMVal(1) )
|
|
kSPMTim1= [1; kSPMTim1]; kSPMVal=[kS(1); kSPMVal ]; %Start must be included as a Max
|
|
end
|
|
if ( kS(end) >= kSPMVal(end))
|
|
kSPMTim1= [kSPMTim1; quntLenS]; kSPMVal=[kSPMVal; kS(end)]; %End must be included as a Max
|
|
end
|
|
end
|
|
if quntMaxCnt==0
|
|
if ( kS(1) > kS(2) )
|
|
kSPMTim1= [1; kSPMTim1]; kSPMVal=[kS(1); kSPMVal ]; %Start must be included as a Max
|
|
end
|
|
if ( kS(end) > kS(end-1))
|
|
kSPMTim1= [kSPMTim1; quntLenS]; kSPMVal=[kSPMVal; kS(end)]; %End must be included as a Max
|
|
end
|
|
end
|
|
if quntMaxCnt<0
|
|
error('rGetPMaxs_s: Invalid MaxCnt value');
|
|
end
|
|
rPMax= sortrows([kSPMTim1, kSPMVal]);
|
|
end
|
|
%---------- make at 17-Jul-07 10:16:59.44
|
|
function rPMin= rGetPMins_s(aS) %Get Parabolic Mins, plateaus out
|
|
% build 20070612001
|
|
kS= aS(:);
|
|
quntLenS=length(kS);
|
|
quntMinCnt=0;
|
|
kSMNdx1= []; kSMVal=[]; %signal S Minima indices and values
|
|
kSPMTim1= []; kSPMVal=[]; %signal S Parabolic Minima times and values
|
|
if (quntLenS>2) %if signal has enough length
|
|
for Cnt=2:(quntLenS-1) %search the Mins
|
|
if ( ((kS(Cnt) < kS(Cnt+1))) && ((kS(Cnt) <= kS(Cnt-1))) || ((kS(Cnt) <= kS(Cnt+1))) && ((kS(Cnt) < kS(Cnt-1))) )
|
|
quntMinCnt=quntMinCnt+1;
|
|
kSMNdx1= [kSMNdx1; Cnt]; kSMVal=[kSMVal; kS(Cnt)];
|
|
end
|
|
end
|
|
end
|
|
% Now we have the Mins, lets get the Parabolic Mins
|
|
oldxv= -Inf; oldyv= -Inf;
|
|
intGapMax= max(kS)-min(kS);
|
|
for jj=1:quntMinCnt %for all Mins
|
|
%xa= -1; xb= 0; xc= 1;
|
|
ya= kS(kSMNdx1(jj)-1); % Sample point before
|
|
yb= kS(kSMNdx1(jj)); % Sample point, == kSMVal(jj)
|
|
yc= kS(kSMNdx1(jj)+1); % Sample point after
|
|
D= (-4*yb+2*ya+2*yc);
|
|
if (D==0), xv= kSMNdx1(jj);
|
|
else xv= kSMNdx1(jj)+(ya-yc)/D; end; % Vertix abscissa
|
|
D= (-16*yb+ 8*ya+ 8*yc);
|
|
if (D==0), yv= yb;
|
|
else yv= yb+ (2*yc*ya- ya*ya- yc*yc)/D; end;
|
|
% Lets check for double minima
|
|
if ( (xv==oldxv)||(abs(yv-oldyv)/abs(xv-oldxv))> (2*intGapMax) )
|
|
xv= (xv+ oldxv)/2; yv= min(yv,oldyv); %Double found
|
|
kSPMTim1(length(kSPMTim1))= xv; kSPMVal(length(kSPMVal))= yv;
|
|
else
|
|
kSPMTim1= [kSPMTim1; xv]; kSPMVal=[kSPMVal; yv];
|
|
end
|
|
oldxv= xv; oldyv= yv;
|
|
end % for jj=1:quntMinCnt
|
|
if quntMinCnt>0
|
|
if ( kS(1) <= kSPMVal(1) )
|
|
kSPMTim1= [1; kSPMTim1]; kSPMVal=[kS(1); kSPMVal ]; %Start must be included as a Min
|
|
end
|
|
if ( kS(end) <= kSPMVal(end))
|
|
kSPMTim1= [kSPMTim1; quntLenS]; kSPMVal=[kSPMVal; kS(end)]; %End must be included as a Min
|
|
end
|
|
end
|
|
if quntMinCnt==0
|
|
if ( kS(1) < kS(2) )
|
|
kSPMTim1= [1; kSPMTim1]; kSPMVal=[kS(1); kSPMVal]; %Start must be included as a Min
|
|
end
|
|
if ( kS(end) < kS(end-1))
|
|
kSPMTim1= [kSPMTim1; quntLenS]; kSPMVal=[kSPMVal; kS(end)]; %End must be included as a Min
|
|
end
|
|
end
|
|
if quntMinCnt<0
|
|
error('rGetPMins_s: Invalid MinCnt value');
|
|
end
|
|
rPMin= sortrows([kSPMTim1, kSPMVal]);
|
|
end
|
|
%---------- make at 17-Jul-07 10:16:59.44 |