- Created on Sunday, 09 June 2019 06:06

As part of my research I deal with blood pressure (BP) waveforms. Often it is useful to ensemble average a continous trace of BP waveforms (spanning multiple cardiac cycles) into a single ensemble average beat (spanning one cardiac cycle). Currently we do this using Matlab, but I have been working on moving this functionality into R in my spare time. This post will detail where I am at.

First, we load the blood pressure trace. This is a single file (txt or csv) with a single column which contains the data.

`BP = read.table("C:/Users/mka4/Desktop/pressure.txt", dec=".") #load text file`

To view the data we could simply plot it out. Always good to check!

```
p = BP$V1
plot(p, type ="l",
ylab = "mmHg",
xlab = "t(ms)",
main = "BP trace")
```

Next comes the trickier part. We need to detect individual cardiac cycles and average the values for each beat to give a single ensembled averaged beat. Below is my function for doing this. Cardiac cycles are detected using the nadir (diastole) of the waveform trace.

```
ensemble = function(P){
P = P$V1 # loading data
p = P # loading data
ns = 1000 # sample rate
nd = 250 # lag factor to pic beats
# Find troughs -----------------------------------------------------------
a = min(p[1:ns]) # find 1st trough
nmin = which(p[1:ns]==a) # time @ first trough
nmin = nmin[1]
j = 1 # nmax= index corresponding to maximum peak
nb = NULL
nb[j] = nmin
# find all troughs
nend = nb[j] + nd + ns
suppressWarnings ( while (nb[j] + nd + ns < length(p)) {
nstart = nb[j] + nd
nend = nstart + ns
a = min(p[nstart:nend])
nmin = which(p[nstart:nend] == a)
j = j + 1
nb[j] = nmin + nstart - 1} )
# Find peaks ---------------------------------------------------------------
v = max(p[1:ns]) # find 1st peak
nmax = which(p[1:ns]==v) # time @ first peak
nmax = nmax[1]
k = 1
nv = NULL
nv[k] = nmax # nmax= index corresponding to minimum trough
# find all peaks
nend2 = nv[k] + nd + ns
suppressWarnings ( while (nv[k] + nd + ns < length(p)) {
nstart2 = nv[k] + nd
nend2 = nstart2 + ns
v = max(p[nstart2:nend2])
nmax = which(p[nstart2:nend2] == v)
k = k + 1
nv[k] = nmax + nstart2 - 1} )
# Ensemble ----------------------------------------------------------------
Nb = length(nb)
nperiod = mean(diff(nb)) # find the average cardiac period
nint = round(1.05 * nperiod) # choose a beat interval slightly longer than nperiod
tens = (0:(nint - 1)) / ns # time variable of the ensemble
pens = array(0, c(Nb - 1, nint))
for (n in 1:(Nb - 1)) {
pens[n, ] = P[nb[n]:(nb[n] + nint - 1)]}
pensavg = apply(pens, 2, mean, na.rm = TRUE) # ensemble presure (2=applys over columns)
pensSD = apply(pens, 2, sd, na.rm = TRUE) # SD for waveform
pensSD = round(mean(pensSD), digits = 2) # mean SD
# Additional calcs --------------------------------------------------------
sysP = round(mean(p[nv]), digits = 2) # Calculate average SBP
diaP = round(mean(p[nb]), digits = 2) # Calculate average DBP
sysSD = round(sd(p[nv]), digits = 2) # Calculate SD for SBP
diaSD = round(sd(p[nb]), digits = 2) # Calculate SD for DBP
# MAP ---------------------------------------------------------------------
# integrate using the trapezoidal rule
x = pensavg/(length(pensavg)-1)
y = x
x = seq(along=x)
m <- length(x)
xp <- c(x, x[m:1])
yp <- c(numeric(m), y[m:1])
n <- 2*m
p1 <- sum(xp[1:(n-1)]*yp[2:n]) + xp[n]*yp[1]
p2 <- sum(xp[2:n]*yp[1:(n-1)]) + xp[1]*yp[n]
MAP = round((0.5*(p1-p2)), digits = 2)
#aritmetic mean map
mapCum = round(mean(pensavg), digits = 2)
#as per Alun's matlab (requires pracma)
#mapAH = round(pracma::trapz(pensavg/(length(pensavg)-1)), digits = 2)
#print(c("MapAH: ", mapAH), quote = F) # Mean arterial AH
# Heart rate --------------------------------------------------------------
#calculate heart rate from sbp peaks
HR = round(
((length(nv)-1)*60) / ((max(nv)-min(nv))/1000),
digits = 2)
# Plots -------------------------------------------------------------------
layout(matrix(c(1,1,2,3), 2, 2, byrow = TRUE))
#plot peaks
plot(P[1:length(P)],
type = "l",
xlab = "t (ms)",
ylab = "mmHg",
main = "Trace",
las = 1)
points(nb,p[nb],col=4, cex = 1.2, pch = 8)
points(nv,p[nv],col=2, cex = 1.2, pch = 8)
#plot ensemble
plot(tens,
pensavg,
type = "l",
ylim = c(min(pensavg)-5, max(pensavg)+10),
ylab = "mmHg",
xlab = "t (s)",
main = "Ensemble",
las = 1)
grid()
for (n in 1:(Nb - 1)) {
lines(tens, pens[n, ], col = n)}
lines(tens, pensavg, col = 1, lwd = 4)
plot(tens,
pensavg,
type = "l",
xlab = "t (s)",
ylab = "mmHg",
main = "Map",
las = 1)
abline(h = MAP, col = "red")
# Print results ---------------------------------------------------------
print(c("Sbp: ", sysP), quote = F) # Systolic BP ave
print(c("Sbp SD:", sysSD), quote = F) # Systoli SD
print(c("Dbp: ", diaP), quote = F) # Diastolic BP ave
print(c("Dbp SD:", diaSD), quote = F) # Diastolic SD
print(c("Map: ", MAP), quote = F) # integrated MAP
print(c("MapC: ", mapCum), quote = F) # arithmatic MAP
print(c("Hr: ", HR), quote = F) # Heart rate
print(c("SD: ", pensSD), quote = F) # Heart rate
return(pensavg) # Ensembled waveform
}
```

