Finding exponential decay in noisy vibration signal

Question

I have to analyse vibrational signals for which the general assumption is that there is one dominant excitation and an exponential decay in amplitude thereafter.

I have created smoothened envelopes from the original signal like is visualised in the following plot:

Now I want to fit the exponential decay in the signal. Before being able to perform a standard fit I have to select the correct subsection of the signal which contains the exponential decay.

I have a kind of working solution where I simply perform a regression on slices of the signal selected by a moving window. Then I select the longest section of the signal where the normalized ${\chi}^2$-value is below a manually adapted threshold.

This leads to a result like in the following plot:

One could then use this section (extending it maybe by the moving average window length) and perform another exponential fit yielding the final values.

I dislike my solution because it seems quite wasteful to perform so many fits before performing another final one. And I think this has to be such a common problem in signal processing that there must be a more elegant and less brute-force solution to this problem.

Another thing I thought about is simply performing a convolution but I expect the exponentials to differ quite widely between individual signals. The only idea for an improvement of the convolution approach I had, was to use multiple different exponential decays for convolution.

I am searching for a standard and well-tested approach for this kind of problem. Ideally also one that is efficient or even a non-iterative one-step solution but I am not sure if that is even mathematically possible.

Answer 1

Consider a noisy exponential decay signal like the following test data, generated in python:

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import simps

dt = 1e-3 # timestep
N = 15000 # number of step
t = np.arange(0,N*dt,dt)


tau = 5 # this is what I want to estimate
np.random.seed(19) # set the random state
noise = [] # list to be filled with noise
n = 14.4 # initial noise value
D = 13.4 # diffusivity
gam = 0.3# damping
n0 = 31
i=0 
while i<N:  # produce an Ornstein uhlenbeck noise (correlated)
    noise.append(n)
    n = n + gam*(n0 - n)*dt + np.sqrt(2*D*dt)*np.random.normal()
    i+=1

q = noise*np.exp(-t/tau)
plt.plot(t,q)
plt.xlabel('t',fontsize=15)
plt.ylabel('f(t)',fontsize=15)

Using the technique of "Schroeder Integration", rather than fitting this noisy curve directly, we can fit the "backward integrated" signal $$ L(t) = \frac{\int_t^T f(u)^2 du}{\int_0^T f(u)^2 du}. $$ This signal produces a line in the $t$ -- $\log L(t)$ plane, provided the original signal $f(t)$ is a noisy exponential decay. The noise gets pushed to the tail. The slope of the line is $2\tau$, where $\tau$ is the original decay rate of the exponential.

For the signal above, schroder integral $L(t)$ appears as the blue curve here:

For the test data above, the predicted decay constant differs by 2% from the one used to generate the test data:

L = np.array([simps(q[i:]**2,t[i:]) for i in range(len(t))])/simps(q**2,t)
L = np.log(L)

mask = ~(np.isnan(L)|np.isinf(L))&(t<12)
a,b = np.polyfit(t[mask],L[mask],1)

plt.semilogy(t,np.exp(L))
plt.semilogy(t,np.exp(b+a*t))
plt.xlabel('t',fontsize=15)
plt.ylabel('L(t)',fontsize=15)


T = -1/a*2
print(round(np.abs(T-tau)/tau*100,2), 'percent error.')
# 1.98 percent error.

This is a nice reference that briefly explains why this technique works and provides more detailed citations https://asa.scitation.org/doi/10.1121/1.3587944