function [h,xs,w] = hsp(data,nyy,t0,t1)
%
%   [nt,t,TT]=hsp(data,ny,t0,t1):   Hilbert spectrum of the SMF-components data(n,k).
%   ny:   the period 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
%    TT:   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;
data=hilbert(data);
a=abs(data);
omg=abs(diff(data)./data(1:npt-1,:)/(2*pi*dt));
%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 -----------------!
tmx=t1-t0;
tmn=1/max(max(omg));
dtt=tmx-tmn;
omg=1../omg;
clear p;
%----- Construct the ploting matrix --------!
h1=zeros(npt-1,nyy+1);
p=round(nyy*(omg-tmn)/dtt)+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(tmn,tmx,ny)';
xs=linspace(t0,t1,n1)';
h=flipud(rot90(abs(h)));