Estimating the Envelope of an Ultrasound Pulse
I have been doing some work with nondestructive evaluation, and one of the things I have been experimenting with is parameter estimation for an ultrasonic pulse model. A pulse might look something like this:
We can characterize this as a gaussian-damped sunusoid, where
\[\begin{aligned} d(t) =& A\cos\left(\theta + 2\pi f_c (t-\tau) + \beta (t-\tau)^2\right)e^{-\alpha(t-\tau)^2}\cr =& \cos\left(\theta + 2\pi f_c (t-\tau) + \beta (t-\tau)^2\right)e^{\ln(A) + \alpha t^2 -2\alpha\tau t + \alpha\tau^2}\end{aligned}.\]
The Gaussian can be estimated in a couple of ways. I started with an exponential curve fit of all maxima of the absolute value of the waveform (because the minimum peaks follow the same shape as the maximums). The curve fit can be formulated as a least-norm optimization problem,
\[\begin{aligned} & \underset{\bf{h}}{\text{minimize}} & & \|\ln(\bf{|d|}) - \bf{X}\bf{h}\| \end{aligned}\]
where \(\ln(|\bf{d}|)\) is a vector of the absolute natural log of the signal, \(\bf{X}\) is a 2nd order Vandermonde matrix of our sample indices, and \(\bf{h}\) is a vector of the coefficients we are estimating, assuming our envelope is of the form
\[f_{\text{env}}(t) = e^{h_0 + h_1t + h_2t^2}.\]
The coefficients are related to our original paramters in the following way:
\[\begin{aligned} h_2 =& \alpha \cr h_1 =& -2\alpha \tau \cr h_0 =& \ln(A) + \alpha\tau^2. \end{aligned}.\]
The least-norm approach, even using the L1 norm, doesn't do a fantastic job of fitting the exponential. The resulting curve is underfitted and misses the true envelope of the signal.
This underestimation is a side effect of regression on log-scaled values: less effort is expended minimizing errors for data points that are further away from zero because larger values become more and more compressed, while the values below 1 become more sensitive to errors. We can mitigate this problem by shrinking our data window so more emphasis is placed on the higher amplitude points near \(\tau\).
It's still not perfect because of the way we are estimating the envelope of the signal. If you notice, this fitted curve does exactly what we are asking it to: it passes through almost all of the peak values in the pulse, but it does not appear to rest on top of the waveform itself. We can do much better by tweaking the way we pose the problem: we will add a constraint that the resulting curve is on or above the data, which should give us a true envelope. It will also simplify preprocessing, getting us around the need to find the peaks of the waveform. The new constrained problem is
\[\begin{aligned} & \underset{\bf{h}}{\text{minimize}} & & \|\ln(\bf{|d|}) - \bf{X}\bf{h}\| \cr & \text{subject to} & & \left(\ln(\bf{|d|}) - \bf{X}\bf{h}\right) \leq \bf{0} \end{aligned}\]
The new estimate fits much better.
Here are the parameter estimates, to full precision, of both methods (with the original window size):
\[\begin{array}{cccc} & \text{Original} & \text{Unconstrained} & \text{Constrained} \cr A & 72 & 63.44982 & 72.06581 \cr \tau & 102 & 103.0854 & 101.973 \cr \alpha & -0.003 & -0.003100467 & -0.002988858. \end{array}\]
They are both close, except the unconstrained method is pretty far off on estimating gain, \(A\). Let's try this on some real data. This is a 5Mhz backsurface echo sampled at 100MHz.
Applying the fitting techniques does not yield great results.
The problem is that there is some dissimilarity between the proposed mathematical model and what the pulse is actually doing. The real pulse is asymmetric and doesn't fit into a true Gaussian shape, and beyond the main signal are small ripples present that distort the estimate.
We can first trim the pulse so we aren't looking outside of the main signal. We can also relax the constraint in such a way that our solution will ignore the outside points and focus more on the center of the signal. The original constraint says the fitted signal needs to lie on or above the original waveform; in other words, the error we are minimizing needs to be nonpositive. If we relax the bound on the error, or allow it to become positive as we move away from the center of the pulse, we can more accurately model the part of the waveform we care about and ignore the rest. Let us reformulate our problem and instead solve the following:
\[\begin{aligned} & \underset{\bf{h}}{\text{minimize}} & & \|\ln(\bf{|d|}) - \bf{X}\bf{h}\| \cr & \text{subject to} & & \left(\ln(\bf{|d|}) - \bf{X}\bf{h}\right) \leq \bf{F} \end{aligned}\]
where \(\bf{F}\) is a vector of constraints that relaxes as we move away from the peak of the function. For this problem I defined the relaxation as follows:
\[\bf{F} = c|\max(\bf{d}) - \bf{d}|^k\]
where \(c\) and \(k\) are tuneable constants. I set \(c = 0.05\) and \(k = .5\) and narrowed the window to eliminate the noise and focus on the most central portion of the pulse, allowing only 7 peaks near the center of the signal. Narrowing that far in gave very similar estimates for both the constrained an unconstrained methods.
Zooming out, we can see the final approximation of the real ultrasound pulse.
It's not perfect but it's close, maybe the best we could do without changing the model. This could have been entirely done by hand, tuning parameters to achieve the same results, but this type of method could also easily get a first estimate to begin working with. If nothing else, this has been a fun little exercise in optimization.
