function [h,xs,w] = nspen(data,nyy,t0,t1)
%
%   [nt,t,f]=nspen(data,ny,t0,t1):   Hilbert spectrum using [f,a]=DESA1m(data,dt);
%
%   Input:
%   data: input data.
%   ny:   the frequency resolution.
%   t0:	  the true start time.
%   t1:   the true end time.
%   
%   Output:
%   nt:   a 2-D matrix of the hilbert spectrum
%    t:	  the time axis
%    f:   the frequency 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   Initial
%       D. Xiang 03-25-2002  Modify
%       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=diff(unwrap(angle(data)))/(2*pi*dt);
%---- force the negative omg to be zero, DX. ---!
%omg=(omg+abs(omg))/2.;
% add the last row to omg so that it has the same dimention
%omg=[omg;omg(npt-1,:)];
[omg,a]=DESA1m(data,dt);

clear data
%----- get local frequency -----------------!
wmx=max(max(omg))
wmn=min(min(omg))
dw=wmx-wmn;
 if wmn<0.
   error('Error: negative frequency appears!');
 end
 clear p;
%----- Construct the ploting matrix --------!
h1=zeros(npt,nyy+1);
p=round(nyy*(omg-wmn)/dw)+1;
for j1=1:npt
   for i1=1:knb
      ii1=p(j1,i1);
      h1(j1,ii1)=h1(j1,ii1)+a(j1,i1);
   end
end
%---- the results ------------------!
w=linspace(wmn,wmx,nyy+1)';
xs=linspace(t0,t1,npt)';
h=flipud(rot90(abs(h1)));