How to remove heartbeat signal from blood pressure signal?

Question

I am working on a medical engineering project that involves processing signals received from a human body. We use a sensor to record the blood pressure, $B(t)$.

However, if you put your finger on many parts of your skin, you can notice your heartbeat. Accordingly, the recorded signal also has some weak heartbeat signal added, $A'(t)$. To obtain accurate blood pressure signal, we have a sensor that records the heartbeat signal, $A(t)$. However, $A'(t)$ is a delayed and weak version of $A(t)$, where the phase lag is unknown.

Knowing $A(t)$, how can we find accurate $B(t)$?

Answer 1

If $A(t)$ is known, it can be zeroed in the synchrosqueezed representation - the remainder is then $B(t)$, recovered by inversion. $A(t)$ need not be known perfectly - just enough to indentify its time-frequency ridges.

Answer 2

Problem

Your model is: $B'(t) = B(t) + A'(t)$ where $B'(t)$ is your observed blood pressure corrupted with the heartbeat signal, $A'(t)$. You want to find $B(t)$ but you don't know $A'(t)$, instead you know $A(t) = \gamma A'(t-d)$.

Option 1

A similar problem led Widrow to invent the Least-mean-squares (LMS) filter.

We want to find a time-varying filter $w(t)$ that when multiplied with the heartbeat signal $A(t)$ approximates the blood pressure signal $B(t)$, i.e., $B(t) = w(t)A(t)$. To do so, we minimize the mean squared error between the measured signal and the model wrt $w$.

The error signal is $\epsilon(t) = B'(t) - B(t) = B'(t) - w(t)A(t)$

The LMS cost function is,

$J(w) = \mathbb{E}(\epsilon^2(t)) \\ = J(w) = \mathbb{E} [(B'(t) - w(t) A(t))^2] \\ = J(w) \mathbb{E} [B'^2(t)] + w(t)^2 \mathbb{E}[A^2(t)] - 2w(t) \mathbb{E} [B(t)A(t)] \\ \frac{\partial J}{\partial w(t)} = 2w(t) \mathbb{E}[A^2(t)] - 2 \mathbb{E}[B'(t)A(t)] $

Now, Widrow suggested an adaptive update to the filter weights to account for changes in the signal over time:

\begin{align} w(t) &= w(t-1) + 0.5 \mu \frac{\partial J}{\partial w(t)} \\ &= w(t-1) + \mu \left(w(t)A^2(t) - B'(t)A(t) \right) \\ &= w(t-1) + \mu \left[A(t) (w(t)A(t) - B'(t))\right] \\ &= w(t-1) + \mu A(t) \epsilon(t). \end{align}

Your best estimate of the blood pressure is $w(t)A(t)$.

Option 2

Let's take the Fourier transform of $B'(t)$ and $A(t)$ and look at them in the frequency domain,

$B'(\omega) = B(\omega) + A'(\omega) \\ A(\omega) = \gamma \exp(-j\omega d) A'(\omega)$

If we look at the cross-power spectrum of the signals, we get

$B'^*(\omega)A(\omega) = \left(B^*(\omega)+A'^*(\omega)\right) \gamma \exp(-j\omega d)A'(\omega) $

We can make the assumption that the blood pressure and the heartbeat signal are uncorrelated (I have a feeling a biologist would not be happy with this). Then, $B^*(\omega)A'(\omega) = 0$ and \begin{align} B'^*(\omega)A(\omega) &= \gamma \exp(-j\omega d) |A'(\omega)|^2 \\ &= \gamma^3 \exp(-j\omega d) |A(\omega)|^2 \\ \frac{B'^*(\omega)A(\omega) } {|A(\omega)|^2} & = \gamma^3 \exp(-j\omega d) \end{align} The inverse Fourier transform of this is some sort of a modified cross-correlation function,

$\hat{R}_{B',A}(t) = \frac{1}{2\pi} \int_{-\pi}^{\pi} \frac{B'^*(\omega)A(\omega)} {|A(\omega)|^2} e^{j\omega t} d\omega = \gamma^3 \delta(t-d).$

From this function, you can find $\gamma$ and $d$ and thus, $A'(t)$ and simply subtract it from $B'(t)$ to get $B(t)$.