What's the fastest way to approximate the peri

问题:

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

回答1:

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

回答2:

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 http://img402.imageshack.us/img402/8775/timedomain.png (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 http://img406.imageshack.us/img406/2696/freqdomain.png 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.