Nothing
R/StationaryFrequencyAnalysis.R
In RHRV: Heart Rate Variability Analysis of ECG Data
Defines functions
selectRegressionRange
checkRegressionRange
EstimatePSDSlope
getEnergyInBand
CalculateEnergyInPSDBands
plotUnderCurve
plotFrequencyBands
buildPeriodogramLegend
PlotPSD
lombAdapter
spec.pgramAdapter
spec.arAdapter
CalculatePSD
Documented in
CalculateEnergyInPSDBands
CalculatePSD
EstimatePSDSlope
PlotPSD
############################## CalculatePSD ##########################
#' Spectral Density Estimation
#' @description Estimate the Power Spectral Density (PSD) of the RR time series.
#' @param HRVData Data structure that stores the beats register and information
#' related to it.
#' @param indexFreqAnalysis An integer referencing the data structure that will contain
#' the frequency analysis.
#' @param method String specifying the method used to estimate the spectral
#' density. Allowed methods are "pgram" (the default), "ar" and "lomb".
#' @param doPlot Plot the periodogram?
#' @param ... Further arguments to specific PSD estimation methods or
#' \code{\link[RHRV]{PlotPSD}}.
#' @details The "pgram" and "ar" methods use the \code{\link[stats]{spec.pgram}}
#' and \code{\link[stats]{spec.ar}} functions. Thus, the same arguments used
#' in \code{\link[stats]{spec.pgram}} or \code{\link[stats]{spec.ar}}
#' can be used when method is "pgram" or "ar", respectively.
#' The "lomb" is based in the \code{\link[lomb]{lsp}} and thus it accepts the
#' same parameters as this function.
#' @return The \emph{CalculatePSD} returns the \emph{HRVData} structure
#' containing a \emph{periodogram} field storing and PSD estimation
#' of the RR time series. When the "pgram" and "ar" methods are used the
#' \emph{periodogram} field is an object of class "spec". If "lomb" is used,
#' the \emph{periodogram} field is just a list. In any case the
#' \emph{periodogram} field will contain:
#' \itemize{
#' \item{freq: vector of frequencies at which the spectral density is estimated.}
#' \item{spec: spectral density estimation }
#' \item{series: name of the series}
#' \item{method: method used to calculate the spectrum}
#' }
#' @examples
#' \dontrun{
#' data(HRVData)
#' HRVData=BuildNIHR(HRVData)
#' HRVData=FilterNIHR(HRVData)
#' # Frequency analysis requires interpolated data (except Lomb)
#' HRVData=InterpolateNIHR(HRVData)
#' # Create a different freqAnalysis for each method
#' HRVData=CreateFreqAnalysis(HRVData)
#' HRVData=CalculatePSD(HRVData,1,"pgram",doPlot = F)
#'
#' HRVData=CreateFreqAnalysis(HRVData)
#' HRVData=CalculatePSD(HRVData,2,"pgram",spans=9, doPlot = F)
#'
#' HRVData=CreateFreqAnalysis(HRVData)
#' HRVData=CalculatePSD(HRVData,3,"ar",doPlot = F)
#'
#' HRVData=CreateFreqAnalysis(HRVData)
#' HRVData=CalculatePSD(HRVData,4,"lomb",doPlot = F)
#' # Plot the results
#' layout(matrix(c(1,2,3,4), 2, 2, byrow = TRUE))
#' PlotPSD(HRVData,1)
#' PlotPSD(HRVData,2)
#' PlotPSD(HRVData,3)
#' PlotPSD(HRVData,4)
#' }
#' @rdname CalculatePSD
#' @seealso \code{\link[stats]{spectrum}}, \code{\link[RHRV]{PlotPSD}}.
CalculatePSD =
function(HRVData, indexFreqAnalysis = length(HRVData$FreqAnalysis),
method=c("pgram","ar","lomb"), doPlot=T, ...){
CheckAnalysisIndex(indexFreqAnalysis, length(HRVData$FreqAnalysis),
"frequency")
switch(match.arg(method),
pgram = spec.pgramAdapter(HRVData,indexFreqAnalysis,doPlot,...),
ar = spec.arAdapter(HRVData,indexFreqAnalysis,doPlot,...),
lomb = lombAdapter(HRVData,indexFreqAnalysis,doPlot,...))
}
# Private Function for CalculatePSD. Performs periodogram
# estimation using AR modelling. It is just an adapter of the
# spec.ar function from the stats package
spec.arAdapter = function(HRVData, indexFreqAnalysis, doPlot = TRUE,
n.freq, order = NULL,
na.action = na.fail, method = "yule-walker",
...){
CheckInterpolation(HRVData)
VerboseMessage(HRVData$Verbose, "Calculating Periodogram using AR modelling")
# Compute the spectrogram using RR data
x = ts(1000.0/(HRVData$HR/60.0),frequency = HRVData$Freq_HR,
start=HRVData$Beat$Time[[1]])
periodogram = stats::spec.ar(x,n.freq,order,plot=F,na.action, method)
HRVData$FreqAnalysis[[indexFreqAnalysis]]$periodogram=periodogram
if (doPlot){
PlotPSD(HRVData,indexFreqAnalysis,...)
}
HRVData
}
# Private Function for CalculatePSD. Performs periodogram
# estimation using Daniell smoothers. It is just an adapter of the
# spec.pgram function from the stats package
spec.pgramAdapter = function(HRVData, indexFreqAnalysis,doPlot = TRUE,
spans = NULL, kernel=NULL, taper = 0.1,
pad = 0, fast = TRUE, demean = FALSE, detrend = TRUE,
na.action = na.fail, ...){
CheckInterpolation(HRVData)
VerboseMessage(HRVData$Verbose,
"Calculating Periodogram using DFT + Daniell smoothers")
# Compute the spectrogram using RR data
x = ts(1000.0/(HRVData$HR/60.0),frequency = HRVData$Freq_HR,
start=HRVData$Beat$Time[[1]])
periodogram = stats::spec.pgram(x, spans, kernel, taper,
pad, fast, demean, detrend,
plot = FALSE, na.action)
HRVData$FreqAnalysis[[indexFreqAnalysis]]$periodogram=periodogram
if (doPlot){
PlotPSD(HRVData,indexFreqAnalysis,...)
}
HRVData
}
lombAdapter = function(HRVData,indexFreqAnalysis,doPlot=TRUE,
from = 0, to = 0.5, type = c("frequency", "period"),
ofac = 1, alpha = 0.01,...){
kCopyNames = c("n","type","ofac","n.out", "alpha",
"peak.at","p.value")
CheckNIHR(HRVData)
VerboseMessage(HRVData$Verbose,
"Calculating Periodogram using Lomb periodogram")
# Compute the spectrogram using RR data
lombper = lomb::lsp(x=HRVData$Beat$RR,times=HRVData$Beat$Time,
from, to, type,
ofac, alpha, plot=FALSE, ...)
periodogram = lombper[kCopyNames]
periodogram$series = "rr"
periodogram$method = "Lomb Periodogram"
periodogram$freq = lombper$scanned
# renormalize spectrum
var.RR = var(HRVData$Beat$RR)
periodogram$spec = lombper$power * var.RR
periodogram$sig.level = lombper$sig.level * var.RR
periodogram$peak = lombper$peak * var.RR
HRVData$FreqAnalysis[[indexFreqAnalysis]]$periodogram=periodogram
if (doPlot){
PlotPSD(HRVData,indexFreqAnalysis,...)
}
HRVData
}
#' Plot Spectral Density Estimation
#' @description Plot the PSD estimate of the RR time series distinguishing
#' the different frequency bands with different colurs.
#' @param HRVData Data structure that stores the beats register and information
#' related to it.
#' @param indexFreqAnalysis An integer referencing the data structure that contains
#' the PSD analysis.
#' @param ULFmin Lower limit ULF band used for distinguish the ULF band.
#' @param ULFmax Upper limit ULF band used for distinguish the ULF band.
#' @param VLFmin Lower limit VLF band.
#' @param VLFmax Upper limit VLF band.
#' @param LFmin Lower limit LF band.
#' @param LFmax Upper limit LF band.
#' @param HFmin Lower limit HF band.
#' @param HFmax Upper limit HF band.
#' @param type 1-character string giving the type of plot desired. See \code{\link[graphics]{plot.default}}.
#' @param xlim the x limits (x1, x2) of the plot. See \code{\link[graphics]{plot.default}}.
#' @param ylim the y limits of the plot.
#' @param log a character string which contains "x" if the x axis is to be
#' logarithmic, "y" if the y axis is to be logarithmic and "xy" or "yx" if
#' both axes are to be logarithmic. Default: "y".
#' @param main a main title for the plot. See \code{\link[graphics]{plot.default}}.
#' @param xlab a label for the x axis. See \code{\link[graphics]{plot.default}}.
#' @param ylab a label for the y axis. See \code{\link[graphics]{plot.default}}.
#' @param addLegend add a simple legend? Default: True.
#' @param addSigLevel Logical value (only used with the lomb method). If true an
#' horizontal line limiting the significance level is included ( Powers >
#' sig.level can be considered significant peaks). See \code{\link[lomb]{lsp}}.
#' @param usePalette A new palette of colors for plotting the frequency bands.
#' @param ... graphical parameters. See \code{\link[graphics]{plot.default}}.
#' @examples
#' \dontrun{
#' data(HRVData)
#' HRVData=BuildNIHR(HRVData)
#' HRVData=FilterNIHR(HRVData)
#' # Frequency analysis requires interpolated data (except Lomb)
#' HRVData=InterpolateNIHR(HRVData)
#' # Create a different freqAnalysis for each method
#' HRVData=CreateFreqAnalysis(HRVData)
#' HRVData=CalculatePSD(HRVData,1,"pgram",doPlot = F)
#'
#' HRVData=CalculatePSD(HRVData,2,"pgram",spans=9,doPlot = F)
#'
#' HRVData=CreateFreqAnalysis(HRVData)
#' HRVData=CalculatePSD(HRVData,3,"ar",doPlot = F)
#'
#' HRVData=CreateFreqAnalysis(HRVData)
#' HRVData=CalculatePSD(HRVData,4,"lomb",doPlot = F)
#' # Plot the results
#' layout(matrix(c(1,2,3,4), 2, 2, byrow = TRUE))
#' PlotPSD(HRVData,1)
#' PlotPSD(HRVData,2)
#' PlotPSD(HRVData,3)
#' PlotPSD(HRVData,4)
#' }
#' @rdname PlotPSD
#' @seealso \code{\link[stats]{spectrum}}, \code{\link[lomb]{lsp}}, \code{\link[RHRV]{CalculatePSD}}.
#'
PlotPSD = function(HRVData,
indexFreqAnalysis=length(HRVData$FreqAnalysis),
ULFmin=0, ULFmax=0.03,
VLFmin=0.03, VLFmax=0.05,
LFmin=0.05, LFmax=0.15,
HFmin=0.15, HFmax=0.4,
log="y", type="l",
xlab="Frequency (Hz) ", ylab="Spectrum",
main=NULL,
xlim=c(min(ULFmin, ULFmax,
VLFmin, VLFmax,
LFmin, LFmax,
HFmin, HFmax),
max(ULFmin, ULFmax,
VLFmin, VLFmax,
LFmin, LFmax,
HFmin, HFmax)),
ylim=NULL,
addLegend=TRUE,
addSigLevel=TRUE,
usePalette=c("#000000", "#E69F00", "#56B4E9", "#009E73",
"#F0E442"),
...){
CheckAnalysisIndex(indexFreqAnalysis, length(HRVData$FreqAnalysis),
"frequency")
# check if the periodogram has been computed
CheckPeriodogram(HRVData, indexFreqAnalysis)
method = HRVData$FreqAnalysis[[indexFreqAnalysis]]$periodogram$method
freq = HRVData$FreqAnalysis[[indexFreqAnalysis]]$periodogram$freq
spect = HRVData$FreqAnalysis[[indexFreqAnalysis]]$periodogram$spec
# if we are plotting a lomb periodogram, check if the plotting range is
# larger than the computed frequencies
if ((method == "Lomb Periodogram") && (max(freq) < max(ULFmin, ULFmax,
VLFmin, VLFmax,
LFmin, LFmax,
HFmin, HFmax))){
warning(paste("Plotting larger frequency range than computed.",
"Increase the 'to' parameter in CalculatePSD()"))
}
# modify main if not supplied
if(is.null(main)) main = method
if(is.null(ylim)) ylim = range(spect[is.finite(spect)])
# avoid using 0 in freq if log="x", log="xy" or log="yx"
if (log== "x" || log == "xy" || log == "yx"){
if (freq[[1]] == 0){
freq = freq[-1]
spect = spect[-1]
}
if (xlim[[1]] == 0){
xlim[[1]] = freq[[1]]
}
}
plot(freq,spect, xlab=xlab, ylab=ylab, main=main,
xlim=xlim, ylim=ylim,
type=type,
log=log,...)
oldPalette = palette()
on.exit(palette(oldPalette))
palette(usePalette)
plotFrequencyBands(freq,spect,
ULFmin, ULFmax,
VLFmin, VLFmax,
LFmin, LFmax,
HFmin, HFmax, ylim)
if(addLegend){
legendVars = buildPeriodogramLegend(ULFmin, ULFmax, VLFmin, VLFmax,
LFmin, LFmax, HFmin, HFmax)
legend("topright",
legendVars$msg,
lwd=3,col=legendVars$col,bty="n")
}
if(addSigLevel && (method == "Lomb Periodogram")){
abline(h=HRVData$FreqAnalysis[[indexFreqAnalysis]]$periodogram$sig.level,
lwd=2,col=1,lty=2)
}
}
buildPeriodogramLegend = function(ULFmin, ULFmax, VLFmin, VLFmax,
LFmin, LFmax, HFmin, HFmax){
msg = c()
col = c()
if ( (!is.null(ULFmin)) && (!is.null(ULFmax)) ){
msg = c(msg, paste( "ULF: [",round(ULFmin,digits=2),", ",round(ULFmax,digits=2),"] Hz",sep=""))
col = c(col,2)
}
if ( (!is.null(VLFmin)) && (!is.null(VLFmax)) ){
msg = c(msg, paste( "VLF: [",round(VLFmin,digits=2),", ",round(VLFmax,digits=2),"] Hz",sep=""))
col = c(col,3)
}
if ( (!is.null(LFmin)) && (!is.null(LFmax)) ){
msg = c(msg, paste( "LF: [",round(LFmin,digits=2),", ",round(LFmax,digits=2),"] Hz",sep=""))
col = c(col,4)
}
if ( (!is.null(HFmin)) && (!is.null(HFmax)) ){
msg = c(msg, paste( "HF: [",round(HFmin,digits=2),", ",round(HFmax,digits=2),"] Hz",sep=""))
col = c(col,5)
}
list(msg=msg,col=col)
}
# Private function for distinguishing the different frequency bands of the
# PSD
plotFrequencyBands = function(freq,spect,
ULFmin, ULFmax,
VLFmin, VLFmax,
LFmin, LFmax,
HFmin, HFmax, ylim){
# modify boundaries to avoid gaps in the plot between
# bands
if( (!is.null(ULFmax)) && (!is.null(VLFmin)) && (ULFmax == VLFmin) )
ULFmax = VLFmin = tail(freq[freq <= ULFmax],1)
if( (!is.null(VLFmax)) && (!is.null(LFmin)) && (VLFmax == LFmin) )
VLFmax = LFmin = tail(freq[freq <= VLFmax],1)
if( (!is.null(LFmax)) && (!is.null(HFmin)) && (LFmax == HFmin) )
LFmax = HFmin = tail(freq[freq <= LFmax],1)
plotUnderCurve(freq, spect, xlim=c(ULFmin,ULFmax), ylim=ylim, col=2)
plotUnderCurve(freq, spect, xlim=c(VLFmin,VLFmax), ylim=ylim, col=3)
plotUnderCurve(freq, spect, xlim=c(LFmin,LFmax), ylim=ylim, col=4)
plotUnderCurve(freq, spect, xlim=c(HFmin,HFmax), ylim=ylim, col=5)
}
# Private function for highlightning the area under a curve
plotUnderCurve = function(x,y,xlim,ylim,col){
# if some member of xlim is null => length(xlim) < 2
if ( length(xlim) == 2) {
index = which(x >= xlim[[1]] & x <= xlim[[2]])
if (length(index) == 0) {
stop("Band out of bounds")
}
coordPolX = x[index]
coordPolY = y[index]
# Close the polygon
coordPolX = c(coordPolX[[1]], coordPolX, tail(coordPolX,1))
coordPolY = c(min(ylim), coordPolY, min(ylim))
polygon(coordPolX ,coordPolY, col = col)
}
}
########################### CalculatePSDBandsEnergy ###########################
#' CalculateSPDBandsEnergy
#' @description Calculates the Energy in the bands of the Power Spectral Density
#' (PSD).
#' @inheritParams PlotPSD
#' @return A vector containing the energy of the ULF, VLF, LF and HF bands in the
#' PSD.
#' @examples
#' \dontrun{
#' data(HRVData)
#' HRVData=BuildNIHR(HRVData)
#' HRVData=FilterNIHR(HRVData)
#' # Frequency analysis requires interpolated data (except Lomb)
#' HRVData=InterpolateNIHR(HRVData)
#' HRVData=CreateFreqAnalysis(HRVData)
#' HRVData=CalculatePSD(HRVData,1,"pgram",doPlot = F)
#' # get Energy in the default ULF, VLF and LF frequency bands.
#' # We modify the limits for the HF band
#' CalculateEnergyInPSDBands(HRVData, 1, HFmin = 0.15, HFmax = 0.3)
#' }
#' @seealso \code{\link{PlotPSD}}, \code{\link{CalculatePSD}}.
CalculateEnergyInPSDBands = function(HRVData,
indexFreqAnalysis = length(HRVData$FreqAnalysis),
ULFmin=0, ULFmax=0.03,
VLFmin=0.03, VLFmax=0.05,
LFmin=0.05, LFmax=0.15,
HFmin=0.15, HFmax=0.4){
CheckAnalysisIndex(indexFreqAnalysis, length(HRVData$FreqAnalysis),
"frequency")
# check if the periodogram has been computed
CheckPeriodogram(HRVData, indexFreqAnalysis)
psd = HRVData$FreqAnalysis[[indexFreqAnalysis]]$periodogram
c(getEnergyInBand(psd,ULFmin, ULFmax),
getEnergyInBand(psd,VLFmin, VLFmax),
getEnergyInBand(psd,LFmin, LFmax),
getEnergyInBand(psd,HFmin, HFmax))
}
getEnergyInBand = function(psd,freq.min, freq.max){
if (freq.min == freq.max){
return(0)
}
if (freq.min > freq.max) {
stop("The band limits are not properly defined")
}
if (freq.min > max(psd$freq)){
stop("The band's lower frequency is greater than the maximum reachable frequency")
}
if (freq.max < min(psd$freq)){
stop("The band's higher frequency is smaller than the minimun reachable frequency")
}
indx = getIndexValuesInRange(psd$freq, c(freq.min,freq.max))
delta = median(diff(psd$freq))
if (length(indx) == 0) {
stop("No values in the specified frequency band")
}else{
return(sum(psd$spec[indx] * delta))
}
}
############################## EstimatePSDSlope ##########################
#' Estimate the slope of the Power Spectral Density (PSD).
#' @description Estimate the slope of the Power Spectral Density (PSD) of the
#' RR time series.
#' @param HRVData Data structure that stores the beats register and information
#' related to it.
#' @param indexFreqAnalysis An integer referencing the periodogram that will be
#' used for estimating the spectral index.
#' @param indexNonLinearAnalysis An integer referencing the structure that will
#' store the resulting estimations.
#' @param regressionRange Range of frequencies in which the regression will be
#' performed. Default is c(1e-4, 1e-2) Hz.
#' @param doPlot Plot the periodogram and the least-squares fit?
#' @param main Title for the plot.
#' @param xlab Title for the x axis.
#' @param ylab Title for the y axis.
#' @param pch Symbol for the plotting points.
#' @param log A character string which contains "x" if the x axis is to be
#' logarithmic, "y" if the y axis is to be logarithmic and "xy" or "yx" if
#' both axes are to be logarithmic (default).
#' @param ... Other arguments for the plotting function.
#' @details The power spectrum of most physiological signals fulfils
#' \eqn{S(f)=Cf^{-\beta}}{S(f)=C*f^-B} (1/f spectrum). This function estimates
#' the \eqn{\beta}{B} exponent, which is usually referred to as the spectral
#' index.
#' @note It should be noted that the PSD must be estimated prior to the
#' use of this function. We do not recommend the use of the AR spectrum when
#' estimating the spectral index.
#' @return The \emph{EstimatePSDSlope} returns the \emph{HRVData} structure
#' containing a \emph{PSDSlope} field storing the spectral index and the
#' proper Hurst exponent.
#' @examples
#' \dontrun{
#' data(HRVProcessedData)
#' # use other name for convenience
#' HRVData=HRVProcessedData
#' # Estimate the periodogram
#' HRVData=CreateFreqAnalysis(HRVData)
#' HRVData=CalculatePSD(HRVData,1,"pgram",doPlot = T,log="xy")
#' HRVData=CreateNonLinearAnalysis(HRVData)
#' HRVData=SetVerbose(HRVData,T)
#' HRVData=EstimatePSDSlope(HRVData,1,1,
#' regressionRange=c(5e-4,1e-2))
#' }
#' @references Voss, Andreas, et al. "Methods derived from nonlinear dynamics
#' for analysing heart rate variability." Philosophical Transactions of the
#' Royal Society A: Mathematical, Physical and Engineering Sciences 367.1887
#' (2009): 277-296.
#'
#' Eke, A., Herman, P., Kocsis, L., & Kozak, L. R. (2002). Fractal
#' characterization of complexity in temporal physiological signals.
#' Physiological measurement, 23(1), R1.
#' @seealso \code{\link[stats]{spectrum}},\code{\link[lomb]{lsp}},
#' \code{\link[RHRV]{CalculatePSD}}.
EstimatePSDSlope = function(HRVData,
indexFreqAnalysis = length(HRVData$FreqAnalysis),
indexNonLinearAnalysis = length(HRVData$NonLinearAnalysis),
regressionRange=NULL,doPlot = T,
main="PSD power law", xlab = "Frequency (Hz)",
ylab="Spectrum", pch = NULL,
log="xy",...){
CheckAnalysisIndex(indexFreqAnalysis,length(HRVData$FreqAnalysis),"frequency")
CheckAnalysisIndex(indexNonLinearAnalysis,length(HRVData$NonLinearAnalysis),
"nonlinear")
CheckPeriodogram(HRVData, indexFreqAnalysis)
VerboseMessage(HRVData$Verbose,"Calculating spectral index")
freq = HRVData$FreqAnalysis[[indexFreqAnalysis]]$periodogram$freq
spec = HRVData$FreqAnalysis[[indexFreqAnalysis]]$periodogram$spec
freqRange = range(freq)
if (is.null(regressionRange)){
regressionRange = selectRegressionRange(freqRange)
}
checkRegressionRange(freqRange,regressionRange)
indx = getIndexValuesInRange(freq, regressionRange)
lfreq = log(freq)[indx]
lspec = log(spec)[indx]
fit = lm(lspec ~ lfreq)
# plot
if (doPlot){
if (is.null(pch)){
pch = rep(1,length(freq))
pch[indx] = 3
}
VerboseMessage(HRVData$Verbose,"Plotting PSD")
plot(freq, spec,log=log,main=main,xlab=xlab,ylab=ylab,pch=pch,...)
lines(exp(lfreq),exp(fit$fitted.values),col=2)
}
# get values of beta and H
beta = -coef(fit)[[2]]
VerboseMessage(HRVData$Verbose,paste("Spectral index =",rhrvFormat(beta)))
if (beta > 1){
HfBm = (beta - 1) / 2
HfGn = NULL
VerboseMessage(HRVData$Verbose,paste("H_fBm =",rhrvFormat(HfBm)))
}else{
HfBm = NULL
HfGn = (beta + 1) / 2
VerboseMessage(HRVData$Verbose,paste("H_fGn =",rhrvFormat(HfGn)))
}
HRVData$NonLinearAnalysis[[indexNonLinearAnalysis]]$PSDSlope =
list(computations=list(indexFreqAnalysis=indexFreqAnalysis,
regressionRange=regressionRange),
statistic=list(beta=beta,HfBm=HfBm,HfGn=HfGn))
HRVData
}
checkRegressionRange = function(freqRange, regressionRange){
if ( (freqRange[[1]] > regressionRange[[1]]) ||
(freqRange[[2]] < regressionRange[[2]])){
stop("Invalid Regression range !!")
}
}
selectRegressionRange = function(freqRange){
# typical range is 1e-4 to 1e-2 Hz (Voss et al)
kTypicalMinFreq = 1e-4
kTypicalMaxFreq = 1e-2
if ( (freqRange[[1]] <= kTypicalMinFreq) && (freqRange[[2]] >= kTypicalMaxFreq) ){
regressionRange = c(1e-4,1e-2)
}else{
minFreq = max(kTypicalMinFreq, freqRange[[1]])
maxFreq = min(kTypicalMaxFreq, freqRange[[2]])
regressionRange = c(minFreq, maxFreq)
}
regressionRange
}
Try the RHRV package in your browser
Any scripts or data that you put into this service are public.
RHRV documentation built on Jan. 16, 2024, 3:05 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions selectRegressionRange checkRegressionRange EstimatePSDSlope getEnergyInBand CalculateEnergyInPSDBands plotUnderCurve plotFrequencyBands buildPeriodogramLegend PlotPSD lombAdapter spec.pgramAdapter spec.arAdapter CalculatePSD
Documented in CalculateEnergyInPSDBands CalculatePSD EstimatePSDSlope PlotPSD
############################## CalculatePSD ##########################
#' Spectral Density Estimation
#' @description Estimate the Power Spectral Density (PSD) of the RR time series.
#' @param HRVData Data structure that stores the beats register and information
#' related to it.
#' @param indexFreqAnalysis An integer referencing the data structure that will contain
#' the frequency analysis.
#' @param method String specifying the method used to estimate the spectral
#' density. Allowed methods are "pgram" (the default), "ar" and "lomb".
#' @param doPlot Plot the periodogram?
#' @param ... Further arguments to specific PSD estimation methods or
#' \code{\link[RHRV]{PlotPSD}}.
#' @details The "pgram" and "ar" methods use the \code{\link[stats]{spec.pgram}}
#' and \code{\link[stats]{spec.ar}} functions. Thus, the same arguments used
#' in \code{\link[stats]{spec.pgram}} or \code{\link[stats]{spec.ar}}
#' can be used when method is "pgram" or "ar", respectively.
#' The "lomb" is based in the \code{\link[lomb]{lsp}} and thus it accepts the
#' same parameters as this function.
#' @return The \emph{CalculatePSD} returns the \emph{HRVData} structure
#' containing a \emph{periodogram} field storing and PSD estimation
#' of the RR time series. When the "pgram" and "ar" methods are used the
#' \emph{periodogram} field is an object of class "spec". If "lomb" is used,
#' the \emph{periodogram} field is just a list. In any case the
#' \emph{periodogram} field will contain:
#' \itemize{
#' \item{freq: vector of frequencies at which the spectral density is estimated.}
#' \item{spec: spectral density estimation }
#' \item{series: name of the series}
#' \item{method: method used to calculate the spectrum}
#' }
#' @examples
#' \dontrun{
#' data(HRVData)
#' HRVData=BuildNIHR(HRVData)
#' HRVData=FilterNIHR(HRVData)
#' # Frequency analysis requires interpolated data (except Lomb)
#' HRVData=InterpolateNIHR(HRVData)
#' # Create a different freqAnalysis for each method
#' HRVData=CreateFreqAnalysis(HRVData)
#' HRVData=CalculatePSD(HRVData,1,"pgram",doPlot = F)
#'
#' HRVData=CreateFreqAnalysis(HRVData)
#' HRVData=CalculatePSD(HRVData,2,"pgram",spans=9, doPlot = F)
#'
#' HRVData=CreateFreqAnalysis(HRVData)
#' HRVData=CalculatePSD(HRVData,3,"ar",doPlot = F)
#'
#' HRVData=CreateFreqAnalysis(HRVData)
#' HRVData=CalculatePSD(HRVData,4,"lomb",doPlot = F)
#' # Plot the results
#' layout(matrix(c(1,2,3,4), 2, 2, byrow = TRUE))
#' PlotPSD(HRVData,1)
#' PlotPSD(HRVData,2)
#' PlotPSD(HRVData,3)
#' PlotPSD(HRVData,4)
#' }
#' @rdname CalculatePSD
#' @seealso \code{\link[stats]{spectrum}}, \code{\link[RHRV]{PlotPSD}}.
CalculatePSD =
function(HRVData, indexFreqAnalysis = length(HRVData$FreqAnalysis),
method=c("pgram","ar","lomb"), doPlot=T, ...){
CheckAnalysisIndex(indexFreqAnalysis, length(HRVData$FreqAnalysis),
"frequency")
switch(match.arg(method),
pgram = spec.pgramAdapter(HRVData,indexFreqAnalysis,doPlot,...),
ar = spec.arAdapter(HRVData,indexFreqAnalysis,doPlot,...),
lomb = lombAdapter(HRVData,indexFreqAnalysis,doPlot,...))
}
# Private Function for CalculatePSD. Performs periodogram
# estimation using AR modelling. It is just an adapter of the
# spec.ar function from the stats package
spec.arAdapter = function(HRVData, indexFreqAnalysis, doPlot = TRUE,
n.freq, order = NULL,
na.action = na.fail, method = "yule-walker",
...){
CheckInterpolation(HRVData)
VerboseMessage(HRVData$Verbose, "Calculating Periodogram using AR modelling")
# Compute the spectrogram using RR data
x = ts(1000.0/(HRVData$HR/60.0),frequency = HRVData$Freq_HR,
start=HRVData$Beat$Time[[1]])
periodogram = stats::spec.ar(x,n.freq,order,plot=F,na.action, method)
HRVData$FreqAnalysis[[indexFreqAnalysis]]$periodogram=periodogram
if (doPlot){
PlotPSD(HRVData,indexFreqAnalysis,...)
}
HRVData
}
# Private Function for CalculatePSD. Performs periodogram
# estimation using Daniell smoothers. It is just an adapter of the
# spec.pgram function from the stats package
spec.pgramAdapter = function(HRVData, indexFreqAnalysis,doPlot = TRUE,
spans = NULL, kernel=NULL, taper = 0.1,
pad = 0, fast = TRUE, demean = FALSE, detrend = TRUE,
na.action = na.fail, ...){
CheckInterpolation(HRVData)
VerboseMessage(HRVData$Verbose,
"Calculating Periodogram using DFT + Daniell smoothers")
# Compute the spectrogram using RR data
x = ts(1000.0/(HRVData$HR/60.0),frequency = HRVData$Freq_HR,
start=HRVData$Beat$Time[[1]])
periodogram = stats::spec.pgram(x, spans, kernel, taper,
pad, fast, demean, detrend,
plot = FALSE, na.action)
HRVData$FreqAnalysis[[indexFreqAnalysis]]$periodogram=periodogram
if (doPlot){
PlotPSD(HRVData,indexFreqAnalysis,...)
}
HRVData
}
lombAdapter = function(HRVData,indexFreqAnalysis,doPlot=TRUE,
from = 0, to = 0.5, type = c("frequency", "period"),
ofac = 1, alpha = 0.01,...){
kCopyNames = c("n","type","ofac","n.out", "alpha",
"peak.at","p.value")
CheckNIHR(HRVData)
VerboseMessage(HRVData$Verbose,
"Calculating Periodogram using Lomb periodogram")
# Compute the spectrogram using RR data
lombper = lomb::lsp(x=HRVData$Beat$RR,times=HRVData$Beat$Time,
from, to, type,
ofac, alpha, plot=FALSE, ...)
periodogram = lombper[kCopyNames]
periodogram$series = "rr"
periodogram$method = "Lomb Periodogram"
periodogram$freq = lombper$scanned
# renormalize spectrum
var.RR = var(HRVData$Beat$RR)
periodogram$spec = lombper$power * var.RR
periodogram$sig.level = lombper$sig.level * var.RR
periodogram$peak = lombper$peak * var.RR
HRVData$FreqAnalysis[[indexFreqAnalysis]]$periodogram=periodogram
if (doPlot){
PlotPSD(HRVData,indexFreqAnalysis,...)
}
HRVData
}
#' Plot Spectral Density Estimation
#' @description Plot the PSD estimate of the RR time series distinguishing
#' the different frequency bands with different colurs.
#' @param HRVData Data structure that stores the beats register and information
#' related to it.
#' @param indexFreqAnalysis An integer referencing the data structure that contains
#' the PSD analysis.
#' @param ULFmin Lower limit ULF band used for distinguish the ULF band.
#' @param ULFmax Upper limit ULF band used for distinguish the ULF band.
#' @param VLFmin Lower limit VLF band.
#' @param VLFmax Upper limit VLF band.
#' @param LFmin Lower limit LF band.
#' @param LFmax Upper limit LF band.
#' @param HFmin Lower limit HF band.
#' @param HFmax Upper limit HF band.
#' @param type 1-character string giving the type of plot desired. See \code{\link[graphics]{plot.default}}.
#' @param xlim the x limits (x1, x2) of the plot. See \code{\link[graphics]{plot.default}}.
#' @param ylim the y limits of the plot.
#' @param log a character string which contains "x" if the x axis is to be
#' logarithmic, "y" if the y axis is to be logarithmic and "xy" or "yx" if
#' both axes are to be logarithmic. Default: "y".
#' @param main a main title for the plot. See \code{\link[graphics]{plot.default}}.
#' @param xlab a label for the x axis. See \code{\link[graphics]{plot.default}}.
#' @param ylab a label for the y axis. See \code{\link[graphics]{plot.default}}.
#' @param addLegend add a simple legend? Default: True.
#' @param addSigLevel Logical value (only used with the lomb method). If true an
#' horizontal line limiting the significance level is included ( Powers >
#' sig.level can be considered significant peaks). See \code{\link[lomb]{lsp}}.
#' @param usePalette A new palette of colors for plotting the frequency bands.
#' @param ... graphical parameters. See \code{\link[graphics]{plot.default}}.
#' @examples
#' \dontrun{
#' data(HRVData)
#' HRVData=BuildNIHR(HRVData)
#' HRVData=FilterNIHR(HRVData)
#' # Frequency analysis requires interpolated data (except Lomb)
#' HRVData=InterpolateNIHR(HRVData)
#' # Create a different freqAnalysis for each method
#' HRVData=CreateFreqAnalysis(HRVData)
#' HRVData=CalculatePSD(HRVData,1,"pgram",doPlot = F)
#'
#' HRVData=CalculatePSD(HRVData,2,"pgram",spans=9,doPlot = F)
#'
#' HRVData=CreateFreqAnalysis(HRVData)
#' HRVData=CalculatePSD(HRVData,3,"ar",doPlot = F)
#'
#' HRVData=CreateFreqAnalysis(HRVData)
#' HRVData=CalculatePSD(HRVData,4,"lomb",doPlot = F)
#' # Plot the results
#' layout(matrix(c(1,2,3,4), 2, 2, byrow = TRUE))
#' PlotPSD(HRVData,1)
#' PlotPSD(HRVData,2)
#' PlotPSD(HRVData,3)
#' PlotPSD(HRVData,4)
#' }
#' @rdname PlotPSD
#' @seealso \code{\link[stats]{spectrum}}, \code{\link[lomb]{lsp}}, \code{\link[RHRV]{CalculatePSD}}.
#'
PlotPSD = function(HRVData,
indexFreqAnalysis=length(HRVData$FreqAnalysis),
ULFmin=0, ULFmax=0.03,
VLFmin=0.03, VLFmax=0.05,
LFmin=0.05, LFmax=0.15,
HFmin=0.15, HFmax=0.4,
log="y", type="l",
xlab="Frequency (Hz) ", ylab="Spectrum",
main=NULL,
xlim=c(min(ULFmin, ULFmax,
VLFmin, VLFmax,
LFmin, LFmax,
HFmin, HFmax),
max(ULFmin, ULFmax,
VLFmin, VLFmax,
LFmin, LFmax,
HFmin, HFmax)),
ylim=NULL,
addLegend=TRUE,
addSigLevel=TRUE,
usePalette=c("#000000", "#E69F00", "#56B4E9", "#009E73",
"#F0E442"),
...){
CheckAnalysisIndex(indexFreqAnalysis, length(HRVData$FreqAnalysis),
"frequency")
# check if the periodogram has been computed
CheckPeriodogram(HRVData, indexFreqAnalysis)
method = HRVData$FreqAnalysis[[indexFreqAnalysis]]$periodogram$method
freq = HRVData$FreqAnalysis[[indexFreqAnalysis]]$periodogram$freq
spect = HRVData$FreqAnalysis[[indexFreqAnalysis]]$periodogram$spec
# if we are plotting a lomb periodogram, check if the plotting range is
# larger than the computed frequencies
if ((method == "Lomb Periodogram") && (max(freq) < max(ULFmin, ULFmax,
VLFmin, VLFmax,
LFmin, LFmax,
HFmin, HFmax))){
warning(paste("Plotting larger frequency range than computed.",
"Increase the 'to' parameter in CalculatePSD()"))
}
# modify main if not supplied
if(is.null(main)) main = method
if(is.null(ylim)) ylim = range(spect[is.finite(spect)])
# avoid using 0 in freq if log="x", log="xy" or log="yx"
if (log== "x" || log == "xy" || log == "yx"){
if (freq[[1]] == 0){
freq = freq[-1]
spect = spect[-1]
}
if (xlim[[1]] == 0){
xlim[[1]] = freq[[1]]
}
}
plot(freq,spect, xlab=xlab, ylab=ylab, main=main,
xlim=xlim, ylim=ylim,
type=type,
log=log,...)
oldPalette = palette()
on.exit(palette(oldPalette))
palette(usePalette)
plotFrequencyBands(freq,spect,
ULFmin, ULFmax,
VLFmin, VLFmax,
LFmin, LFmax,
HFmin, HFmax, ylim)
if(addLegend){
legendVars = buildPeriodogramLegend(ULFmin, ULFmax, VLFmin, VLFmax,
LFmin, LFmax, HFmin, HFmax)
legend("topright",
legendVars$msg,
lwd=3,col=legendVars$col,bty="n")
}
if(addSigLevel && (method == "Lomb Periodogram")){
abline(h=HRVData$FreqAnalysis[[indexFreqAnalysis]]$periodogram$sig.level,
lwd=2,col=1,lty=2)
}
}
buildPeriodogramLegend = function(ULFmin, ULFmax, VLFmin, VLFmax,
LFmin, LFmax, HFmin, HFmax){
msg = c()
col = c()
if ( (!is.null(ULFmin)) && (!is.null(ULFmax)) ){
msg = c(msg, paste( "ULF: [",round(ULFmin,digits=2),", ",round(ULFmax,digits=2),"] Hz",sep=""))
col = c(col,2)
}
if ( (!is.null(VLFmin)) && (!is.null(VLFmax)) ){
msg = c(msg, paste( "VLF: [",round(VLFmin,digits=2),", ",round(VLFmax,digits=2),"] Hz",sep=""))
col = c(col,3)
}
if ( (!is.null(LFmin)) && (!is.null(LFmax)) ){
msg = c(msg, paste( "LF: [",round(LFmin,digits=2),", ",round(LFmax,digits=2),"] Hz",sep=""))
col = c(col,4)
}
if ( (!is.null(HFmin)) && (!is.null(HFmax)) ){
msg = c(msg, paste( "HF: [",round(HFmin,digits=2),", ",round(HFmax,digits=2),"] Hz",sep=""))
col = c(col,5)
}
list(msg=msg,col=col)
}
# Private function for distinguishing the different frequency bands of the
# PSD
plotFrequencyBands = function(freq,spect,
ULFmin, ULFmax,
VLFmin, VLFmax,
LFmin, LFmax,
HFmin, HFmax, ylim){
# modify boundaries to avoid gaps in the plot between
# bands
if( (!is.null(ULFmax)) && (!is.null(VLFmin)) && (ULFmax == VLFmin) )
ULFmax = VLFmin = tail(freq[freq <= ULFmax],1)
if( (!is.null(VLFmax)) && (!is.null(LFmin)) && (VLFmax == LFmin) )
VLFmax = LFmin = tail(freq[freq <= VLFmax],1)
if( (!is.null(LFmax)) && (!is.null(HFmin)) && (LFmax == HFmin) )
LFmax = HFmin = tail(freq[freq <= LFmax],1)
plotUnderCurve(freq, spect, xlim=c(ULFmin,ULFmax), ylim=ylim, col=2)
plotUnderCurve(freq, spect, xlim=c(VLFmin,VLFmax), ylim=ylim, col=3)
plotUnderCurve(freq, spect, xlim=c(LFmin,LFmax), ylim=ylim, col=4)
plotUnderCurve(freq, spect, xlim=c(HFmin,HFmax), ylim=ylim, col=5)
}
# Private function for highlightning the area under a curve
plotUnderCurve = function(x,y,xlim,ylim,col){
# if some member of xlim is null => length(xlim) < 2
if ( length(xlim) == 2) {
index = which(x >= xlim[[1]] & x <= xlim[[2]])
if (length(index) == 0) {
stop("Band out of bounds")
}
coordPolX = x[index]
coordPolY = y[index]
# Close the polygon
coordPolX = c(coordPolX[[1]], coordPolX, tail(coordPolX,1))
coordPolY = c(min(ylim), coordPolY, min(ylim))
polygon(coordPolX ,coordPolY, col = col)
}
}
########################### CalculatePSDBandsEnergy ###########################
#' CalculateSPDBandsEnergy
#' @description Calculates the Energy in the bands of the Power Spectral Density
#' (PSD).
#' @inheritParams PlotPSD
#' @return A vector containing the energy of the ULF, VLF, LF and HF bands in the
#' PSD.
#' @examples
#' \dontrun{
#' data(HRVData)
#' HRVData=BuildNIHR(HRVData)
#' HRVData=FilterNIHR(HRVData)
#' # Frequency analysis requires interpolated data (except Lomb)
#' HRVData=InterpolateNIHR(HRVData)
#' HRVData=CreateFreqAnalysis(HRVData)
#' HRVData=CalculatePSD(HRVData,1,"pgram",doPlot = F)
#' # get Energy in the default ULF, VLF and LF frequency bands.
#' # We modify the limits for the HF band
#' CalculateEnergyInPSDBands(HRVData, 1, HFmin = 0.15, HFmax = 0.3)
#' }
#' @seealso \code{\link{PlotPSD}}, \code{\link{CalculatePSD}}.
CalculateEnergyInPSDBands = function(HRVData,
indexFreqAnalysis = length(HRVData$FreqAnalysis),
ULFmin=0, ULFmax=0.03,
VLFmin=0.03, VLFmax=0.05,
LFmin=0.05, LFmax=0.15,
HFmin=0.15, HFmax=0.4){
CheckAnalysisIndex(indexFreqAnalysis, length(HRVData$FreqAnalysis),
"frequency")
# check if the periodogram has been computed
CheckPeriodogram(HRVData, indexFreqAnalysis)
psd = HRVData$FreqAnalysis[[indexFreqAnalysis]]$periodogram
c(getEnergyInBand(psd,ULFmin, ULFmax),
getEnergyInBand(psd,VLFmin, VLFmax),
getEnergyInBand(psd,LFmin, LFmax),
getEnergyInBand(psd,HFmin, HFmax))
}
getEnergyInBand = function(psd,freq.min, freq.max){
if (freq.min == freq.max){
return(0)
}
if (freq.min > freq.max) {
stop("The band limits are not properly defined")
}
if (freq.min > max(psd$freq)){
stop("The band's lower frequency is greater than the maximum reachable frequency")
}
if (freq.max < min(psd$freq)){
stop("The band's higher frequency is smaller than the minimun reachable frequency")
}
indx = getIndexValuesInRange(psd$freq, c(freq.min,freq.max))
delta = median(diff(psd$freq))
if (length(indx) == 0) {
stop("No values in the specified frequency band")
}else{
return(sum(psd$spec[indx] * delta))
}
}
############################## EstimatePSDSlope ##########################
#' Estimate the slope of the Power Spectral Density (PSD).
#' @description Estimate the slope of the Power Spectral Density (PSD) of the
#' RR time series.
#' @param HRVData Data structure that stores the beats register and information
#' related to it.
#' @param indexFreqAnalysis An integer referencing the periodogram that will be
#' used for estimating the spectral index.
#' @param indexNonLinearAnalysis An integer referencing the structure that will
#' store the resulting estimations.
#' @param regressionRange Range of frequencies in which the regression will be
#' performed. Default is c(1e-4, 1e-2) Hz.
#' @param doPlot Plot the periodogram and the least-squares fit?
#' @param main Title for the plot.
#' @param xlab Title for the x axis.
#' @param ylab Title for the y axis.
#' @param pch Symbol for the plotting points.
#' @param log A character string which contains "x" if the x axis is to be
#' logarithmic, "y" if the y axis is to be logarithmic and "xy" or "yx" if
#' both axes are to be logarithmic (default).
#' @param ... Other arguments for the plotting function.
#' @details The power spectrum of most physiological signals fulfils
#' \eqn{S(f)=Cf^{-\beta}}{S(f)=C*f^-B} (1/f spectrum). This function estimates
#' the \eqn{\beta}{B} exponent, which is usually referred to as the spectral
#' index.
#' @note It should be noted that the PSD must be estimated prior to the
#' use of this function. We do not recommend the use of the AR spectrum when
#' estimating the spectral index.
#' @return The \emph{EstimatePSDSlope} returns the \emph{HRVData} structure
#' containing a \emph{PSDSlope} field storing the spectral index and the
#' proper Hurst exponent.
#' @examples
#' \dontrun{
#' data(HRVProcessedData)
#' # use other name for convenience
#' HRVData=HRVProcessedData
#' # Estimate the periodogram
#' HRVData=CreateFreqAnalysis(HRVData)
#' HRVData=CalculatePSD(HRVData,1,"pgram",doPlot = T,log="xy")
#' HRVData=CreateNonLinearAnalysis(HRVData)
#' HRVData=SetVerbose(HRVData,T)
#' HRVData=EstimatePSDSlope(HRVData,1,1,
#' regressionRange=c(5e-4,1e-2))
#' }
#' @references Voss, Andreas, et al. "Methods derived from nonlinear dynamics
#' for analysing heart rate variability." Philosophical Transactions of the
#' Royal Society A: Mathematical, Physical and Engineering Sciences 367.1887
#' (2009): 277-296.
#'
#' Eke, A., Herman, P., Kocsis, L., & Kozak, L. R. (2002). Fractal
#' characterization of complexity in temporal physiological signals.
#' Physiological measurement, 23(1), R1.
#' @seealso \code{\link[stats]{spectrum}},\code{\link[lomb]{lsp}},
#' \code{\link[RHRV]{CalculatePSD}}.
EstimatePSDSlope = function(HRVData,
indexFreqAnalysis = length(HRVData$FreqAnalysis),
indexNonLinearAnalysis = length(HRVData$NonLinearAnalysis),
regressionRange=NULL,doPlot = T,
main="PSD power law", xlab = "Frequency (Hz)",
ylab="Spectrum", pch = NULL,
log="xy",...){
CheckAnalysisIndex(indexFreqAnalysis,length(HRVData$FreqAnalysis),"frequency")
CheckAnalysisIndex(indexNonLinearAnalysis,length(HRVData$NonLinearAnalysis),
"nonlinear")
CheckPeriodogram(HRVData, indexFreqAnalysis)
VerboseMessage(HRVData$Verbose,"Calculating spectral index")
freq = HRVData$FreqAnalysis[[indexFreqAnalysis]]$periodogram$freq
spec = HRVData$FreqAnalysis[[indexFreqAnalysis]]$periodogram$spec
freqRange = range(freq)
if (is.null(regressionRange)){
regressionRange = selectRegressionRange(freqRange)
}
checkRegressionRange(freqRange,regressionRange)
indx = getIndexValuesInRange(freq, regressionRange)
lfreq = log(freq)[indx]
lspec = log(spec)[indx]
fit = lm(lspec ~ lfreq)
# plot
if (doPlot){
if (is.null(pch)){
pch = rep(1,length(freq))
pch[indx] = 3
}
VerboseMessage(HRVData$Verbose,"Plotting PSD")
plot(freq, spec,log=log,main=main,xlab=xlab,ylab=ylab,pch=pch,...)
lines(exp(lfreq),exp(fit$fitted.values),col=2)
}
# get values of beta and H
beta = -coef(fit)[[2]]
VerboseMessage(HRVData$Verbose,paste("Spectral index =",rhrvFormat(beta)))
if (beta > 1){
HfBm = (beta - 1) / 2
HfGn = NULL
VerboseMessage(HRVData$Verbose,paste("H_fBm =",rhrvFormat(HfBm)))
}else{
HfBm = NULL
HfGn = (beta + 1) / 2
VerboseMessage(HRVData$Verbose,paste("H_fGn =",rhrvFormat(HfGn)))
}
HRVData$NonLinearAnalysis[[indexNonLinearAnalysis]]$PSDSlope =
list(computations=list(indexFreqAnalysis=indexFreqAnalysis,
regressionRange=regressionRange),
statistic=list(beta=beta,HfBm=HfBm,HfGn=HfGn))
HRVData
}
checkRegressionRange = function(freqRange, regressionRange){
if ( (freqRange[[1]] > regressionRange[[1]]) ||
(freqRange[[2]] < regressionRange[[2]])){
stop("Invalid Regression range !!")
}
}
selectRegressionRange = function(freqRange){
# typical range is 1e-4 to 1e-2 Hz (Voss et al)
kTypicalMinFreq = 1e-4
kTypicalMaxFreq = 1e-2
if ( (freqRange[[1]] <= kTypicalMinFreq) && (freqRange[[2]] >= kTypicalMaxFreq) ){
regressionRange = c(1e-4,1e-2)
}else{
minFreq = max(kTypicalMinFreq, freqRange[[1]])
maxFreq = min(kTypicalMaxFreq, freqRange[[2]])
regressionRange = c(minFreq, maxFreq)
}
regressionRange
}
Try the RHRV package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.