How to recover frequency and damping coefficient from superposition of damped oscillators? 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.
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: