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