This function does all the work for you. Now you only need apply it to the BP trace that was loaded in earlier. The function will output some figures to the plot tab and values to the console.

`ensemble_waveform <- ensemble(BP)`

```
## [1] Sbp: 144.92
## [1] Sbp SD: 2.66
## [1] Dbp: 79.09
## [1] Dbp SD: 1.22
## [1] Map: 105.55
## [1] MapC: 105.53
## [1] Hr: 62.14
## [1] SD: 1.97
```

The top plot is the full BP trace with red and blue points intersecting where the code has identified the peak (systolic) and nadir (diastolic) on the waveform. The bottom left figure is all the individual beats plotted on top of each other with the ensemble averaged waveform plotted in black. The bottom right figure is the ensemble averaged waveform with the red line intersecting at the mean arterial pressure (integrated MAP).

**Console output:**

**Sbp** = Systolic BP

**Sbp SD** = Standard deviation of systolic BP

**Dbp** = Diastolic BP

**Dbp SD** = Standard deviation of diastolic BP

**Map** = Mean arterial pressure (via integration of the ensemble averaged waveform)

**MapC** = Mean arterial pressure (via arithemtic mean)

**Hr** = Heart rate

**SD** = Standard deviation of the ensemble averaged waveform

**Notes:**

- Any sample rate can be used, just be sure to change the value “ns” within the function to match the sampling frequency.

- The function may not pick the right peaks and nadirs, small tweaks to the value “nd” within the function should fix this.

**Future additions:**

- Calculate MAP for each individual beat.

- Calculate augmentation index using the zero crossing point of the forth derivative (as per Kelly et al., Circ 1989). Ideally this will be calculated for each individual beat and plotted and printed in the output.

- Ensemble average the waveform using the maximum slope of systole. This is simply calculated as the peak of the first derivative.

- Tweak the function so that it can be easily applied over multiple files using lapply() from the dplyr package. This includes saving all outputs to a single spread sheet.

For now it remains a work in progress, but it works and is free to use if you might find it useful. I would ask if you make improvments, additions or find any errors that you share them here so that others may benifit, you will be credited with any additions made.

Copyright M.K.Armstrong 2019