%....................Perform HHT analysis on wave data.....................
%Author: James Bramante
%Date: November 28, 2018

function [hi,lo,res,HHT,f] = wave_hht(wavDat,fs,options)
% Computes high and low freq components and HHT of discrete sampled signal
%
%INPUTS
%wavDat = discretely sampled signal (e.g. water level, depth-averaged vel)
%fs = sampling frequency
%options = cell array of cells containing {name,value} pairs of parameters
%to be used as input
%
%OUTPUTS
%hi = high frequency component of wavDat
%lo = low frequency component of wavDat
%res = residual component of wavDat
%HHT = frequency, time, IMF -space matrix of Hilbert-Huang transform of IMF
%f = vector containing frequencies from the HHT
%
% High and low frequency components are constructed from intrinsic mode 
% functions (IMF) generated from ensemble empirical mode decomposition 
% (EEMD). Then the IMFs are run through a Hilbert-Huang Transformation, to
% get the time-dependent spectra for each IMF.

%% Final variables
params = struct(...
    'freqBound',0.05,...
    'Interpolation','pchip',...
    'MaxNumIMF',10,...
    'WhiteNoise',0.1,...
    'NumEnsemble',10,...
    'fs',1 ...
    );

%% User-input 
 if (nargin > 1) & ~iscell(fs)
     params.fs = fs;
 end
 if (nargin > 2) | (iscell(fs) & (nargin>1))
    if iscell(fs)
        options = fs;
    end
     if iscell(options) 
         if iscell(options{1}) 
             if (length(options) > 1)
                for x = 1 : length(options)-1
                    if ismember(options{x}{1},fieldnames(params))
                        params.(options{x}{1}) = options{x}{2};
                    else
                        error('Unknown parameter',sprintf('The input option "%s" is not a valid parameter.',options{x}{1}));
                    end
                end
             end
             options = options{end};
         end
         if ismember(options{1},fieldnames(params))
            params.(options{1}) = options{2};
         else
            error('Unknown parameter',sprintf('The input option "%s" is not a valid parameter.',options{1}));
         end
     end
 end
 
 %% Compute the IMFs, then the HHT
 [imf,res] = eemd(wavDat,params.MaxNumIMF,params.NumEnsemble,params.WhiteNoise,params.Interpolation);
 
 % Initialize the output HHT matrix
 [hs,f] = hht(imf(:,1),params.fs,'FrequencyResolution',(params.fs/2-0)/(100*params.fs));
 HHT = zeros(length(f),size(imf,1),size(imf,2));
 HHT(:,:,1) = double(hs);
 meanFreqs = zeros(size(imf,1),size(imf,2));
 meanFreqs(:,1) = nansum(double(hs)./repmat(nansum(double(hs)),size(hs,1),1) .* repmat(f,1,size(hs,2)));
 
 for ii = 2 : size(imf,2)
     [hs,f] = hht(imf(:,ii),params.fs,'FrequencyResolution',(params.fs/2-0)/(100*params.fs));
     HHT(:,:,ii) = double(hs);
     meanFreqs(:,ii) = nansum(double(hs)./repmat(nansum(double(hs)),size(hs,1),1) .* repmat(f,1,size(hs,2)));
 end
 
 hi = sum(imf(:,nanmedian(meanFreqs,1)>params.freqBound),2);
 lo = sum(imf(:,nanmedian(meanFreqs,1)<=params.freqBound),2);
 
end
 
    





