hht时频图三维_Matlab怎样画HHT三维时频图？

关于三维图我看有人用nspab来画的，具体程序：

function [ output_args ] = Untitled1( input_args )

%UNTITLED1 Summary of this function goes here

%  Detailed explanation goes here

%nspab:

function h1=nspab(data,nyy,minw,maxw,dt)%

% The function NSPAB generates a smoothed HHT spectrum of data(n,k)

% % in time-frequency space, where

% n specifies the length of time series, and

% k is the number of IMF components.

% The frequency-axis range is prefixed.

% Negative frequency sign is reversed.

%

% MATLAB Library function HILBERT is used to calculate the Hilbert transform.

%

% Example, [h,xs,w] = nspab(lod78_p',200,0,0.12,1,3224).

%

% Functions CONTOUR or IMG can be used to view the spectrum,

%    for example contour(xs,w,h) or img(xs,w,h).

%

% Calling sequence-

% [h,xs,w] = nspab(data,nyy,minw,maxw,t0,t1)

%

% Input-

%        data        - 2-D matrix data(n,k) of IMF components%        nyy        - the frequency resolution

%        minw        - the minimum frequency

%        maxw        - the maximum frequency

%        t0        - the start time

%        t1        - the end time

% Output-

%        h        - 2-D matrix of the HHT spectrum, where

%                  the 1st dimension specifies the number of frequencies,

%                  the 2nd dimension specifies the number of time values

%        xs        - vector that specifies the time-axis values

%        w        - vector that specifies the frequency-axis values

% Z. Shen (JHU)                July 2, 1995 Initial

%----- Get dimensions (number of time points and components)

[npt,knb] = size(data);%knb表示时间点

%----- Get time interval

%----- Apply Hilbert Transform

data=hilbert(data);

a=abs(data);%

omg=abs(diff(unwrap(angle(data))))/(2*pi*dt)

%----- Smooth amplitude and frequency

filtr=fir1(8,.1);

%b = fir1(n,Wn) returns row vector b containing the n+1 coefficients of an

%order n lowpass FIR filter. This is a Hamming-window based, linear-phase

%filter with normalized cutoff frequency Wn. The output filter

%coefficients, b, are ordered in descending powers of z.

for i=1:knb

a(:,i)=filtfilt(filtr,1,a(:,i));%

omg(:,i)=filtfilt(filtr,1,omg(:,i));

end

%----- Limit frequency and amplitude

for i=1:knb

for i1=1:npt-1%针对npt-1个imf而言

if omg(i1,i) >=maxw,%maxw表示最大的频率

omg(i1,i)=maxw;%在(i1,i)处瞬时频率值为最大的值

a(i1,i)=0;%在该点处瞬时振幅值为0

elseif omg(i1,i)<=minw,

omg(i1,i)=minw;

a(i1,i)=0;

else

end

end

end

clear filtr data

%va=var(omg(200:1200))

%----- Get local frequency

dw=maxw - minw;

wmx=maxw;

wmn=minw;

%----- Construct the ploting matrix

clear p;

h1=zeros(npt-1,nyy+1);

p=round(nyy*(omg-wmn)/dw)+1;

for j1=1:npt-1

for i1=1:knb

ii1=p(j1,i1);

h1(j1,ii1)=h1(j1,ii1)+a(j1,i1);

end

end

%----- Do 3-point to 1-point averaging

[nx,ny]=size(h1);

%n1=fix(nx/3);

%h=zeros(n1,ny);

%for i1=1:n1

%h(i1,:)=(h1(3*i1,:)+h1(3*i1-1,:)+h1(3*i1-2,:));

%end

%clear h1;

%----- Do 3-points smoothing in x-direction

fltr=1./3*ones(3,1);

for j1=1:ny

h1(:,j1)=filtfilt(fltr,1,h1(:,j1));

end

clear fltr;

%----- Define the results

%w=linspace(wmn,wmx,ny-1)';

%xs=linspace(t0,t1,nx)';

h1=flipud(rot90(h1));

h1=h1(1:ny-1,:);

运用这个程序画很简单的几个频率相差很大的正余弦信号还行，但来画具体的实际信号就不好了，还有一点一直不明白，就是这个函数中的各个参数具体代表什么意思？感觉这里解释的不太清楚。