scipy.signal.

csd#

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

Schätzt die Kreuzleistungsspektraldichte, Pxy, mit der Welch-Methode.

Parameter:

xarray_like: Zeitreihen von Messwerten
yarray_like: Zeitreihen von Messwerten
fsfloat, optional: Abtastfrequenz der Zeitreihen x und y. Standard 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.
noverlap: int, optional: Anzahl der Punkte, die zwischen den Segmenten überlappen. Wenn None, ist noverlap = nperseg // 2. Standard ist None und darf nicht größer als nperseg sein.
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 der Kreuzspektraldichte ('density'), bei der Pxy die Einheiten V**2/Hz hat, und der Berechnung des Kreuzspektrums ('spectrum'), bei dem Pxy die Einheiten V**2 hat, wenn x und y in V und fs in Hz gemessen werden. Standard ist 'density'.
axisint, optional: Achse, entlang derer die CSD für beide Eingaben berechnet wird; Standard ist die letzte Achse (d.h. axis=-1).
average{ ‘mean’, ‘median’ }, optional: Methode, die beim Mitteln von Periodogrammen verwendet wird. Wenn das Spektrum komplex ist, wird der Mittelwert getrennt für den Real- und Imaginärteil berechnet. Standard ist 'mean'.

Hinzugefügt in Version 1.2.0.

Rückgabe:

fndarray: Array von Abtastfrequenzen.
Pxyndarray: Kreuzspektraldichte oder Kreuzleistungsspektrum von x,y.

Siehe auch

periodogram: Einfaches, optional modifiziertes Periodogramm
lombscargle: Lomb-Scargle-Periodogramm für ungleichmäßig abgetastete Daten
welch: Leistungsspektraldichte nach der Welch-Methode. [Äquivalent zu csd(x,x)]
coherence: Betragsquadrat-Kohärenz nach der Welch-Methode.

Hinweise

Per Konvention wird Pxy mit der konjugierten FFT von X multipliziert mit der FFT von Y berechnet.

Wenn die Eingangsreihen unterschiedliche Längen haben, wird die kürzere Reihe mit Nullen aufgefüllt, um sie anzupassen.

Ein angemessener Überlappungsgrad hängt von der Wahl des Fensters und von 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 einer Überzählung von Daten. Engere Fenster erfordern möglicherweise eine größere Überlappung.

Das Verhältnis des Kreuzspektrums (scaling='spectrum') zum Kreuzleistungsspektrum (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 addiert.

Konsultieren Sie den Abschnitt Spektralanalyse des SciPy-Benutzerhandbuchs für eine Diskussion über die Skalierungen einer Spektraldichte und eines (Amplituden-)Spektrums.

Die Welch-Methode kann als Mittelwertbildung über die Zeitschlitze eines (Kreuz-)Spektrogramms interpretiert werden. Intern verwendet diese Funktion die ShortTimeFFT, um das erforderliche (Kreuz-)Spektrogramm zu bestimmen. Ein Beispiel unten zeigt, dass es einfach ist, Pxy direkt mit der ShortTimeFFT zu berechnen. Es gibt jedoch einige bemerkenswerte Unterschiede im Verhalten der ShortTimeFFT

Es gibt kein direktes ShortTimeFFT-Äquivalent für die csd-Parameterkombination return_onesided=True, scaling='density', da fft_mode='onesided2X' die Skalierung 'psd' erfordert. Dies liegt daran, dass csd in diesem Fall das verdoppelte quadrierte Betragsquadrat zurückgibt, was keine sinnvolle Interpretation hat.
ShortTimeFFT verwendet intern float64 / complex128, was auf das Verhalten des verwendeten fft-Moduls zurückzuführen ist. Daher sind dies die zurückgegebenen Dtypes. Die Funktion csd wandelt die Rückgabewerte in float32 / complex64 um, wenn die Eingabe ebenfalls float32 / complex64 ist.
Die Funktion csd berechnet np.conj(Sx[q,p]) * Sy[q,p], während spectrogram Sx[q,p] * np.conj(Sy[q,p]) berechnet, wobei Sx[q,p] und Sy[q,p] die STFTs von x und y sind. Außerdem ist die Fensterpositionierung unterschiedlich.

Hinzugefügt in Version 0.16.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]

Rabiner, Lawrence R., und B. Gold. „Theory and Application of Digital Signal Processing“ Prentice-Hall, S. 414-419, 1975

Beispiele

Das folgende Beispiel plottet die Kreuzleistungsspektraldichte zweier Signale mit einigen gemeinsamen Merkmalen

>>> import numpy as np
>>> from scipy import signal
>>> import matplotlib.pyplot as plt
>>> rng = np.random.default_rng()
...
... # Generate two test signals with some common features:
>>> N, fs = 100_000, 10e3  # number of samples and sampling frequency
>>> amp, freq = 20, 1234.0  # amplitude and frequency of utilized sine signal
>>> noise_power = 0.001 * fs / 2
>>> time = np.arange(N) / fs
>>> b, a = signal.butter(2, 0.25, 'low')
>>> x = rng.normal(scale=np.sqrt(noise_power), size=time.shape)
>>> y = signal.lfilter(b, a, x)
>>> x += amp*np.sin(2*np.pi*freq*time)
>>> y += rng.normal(scale=0.1*np.sqrt(noise_power), size=time.shape)
...
... # Compute and plot the magnitude of the cross spectral density:
>>> nperseg, noverlap, win = 1024, 512, 'hann'
>>> f, Pxy = signal.csd(x, y, fs, win, nperseg, noverlap)
>>> fig0, ax0 = plt.subplots(tight_layout=True)
>>> ax0.set_title(f"CSD ({win.title()}-window, {nperseg=}, {noverlap=})")
>>> ax0.set(xlabel="Frequency $f$ in kHz", ylabel="CSD Magnitude in V²/Hz")
>>> ax0.semilogy(f/1e3, np.abs(Pxy))
>>> ax0.grid()
>>> plt.show()

../../_images/scipy-signal-csd-1_00_00.png

Die Kreuzleistungsspektraldichte wird durch Mittelung über die Zeitschlitze eines Spektrogramms berechnet

>>> SFT = signal.ShortTimeFFT.from_window('hann', fs, nperseg, noverlap,
...                                       scale_to='psd', fft_mode='onesided2X',
...                                       phase_shift=None)
>>> Sxy1 = SFT.spectrogram(y, x, detr='constant', k_offset=nperseg//2,
...                        p0=0, p1=(N-noverlap) // SFT.hop)
>>> Pxy1 = Sxy1.mean(axis=-1)
>>> np.allclose(Pxy, Pxy1)  # same result as with csd()
True

Wie im Abschnitt Hinweise erläutert, können die Ergebnisse der Verwendung eines Ansatzes, der analog zum obigen Code-Snippet ist, und der Funktion csd aufgrund von Implementierungsdetails abweichen.

Beachten Sie, dass das obige Code-Snippet leicht angepasst werden kann, um andere statistische Eigenschaften als den Mittelwert zu bestimmen.