function [p,f,t] = eemd_hht(s,rnoise,nensemble,fs,fbin,minf,maxf,graph_option)
% Time-frequency analysis using Hilbert-Huang Transform and Empirical Mode
% Decomposition (EEMD)
%
% Syntax: [p,f,t] = shht(s,rnoise,nensemble,fs,fbin,minf,maxf,graph_option)
%
% Input parameters:
% s: one-dimension time series
% rnoise: the level of noise used in ensemble empirical mode decomposition
%         the rnoise value represents the fraction of standard deviation of
%         time series (e.g., 0.1)
% nensemble: the number of ensemble to average the results of noise-assisted
%            empirical mode decomposition. 
%            The error of IMFs equals to rnoise/sqrt(nensemble)        
% fs: sampling frequency (Hz) 
% fbin: number of bins in dividing frequency axis
% minf: minimum of observed frequency
% maxf: maximum of observed frequency
% graph_option: option for visualizatoin of time-frequency HHT mapping
%
% Output parameters:
% p: two-dimensional power spectrum matrix
% f: frequency axis between minimum and maximum of observed frequency
% t: time axis of the given signal
%
% Example: fs=100;
%          t=1/fs:1/fs:10;
%          s=sin(2*pi*10*t);
%          [p,f,t] = shht(s,0.05,100,fs,200,0,20,1);
%
% Ver 1.0: Albert C. Yang, MD, PhD 9/9/2019
%          Revised 10/2: adjust c-axis lim in graph option
%
% Referece
% 1. Huang, M. L. Wu, S. R. Long, S. S. Shen, W. D. Qu, P. Gloersen, and 
%    K. L. Fan (1998) The empirical mode decomposition and the Hilbert 
%    spectrum for nonlinear and non-stationary time series analysis. 
%    Proc. Roy. Soc. Lond., 454A, 903-993.
% 2. Wu, Z., Huang, N. E, S. R. Long, and C.-K. Peng (2007) On the trend, 
%    detrending, and the variability of nonlinear and non-stationary time 
%    series. Proc. Natl. Acad. Sci. USA., 104, 14889-14894.

% Empirical mode decomposition (EEMD)
imf = eemd(s,rnoise,nensemble,0)';

% Determine parameters
[nt nimf] = size(imf);
dt = 1/fs;
p = zeros(nt-1,fbin);
df = (maxf-minf)/fbin;
t = (0:1/fs:nt/fs)';
f = (minf:df:maxf)';

% Hilbert Transform
data = hilbert(imf);
a = abs(data);
omg = abs(diff(unwrap(angle(data))))/(2*pi*dt);

%----- Smooth amplitude and frequency
filtr=fir1(12,0.1);
for i=1:nimf
   a(:,i)=filtfilt(filtr,1,a(:,i));
   omg(:,i)=filtfilt(filtr,1,omg(:,i));
end

% Fill time frequency spectral density matrix
findex = round((omg-minf)/df)+1;
findex(findex > fbin) = 0;
for i=1:nt-1
    for j=1:nimf
        if findex(i,j)>0
            p(i,findex(i,j)) = p(i,findex(i,j)) + a(i,j);
        end
    end
end

if graph_option == 1
    meanp=mean(reshape(p,size(p,1)*size(p,2),1));
    stdp=std(reshape(p,size(p,1)*size(p,2),1));
    img(t,f,p');
    ax=gca;
    ax.CLim = [0 meanp+stdp];
    xlabel('Time');
    ylabel('Hz')
end