nsptab.m

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

function [h,xs,w] = nsptab(data,nyy,min_w,max_w,t0,t1)
%   [nt,t,p]=nsptab(data,ny,min_p,max_p,t0,t1)
%   Hilbert spectrum on data(n,k), with the frequency-axis range prefixed.
%      ny:   the frequency resolution
%   min_p:   the minimum period
%   max_w:   the maximum period
%      t0:   the start time
%      t1:   the end time
%   The outputs ar2:
%      nt:   2-D matrix of the Hilbert transform
%       t:   the time-axis
%       p:   the period-axis

%       Z. SHEN 07-2-1995
%       At The Johns Hopkins University.

[npt,knb] = size(data);            %read the dimensions
dt=(t1-t0)/(npt-1);
%-----Hilbert Transform --------------------!
data=hilbert(data);
a=abs(data);
%omg=abs(diff(data)./data(1:npt-1,:)/(2*pi*dt));
%omg=1./omg;
omg=abs(diff(unwrap(angle(data))))/(2*pi*dt);
%----- 5-points smoothing of a, p ----------!
filtr=fir1(8,.1);
for i=1:knb
   a(:,i)=filtfilt(filtr,1,a(:,i));
   omg(:,i)=filtfilt(filtr,1,omg(:,i));
end
for i=1:knb
   a(:,i)=filtfilt(filtr,1,a(:,i));
   omg(:,i)=filtfilt(filtr,1,omg(:,i));
end
for i=1:knb
   for i1=1:npt-1
      if omg(i1,i) >=max_w,
         omg(i1,i)=max_w;
	 a(i1,i)=0;
      elseif omg(i1,i)<=min_w,
         omg(i1,i)=min_w;
	 a(i1,i)=0;
      else
      end
   end
end

clear filtr data
%----- get local frequency -----------------!
dw=max_w - min_w;
wmx=max_w;
wmn=min_w;
clear p;
%----- Construct the ploting matrix --------!
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
%---- Smoothing in x-direction ------------!
[nx,ny]=size(h1);               %3-points to 1-point.
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,:))/3.;
end
clear h1;
fltr=1./3*ones(3,1);  %3-points smoothing in x-direction
for j1=1:ny
   h(:,j1)=filtfilt(fltr,1,h(:,j1));
end
clear fltr;

w=linspace(wmn,wmx,ny-1)';
xs=linspace(t0,t1,n1)';
h=rot90(h);
h=h(1:ny-1,:);