p147685/CardiacPhase
1
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity
8f113d279a
CardiacPhase/Code/rParabEmd__L.m
João Narciso 431ff5f7d4 Initial class construction
2019-05-06 16:34:28 +02:00

306 lines
13 KiB
Matlab
Raw Blame History

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
Reference in New Issue View Git Blame Copy Permalink