• wendy222
  • matlab
  • 3KB
  • rar
  • 0
  • 10 积分
  • 26
  • 2017-06-13 11:13
  • freqDomainHRV.m
function output = freqDomainHRV(ibi,VLF,LF,HF,AR_order,window,noverlap,nfft,fs,methods,flagPlot) %freqDomainHRV - calculates freq domain HRV using FFT, AR, and Lomb-Scargle %methods % %Inputs: ibi = 2Dim array of time (s) and inter-beat interval (s) % AR_order = order of AR model % window = # of samples in window % noverlap = # of samples to overlap % fs = cubic spline interpolation rate / resample rate (Hz) % nfft = # of points in the frequency axis % methods = methods of calculating freqDomain. Default is all 3 % methhods %Outputs: output is a structure containg all HRV. One field for each PSD method % Output units include: % peakHF,LF,VLF (Hz) % aHF,aLF,aVLF (ms^2) % pHF,pLF,pVLF (%) % nHF,nLF,nVLF (%) % PSD (ms^2/Hz) % F (Hz) %Usage: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Copyright (C) 2010, John T. Ramshur, jramshur@gmail.com % % This file is part of HRVAS % % HRVAS is free software: you can redistribute it and/or modify % it under the terms of the GNU General Public License as published by % the Free Software Foundation, either version 3 of the License, or % (at your option) any later version. % % HRVAS is distributed in the hope that it will be useful, % but WITHOUT ANY WARRANTY; without even the implied warranty of % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the % GNU General Public License for more details. % % You should have received a copy of the GNU General Public License % along with Foobar. If not, see <http://www.gnu.org/licenses/>. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %check input if nargin<9 error('Not enough input arguments!') elseif nargin<10 methods={'welch','ar','lomb'}; flagPlot=false; elseif nargin<11 flagPlot=false; end flagWelch=false; flagAR=false; flagLomb=false; for m=1:length(methods) if strcmpi(methods{m},'welch') flagWelch=true; elseif strcmpi(methods{m},'ar') flagAR=true; elseif strcmpi(methods{m},'lomb') flagLomb=true; end end t=ibi(:,1); %time (s) y=ibi(:,2); %ibi (s) y=y.*1000; %convert ibi to ms %assumes ibi units are seconds maxF=fs/2; %prepare y y=detrend(y,'linear'); y=y-mean(y); %Welch FFT if flagWelch [output.welch.psd,output.welch.f]=calcWelch(t,y,window,noverlap,nfft,fs); output.welch.hrv=calcAreas(output.welch.f,output.welch.psd,VLF,LF,HF); else output.welch=emptyData(nfft,maxF); end %AR if flagAR [output.ar.psd,output.ar.f]=calcAR(t,y,fs,nfft,AR_order); output.ar.hrv=calcAreas(output.ar.f,output.ar.psd,VLF,LF,HF); else output.ar=emptyData(nfft,maxF); end %Lomb if flagLomb [output.lomb.psd,output.lomb.f]=calcLomb(t,y,nfft,maxF); output.lomb.hrv=calcAreas(output.lomb.f,output.lomb.psd,VLF,LF,HF,true); else output.lomb=emptyData(nfft,maxF); end %plot all three psd if flagPlot figure; h1=subplot(3,1,1); plotPSD(h1,output.welch.f,output.welch.psd,VLF,LF,HF,[0 0.6],[]); legend('welch') h2=subplot(3,1,2); plotPSD(h2,output.ar.f,output.ar.psd,VLF,LF,HF,[0 0.6],[]); legend('AR') h3=subplot(3,1,3); plotPSD(h3,output.lomb.f,output.lomb.psd,VLF,LF,HF,[0 0.6],[]); legend('Lomb-Scargle') end end function [PSD,F]=calcWelch(t,y,window,noverlap,nfft,fs) %calFFT - Calculates the PSD using Welch method. % %Inputs: %Outputs: %Prepare y t2 = t(1):1/fs:t(length(t));%time values for interp. y=interp1(t,y,t2','spline')'; %cubic spline interpolation y=y-mean(y); %remove mean %Calculate PSD [PSD,F] = pwelch(y,window,noverlap,(nfft*2)-1,fs,'onesided'); %uses a hamming window end function [PSD,F]=calcAR(t,y,fs,nfft,AR_order) %calAR - Calculates the PSD using Auto Regression model. % %Inputs: %Outputs: %Prepare y t2 = t(1):1/fs:t(length(t)); %time values for interp. y=interp1(t,y,t2,'spline')'; %cubic spline interpolation y=y-mean(y); %remove mean y = y.*hamming(length(y)); %hamming window %Calculate PSD %Method 1 % [A, variance] = arburg(y,AR_order); %AR using Burg method % [H,F] = freqz(1,A,nfft,fs); % PSD=(abs(H).^2).*(variance/fs); %malik, p.67 %Method 2 [PSD,F]=pburg(y,AR_order,(nfft*2)-1,fs,'onesided'); %Method 3 % h=spectrum.burg; % hpsd = psd(h, y, 'NFFT', nfft, 'Fs', 2); % F=hpsd.Frequencies; % PSD=hpsd.Data; end function [PSD,F]=calcLomb(t,y,nfft,maxF) %calLomb - Calculates the PSD using Lomb-Scargle method. % %Inputs: %Outputs: %Calculate PSD deltaF=maxF/nfft; F = linspace(0.0,maxF-deltaF,nfft); PSD=lomb2(y,t,F,false); %calc lomb psd end function output=calcAreas(F,PSD,VLF,LF,HF,flagNorm) %calcAreas - Calulates areas/energy under the PSD curve within the freq %bands defined by VLF, LF, and HF. Returns areas/energies as ms^2, %percentage, and normalized units. Also returns LF/HF ratio. % %Inputs: % PSD: PSD vector % F: Freq vector % VLF, LF, HF: array containing VLF, LF, and HF freq limits % flagNormalize: option to normalize PSD to max(PSD) %Output: % %Usage: % % % Modified from Gary Clifford's ECG Toolbox: calc_lfhf.m if nargin<6 flagNorm=false; end %normalize PSD if needed if flagNorm PSD=PSD/max(PSD); end % find the indexes corresponding to the VLF, LF, and HF bands iVLF= (F>=VLF(1)) & (F<=VLF(2)); iLF = (F>=LF(1)) & (F<=LF(2)); iHF = (F>=HF(1)) & (F<=HF(2)); %Find peaks %VLF Peak tmpF=F(iVLF); tmppsd=PSD(iVLF); [pks,ipks] = zipeaks(tmppsd); if ~isempty(pks) [tmpMax i]=max(pks); peakVLF=tmpF(ipks(i)); else [tmpMax i]=max(tmppsd); peakVLF=tmpF(i); end %LF Peak tmpF=F(iLF); tmppsd=PSD(iLF); [pks,ipks] = zipeaks(tmppsd); if ~isempty(pks) [tmpMax i]=max(pks); peakLF=tmpF(ipks(i)); else [tmpMax i]=max(tmppsd); peakLF=tmpF(i); end %HF Peak tmpF=F(iHF); tmppsd=PSD(iHF); [pks,ipks] = zipeaks(tmppsd); if ~isempty(pks) [tmpMax i]=max(pks); peakHF=tmpF(ipks(i)); else [tmpMax i]=max(tmppsd); peakHF=tmpF(i); end % calculate raw areas (power under curve), within the freq bands (ms^2) aVLF=trapz(F(iVLF),PSD(iVLF)); aLF=trapz(F(iLF),PSD(iLF)); aHF=trapz(F(iHF),PSD(iHF)); aTotal=aVLF+aLF+aHF; %calculate areas relative to the total area (%) pVLF=(aVLF/aTotal)*100; pLF=(aLF/aTotal)*100; pHF=(aHF/aTotal)*100; %calculate normalized areas (relative to HF+LF, n.u.) nLF=aLF/(aLF+aHF); nHF=aHF/(aLF+aHF); %calculate LF/HF ratio lfhf =aLF/aHF; %create output structure if flagNorm output.aVLF=round(aVLF*1000)/1000; output.aLF=round(aLF*1000)/1000; output.aHF=round(aHF*1000)/1000; output.aTotal=round(aTotal*1000)/1000; else output.aVLF=round(aVLF*100)/100; % round output.aLF=round(aLF*100)/100; output.aHF=round(aHF*100)/100; output.aTotal=round(aTotal*100)/100; end output.pVLF=round(pVLF*10)/10; output.pLF=round(pLF*10)/10; output.pHF=round(pHF*10)/10;