hilbtm1.asv

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

function h=hilbtm1(x)
% hilbtm1(x): Improved Hilbert transform on x(n,k)
%        k: the sifted components; n: length of the time series.

%	Z. SHEN    Jan. 1996   Initial
%  D. XIANG   04-04-2002  Modified
%  At the Johns Hopkins University

[n,k]=size(x);
for j=1:k
    j2=2; %-------------------Treat the head ----------------------------------
    indx1=0; indx2=0;
    if (diff(x(:,j))<0)
        n_mx=1;
        x_mx=x(1,j);
        n_mn=n;
        x_mn=x(n,j);
    end
    if (diff(x(:,j))>0)
        n_mx=n;
        x_mx=x(n,j);
        n_mn=1;
        x_mn=x(1,j);
    end
    while j2<=n-1 ,
        if (x(j2-1,j)<x(j2,j))&(x(j2,j)>x(j2+1,j))	% max point
            n_mx=j2;
            x_mx=x(j2,j);
            indx1=1;
        elseif (x(j2-1,j)>x(j2,j))&(x(j2,j)<x(j2+1,j))	% min point
            n_mn=j2;
            x_mn=x(j2,j);
            indx2=1;
        end
        if indx1>=.9 & indx2>=.9
            break
        end
        j2=j2+1;
     end
     
    if n_mn<n_mx,
		zz=x(n_mn:n_mx,j);
		mm1=size(zz);mm=mm1(1);
		ia=fix(n_mn/(2*(mm-1)))+1;
		iaa=max(ia,2); %DX, modify constant from 1 to 2.
		zz1=flipud(zz);
		x1=zz1;
		for jj=1:iaa,
	  	  x1=[x1;zz(2:mm);zz1(2:mm)];
        end
        x1(:,j)=x1;
      %DX, find the first zero-cross and the slop
      [r,c]=size(x1);
      for kk=1:r-1
         if((x1(kk,j)*x1(kk+1,j))<=0)
            if(abs(x1(kk,j)) > abs(x1(kk+1,j)))
               n0=kk;
            else
               n0=kk+1;
            end
            s0=x1(kk+1,j)-x1(kk,j);
            break;
         end
      end
      x1=x1(n0:r,j);
      
      x1=[x1;x(n_mn+1:n,j)];
		sz=max(size(x1));
		np1=sz-n;
    else
		zz=x(n_mx:n_mn,j);
		mm1=size(zz);mm=mm1(1);
        ia=fix(n_mx/(2*(mm-1)))+1;
		iaa=max(ia,2); %DX, modify constant from 1 to 2.
		zz1=flipud(zz);
        x1=[zz;zz1(2:mm)];
		for jj=1:iaa,
	  	  x1=[x1;zz(2:mm);zz1(2:mm)];
        end
       
      %DX, find the first zero-cross and the slope
      [r,c]=size(x1);
      for kk=1:r-1
         if((x1(kk,j)*x1(kk+1,j))<=0)
            if(abs(x1(kk,j)) > abs(x1(kk+1,j)))
               n0=kk;
            else
               n0=kk+1;
            end
            s0=x1(kk+1,j)-x1(kk,j);
            break;
         end
      end
      x1=x1(n0:r,j);

		x1=[x1;x(n_mx+1:n,j)];
		sz=max(size(x1));
		np1=sz-n;
    end

    j2=sz-1;        %---------------------Treat the tail----------------------------
    indx1=0; indx2=0;	
    while j2>=2 ,
        if (x1(j2-1)<x1(j2))&(x1(j2)>=x1(j2+1))	% max point
            n_mx=j2;
            x_mx=x1(j2);
            indx1=1;
        elseif (x1(j2-1)>=x1(j2))&(x1(j2)<x1(j2+1))	% min point
            n_mn=j2;
            x_mn=x1(j2);
            indx2=1;
        end
        if indx1>=.9 & indx2>=.9
            break
        end
        j2=j2-1;
    end
    if n_mn<n_mx,
		zz=x1(n_mn:n_mx);
		mm=max(size(zz));
		zz1=flipud(zz);	
		ia=fix((sz-n_mx)/(2*(mm-1)))+1;
		iaa=max(ia,2); %DX, modify constant from 1 to 2.
		iaa1=iaa*2;	
		x1=[x1(1:n_mx);zz1(2:mm)];
		for jj=1:iaa
	   	 x1=[x1;zz(2:mm);zz1(2:mm)];
		end
      x1=[x1;zz(2:mm-1)];
      
       %DX, find the first zero-cross and the slop
      [r,c]=size(x1);
      for k=1:r-1
         kk=r-k+1;
         if((x1(kk,j)*x1(kk-1,j))<=0)
            if(abs(x1(kk,j)) > abs(x1(kk-1,j)))
               n0=kk;
            else
               n0=kk-1;
            end
            
            if(((x1(kk,j)-x1(kk-1,j))*s0)>0)   %same sign as s0
               break;
            end
         end
      end
      x1=x1(1:n0,j);

    else
		zz=x1(n_mx:n_mn);
		mm=max(size(zz));
		zz1=flipud(zz);
		ia=fix((sz-n_mn)/(2*(mm-1)))+1;
		iaa=max(ia,2); %DX, modify constant from 1 to 2.
		iaa1=iaa*2;
		x1=x1(1:n_mn);
		for jj=1:iaa
	    	x1=[x1;zz1(2:mm);zz(2:mm)];
		end
      x1=[x1;zz1(2:mm-1)];
      x1(:,j)=x1;
      
      %DX, find the first zero-cross and the slop
      [r,c]=size(x1);
      for k=1:r-1
         kk=r-k+1;
         if((x1(kk,j)*x1(kk-1,j))<=0)
            if(abs(x1(kk,j)) > abs(x1(kk-1,j)))
               n0=kk;
            else
               n0=kk-1;
            end
            
            if(((x1(kk,j)-x1(kk-1,j))*s0)>0)   %same sign as s0
               break;
            end
         end
      end
      x1=x1(1:n0,j);

	end
	x1=hilbert(x1);
    if ((np1+n-1)<=n0)
	    x(:,j)=x1(np1:np1+n-1);
    else
        np2=n0-np1;
        x(1:np2,j)=x1(np1:n0-1);
    end
end
h=x;
clear x