Quantified Self Experiments / Toolbox - Cooking ECG / RR for time-varying HRV and Respiration

Daily HRV measurements seems to be an important part of my QS adventures, mostly being used to assess ANS status and describe some states like overtraining, sickness, recovery. I'm utilizing contextual information + ANS balance to fine-tune my daily decisions.

Devices & Apps

There are some devices and apps to collect HRV, i've experimented with these setups

  1. Wearables like Oura or Fitbit with nightly average seems to have low sensitivity to short-term daily stressors (you may read Andrew Flatt twitter posts about nightly HRV or read my investigations 1, 2, 3). I look at nightly averages as marker of long-term health and dont use them for daily decisions, since i have more sensitive tool.
  2. Polar H10 + water +  Kubios. The best choice for someone who want simplicity
  3. Polar H10 + water + Polar Equine for ECG Raw data. Raw data processed in Kubios HRV Standard. It's more complicated then just looking at 1 number from Kubios app, but you get your raw data and may cook it the way you want. Desktop Kubios provides more HRV parameters to analyse.
  4. Polar H10 + water + Polar Vantage V2 for Orthostatic test. No raw data and test takes 4 minutes. Lying 2 mins seems to be a loss of time (more details on topic) and after doing some time of 4 mins orthostatic test i've switched to only standing test, because it take 2 minutes (60s stabilization and 60 sec measurement) and seems to have good sensitivity to short-term stressors.
  5. Polar Verity Sense + Polar Sensor Logger for PPG raw data. Not practical since PPG too sensitive to movements. I wasnt able to get same values as H10 (worn simultaneously)
  6. Movesense MD seems to be the best option for EEG raw data (512 Hz sampling rate). I've found it by exploring Bruce Rogers blog. The signal quality is better than H10 and i've getting zero artifacts for measurements at rest (H10 sometimes have artifacts). Raw data can be processed in Kubios HRV or cooked with own processing (my choise).

Raw data cooking

My current preference is Movesense MD + Raw data at 512Hz. I have 2 minutes of ECG every day, artifact free and with nice ECG shape. To cook raw data a few steps will be needed:

  1. ECG: preprocess, clean, asess quality, detect R-peaks and calculate RR intervals
  2. RRi: find artifacts, correct them, detrend
  3. HRV: calculate some metrics (like rmssd) from trended / detrended data

There are ready to use libraries like Neurokit to do all of these steps. But some methodological choices should be made. I mostly follow Kubios approaches, since they one of the current HRV software leaders. Their methodology explained here and in their software manual. First of all we have to get raw ECG signal in array and then pass it to Neurokit functions.

