matplotlib.mlab.psd(...)
und scipy.signal.welch(...)
implementieren beide die Welch-Methode des durchschnittlichen Periodogramms, um die spektrale Leistungsdichte (PSD) eines Signals abzuschätzen. Unter der Annahme, dass äquivalente Parameter an jede Funktion übergeben werden, sollte die zurückgegebene PSD äquivalent sein.
Bei Verwendung einer einfachen sinusförmigen Testfunktion gibt es jedoch systematische Unterschiede zwischen den beiden Methoden, wenn die Anzahl der pro Fenster verwendeten Punkte (das nperseg
-Schlüsselwort für scipy.signal.welch
; das NFFT
-Schlüsselwort für mlab.psd
) ist even , wie unten für den Fall von 64 Punkten pro Fenster gezeigt
Das obere Diagramm zeigt die PSD, die mit beiden Methoden berechnet wurde, während das untere Diagramm ihren relativen Fehler anzeigt (d. h. die absolute Differenz der beiden PSDs dividiert durch ihren absoluten Durchschnitt). Ein größerer relativer Fehler weist auf größere Unstimmigkeiten zwischen den beiden Methoden hin.
Die beiden Funktionen haben eine viel bessere Übereinstimmung, wenn die Anzahl der pro Fenster verwendeten Punkte ungerade ist, wie unten für das gleiche Signal gezeigt, aber mit 65 Punkten pro Fenster Das Hinzufügen anderer Merkmale zum Signal (z. B. Rauschen) verringert diesen Effekt etwas, ist aber immer noch vorhanden, mit relativen Fehlern von ~ 10% zwischen den beiden Methoden, wenn eine gerade Anzahl von Punkten pro Fenster verwendet wird. (Mir ist klar, dass der relative Fehler bei der höchsten Frequenz für die mit 65 Punkten pro Fenster berechneten PSDs relativ groß ist. Dies ist jedoch bei der Nyquist-Frequenz des Signals, und ich bin nicht so besorgt über Features bei so hohen Frequenzen bin mehr besorgt über den großen und systematischen relativen Fehler über den Großteil der Bandbreite des Signals, wenn die Anzahl der Punkte pro Fenster gerade ist. Gibt es einen Grund für diese Diskrepanz? Ist eine Methode in Bezug auf die Genauigkeit der anderen vorzuziehen? Ich benutze scipy Version 0.16.0 und matplotlib Version 1.4.3. Der Code, der zum Erzeugen der obigen Plots verwendet wurde, ist unten enthalten. Stellen Sie für ein reines Sinussignal A_noise = 0
; Stellen Sie für ein verrauschtes Signal A_noise
auf einen endlichen Wert ein.
Obwohl die Parameter möglicherweise gleichwertig sind, kann sich der Fensterparameter für ein Fenster mit gerader Größe leicht unterscheiden. Genauer gesagt wird das von scipys welch
-Funktion verwendete Fenster mit
, das den Standardparameter fftbins=True
verwendet und gemäß scipy Dokumentation :
Falls True, erstellen Sie ein "periodisches" Fenster, das mit ifftshift verwendet werden kann und mit dem Ergebnis eines fft multipliziert wird (SIEHE AUCH fftfreq).
Dies führt zu einem anderen generierten Fenster für gerade Längen. Von diesem Abschnitt des Eintrags der Fensterfunktion in Wikipedia , könnte dies einen leichten Leistungsvorteil gegenüber Matplotlibs window_hanning
ergeben. was immer die symmetrische Version zurückgibt.
Um das gleiche Fenster zu verwenden, können Sie den Fenstervektor für beide PSD-Schätzfunktionen explizit angeben. Sie könnten zum Beispiel dieses Fenster mit:
berechnen %Vor% Wenn Sie dieses Fenster als Parameter verwenden (mit window=win
in beiden Funktionen), erhalten Sie das folgende Diagramm, in dem Sie eine viel engere Übereinstimmung zwischen den beiden PSD-Schätzfunktionen feststellen können:
Wenn Sie die Eigenschaft periodische Fenster nicht verwenden möchten, können Sie alternativ auch Folgendes verwenden:
%Vor%Tags und Links python-2.7 scipy signal-processing discrete-mathematics mlab