Plotting the decay envelope for dampening function without the usual “fit with a decaying sinusoidal first” method?

To be more clear, by the decay rate of a dampening function I mean fitting the exponential curve to the amplitude of the damped oscillations and see how it decays.

Something like the following :

Damped Sine wave

but for my graph which is certainly not smooth like the above and whose picture I still cannot upload due to my low reputation ):

Now I can think of two ways to do it to:

1. First fit a damped sinusoidal to my data points and then find the decay rate using this damped sinusoidal function.

2. Break the function into a finite number of intervals and find the maximum value of the function in that interval. Then, treat it as the amplitude and fit these with an exponential.

I though think both the above implementations will be very imprecise. The damped sinusoidal won’t be a good fit for my function as even though ti would probably capture the decay in my function it cannot keep up with it’s non-uniform oscillations.
For the second case I am guessing the accuracy depends on the length I choose to discretize the function with and so it’s not clear what should be the appropriate length to choose.

Is there a method to do this that is better than the above two or maybe that uses the in-built Python functions?

Edit:

Here is my plot:

Note: The erratic fast local oscillations are what made me think why the damped sine method might be less accurate for this (even though it seems that it might actually be more accurate as shown in K.CI’s answer)

Answer

I’ve made some toy functions that could help you decide on your method. I’ve found the exponential factors using the method I described (only peaks, found naively) and also by fitting the whole curve. If you could provide an example data then we might approach a problem directly.

import matplotlib.pyplot as plt
import numpy as np
from scipy.optimize import curve_fit

def exp(t, A, lbda):
    r"""y(t) = A cdot exp(-lambda t)"""
    return A * np.exp(-lbda * t)

def sine(t, omega, phi):
    r"""y(t) = sin(omega cdot t + phi)"""
    return np.sin(omega * t + phi)

def damped_sine(t, A, lbda, omega, phi):
    r"""y(t) = A cdot exp(-lambda t) cdot left( sin left( omega t + phi ) right)"""
    return exp(t, A, lbda) * sine(t, omega, phi)

def add_noise(y, loc=0, scale=0.01):
    noise = np.random.normal(loc=loc, scale=scale, size=y.shape)
    return y + noise

You can increase the noise with the scale value in add_noise. Let’s get some data:

# Real vals
t = np.linspace(0, 5, num=300)
A = 1
lbda = 1
omega = 8 * np.pi
phi = 0

real_y = damped_sine(t, A, lbda, omega, phi)
np.random.seed(42)
noisy_y = add_noise(real_y, scale=0.05)

The plot:

Naive method to find peaks:

def find_peaks(x, y):
    peak_x = []
    peak_vals = []
    for i in range(len(y)):
        if i == 0:
            continue
        if i == len(y) - 1:
            continue
        if (y[i-1] < y[i]) and (y[i+1] < y[i]):
            peak_x.append(x[i])
            peak_vals.append(y[i])
    return np.array(peak_x), np.array(peak_vals)

Yields this:

Fitting only the peaks with the exponential:

popt, pcov = curve_fit(exp, noisy_peak_x, noisy_peak_y)
print(*[f"{val:.2f}+/-{err:.2f}" for val, err in zip(popt, np.sqrt(np.diag(pcov)))])
print(A, lbda)

Yields: 1.06+/-0.09; 1.39+/-0.13. Pretty bad, but the noise is quite big and the peaks are badly found. Here’s the plot:

If you fit the entire sinusoidal curve, the result is better:

popt, pcov = curve_fit(damped_sine, t, noisy_y)
print(*[f"{val:.2f}+/-{err:.2f}" for val, err in zip(popt, np.sqrt(np.diag(pcov)))])
print(A, lbda, omega, phi)

Yields this 1.02+/-0.02 1.03+/-0.03 25.16+/-0.03 -0.02+/-0.02 compared to 1 1 25.13 0, which is quite better, after all, there’s many more points. Here’s the graph:

If you expand the peak finding method to consider a window:

def find_peaks_window(x, y, window=1):
    peak_x = []
    peak_vals = []
    for i in range(len(y)):
        if i < window:
            continue
        if i > len(y) - window:
            continue
        left_side = y[i-window:i]
        right_side = y[i+1:i+1+window]
        if all(left_side < y[i]) and all(right_side < y[i]):
            peak_x.append(x[i])
            peak_vals.append(y[i])
    return np.array(peak_x), np.array(peak_vals)

You get this 0.93+/-0.04 0.78+/-0.05 compared to 1 1: