envlp.m

一种新的时频分析方法的matlab源程序。

function [envx, envm]=envlp(xx)
  
% The function ENVLP generates an envelope of maximum and minimum points
% of improved data x(n,k), where n is the length of the time series
% and k is the number of IMF components.
% The procedure extends both ends of data for each component
% with 2 waves, and then uses cubic spline to find the envelope for
% maximum points.
%
% Calling sequence-
% [envx,envm]=envlp(xx) 
%
% Input-
%	xx		    - 2-D matrix xx(n,k) of IMF components
% Output:
%	envx		- 2-D matrix envx(n,k) of envelope of max points
%	envm		- 2-D matrix envm(n,k) of envelope of min points
%
% Used by-
%   NFAM5X
   
% D. Xiang (JHU)		June, 10 2002  Initial
%	(extended both ends with 2 waves,
%	then used cubic spline to find the envelope for
%	max and min points). 

%----- Get dimensions
[n,kpt]=size(xx);

%----- Initialize
yy=zeros(n,kpt);

%----- Process each IMF
for j=1:kpt

%----- Find extrema
    maxx=0;
    maxy=0;
    minx=0;
    miny=0;
    j2=2;  
    while j2<=n-1 ,    
        if (xx(j2-1,j)<xx(j2,j))&(xx(j2,j)>=xx(j2+1,j))	% max point
            maxx=[maxx;j2];
            maxy=[maxy;xx(j2,j)];
        elseif (xx(j2-1,j)>xx(j2,j))&(xx(j2,j)<=xx(j2+1,j)) % min point
            minx=[minx;j2];
            miny=[miny;xx(j2,j)];
        end
        j2=j2+1;
    end
   
    maxx=maxx(2:end);
    maxy=maxy(2:end);
    minx=minx(2:end);
    miny=miny(2:end);
    
%----- Treat the head
    n_mx=maxx(1);
    n_mn=minx(1);
    e_mx=maxx(end);
    e_mn=minx(end);
    
    if (n_mx==-1) | (n_mn==-1) % no max or min
        disp('At least one component does not have Max or Min!');
        break;
    elseif n_mn<n_mx,
        dlt=n_mx-n_mn;
    else
        dlt=n_mn-n_mx;
    end
    
    maxx=[maxx(1)-2*dlt; maxx];
    maxx=[maxx(1)-2*dlt; maxx];
    maxy=[maxy(1);maxy];
    maxy=[maxy(1);maxy];
    minx=[minx(1)-2*dlt; minx];
    minx=[minx(1)-2*dlt; minx];
    miny=[miny(1);miny];
    miny=[miny(1);miny];

    %----- Treat the tail
    if e_mn<e_mx,
        dlt=e_mx-e_mn;
    else
        dlt=e_mn-e_mx;
    end
    
    maxx=[maxx; maxx(end)+2*dlt];
    maxx=[maxx; maxx(end)+2*dlt];
    maxy=[maxy; maxy(end)];
    maxy=[maxy; maxy(end)];
    minx=[minx; minx(end)+2*dlt];
    minx=[minx; minx(end)+2*dlt];
    miny=[miny; miny(end)];
    miny=[miny; miny(end)];
       
    if maxx(1)>0|maxx(end)<n|minx(1)>0|minx(end)<n
        disp('Extending end fail!');
        break;
    end
    
    mx=(maxx(1):maxx(end))';
    my=spline(maxx, maxy, mx);
	nx=(minx(1):minx(end))';
    ny=spline(minx, miny, nx);
    
    np1=abs(maxx(1))+2;
    np2=abs(minx(1))+2;
    
    xx(:,j)=my(np1:np1+n-1);
    yy(:,j)=ny(np2:np2+n-1);
end
%h=x;
envx=xx;
envm=yy;
clear xx yy