# How to decompose a signal into exponentally decaying sinusoids?

I have a system of coupled oscillators. One is a pendulum acting as a tuned mass damper and the other is a wooden structure. When I loosen the screw that allows the pendulum to oscillate I get a signal of the acceleration of the system like the following:

I know for sure that the amplitude of the waves composing this signal decay exponentially, however, I would like to know what their decay rate is. If I perform FFT on the signal I can only get the frequencies of the waves composing the signal. Is there a procedure similar to FFT I can implement in Python/Matlab with which I can find the frequencies, decay constants and initial amplitude of the waves composing this signal?

I found the following library in github (https://github.com/NanoComp/harminv), which would exactly do what I want, but it only runs on Linux systems (I have a Windows PC). Do you know of any programme that I can use to do this?

• This is exactly what the Laplace transform is classically for! – Envidia Feb 22 at 3:28

For a continuous-time signal this decomposition is precisely the Laplace Transform, and for a discrete time signal, the z-transform. Both decompose a time domain waveform into decaying and growing exponentials.

This is clear from observing the formulas and seeing the underlying operation similar to generalized correlation which is a sum of products:

Here is the general form for the Unilateral Laplace Transform:

$$\mathscr{L\{x(t)\}} = X(s) = \int_0^\infty x(t)e^{-st}dt$$

Where $$s= \sigma + j\omega$$

Note that $$e^{-st}$$ represents a rotating phasor that is growing or decaying with time!

If it isn't yet clear, $$e^{-st} = e^{-(\sigma + j \omega) t}= e^{-\sigma t}e^{-j\omega t}$$, and thus we have a rotating phasor ($$e^{-j\omega t}$$) that is growing or decaying $$e^{-\sigma t}$$ both as a function of time. If that is confusing, it may further help to know that the general form of the complex exponential $$Ke^{j\theta}$$ describes a vector (phasor) of magnitude K and angle $$\theta$$. Rotating phasors are unique frequency tones and fundamental elements (basis functions) of the Fourier Transform. Such a phasor rotating counter-clockwise on a complex plane is a positive frequency (the phase of the phasor is increasing linearly with time, and frequency is $$d\theta/dt$$), and one rotating clockwise would be a negative frequency. If we observed the real axis of any such rotating phasor, we would observe a sinusoidal waveform which is how we are typically first introduced to "frequencies". Real sinusoids consist of two such tones, a positive and negative frequency, as given by Euler's equation: $$2cos(\omega t) = e^{j\omega t} + e^{-j\omega t}$$. Thus complex tones of the form $$e^{j\omega t}$$ are even simpler subcomponents of real tones and once that is well understood, thinking of "frequencies" in terms of complex tones rather than real sinusoids greatly simplifies signal processing and leads to much greater possibilities! Complex need not be "so complex"!!

Let's call $$e^{st}$$ our hypothesis of what may be a component of $$x(t)$$ for any given value of $$s$$. We vary $$s$$ over all possible values to establish the result $$X(s)$$. In the Unilateral Laplace Transform, as given by the formula above, we multiply our signal $$x(t)$$ in time with $$e^{-st}$$ and integrate for all time, $$t\ge 0$$. Note the product is negated in the exponent: We spin exactly opposite to our hypothesis, and decay if it is growing and vice versa at the exact same rate as the mechanism to locate all these fundamental components! When a component is found, the rotation is stopped and the growth or decay is leveled out!

For any particular value (point) $$s_1$$ of $$e^{s_1 t}$$ that may be a fundamental component of our signal $$x(t)$$, meaning $$x(t)$$ contains as part of it one such fundamental decaying, rotating complex exponential; then for that component, it's product with $$e^{-s_1 t}$$ will be $$e^{s_1 t}e^{-s_1 t} = 1$$ for all time (for t>0)!... which then when integrated to infinity as given by the formula for the Unilateral Laplace Transform would go to infinity-- but notice any other value for $$s$$ where $$\sigma > \sigma _1$$ (for $$s_1 = \sigma _1 + j\omega_1$$) would necessarily decay and therefore converge to a finite value. Thus we see how we get distinct points that go to infinity (poles) and a vertical line as the boundary with everything to the right in the "Region of Convergence". Anything to the left of this vertical line (as given by the left most poles in the solution to the Laplace Transform) will not converge and therefore is an invalid region where no solution exists (and therefore our signal has no components there).

If plotted the Laplace Transform would be a surface of complex results for every combination on the 2 dimensional surface given by $$(\sigma, j\omega)$$. We often just depict this as poles (points going to infinity) and zeros (points equal to zero) since every other point on that surface is defined and unique by those pole and zero locations (but we may then miss understanding the full picture).

If we restrict $$s$$ to be only the $$j\omega$$ axis, then the Laplace Transform reduces to the Fourier Transform, which is a correlation with every possible rotating phasor that doesn't decay or grow (stays constant for all time). Specifically for the Unilateral Laplace Transform where $$t\ge 0$$, this would also involve multiplying the signal by a unit step ($$u(t)$$):

$$\mathscr{F\{u(t)x(t)\}} = \int_0^\infty x(t)e^{-j\omega t}dt$$

Notice how similar that equation is to the Unilateral Laplace Transform first given, and observe the surface plot of the magnitude of an example Laplace Transform below (the results of the Laplace Transform are themselves complex values, so each result on the surface has a magnitude and phase) for a low pass filter, and if we slice that 3 dimensional shape along the $$j\omega$$ axis we end up with the frequency response (Fourier Transform of the filter's impulse response, where $$X(s)$$ is the Laplace transform of the filter's impulse response), as given in the graphic below in the upper right hand corner (a low pass filter), and we see from the low pass filter that the impulse response consists of a decaying sinusoid (if the poles were on the $$j\omega$$ axis, this would be a constant level sinusoid, but given they are shifted to the left, that sinusoid is decaying with time.)

For discrete-time systems, the z-transform does the exact same thing, just with discrete phasors rotating in discrete time that are growing and decaying and instead of an s-plane we have a z-plane where the unit circle corresponds to the $$j\omega$$ axis. With discrete processing we use algorithms that compute the Discrete Fourier Transform (DFT) such as the ubiquitous FFT. The corollary for computing the Z-Transform is the chirp Z-transform (CZT), and for efficient implementations of this see Bluestein's algorithm.

You can use a time-frequency transform like CWT or STFT and focus it via reassignment for component extraction - example, synchrosqueezing (full post):

From there you run basic math on the relevant row(s) (plot on right) to find decay rate. Advantage over other methods is you can manually inspect the representation to separate noise from signal and focus on the non-noise. Library's ssqueezepy.