views:

123

answers:

2

I have a set of data that is periodic (but not sinusoidal). I have a set of time values in one vector and a set of amplitudes in a second vector. I'd like to quickly approximate the period of the function. Any suggestions?

Specifically, here's my current code. I'd like to approximate the period of the vector x(:,2) against the vector t. Ultimately, I'd like to do this for lots of initial conditions and calculate the period of each and plot the result.

function xdot = f (x,t)
         xdot(1) =x(2);
         xdot(2) =-sin(x(1));
endfunction

x0=[1;1.75];     #eventually, I'd like to try lots of values for x0(2)
t = linspace (0, 50, 200);


x = lsode ("f", x0, t)

plot(x(:,1),x(:,2));

Thank you!

John

+2  A: 

The Discrete Fourier Transform can give you the periodicity. A longer time window gives you more frequency resolution so I changed your t definition to t = linspace(0, 500, 2000). time domain (here's a link to the plot, it looks better on the hosting site). You could do:

h = hann(length(x), 'periodic'); %# use a Hann window to reduce leakage
y = fft(x .* [h h]); %# window each time signal and calculate FFT
df = 1/t(end); %# if t is in seconds, df is in Hz
ym = abs(y(1:(length(y)/2), :)); %# we just want amplitude of 0..pi frequency components
semilogy(((1:length(ym))-1)*df, ym);

frequency domain Plot link.

Looking at the graph, the first peak is at around 0.06 Hz, corresponding to the 16 second period seen in plot(t,x).

This isn't computationally that fast though. The FFT is N*log(N) operations.

mtrw
The FFT won't work well for determining general periodicity. For example, if the signal is the sum of two components, one with period 5T and the other with period 7T, the signal's period will be 35T, but this won't show in the FFT. The autocorrelation is a better choice for this task.
tom10
@tom10 - The FFT method depends on having several periods of the signal in the observation, but I think that would be true for any non-analytical method, wouldn't it?
mtrw
@mtrw The sample length is not at what I'm talking about here (though you're correct about it). Try plotting sin(t/2) + sin(t/3) and look at the periodicity, and then compare this to the FFT, and you can see that the period and FFT aren't so obviously related.
tom10
@tom10 - Really good point, I didn't think of that.
mtrw
ummm... O(N log N) is pretty fast. If you suspect a particular frequency or band of frequencies, there are faster algorithms, but you're definitely not going to do better than O(N).
Jason S
+5  A: 

Take a look at the auto correlation function.

From Wikipedia

Autocorrelation is the cross-correlation of a signal with itself. Informally, it is the similarity between observations as a function of the time separation between them. It is a mathematical tool for finding repeating patterns, such as the presence of a periodic signal which has been buried under noise, or identifying the missing fundamental frequency in a signal implied by its harmonic frequencies. It is often used in signal processing for analyzing functions or series of values, such as time domain signals.

Paul Bourke has a description of how to calculate the autocorrelation function effectively based on the fast fourier transform (link).

midtiby
This is much better than my suggestion. In Octave, you can do the following: `[xx, lags] = xcorr(x); plot(lags, xx(:,[1 4]));` (the 2nd and 3rd columns of xx are the cross-correlations) and look at the separation between peaks of the two autocorrelations.
mtrw