scipy.signal.

welch

scipy.signal.welch(x, fs=1.0, window='hann', nperseg=None, noverlap=None, nfft=None, detrend='constant', return_onesided=True, scaling='density', axis=-1, average='mean')

Schätzt die Leistungsdichtespektrum mittels der Welch-Methode.

Die Welch-Methode [1] berechnet eine Schätzung des Leistungsdichtespektrums, indem die Daten in überlappende Segmente unterteilt werden, für jedes Segment ein modifiziertes Periodogramm berechnet und die Periodogramme gemittelt werden.

Parameter:

xarray_like

Zeitreihen von Messwerten

fsfloat, optional

Abtastfrequenz der Zeitreihe x. Standardwert ist 1.0.

windowstr oder Tupel oder Array_like, optional

Gewünschtes Fenster zur Verwendung. Wenn window eine Zeichenkette oder ein Tupel ist, wird es an get_window übergeben, um die Fensterwerte zu generieren, die standardmäßig DFT-gerade sind. Siehe get_window für eine Liste von Fenstern und erforderliche Parameter. Wenn window Array_like ist, wird es direkt als Fenster verwendet und seine Länge muss nperseg sein. Standard ist ein Hann-Fenster.

npersegint, optional

Länge jedes Segments. Standard ist None, aber wenn window str oder tuple ist, wird es auf 256 gesetzt, und wenn window array_like ist, wird es auf die Länge des Fensters gesetzt.

noverlapint, optional

Anzahl der Punkte, die zwischen Segmenten überlappen. Wenn None, dann ist noverlap = nperseg // 2. Standardwert ist None.

nfftint, optional

Länge der FFT, die verwendet wird, wenn eine Null-gepolsterte FFT gewünscht wird. Wenn None, ist die FFT-Länge nperseg. Standard ist None.

detrendstr oder Funktion oder False, optional

Gibt an, wie jedes Segment entzerrt werden soll. Wenn detrend eine Zeichenkette ist, wird sie als Argument type an die Funktion detrend übergeben. Wenn es sich um eine Funktion handelt, nimmt sie ein Segment entgegen und gibt ein entzerrtes Segment zurück. Wenn detrend False ist, wird keine Entzerrung durchgeführt. Standard ist 'constant'.

return_onesidedbool, optional

Wenn True, wird ein einseitiges Spektrum für reale Daten zurückgegeben. Wenn False, wird ein zweiseitiges Spektrum zurückgegeben. Standard ist True, aber für komplexe Daten wird immer ein zweiseitiges Spektrum zurückgegeben.

scaling{ ‘density’, ‘spectrum’ }, optional

Wählt zwischen der Berechnung des Leistungsdichtespektrums ('density'), bei der Pxx die Einheiten V**2/Hz hat, und der Berechnung des quadrierten Amplitudenspektrums ('spectrum'), bei der Pxx die Einheiten V**2 hat, wenn x in V und fs in Hz gemessen wird. Standardwert ist 'density'.

axisint, optional

Achse, entlang der das Periodogramm berechnet wird; Standard ist die letzte Achse (d.h. axis=-1).

average{ ‘mean’, ‘median’ }, optional

Methode zur Mittelung der Periodogramme. Standardwert ist 'mean'.

Hinzugefügt in Version 1.2.0.

Rückgabe:

fndarray

Array von Abtastfrequenzen.

Pxxndarray

Leistungsdichtespektrum oder Leistungsspektrum von x.

Siehe auch

csd

Kreuz-Leistungsdichtespektrum mittels der Welch-Methode

periodogram

Einfaches, optional modifiziertes Periodogramm

lombscargle

Lomb-Scargle-Periodogramm für ungleichmäßig abgetastete Daten

Hinweise

Eine geeignete Überlappungsmenge hängt von der Wahl des Fensters und Ihren Anforderungen ab. Für das Standard-Hann-Fenster ist eine Überlappung von 50 % ein vernünftiger Kompromiss zwischen der genauen Schätzung der Signalenergie und der Vermeidung von Doppelzählungen von Daten. Engere Fenster erfordern möglicherweise eine größere Überlappung. Wenn noverlap 0 ist, ist diese Methode äquivalent zur Bartlett-Methode [2].

