0
$\begingroup$

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)

Signal

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.

enter image description here

| improve this question | |
$\endgroup$
  • 1
    $\begingroup$ Can I please ask if this was resolved? $\endgroup$ – A_A Jun 11 '18 at 5:49
  • $\begingroup$ @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. $\endgroup$ – Alex I Dec 18 '19 at 23:11
  • $\begingroup$ @A_A Let me know if you have any other ideas - I think this is a pretty useful problem to solve in general. $\endgroup$ – Alex I Dec 18 '19 at 23:11
  • $\begingroup$ Can you please post a representative signal that you are working with? $\endgroup$ – A_A Dec 19 '19 at 11:36
2
$\begingroup$

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:

enter image description here

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

enter image description here

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:

enter image description here

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:

enter image description here

| improve this answer | |
$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – Alex I May 9 '18 at 6:36
  • $\begingroup$ @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? $\endgroup$ – A_A May 9 '18 at 6:52
  • $\begingroup$ I uploaded an example of a synthetic signal (very similar to real one) here: pastebin.com/jUFYXyrA $\endgroup$ – Alex I Jun 2 '18 at 5:23
  • $\begingroup$ 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). $\endgroup$ – Alex I Jun 2 '18 at 5:35
  • $\begingroup$ @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 (?) $\endgroup$ – A_A Jun 2 '18 at 5:43
1
$\begingroup$

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

| improve this answer | |
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.