# 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. – A_A
Jun 11, 2018 at 5:49
• @A_A Sorry for the slow reply. It is not completely resolved - the method I posted above works for recovering the lowest frequency/largest amplitude component if a good starting point is picked by hand. I don't have a reliable way of getting the other components. I also couldn't get STFT to work at all for these types of signals. Dec 18, 2019 at 23:11
• @A_A Let me know if you have any other ideas - I think this is a pretty useful problem to solve in general. Dec 18, 2019 at 23:11
• Can you please post a representative signal that you are working with?
– A_A
Dec 19, 2019 at 11:36

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: • Hi - Thanks for the answer! I tried STFT, but I didn't get really good results so far. I could give it another try. I think one difference is how fast the damping is relative to the period of the sine waves. In my case the damping time constants are 1 to 30 periods or so, while in the synthetic example you showed it is 80 and 440 periods. It seems like STFT works best when damping has a time constant of tens of periods at least. May 9, 2018 at 6:36
• @AlexI Can you please post the signals you are dealing with? 1 period is too short anyway. 30 periods is recoverable. There are a few more options in this space. Are these simulation outputs or real-world signals?
– A_A
May 9, 2018 at 6:52
• I uploaded an example of a synthetic signal (very similar to real one) here: pastebin.com/jUFYXyrA Jun 2, 2018 at 5:23
• Please have a look at the updated question, I added some python code. Even 1 period works okay using curve fitting (but depends on picking good starting values). Jun 2, 2018 at 5:35
• @AlexI Thanks for providing a signal, I will have a quick look at it but I do not understand what am I supposed to be looking for in the updated question (?). Is it about clarification? Further question? Something else? If you have found a satisfactory solution to the question that is not covered by the proposed ones, then you should write it as a self-answer to this question and accept it (?)
– A_A
Jun 2, 2018 at 5:43

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).