function [nhf, nha, hex, hE, hN] = nhfa(data,dt)
%
%    [nhf,nha,hex,hE,hN]=nhfa(data,dt) :
%  
%    Function to calculate the modified Hilbert frequency and amplitude
%    of data(n,k), where n specifies the length of time series, and 
%    k is the IMF component number.
%    Except frequency and amplitude the function returns other fields.
%    It is based on normalized data with modified Hilbert transformed envelope.
%    Final frequency and amplitude values are not smoothed with any filter.
%    Modified Hilbert transform function 'hilbtm()' with the end effect
%    eliminated by extension of data is used to perform a Hilbert transform.

%    Input-
%	data	- 2-D matrix data(n,k) of one IMF component with
%		  2-nd dimension being equal to 1
%	dt	- sampling period in seconds
%    Output-
%	nhf	- 2-D matrix nhf(n,k) that specifies the Hilbert frequency in Hz
%	nha	- 2-D matrix nha(n,k) that specifies the Hilbert amplitude
%		  of normalized data
%	hex	- 2-D matrix hex(n,k) that specifies the Hilbert envelope data
%	hE	- 2-D matrix he(n,k) that specifies the normalized data
%	hN	- 2-D matrix hN(n,k) that specifies the nonlinearity index
%
%    Norden Huang (NASA GSFC)	June 2, 2002 Initial
%    Jelena Marshak (NASA GSFC)	November 8, 2003 Edited
%
%    Notes-
%    Non MATLAB Library routine used in the function is:
%	'hilbtm.m'.
%
%    Temporary comments-
%    nhfa is the same as nhfam: even though in the description it first
%    specified the use of "simple Hilbert Transformed envelope",
%    instead the modified Hilbert Transform is used in both functions
%    to calculate an envelope.
%    To be consistent with the function naming it is suggested to apply 
%    the MATLAB Hilbert transform to the data. 

%----- Get dimensions
[n,m] = size(data); 
te=1:n; te=te';

%----- Apply modified Hilbert transform to get an envelope
hex=abs(hilbtm(data));

%----- Normalize the data by Hilbert envelope
hE=data./hex;

%----- Apply modified Hilbert transform
H=hilbtm(hE);

nhf=diff(unwrap(angle(H)))./(2*pi*dt);
nhf=[nhf;nhf(n-1)];
nha=abs(H);
    
%----- Compute the Nonlinear Index
hN=(nha-1).^2;

%----- Re-constitute na
%  na=na.*hex;

%----- Plot the results
plot(te,data,te,hex,te,hE,te,nhf,te,nha,te,hN, 'LineWidth', 1.5);
legend('Data','hex','hE','nhf','nha','hN');