Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

Periodogram in Octave/Matlab vs Scipy

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

Matlab/Octave code

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

Python code

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

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!


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

like image 742
mFoxRU Avatar asked Jun 16 '13 06:06

mFoxRU


2 Answers

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.

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.

like image 124
mFoxRU Avatar answered Sep 20 '22 13:09

mFoxRU


I faced the same problem but then I came across the documentation of scipy's periodogram

As you would see there that detrend='constant' is the default argument. This means that python automatically subtracts the mean of the input data from each point. (Read here). While Matlab/Octave do no such thing. I believe that is the reason why the outputs are different. Try specifying detrend=False, while calling scipy's periodogram you should get the same output as Matlab.

like image 34
Utkarsh Dixit Avatar answered Sep 22 '22 13:09

Utkarsh Dixit