function [h,xs,w] = hhtp_m(data,nyy,t0,t1,ni0)
%
%   [ht,t,f]=hhtp_m(data,ny,t0,t1,ni):   Hilbert spectrum -- phase-image.
%   ny:   the frequency resolution.
%   t0:	  the true strat time.
%   t1:   the true end time.
%   ni:   vertical smouth points (Ex. ni=5)
%The outputs are:
%   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-26-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=angle(exp(j*pi)*data);
omg=diff(data)./data(1:npt-1,:)/(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
%----- 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);
      if (ii1 > 1)&(ii1 <= nyy)
         h1(j1,ii1-1)=h1(j1,ii1);
	 h1(j1,ii1+1)=h1(j1,ii1);
      end
   end
end

%---- the results ------------------!
w=linspace(wmn,wmx,nny+1)';
xs=linspace(t0,t1,npt)';
h=flipud(rot90(abs(h1)));