What is the fastest way to get closer to a data period with Octave?

Question

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

In particular, here is my current code. I would like to approximate the period of the vector x (:, 2) against the vector t. Ultimately, I would like to do this for a variety of initial conditions and calculate the period of each of them and draw a 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)); 

Thanks!

John

Answer 1

Take a look at the auto-correlation feature.

From Wikipedia

Autocorrelation is a mutual correlation of a signal with itself. Informally, this similarity between observations is a function of the time separation between them. This is a mathematical tool for finding repetitive patterns, such as the presence of a periodic signal that was buried by noise or identifying the missing main frequency in the signal implied by its harmonic frequencies. It is often used to process signals to analyze functions or series of values, such as time-domain signals.

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

Answer 2

The discrete Fourier transform can give you periodicity. A longer time window gives you a higher frequency resolution, so I changed the definition of t to t = linspace(0, 500, 2000) . time domain http://img402.imageshack.us/img402/8775/timedomain.png (here is a link to the chart , it looks better on the hosting site). You can do:

 h = hann(length(x), 'periodic'); %# use a Hann window to reduce leakage y = fft(x .* [hh]); %# 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 http://img406.imageshack.us/img406/2696/freqdomain.png Link to the chart.

Looking at the graph, the first peak is about 0.06 Hz, which corresponds to the period of 16 seconds observed in plot(t,x) .

It is not computational, but fast. FFTs are N * log (N) operations.