Finding exponential decay in noisy vibration signal

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.

• Instead of smoothing the data, I would try using only the adjacent peaks, or the maximum values found within some narrow sliding window. Commented Jan 24, 2020 at 21:04
• I thought about the peak idea too but couldn't think of a simple algorithm to achieve it. Do you have a suggestion on that? The maximum values approach is of course also a good idea I will definitely try that, thank you! Do you have any ideas on the core issue of my question? Commented Feb 7, 2020 at 9:15
• Did you ever find a standard approach lineInk? I have a similar problem Commented Jul 19, 2022 at 3:17
• @kevinkayaks No, I ended up using an entirely different approach. Given my current state of knowledge, I would probably try to use a Bayesian model instead of a non-linear least squares regression model. Such a model should contain a dampened sinusoidal that starts with a certain offset t_0 with its frequency and the dampening coefficient also being free parameters. Then you could also add a constant noise term. Given the data shown, maybe adding additional components might be necessary. Commented Jul 20, 2022 at 11:14
• Thanks @lineInk. I will post a tentative solution from my recent research on the problem Commented Jul 20, 2022 at 17:22

1 Answer

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