function [h,xs,w] = nspaben(data,nyy,min_w,max_w,t0,t1)
%   [nt,t,f]=nspaben(data,ny,min_w,max_w,t0,t1):  Hilbert spectrum of data(n,k)
%                                                with [f,a]=DESA1m(data,dt);
%   Input:
%      ny:   the frequency resolution
%   min_w:   the minimum frequency
%   max_w:   the maximum frequency
%      t0:   the start time
%      t1:   the end time
%   Output:
%      nt:   2-D matrix of the Hilbert transform
%       t:   the time-axis
%       f:   the frequency-axis

%       Z. SHEN 07-2-1995	Initial
%       D. XIANG 03-27-2002 Modify
%       At The Johns Hopkins University.

[npt,knb] = size(data);            %read the dimensions
dt=(t1-t0)/(npt-1);
%-----Hilbert Transform --------------------!
[omg,a]=DESA1m(data,dt);
%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,:)];
   
clear data
for i=1:knb
   for i1=1:npt
      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

%----- get local frequency -----------------!
dw=max_w - min_w;
wmx=max_w;
wmn=min_w;
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

w=linspace(wmn,wmx,nyy+1)';
xs=linspace(t0,t1,npt)';
h=flipud(rot90(h1));
h=h(1:nyy,:);