function out = EvalFrequencyDomainHRVstats(NN, tNN, sqi, HRVparams, windows_all)
%
% out = EvalFrequencyDomainHRVstats (NN, tNN, , sqi, settings)
%   
%   OVERVIEW:   This function returns frequency domain HRV metrics 
%               calculated on input NN intervals.
%
%   INPUT:      MANDATORY:
%               NN          : a single row of NN (normal normal) interval
%                             data in seconds
%               tNN         : a single row of time indices of the rr interval 
%                             data (seconds)
%               sqi         : (Optional) Signal Quality Index; Requires a 
%                             matrix with at least two columns. Column 1 
%                             should be timestamps of each sqi measure, and 
%                             Column 2 should be SQI on a scale from 0 to 1.
%                             Additional columns can be included with
%                             additional sqi at the same timestamps
%               HRVparams   : struct of settings for hrv_toolbox analysis
%               windows_all : vector containing the starting time of each
%                             windows (in seconds) 
%
%   OUTPUT:     out.ulf        : (ms^2) Power in the ultra low frequency range (default < 0.003 Hz)
%               out.vlf        : (ms^2) Power in very low frequency range (default 0.003 <= vlf < 0.04 Hz)
%               out.lf         : (ms^2) Power in low frequency range (default 0.04Hz  <= lf < 0.15 Hz)
%               out.hf         : (ms^2) Power in high frequency range (default 0.15 <= hf < 0.4 Hz)
%               out.lfhf       : Ratio LF [ms^2]/HF [ms^2]
%               out.ttlpwr     : (ms^2) Total spectral power (approximately <0.4 Hz)
%               out.fdflag     : 1 - Lomb Periodogram or other method failed
%                                2 - Not enough high SQI data
%                                3 - Not enough data in the window to analyze
%                                4 - Window is missing too much data
%                                5 - Success
%
%	REPO:       
%       https://github.com/cliffordlab/PhysioNet-Cardiovascular-Signal-Toolbox
%   ORIGINAL SOURCE AND AUTHORS:     
%       Gari Clifford HRV Tools
%                   G. Clifford 2001 gari@mit.edu, calc_lfhf.m 
%                   http://www.robots.ox.ac.uk/~gari/CODE/HRV/
%       ChengYu Lui
%       Adriana Vest
%       Dependent scripts written by various authors 
%       (see functions for details)       
%	COPYRIGHT (C) 2016 
%   LICENSE:    
%       This software is offered freely and without warranty under 
%       the GNU (v3 or later) public license. See license file for
%       more information  
%
%   10-28-2017 Modified by Giulia Da Poian: 
%   Initlialiezed output as NaN
%   Removed fft, Welch and Burg method without resampling. 
%   Added function for resampling when required by the method 
%   Acceppting only one method at time to make the function output more
%   clear
%   Removed option to write results on file directly from this function
%   Converted output to ms^2



% Verify input arguments

if nargin < 5
    error('Eval_FrequencydomainHRVstats: wrong number of input arguments!')
end

if isempty(sqi) 
     sqi(:,1) = tNN;
     sqi(:,2) = ones(length(tNN),1);
end

% Set Defaults
windowlength = HRVparams.windowlength;
fd_threshold1 = HRVparams.sqi.LowQualityThreshold;
fd_threshold2 = HRVparams.RejectionThreshold;
limits = HRVparams.freq.limits;
method = HRVparams.freq.method;
plot_on = HRVparams.freq.plot_on;
debug_sine = HRVparams.freq.debug_sine;   % if use the sine wave to debug
f_sine = HRVparams.freq.debug_freq;       % the frequency of the added sine wave
weight = HRVparams.freq.debug_weight;
sf = HRVparams.freq.resampling_freq;


% Preallocate arrays before entering the loop

out.ulf = nan(1,length(windows_all));
out.vlf = nan(1,length(windows_all));
out.lf = nan(1,length(windows_all));
out.hf = nan(1,length(windows_all));
out.lfhf = nan(1,length(windows_all));
out.ttlpwr = nan(1,length(windows_all));
out.fdflag = nan(1,length(windows_all));
% Window by Window Analysis

