关于三维图我看有人用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,:);
运用这个程序画很简单的几个频率相差很大的正余弦信号还行,但来画具体的实际信号就不好了,还有一点一直不明白,就是这个函数中的各个参数具体代表什么意思?感觉这里解释的不太清楚。