3
$\begingroup$

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:

Absolute of vibrational signal with over-plotted RMS envelope and smoothened envelope.

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:

enter image description here

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.

$\endgroup$
5
  • 1
    $\begingroup$ Instead of smoothing the data, I would try using only the adjacent peaks, or the maximum values found within some narrow sliding window. $\endgroup$ Jan 24, 2020 at 21:04
  • $\begingroup$ 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? $\endgroup$
    – lineInk
    Feb 7, 2020 at 9:15
  • $\begingroup$ Did you ever find a standard approach lineInk? I have a similar problem $\endgroup$ Jul 19, 2022 at 3:17
  • $\begingroup$ @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. $\endgroup$
    – lineInk
    Jul 20, 2022 at 11:14
  • $\begingroup$ Thanks @lineInk. I will post a tentative solution from my recent research on the problem $\endgroup$ Jul 20, 2022 at 17:22

1 Answer 1

1
$\begingroup$

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)

enter image description here

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: enter image description 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

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.