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

【问题标题】：What's the fastest way to approximate the period of data using Octave?使用 Octave 估算数据周期的最快方法是什么？
【发布时间】：2010-04-03 19:43:28
【问题描述】：

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

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

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

谢谢！

约翰

【问题讨论】：

标签： matlab signal-processing fft octave

【解决方案1】：

看一下自相关函数。

来自Wikipedia

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

Paul Bourke 描述了如何基于快速傅立叶变换有效地计算自相关函数 (link)。

【讨论】：

这比我的建议要好得多。在 Octave 中，您可以执行以下操作：[xx, lags] = xcorr(x); plot(lags, xx(:,[1 4]));（xx 的第 2 列和第 3 列是互相关）并查看两个自相关峰之间的分离。

【解决方案2】：

离散傅里叶变换可以为您提供周期性。更长的时间窗口为您提供更高的频率分辨率，因此我将您的 t 定义更改为 t = linspace(0, 500, 2000)。 time domain http://img402.imageshack.us/img402/8775/timedomain.png（这是link to the plot，它在托管站点上看起来更好）。你可以这样做：

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.

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

不过，这在计算上并没有那么快。 FFT 是 N*log(N) 次操作。

【讨论】：

FFT 不能很好地确定一般周期性。例如，如果信号是两个分量的总和，一个具有周期 5T，另一个具有周期 7T，则信号的周期将为 35T，但这不会显示在 FFT 中。自相关是此任务的更好选择。
@tom10 - FFT 方法取决于在观察中具有多个信号周期，但我认为对于任何非分析方法都是如此，不是吗？
样本长度不是我在这里所说的（尽管你是对的）。尝试绘制 sin(t/2) + sin(t/3) 并查看周期性，然后将其与 FFT 进行比较，您会发现周期和 FFT 并没有那么明显的相关性。
嗯... O(N log N) 相当快。如果您怀疑某个特定频率或频带，有更快的算法，但您肯定不会比 O(N) 做得更好。