# How to remove heartbeat signal from blood pressure signal?

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)$$?

• @rodrigo-de-azevedo , as mentioned we have recorded h(t) for a period of time and have some recorded signal like B'(t) in the same period, which has B(t) and h'(t). h'(t) has lower amplitudes than h(t) (as the sensor captures it weakly). However, it is not something like multiplying h(t) by some constant. Usually, h(t) has higher frequencies than B(t) so a low pass filter may be used. However, as h'(t) may change from person to person and B(t) may have high frequency components, I am looking for a better solution. for instance some sort of decomposing B'(t) and finding components related to.. Sep 12, 2021 at 14:28
• h'(t) is unknown and we wish to find it, If you mean B'(t), I dont know if cross correlation would help for finding delay. besides, while h' is much similar to h at each point t we have something like h'(t)=ah(t-v) for some noisy a. so finding only the delay would not solve the problem Sep 12, 2021 at 14:48
• Are you trying to record the blood pressure averaged over time, or are you trying to record systolic and diastolic blood pressure? If the latter, I'm not sure how you're going to separate these out, because the instantaneous blood pressure varies in sync with the heartbeat. Sep 12, 2021 at 14:55
• Just as a note -- you may want to review your variable names. In signal processing, $h(t)$ is almost universally used to denote a system impulse response, so unless $h(t)$ is already used to denote variations due to heartbeat in biomechanical circles, you'll probably confuse more than you clarify by that choice of signal name. I'd see if I can find what people tend to use in papers; if that doesn't work I'm not sure what to suggest -- $x_h(t)$, or $p_h(t)$, maybe. Sep 12, 2021 at 14:57
• @TimWescott , I am not only recording the blood pressure, It was just one sample. There are other types of signals which are mixed with heart beat. The choice of h(t) was poor and I am changing it. Sep 12, 2021 at 15:53

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.

• h(t) is not the impulse response it was just a poor naming, I changed it to A(t). thanks for the program I am working on it. Sep 12, 2021 at 16:01

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)$$.