% Loop through each window of RR data
for iWin = 1:length(windows_all)
    % Check window for sufficient data
    if ~isnan(windows_all(iWin))    
        % Isolate data in this window
        idx_NN_in_win = find(tNN >= windows_all(iWin) & tNN < windows_all(iWin) + windowlength);
        idx_sqi_win = find(sqi(:,1) >= windows_all(iWin) & sqi(:,1) < windows_all(iWin) + windowlength);

        sqi_win = sqi(idx_sqi_win,:);
        t_win = tNN(idx_NN_in_win);
        nn_win = NN(idx_NN_in_win);

        % Analysis of SQI for the window
        lowqual_idx = find(sqi_win(:,2) < fd_threshold1);

        % If enough data has an adequate SQI, perform the calculations
        if numel(lowqual_idx)/length(sqi_win(:,2)) < fd_threshold2

            % Initialize variables
            % maxF=fs/2; % This calculation works for regularly sampled data
            N     =  length(nn_win);       % RR interval series length
            % m_fs  = 1/mean(nn_win);     % mean frequency of heart rate, i.e., the mean sample rate (fs) of RR sereis  
            % max_f = .5/(min(nn_win));   % max frequency of RR interval series
            % nfft  = 2^nextpow2(N);      % Next power of 2 from N
            nfft = 1024;
            %F = [1/nfft:1/nfft:m_fs];  % setting up frequency vector
            F = [1/nfft:1/nfft:.5];  % setting up frequency vector

            % add sine wave to RR signal
            if debug_sine
                s_sin  = weight*sin(2*pi*f_sine*t);
                nn_win = nn_win + s_sin;
            end

            % subtract mean of segment
            if HRVparams.freq.zero_mean
                rr_0 = nn_win - mean(nn_win); %rudimentary detrending
            else
                rr_0 = nn_win;
            end
            
            switch method
                % 1. Lomb-Scargle Periodogram (no resampling)
                case 'lomb'
                    try
                        [PSDlomb,Flomb] = CalcLomb(t_win,rr_0,F,nfft,HRVparams.freq.normalize_lomb);
                        % plomb equivalent to CalcLomb when normalized 
                        % lomb-scargle periodogram
                        %[PSDmatlablomb,fplombout] = plomb(vv,tt,F,'normalized');
                        [out.ulf(iWin), out.vlf(iWin), out.lf(iWin), out.hf(iWin), out.lfhf(iWin),...
                            out.ttlpwr(iWin)] = CalcLfHfParams(PSDlomb, Flomb, limits, plot_on);
                        
                        out.fdflag(iWin) = 5; %'sucess';

                    catch
                        out.fdflag(iWin) = 1; %'lomb_failed'
                    end

                % 2. Pwelch (after resampling)
                case 'welch'                    
                    % Resampling
                    rr_int = ApplyResamplingToRR(t_win,rr_0,HRVparams);
                    try
                        [PSDwelch,Fwelch] = pwelch(rr_int,[],[],2^nextpow2(length(rr_int)),sf);
                        [out.ulf(iWin), out.vlf(iWin), out.lf(iWin), out.hf(iWin), out.lfhf(iWin),...
                                out.ttlpwr(iWin)] = CalcLfHfParams(PSDwelch, Fwelch, limits, plot_on) ;  
                        out.fdflag(iWin) = 5; %'sucess';
                    catch
                        out.fdflag(iWin) = 1; %'Pwelch faild';
                    end
                    
                % 3. FFT (after resampling)
                case 'fft'
                    % Resampling
                    rr_int = ApplyResamplingToRR(t_win,rr_0,HRVparams);
                    try
                        Y_FFT2 = fft(rr_int);  % fft using the original sample length
                        PSDfft = Y_FFT2.*conj(Y_FFT2)/length(rr_int);
                        Ffft = sf*(0:floor(length(rr_int)/2)+1)/length(rr_int);
                        [out.ulf(iWin), out.vlf(iWin), out.lf(iWin), out.hf(iWin), out.lfhf(iWin),...
                                out.ttlpwr(iWin)] = CalcLfHfParams(PSDfft(1:length(Ffft)), Ffft, limits, plot_on);
                        out.fdflag(iWin) = 5; %'sucess';
                    catch
                        out.fdflag(iWin) = 1; %'fft faild';
                    end

                % 4. Burg (after resampling)
                case 'burg'
                    % Resampling
                    rr_int = ApplyResamplingToRR(t_win,rr_0,HRVparams);
                    try
                    % ChengYu's Burg Method with resampled data
                    pb = HRVparams.freq.resampled_burg_poles; % pole setting
                    [PSD_Burg1,f_Burg1] = pburg(rr_int,pb,2^nextpow2(length(rr_int)),sf);
                    [out.ulf(iWin), out.vlf(iWin), out.lf(iWin), out.hf(iWin), out.lfhf(iWin),...
                            out.ttlpwr(iWin)] = CalcLfHfParams(PSD_Burg1,f_Burg1, limits, plot_on);
                        out.fdflag(iWin) = 5; %'sucess';
                    catch
                        out.fdflag(iWin) = 1; %'Burg faild';
                    end
            end
        else
            out.fdflag(iWin) = 2; %'nt_enuf_hi_sqi_data'
        end % end conditional statements that run only if SQI is adequate

    else
        out.fdflag(iWin) = 3; %'nt_enuf_data_in_win';   
    end % end check for sufficient data
end % end of loop through windows


%   References:
%               5 minutes is the shortest period that HRV spectral metrics
%               should be calculated according to 184(Clifford Thesis). With a 5 min
%               window, the lowest frequency that can be theoretically resolved is
%               1/300 (.003 Hz).