Octave/Matlab与Scipy中的周期图 [英] Periodogram in Octave/Matlab vs Scipy

问题描述

我正在使用scipy将一些matlab代码移植到python，并陷入以下行:

I am porting some matlab code to python using scipy and got stuck with the following line:

Matlab/八度代码

Matlab/Octave code

[Pxx, f] = periodogram(x, [], 512, 5)

Python代码

f, Pxx = signal.periodogram(x, 5, nfft=512)

问题是我在同一数据上得到不同的输出.更具体地说，Pxx向量是不同的.我尝试了不同的signal.periodogram窗口，但是没有运气(似乎默认scypy的boxcar窗口与默认matlab的矩形窗口相同)另一个奇怪的行为是在python中，无论什么数据，Pxx的第一个元素始终为0输入是.

The problem is that I get different output on the same data. More specifically, Pxx vectors are different. I tried different windows for signal.periodogram, yet no luck (and it seems that default scypy's boxcar window is the same as default matlab's rectangular window) Another strange behavior is that in python, first element of Pxx is always 0, no matter what data input is.

我想念什么吗?任何建议将不胜感激！

Am i missing something? Any advice would be greatly appreciated!

带有实际数据的简单Matlab/八度代码: http://pastebin.com/czNeyUjs
具有实际数据的简单Python + scipy代码: http://pastebin.com/zPLGBTpn

Simple Matlab/Octave code with actual data: http://pastebin.com/czNeyUjs
Simple Python+scipy code with actual data: http://pastebin.com/zPLGBTpn

推荐答案

研究八度音阶和scipy的周期图源代码后，我发现它们使用不同的算法来计算功率谱密度估计值.八度(和MATLAB)使用FFT ，而scipy的周期图使用韦尔奇方法.

After researching octave's and scipy's periodogram source code I found that they use different algorithm to calculate power spectral density estimate. Octave (and MATLAB) use FFT, whereas scipy's periodogram use the Welch method.

正如@georgesl所提到的，输出看起来非常相似，但是仍然有所不同.出于移植的原因，这一点至关重要.最后，我只是编写了一个小函数来使用FFT计算PSD估算值，现在输出是相同的.根据timeit测试，它的运行速度提高了约50％(在1万次迭代的循环中，速度为1.9006s与2.9176s).我认为这是由于FFT在scipy的实现中比Welch更快，而只是更快.

As @georgesl has mentioned, the output looks quite alike, but still, it differs. And for porting reason it was critical. In the end, I simply wrote a small function to calculate PSD estimate using FFT, and now output is the same. According to timeit testing, it works ~50% faster (1.9006s vs 2.9176s on a loop with 10.000 iterations). I think it's due to the FFT being faster than Welch in scipy's implementation, of just being faster.

感谢所有表现出兴趣的人.

Thanks to everyone who showed interest.

这篇关于Octave/Matlab与Scipy中的周期图的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！