function [ulf, vlf, lf, hf, lfhf, ttlpwr] = CalcLfHfParams(PSD, F, limits,plot_on)
% [ulf, vlf, lf, hf, lfhf, ttlpwr] = CalcLfHfParams(PSD, F, limits,plot_on)
%
% OVERVIEW: Compute the frequency domain features for a given PSD and
% frequency bans limits
%
% INPUT:
% PSD - power spectral density
% F - frequency vector
% limits - frequency domain analysis limits
% plot_on -
%
% OUTPUT:
% - ulf : (ms^2) Power in the ultra low frequency range (default < 0.003 Hz)
% - vlf : (ms^2) Power in very low frequency range (default 0.003 <= vlf < 0.04 Hz)
% - lf : (ms^2) Power in low frequency range (default 0.04Hz <= lf < 0.15 Hz)
% - hf : (ms^2) Power in high frequency range (default 0.15 <= hf < 0.4 Hz)
% - lfhf : Ratio LF [ms^2]/HF [ms^2]
% - ttlpwr : (ms^2) Total spectral power (approximately <0.4 Hz)
%
% REPO:
% https://github.com/cliffordlab/PhysioNet-Cardiovascular-Signal-Toolbox
% ORIGINAL SOURCE AND AUTHORS:
% Main script written by Adriana N. 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
%%
if nargin <3
ULF = [0 .003];
VLF = [0.003 .04];
LF = [.04 .15];
HF = [0.15 0.4];
limits = [ULF; VLF; LF; HF];
end
if nargin < 4
plot_on =1;
end
Indx_ULF = find( (limits(1,1) <= F) & (F <= limits(1,2)) );
Indx_VLF = find( (limits(2,1) <= F) & (F <= limits(2,2)) );
Indx_LF = find( (limits(3,1) <= F) & (F <= limits(3,2)) );
Indx_HF = find( (limits(4,1) <= F) & (F <= limits(4,2)) );
space = F(2)-F(1);
ulf = sum(PSD(Indx_ULF)*space) * 1e6; % convert to ms^2
vlf = sum(PSD(Indx_VLF)*space) * 1e6; % convert to ms^2
lf = sum(PSD(Indx_LF)*space) * 1e6; % convert to ms^2
hf = sum(PSD(Indx_HF)*space) * 1e6; % convert to ms^2
ttlpwr = sum([ulf vlf lf hf]);
lf_n = lf/ttlpwr; % normalized
hf_n = hf/ttlpwr;
lfhf = round(lf_n/hf_n*100)/100; % lf/hf ratio
if plot_on
figure
% plot PSD
plot(F,10*log10(PSD),'b','linewidth',2)
hold on
% plot limits on graph for lf and hf
plot([F(Indx_LF(1)) F(Indx_LF(1))],[-80 40],'k:')
hold on
plot([F(Indx_LF(end)) F(Indx_LF(end))],[-80 40],'k:')
hold on
plot([F(Indx_HF(end)) F(Indx_HF(end))],[-80 40],'k:')
% labelsc
text(0.07,30,'LF','Fontname','Times New Roman','Fontsize',10)
text(0.25,30,'HF','Fontname','Times New Roman','Fontsize',10)
%text(0.15, 35, 'Power Spectral Density','Fontname','Times New Roman','Fontsize',10)
text(0.3, -60, strcat('LF/HF=',num2str(lfhf)),'Fontname','Times New Roman','Fontsize',10)
ylabel('Normalized PSD (db/Hz)','Fontname','Times New Roman','fontsize',10)
xlabel('Frequency (Hz)','Fontname','Times New Roman','fontsize',10)
axis([0 .45 -80 40]);
box off
end % end plot
end % end function