dst.m

Hilbert 黄变换代码,主要实现黄变换的分解和谱分析

function z=dst(nt,dt,nzp)

% The function DST analyzes the degree of stationarity of data nt(n,m) by
% calculating the variance  for each n, where n specifies the number 
% of frequency values, and the length of time series.
%
% The first parameter is required. If not passed dt=0 and nzp=0 by default.
% If mean value is zero, its components are zero.
%
% Note: an input matrix can be generated by calling a function that
% generates the spectrum, for example NSP or NSPAB.
%
% Calling sequence-
% z=dst(nt[,dt][,nzp])
%
% Input-
%   nt	- 2-D matrix nt(n,m) that specifies the HHT spectrum
%   dt	- positive integer number specifying the smoothing length 
%   nzp	- number of non zero value points that
%		  can be ignored for each frequency    

% Output-
%   z	- vector that specifies the degree of stationarity 
%		  for each frequency value
 
% Norden Huang (NASA GSFC)	Initial
% Jelena Marshak (NASA GSFC)	Nov.7, 2003 Modified 
%   (corrected the case when mean value is equal to zero
%   and allowed the variable number of parameters).


%----- Define default parameters
if nargin < 1
   error('dst: first parameter must be passed');
end
if nargin < 3
    nzp=0;
end
if nargin < 2
    dt=0;
end

%---- Get dimensions (number of frequency and time values)
ii=size(nt);
jf1=ii(1);

%---- Calculate the mean for each frequency value
ms=sum(nt,2);
ms=ms/ii(2);

%---- Filter, normalize and calculate the variance
if (dt > 0) 
    dt1=ones(dt)/dt;
end

for i=1:jf1
    if (dt > 0) 
        nt(i,:)=filtfilt(dt1,1,nt(i,:));
    end
    nnt(i)=0;
    for j=1:ii(2)
       %----- Count the number of non zero values
       if(nt(i,j)~= 0)
          nnt(i)=nnt(i)+1;
       end
    end
    if (ms(i) ~= 0)
       %----- Normalize non zero values
       if(nnt(i) > nzp )
          nt(i,:)=nt(i,:)/ms(i);
       else
          %----- Ignore requested number of points
          nt(i,:)=0;
          % nt(i,:)=1-nt(i,:);
       end
    else
        %line=[ms(i)];
        fprintf(1,'WARNING: Zero mean at i= %6d : %6.2f\n',i,ms(i))
        nt(i,:)=0.;
    end
end
nt=1-nt;
nt=nt.*nt;
z=sum(nt,2);
z=z/ii(2);
%z(1)=0.;

%----- Plot number of non zero values and variance for each frequency
x=1:ii(1);
%set(gcf,'DefaultAxesColorOrder',[0 0 1;0 1 0;1 0 0])
%set(gcf,'DefaultAxesLineStyleOrder','-|:|--|-.')
plot(x,z,x,nnt);
legend('variance','# of non zero points');
disp('Done.');