--- title: "MRI physiological signal processing" output: rmarkdown::html_vignette vignette: > %\VignetteIndexEntry{MRI physiological signal processing} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r setup, include = FALSE} knitr::opts_chunk$set( collapse = TRUE, comment = "#>", fig.width = 7, fig.height = 3.5 ) ``` Philips MRI scanners record cardiac, respiratory, and motion signals alongside the BOLD acquisition via the Physiological Monitoring Unit (PMU). mrpheus can read these `.log` files, detect QRS complexes with the Pan-Tompkins algorithm, and compute instantaneous heart rate — producing cardiac regressors ready for RETROICOR or HRV analysis. ```{r libs, message = FALSE} library(mrpheus) library(ggplot2) library(gsignal) ``` ```{r theme, include = FALSE} theme_mri <- function() { theme_minimal(base_size = 13) + theme( panel.grid = element_blank(), axis.line = element_line(colour = "black"), axis.line.x.top = element_blank(), axis.line.y.right = element_blank() ) } ``` --- ## Reading a Philips PMU log `read_philips_physlog()` parses any Philips PMU `.log` file and returns a `mrpheus_physlog` object containing the signal matrix, event markers, and recording metadata. ```{r read} path <- system.file("extdata", "example_physlog.log", package = "mrpheus") rec <- read_philips_physlog(path) rec ``` The `system` argument sets the hardware preset: ```{r system, eval = FALSE} # Wireless VCG (Achieva/Ingenia dStream) — default rec <- read_philips_physlog(path, system = "wBTU") # Older wired ECG system (500 Hz) rec <- read_philips_physlog(path, system = "wired") # Site-specific sampling rate rec <- read_philips_physlog(path, system = "custom", sfreq = 400) ``` The resolved sampling frequency is always at `rec$HDR$sfreq` so nothing downstream needs to hardcode it. ### Signal channels ```{r channels} ecg <- as.double(rec$C[, "v1raw"]) ppu <- as.double(rec$C[, "ppu"]) resp <- as.double(rec$C[, "resp"]) ``` The wBTU layout provides ten channels: | Channel | Contents | |---|---| | `v1raw`, `v2raw` | Raw VCG electrode voltages | | `v1`, `v2` | Processed VCG signals | | `ppu` | Peripheral pulse unit (PPG, foot) | | `resp` | Respiratory belt | | `gx`, `gy`, `gz` | 3-axis accelerometer | | `mark` | Bit-encoded event marker | ```{r plot-signals, fig.height = 5, echo = FALSE} sfreq <- rec$HDR$sfreq t <- seq_along(ecg) / sfreq df_signals <- data.frame( t = rep(t, 3), value = c(ecg, ppu, resp), channel = rep(c("ECG (v1raw)", "PPU", "RESP"), each = length(t)) ) df_signals$channel <- factor(df_signals$channel, levels = c("ECG (v1raw)", "PPU", "RESP")) ggplot(df_signals, aes(t, value)) + geom_line(linewidth = 0.3, colour = "#014370") + facet_wrap(~channel, ncol = 1, scales = "free_y") + labs(x = "Time (s)", y = "Amplitude (a.u.)") + theme_mri() + theme(strip.text = element_text(face = "bold", size = 11)) ``` ### Event markers ```{r markers} rec$I$ScannerStart rec$I$ScannerStop head(rec$I$VcgOnset) ``` --- ## QRS detection `detect_qrs()` implements the Pan-Tompkins (1985) algorithm and returns a `mrpheus_qrs` object. ```{r qrs} qrs <- detect_qrs(ecg, fs = rec$HDR$sfreq) qrs ``` R-peak indices are in `qrs$qrs_i`. Account for `qrs$delay` if you need indices aligned to the raw (pre-filter) ECG. ```{r plot-peaks, echo = FALSE} bp <- gsignal::butter(3L, c(5, 15) * 2 / sfreq, type = "pass") ecg_bp <- as.double(gsignal::filtfilt(bp, ecg)) df_bp <- data.frame(t = t, ecg = ecg_bp) df_rpeaks <- data.frame(t = qrs$qrs_i / sfreq, amp = qrs$qrs_amp) ggplot(df_bp, aes(t, ecg)) + geom_line(linewidth = 0.3, colour = "#1B6799") + geom_point(data = df_rpeaks, aes(x = t, y = amp), colour = "#FC544A", size = 2, shape = 21, fill = "#FC544A", alpha = 0.85) + labs(x = "Time (s)", y = "Amplitude (normalised)", title = "Bandpass ECG with detected R-peaks") + theme_mri() ``` ```{r rr} rr_ms <- diff(qrs$qrs_i) / rec$HDR$sfreq * 1000 cat("Mean RR :", round(mean(rr_ms)), "ms\n") cat("SDNN :", round(sd(rr_ms)), "ms\n") ``` ### PPU as an alternative In neonates, PPU often gives a cleaner signal than the ECG electrodes. `detect_qrs()` accepts any pulsatile channel: ```{r ppu, eval = FALSE} qrs_ppu <- detect_qrs(as.double(rec$C[, "ppu"]), fs = rec$HDR$sfreq) ``` --- ## Heart rate signal `compute_hr_signal()` converts R-peak indices into a sample-by-sample instantaneous HR trace (bpm). It reads `fs` from the `mrpheus_qrs` object directly: ```{r hr} hr <- compute_hr_signal(qrs) ``` ```{r plot-hr, echo = FALSE} df_hr <- data.frame(t = seq_along(hr) / sfreq, hr = hr) ggplot(subset(df_hr, hr > 0), aes(t, hr)) + geom_line(linewidth = 0.4, colour = "#FFA75D") + labs(x = "Time (s)", y = "Heart rate (bpm)", title = "Instantaneous heart rate") + theme_mri() ``` --- ## Epoch extraction To extract the physiology aligned to the fMRI scan, use the scanner-stop marker. This is more reliable than the start marker because the exact recording onset varies but the stop is always clean: ```{r epoch} sfreq <- rec$HDR$sfreq enddelay <- 42L tr <- 0.392 vols <- 47L stopidx <- tail(rec$I$ScannerStop, 1L) - enddelay startidx <- max(1L, stopidx - round(vols * tr * sfreq)) epoch <- rec$C[startidx:stopidx, ] cat("Epoch:", nrow(epoch), "samples /", round(nrow(epoch) / sfreq, 1), "s\n") ``` ### RR intervals per TR ```{r rr-tr} qrs_epoch <- detect_qrs(as.double(epoch[, "v1raw"]), fs = sfreq) rr_t <- qrs_epoch$qrs_i[-1] / sfreq rr_s <- diff(qrs_epoch$qrs_i) / sfreq tr_times <- seq(0, vols * tr, by = tr) rr_tr <- approx(rr_t, rr_s, xout = tr_times, method = "constant")$y cat("RR per TR (first 5):", round(rr_tr[1:5] * 1000), "ms\n") ``` `rr_tr` is a length-`vols` vector ready for RETROICOR phase computation or inclusion as a nuisance regressor in a GLM. --- ## Full pipeline ```{r pipeline, eval = FALSE} library(mrpheus) rec <- read_philips_physlog("sub-01_ses-01_physlog.log") sfreq <- rec$HDR$sfreq stopidx <- tail(rec$I$ScannerStop, 1L) - 42L startidx <- stopidx - round(2300L * 0.392 * sfreq) epoch <- rec$C[startidx:stopidx, ] qrs <- detect_qrs(as.double(epoch[, "v1raw"]), fs = sfreq) hr <- compute_hr_signal(qrs) rr_tr <- approx( x = qrs$qrs_i[-1] / sfreq, y = diff(qrs$qrs_i) / sfreq, xout = seq(0, 2300 * 0.392, by = 0.392), method = "constant" )$y ``` --- ## References Pan J, Tompkins WJ. A real-time QRS detection algorithm. *IEEE Trans Biomed Eng.* 1985;32(3):230–236. Groot PFC. *ReadPhilipsScanPhysLog* (MATLAB). AMC Amsterdam. 2010. Glover GH, Li TQ, Ress D. Image-based method for retrospective correction of physiological motion effects in fMRI: RETROICOR. *Magn Reson Med.* 2000;44(1):162–167.