使用 Octave 估算数据周期的最快方法是什么？

Question

我有一组周期性的数据（但不是正弦曲线）。我在一个向量中有一组时间值，在第二个向量中有一组幅度。我想快速估算出函数的周期。有什么建议吗？

具体来说，这是我当前的代码。我想近似向量 x(:,2) 相对于向量 t 的周期。最终，我想对许多初始条件执行此操作，并计算每个条件的周期并绘制结果。

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

谢谢你！

约翰

Answer 1

看一下自相关函数。

来自 维基百科

自相关是
信号与的互相关
本身。非正式地，它是
观测值之间的相似性
时间间隔函数
他们之间。它是一个数学
寻找重复模式的工具，
例如周期性的存在
信号已被掩埋
噪音，或识别失踪的
信号的基频
由其谐波频率暗示。
常用于信号处理
用于分析函数或系列
值，例如时域信号。

Paul Bourke 描述了如何基于快速傅里叶变换有效计算自相关函数（链接 ）。

Answer 2

离散傅立叶变换可以给出周期性。较长的时间窗口可以提供更高的频率分辨率，因此我将您的 t 定义更改为 t = linspace(0, 500, 2000)。
时域 http://img402.imageshack.us/img402/8775/timedomain.png< /a> （这里有一个绘图链接，它在托管站点）。
你可以这样做：

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

频域 http://img406.imageshack.us/img406/2696/ freqdomain.png 绘图链接。

查看图表，第一个峰值位于 0.06 Hz 左右，对应于 plot(t,x) 中看到的 16 秒周期。

但这在计算上并没有那么快。 FFT是N*log(N)运算。