Das Verhältnis des quadrierten Betrags (scaling='spectrum') zum Leistungsdichtespektrum (scaling='density') ist der konstante Faktor sum(abs(window)**2)*fs / abs(sum(window))**2. Wenn return_onesided True ist, werden die Werte der negativen Frequenzen zu den Werten der entsprechenden positiven Frequenzen addiert.

Konsultieren Sie den Abschnitt Spektralanalyse im SciPy-Benutzerhandbuch für eine Diskussion der Skalierungen des Leistungsdichtespektrums und des (quadrierten) Betragsspektrums.

Hinzugefügt in Version 0.12.0.

Referenzen

[1]

P. Welch, „The use of the fast Fourier transform for the estimation of power spectra: A method based on time averaging over short, modified periodograms“, IEEE Trans. Audio Electroacoust. Bd. 15, S. 70-73, 1967.

[2]

M.S. Bartlett, „Periodogram Analysis and Continuous Spectra“, Biometrika, vol. 37, S. 1-16, 1950.

Beispiele

Probieren Sie es in Ihrem Browser aus!

>>> import numpy as np
>>> from scipy import signal
>>> import matplotlib.pyplot as plt
>>> rng = np.random.default_rng()

Generieren Sie ein Testsignal, eine Sinuswelle mit 2 Veff bei 1234 Hz, verrauscht mit 0,001 V**2/Hz weißem Rauschen, abgetastet mit 10 kHz.

>>> fs = 10e3
>>> N = 1e5
>>> amp = 2*np.sqrt(2)
>>> freq = 1234.0
>>> noise_power = 0.001 * fs / 2
>>> time = np.arange(N) / fs
>>> x = amp*np.sin(2*np.pi*freq*time)
>>> x += rng.normal(scale=np.sqrt(noise_power), size=time.shape)

Berechnen und plotten Sie das Leistungsdichtespektrum.

>>> f, Pxx_den = signal.welch(x, fs, nperseg=1024)
>>> plt.semilogy(f, Pxx_den)
>>> plt.ylim([0.5e-3, 1])
>>> plt.xlabel('frequency [Hz]')
>>> plt.ylabel('PSD [V**2/Hz]')
>>> plt.show()

Wenn wir die letzte Hälfte des Leistungsspektrums mitteln, um den Peak auszuschließen, können wir die Rauschleistung des Signals wiederherstellen.

>>> np.mean(Pxx_den[256:])
0.0009924865443739191

Berechnen und plotten Sie nun das Leistungsspektrum.

>>> f, Pxx_spec = signal.welch(x, fs, 'flattop', 1024, scaling='spectrum')
>>> plt.figure()
>>> plt.semilogy(f, np.sqrt(Pxx_spec))
>>> plt.xlabel('frequency [Hz]')
>>> plt.ylabel('Linear spectrum [V RMS]')
>>> plt.show()

Die Peak-Höhe im Leistungsspektrum ist eine Schätzung der Effektivwertamplitude.

>>> np.sqrt(Pxx_spec.max())
2.0077340678640727

Wenn wir nun eine Diskontinuität im Signal einführen, indem wir die Amplitude eines kleinen Teils des Signals um 50 erhöhen, können wir die Verfälschung des durchschnittlichen Leistungsdichtespektrums sehen, während eine mediane Mittelung das normale Verhalten besser schätzt.

>>> x[int(N//2):int(N//2)+10] *= 50.
>>> f, Pxx_den = signal.welch(x, fs, nperseg=1024)
>>> f_med, Pxx_den_med = signal.welch(x, fs, nperseg=1024, average='median')
>>> plt.semilogy(f, Pxx_den, label='mean')
>>> plt.semilogy(f_med, Pxx_den_med, label='median')
>>> plt.ylim([0.5e-3, 1])
>>> plt.xlabel('frequency [Hz]')
>>> plt.ylabel('PSD [V**2/Hz]')
>>> plt.legend()
>>> plt.show()

ZurückIn Tab öffnen