ECG cleaning is done by using Pan-Tompkins algo (1985). Find peaks done by using Neurokit method. Fix artifacts to get clean RR signal is done by Kubios method. Sample frequency is 512 for Movesense MD and 130 for Polar H10. 

  1. ecg_clean(ecg, sampling_rate=sample.freq, method="pantompkins1985")
  2. ecg_findpeaks(ecg.clean, sampling_rate=sample.freq, method = "neurokit)
  3. signal_fixpeaks(peaks, sampling_rate=sample.freq, method = "kubios", iterative=T, show=F)

If we only have RR intervals and dont have raw ECG, that's fine - just skip ecg_clean and ecg_findpeaks and pass RRi to signal_fixpeaks. H10 RR sample frequency should be set to 1000

Now we have cleaned RRi intervals and next step is to remove slow stationary trends (detrending). I prefer Smothness Priors with lambda = 500 as in Kubios. It's broken in Neurokit and i've found a fix for it:

  1. def signal_detrend_tarvainen2002(signal, regularization=500): 
  2.   N = len(signal)
  3.   identity = np.eye(N)
  4.   B = np.dot(np.ones((N, 1)), np.array([[1, -2, 1]]))
  5.   D_2 = scipy.sparse.dia_matrix((B.T, [0, 1, 2]), shape=(N - 2, N))  # pylint: disable=E1101
  6.   inv = np.linalg.inv(identity + regularization**2 * D_2.T @ D_2)
  7.   z_stat = ((identity - inv)) @ signal
  8.   trend = np.squeeze(np.asarray(signal - z_stat))
  9.   return signal - trend

This is a Python code, i call it from RStudio because for now i'm R guy :) In a future i will learn python and may move some part of my R code to Python. Reticulate package allows to import python packages and call them from R.

After detrending being done, we can calculate some HRV parameters. HR should be calculated from trended data and most of other metrics should be calculated from detrended data. Neurokit offers most of HRV metrics being calculated with single function: hrv(peaks, sampling_rate=sample.freq)

At that step we have most of HRV metrics for our measurement. We can loop through each timestamp and call hrv() function with part of signal (60 secs window) and get time-varying HRV. So for now we got most of what Kubios does and even more (time-varying feature is for paid version, which costs ~$380 per year).

Exercise

We can use Movesense MD or Polar H10 during exercise to record data and analyze it after. During exercise we will get more artifacts and might need more aggressive artifacts detection. I've made a simple thresholding function which works in similar way as Kubios thresholds. You can read my comments or read Kubios documentation.

  1. hrv.peaks.threshold <- function(peaks, sample.freq, filter.threshold = .25, hr.rr.slope = -.005, median.filter.window = 7, f.min.hr = 40, f.max.hr = 200) {
  2.   # adjust threshold - when hr increases threshold should be decreased and vise versa
  3.   # http://www.kubios.com/downloads/Kubios_HRV_Users_Guide_3_1_0.pdf 
  4.   # 5.2.1 Artefact correction methods
  5.   # That is, thresholds shown above are for 60 beats/min heart rate, but for higher heart rates the thresholds are smaller
  6.   # (because the variability is expected to decrease when HR increases) and vice versa for lower heart rates.
  7.   # glm <- lm(rr ~ hr,hrv.p); hr.rr.slope <- as.numeric(glm$coefficients["hr"])/1000; summary(glm)
  8.   rrs <- data.frame(peaks = peaks, rr = c(diff(c(0,peaks))))
  9.   rrs$rr.mf <- runmed(rrs$rr, median.filter.window) #calc rolling median RR
  10.   rrs$hr <- sample.freq*60/(rrs$rr) #calc hr
  11.   rrs$hr.mf <- sample.freq*60/(rrs$rr.mf) #calc rolling median HR
  12.   rrs$rr.diff.rr.mf <- abs(rrs$rr - rrs$rr.mf) #abs diff for each RR between rolling median
  13.   # adjusted threshold for bad rrs: K adj = (median RR + (60 - median HR) * slope)
  14.   # the lower threshold means more aggressive correction
  15.   # for example HR = 195: (60/195 + (195-60)*-.007) = -.15 => .25
  16.   # for example HR = 160: (60/160 + (160-60)*-.007) = -.325 => .25
  17.   # for example HR = 100: (60/100 + (100-60)*-.007) = 0.32
  18.   # for example HR = 40: (60/40 + (40-60)*-.007) = 1.64
  19.   # .25 is a medium threshold in kubios
  20.   rrs$mf.th <- (rrs$rr.mf + (rrs$hr.mf-sample.freq*60)*hr.rr.slope) * filter.threshold
  21.   # find HRs which bigger than specified
  22.   hr.bad.i <- as.numeric(rownames(rrs[(rrs$hr <= f.min.hr) | (rrs$hr >= f.max.hr),]))
  23.   # find RR which differ vs median more than threshold
  24.   rr.bad.i <- as.numeric(rownames(rrs[rrs$rr.diff.rr.mf >= rrs$mf.th,]))
  25.   bad.i <- unique(c(hr.bad.i, rr.bad.i))
  26.   # replace bad HR / RRs with NA and interpolate
  27.   rrs$rr.f <- rrs$rr; if(length(bad.i) > 0) rrs[c(bad.i),]$rr.f <- NA
  28.   rrs$rr.fs <- na.spline(rrs$rr.f, method = "natural")
  29.   rrs$peaks.clean <- rrs[1,"peaks"] - rrs[1,"rr.fs"] + cumsum(rrs$rr.fs)
  30.   return(list(peaks.clean = rrs$peaks.clean, bad.i = bad.i, hr.bad.i = hr.bad.i, rr.bad.i = rr.bad.i))
  31. }

So, for exercise additional step is added at the end of processing. hrv.peaks.threshold(peaks, sample.freq) will do the job and filter more exercise induced artifacts.

Breathing Rate

We can also use ECG/RR intervals to calculate breathing rate. It's called ECG Derived Respiration (EDR) and it base on Respiratory Sinus Arrythmia (RSA) which influence ECG Amplitude / RR peaks. Kubios have this detection in their software and i find EDR as a useful metric, especially during exercise. Neurokit have a function to detect EDR so it's easy to add it at the end of processing framework. I made 3 functions - first is getting respiration signal from RR peaks, second is finding peaks and lows from respiration signal, third calculates breathing rate = (n.peaks + n.lows) / 2 for a period of 60 secs:

  1. hrv.edr.rsp.from.peaks <- function(peaks, len, sample.freq, algo = 'vangent2019') {
  2.   ecg_rate = nk$ecg_rate(peaks, sampling_rate=sample.freq, desired_length=len, interpolation_method='monotone_cubic')
  3.   ecg_rates <- data.frame(i = 1:length(ecg_rate), rate = ecg_rate)
  4.   rsp <- nk$ecg_rsp(ecg_rate, sampling_rate=sample.freq, method=algo)
  5.   rsp_rates <- data.frame(i = (1:length(rsp)), r = rsp)
  6.   return(rsp_rates)
  7. }
  9. hrv.edr.rsp.find.peaks <- function(rsp_rates, sample.freq) {
  10.   dist <- 60/100*sample.freq #assume respiration is not fastest ins/exp 0.5s or 100br/min
  11.   rsp_rates$r <- detrend(rsp_rates$r)
  12.   rsp_rates$r <- as.numeric(scale(rsp_rates$r))
  13.   peaks <- pracma::findpeaks(rsp_rates$r, minpeakdistance = dist)
  14.   mins <- pracma::findpeaks(-1*rsp_rates$r, minpeakdistance = dist)
  15.   rrp <- rsp_rates[peaks[, 2], ]; rrm <- rsp_rates[mins[, 2], ]
  16.   return(list(rrp = rrp, rrm = rrm))
  17. }
  19. hrv.edr.respiration.from.rsp <- function(rrp, rrm, rrs, sample.freq) {
  20.   rrpm <- c(rrp,rrm)
  21.   rsa <- length(rrpm)/2
  22.   window <- rrs[which.max(rrs)] - rrs[which.min(rrs)]
  23.   return(round(rsa * (60/(window/sample.freq)),2))
  24. }

Finally we get a breathing rate:

  1. rsp_rates <- hrv.edr.rsp.from.peaks(ecg.peaks.i,length(ecg.peaks.i$ECG_R_Peaks),sample.freq,'sarkar2015')
  2. rsp_peaks <- hrv.edr.rsp.find.peaks(rsp_rates,sample.freq)
  3. edr.rrp <- rsp_peaks[[1]]; edr.rrm <- rsp_peaks[[2]]; edr.rrs <- rsp_rates
  4. respiratory.rate <- hrv.edr.respiration.from.rsp(edr.rrp$i,edr.rrm$i,edr.rrs$i, sample.freq)

Ok, now we have worked out how to cook raw ECG into HRV metrics salted with ECG Derived Respiration. We can do that even if we only have RR intervals without raw ECG.

I will improve post in a future by adding github repo with code and some plotting.