How to recover frequency and damping coefficient from superposition of damped oscillators?

Question

I have a signal which is (more or less) the superposition of multiple damped oscillators. There are relatively few oscillators involved, 2-3 would account for most of the signal. They have periods which are not multiples of each other (ie not harmonics), and they all have different damping coefficients / time constants. They also are not necessarily in phase initially.

I would like to recover the periods, damping coefficients and phases of the oscillators from the signal.

A synthetic example of a signal like this can be produced by the following python code:

from pylab import *
x = linspace(0,1000,num=1000)
y1 = 1 * sin(2*pi*x/11 + 0) * exp(-x/300)
y2 = 0.5 * sin(2*pi*x/17 + 1.2) * exp(-x/300)
y3 = 2 * sin(2*pi*x/50 - 0.25) * exp(-x/50)
y = y1 + y2 + y3
figure(figsize=(20,10))
plot(y)

This looks pretty similar to the real signal qualitatively, including having some beat frequencies, and a pretty large initial low-frequency component.

This seems like a pretty common problem when analyzing a physical system with multiple coupled oscillators, eg a vibrating object with multiple resonances.

UPDATE Using direct curve-fitting works pretty well, even for strongly damped signals - depending on the choice of initial conditions, that is:

from scipy.optimize import curve_fit

def func(x, a, b, c, d):
     return a * sin(2*pi*x/b - c) * exp(-x/d)

p, pcov = curve_fit(func, x, y, p0=(1,100,0,1000))
print(p)

# [  1.78453351  51.53572084   0.05566558  55.10643076]

This recovers the largest-amplitude (but very strongly damped: time constant = 1 period) part of the signal.

Answer 1

If you are absolutely sure that you have a signal that is the superposition of damped oscillators then you can simply track their evolution in time.

Fs = 2000; % Sampling frequency (Hz)
T = 4;     % Time duration (seconds) 
t = 0:(1./Fs):(T-(1./Fs)); % Time vector (in seconds)
p = 2.0.*pi.*t; % Phase vector
S = 1.0*sin(440.*p).*exp(-t) + 0.25 .* sin(800.*p).*exp(-t*10); % Signal

This signal looks like:

And via a Short Time Fourier Transform (STFT) (surf(stft(S))) it looks like this:

Granted, if you have beating frequencies, the task is getting a bit more challenging. For instance:

S = 1.0*sin(440.*p).*exp(-t) + 0.25 .* sin(442.*p).*exp(-t*2) + sin(800.*p).*exp(-t*10);

Would show up as:

But once you see oscillations on the damped decay, you know you have beating and you can recover how close the beating pair is by measuring the beating frequency and how quickly it fades (to derive its decay).

In a similar way (i.e. using the STFT) you can also recover the phases of the sinusoids.

Hope this helps.

UPDATE:

The STFT of the signal that was shared via pastebin is capable of distinguishing three components in this signal. This is as follows:

%Given a signal in s and that the signal package is installed (in Octave)
[S,F,T] = specgram(s);
surf(abs(S));

This returns:

Answer 2

There is an actual group of methods that deal exactly with this: the so-called High Resolution Methods. They are based in covariance of the signal and are very efficient separing noise from sinusoidal signal recovering frequencya and damping factor (later, you can recover amplitude and phase with a least mean squares algorithm).

In particular, I recommend to you the Subspace Methods. The theory is pretty harsh, deeply based in Linear Algebra, but the results are spectacular. Two algorithms are particularly efficient: MUSIC and ESPRIT. In this link you can download MATLAB routines. You can search for Python implementations (I have googled Python ESPRIT MUSIC and found this repository).

I really encourage you to use these tools; the problem you can encounter when using STFT (as suggested in another answer) is that STFT cannot separate properly near frequencies; you would be forced to use big windows and you can face problems of modulations of amplitudes (which fake the damping factor).

I learned about these methods from Roland Badeau who did his thesis about this topic in the case of Music, but there are a lot of works about this methods.

Answer 3

Google

kumaresan-tufts method.

which is based on Prony’s method. Matlab also has prony

Some people also say Tufts-Kumaresan.

This is a direct solution while the STFT suggestion is probably better is the noise is very colored