2013-06-16 20 views
8

jestem przenoszenie kodu Matlab do Python przy użyciu scipy i utknął z następującej linii:Periodogram w Octave/Matlab vs scipy

Matlab/Octave kod

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

kod Python

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

Problem polega na tym, że otrzymuję różne dane wyjściowe na tych samych danych. Dokładniej, wektory Pxx są różne. Próbowałem różnych okien dla signal.periodogram, ale bez powodzenia (i wygląda na to, że domyślne okno scypijskiego boxa jest takie samo, jak domyślne okno matlaba). Innym dziwnym zachowaniem jest to, że w pythonie, pierwszym elementem Pxx jest zawsze 0, bez względu na dane wejście to.

Czy brakuje mi czegoś? Każda rada byłaby bardzo doceniona!


Prosty kod Matlab/Octave z rzeczywistych danych: http://pastebin.com/czNeyUjs
Proste Python + kod scipy z rzeczywistymi danymi: http://pastebin.com/zPLGBTpn

+0

Spisek byłby bardzo pomocny. –

Odpowiedz

5

Po zbadaniu kodu źródłowego periodogram oktawę i scipy za Okazało się, że używają innego algorytmu do obliczania widmo gęstości widmowej mocy. Octave (i MATLAB) use FFT, podczas gdy periodogram scipy używa Welch method.

Jak wspomniał @georgesl, wyjście wygląda dość podobnie, ale mimo to różni się. A dla przeniesienia rozsądku było to krytyczne. W końcu napisałem małą funkcję do obliczenia oszacowania PSD za pomocą FFT, a teraz wynik jest taki sam. Według testów timeit działa on o około 50% szybciej (1.9006 s vs 2.9176s w pętli z 10.000 iteracjami). Myślę, że to dlatego, że FFT jest szybszy niż Welch w implementacji scipy, a po prostu jest szybszy.

Dziękuję wszystkim, którzy okazali zainteresowanie.

0

Czy spojrzeć na wyniki?

Outputs

Niewielkie różnice między tymi dwoma wynikami może pochodzi z optymalizacje/default Windows/wdrożeń/cokolwiek itp

1

I w obliczu tego samego problemu, ale wtedy natknąłem dokumentacji scipy's periodogram

Jak można zobaczyć tam, że detrend = „stałej” jest domyślnym argumentem. Oznacza to, że python automatycznie odejmuje średnią danych wejściowych z każdego punktu. (Read here). Podczas gdy Matlab/Octave nie robią czegoś takiego. Uważam, że jest to powód, dla którego wyniki są różne. Spróbuj podać detrend = False, podczas wywoływania periodogramu scipy powinieneś otrzymać takie samo wyjście jak Matlab.

2

Po przeczytaniu dokumentacji Matlaba i Scipy, innym wkładem do różnych wartości może być to, że używają innej domyślnej funkcji okna. Matlab używa okna Hamminga, a Scipy używa Hanninga. Dwie funkcje okna i podobne, ale nie identyczne.