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')[源]#

使用Welch方法估计互功率谱密度 Pxy。

参数:
xarray_like

测量值的时间序列

yarray_like

测量值的时间序列

fsfloat, optional

`x` 和 `y` 时间序列的采样频率。默认为1.0。

windowstr or tuple or array_like, optional

要使用的期望窗口。如果 `window` 是字符串或元组，它将被传递给 get_window 以生成窗口值，这些值默认是DFT偶对称的。有关窗口列表和所需参数，请参阅 get_window。如果 `window` 是 array_like 类型，它将直接用作窗口，并且其长度必须为 nperseg。默认为汉宁（Hann）窗口。

npersegint, optional

每个段的长度。默认为 None，但如果 `window` 是字符串或元组，则设置为256；如果 `window` 是 array_like，则设置为 `window` 的长度。

noverlap: int, optional

段之间的重叠点数。如果为 `None`，则 noverlap = nperseg // 2。默认为 `None`，且不得大于 `nperseg`。

nfftint, optional

使用的FFT长度，如果需要零填充FFT。如果为 `None`，则 FFT 长度为 `nperseg`。默认为 `None`。

detrendstr or function or False, optional

指定如何对每个段进行去趋势处理。如果 detrend 是字符串，它将作为 type 参数传递给 detrend 函数。如果它是一个函数，它会接受一个段并返回一个去趋势处理后的段。如果 detrend 为 False，则不执行去趋势处理。默认为 'constant'。

return_onesidedbool, optional

如果为 True，则返回实数数据的一侧频谱。如果为 False，则返回双侧频谱。默认为 True，但对于复数数据，始终返回双侧频谱。

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

选择计算互谱密度（‘density’），此时 Pxy 的单位是 V**2/Hz；或者计算互谱（‘spectrum’），此时 Pxy 的单位是 V**2，前提是 x 和 y 以伏特（V）测量，fs 以赫兹（Hz）测量。默认为 'density'

axisint, optional

计算CSD的轴，适用于两个输入；默认为最后一个轴（即 axis=-1）。

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

平均周期图时使用的方法。如果频谱是复数，则分别计算实部和虚部的平均值。默认为 'mean'。

1.2.0版本新增。

返回:
fndarray

采样频率数组。

Pxyndarray

`x`,`y`的互谱密度或互功率谱。

另请参阅

periodogram

简单，可选修改的周期图

lombscargle

用于不均匀采样数据的Lomb-Scargle周期图

welch

使用Welch方法计算的功率谱密度。[等同于 csd(x,x)]

coherence

使用Welch方法计算的幅度平方相干性。

备注

按照惯例，Pxy 是通过 X 的共轭FFT乘以 Y 的FFT来计算的。

如果输入序列长度不同，较短的序列将进行零填充以匹配。

适当的重叠量将取决于窗口的选择和您的要求。对于默认的汉宁（Hann）窗口，50%的重叠是准确估计信号功率与不过度计算任何数据之间的合理折衷。更窄的窗口可能需要更大的重叠。

互谱（scaling='spectrum'）除以互谱密度（scaling='density'）的比率是常数因子 sum(abs(window)**2)*fs / abs(sum(window))**2。如果 return_onesided 为 True，则负频率的值将添加到相应正频率的值中。

请查阅 频谱分析 部分的 SciPy 用户指南 以获取关于谱密度和（幅度）谱的缩放比例的讨论。

Welch方法可以解释为对（互）谱图的时间切片进行平均。在内部，此函数利用 ShortTimeFFT 来确定所需的（互）谱图。下面的示例说明，直接使用 ShortTimeFFT 计算 Pxy 是很简单的。然而，ShortTimeFFT 的行为存在一些显著差异：

  • 对于 csd 参数组合 return_onesided=True, scaling='density'，没有直接的 ShortTimeFFT 等效项，因为 fft_mode='onesided2X' 需要 `'psd'` 缩放。这是因为 csd 在这种情况下返回的是两倍的平方幅度，这没有合理的解释。

  • ShortTimeFFT 内部使用 float64 / complex128，这是由于所使用的 fft 模块的行为所致。因此，返回的数据类型是这些。如果输入也是 float32 / complex64，csd 函数会将返回值转换为 float32 / complex64。

  • csd 函数计算 np.conj(Sx[q,p]) * Sy[q,p]，而 spectrogram 计算 Sx[q,p] * np.conj(Sy[q,p])，其中 Sx[q,p]、Sy[q,p] 是 x 和 y 的STFT。此外，窗口定位也不同。

0.16.0版本新增。

参考文献

[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. vol. 15, pp. 70-73, 1967.

[2]

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

示例

以下示例绘制了两个具有某些共同特征的信号的互功率谱密度

>>> 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

互谱密度是通过对频谱图的时间切片进行平均来计算的

>>> 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

如“备注”部分所述，使用与上述代码片段类似的方法和 csd 函数的结果可能会因实现细节而有所不同。

请注意，上述代码片段可以轻松调整以确定除平均值之外的其他统计属性。

本页内容