function [h,xs,w] = nspt(data,nyy,t0,t1)
%
%   [nt,t,p]=nspt(data,ny,t0,t1):   Hilbert spectrum of the SMF-components data(n,k).
%   ny:   the frequency resolution.
%   t0:	  the true strat time.
%   t1:   the true end time.
%The outputs are:
%   nt:   a 2-D matrix of the hilbert spectrum
%    t:	  the time axis
%    p:   the period axis

%NOTE: (1) On the screen it will appear 'max-frequency' and 'min-frequency'
%      (2) contour(t,f,nt) will plot the contour of the hilbert spectrum. 

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

[npt,knb] = size(data);            %read the dimensions
%-----Hilbert Transform --------------------!
dt=(t1-t0)/(npt-1);
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
clear filtr data
%----- get local frequency -----------------!
wmx=max(max(omg));
wmn=min(min(omg));
tmx=wmx
tmn=wmn

dw=wmx-wmn;
if wmn<0.
   error('Error: negative frequency appears!');
end
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;
%---- the results ------------------!
w=linspace(wmn,wmx,ny)';
xs=linspace(t0,t1,n1)';
h=rot90(abs(h));