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，可选: x 和 y 时间序列的采样频率。默认为 1.0。
windowstr 或 tuple 或 array_like，可选: 要使用的期望窗口。如果 window 是字符串或元组，则将其传递给 get_window 以生成窗口值，默认情况下这些值是 DFT 偶数的。有关窗口和所需参数的列表，请参见 get_window。如果 window 是 array_like，则它将直接用作窗口，并且其长度必须为 nperseg。默认为汉宁窗口。
npersegint，可选: 每个段的长度。默认为 None，但是如果 window 是 str 或 tuple，则设置为 256，如果 window 是 array_like，则设置为窗口的长度。
noverlapint，可选: 段之间重叠的点数。如果 None，则 noverlap = nperseg // 2。默认为 None。
nfftint，可选: 如果需要零填充 FFT，则使用的 FFT 长度。如果 None，则 FFT 长度为 nperseg。默认为 None。
detrendstr 或函数或 False，可选: 指定如何对每个段进行去趋势。如果 detrend 是字符串，则将其作为 type 参数传递给 detrend 函数。如果它是一个函数，则它接受一个段并返回一个去趋势的段。如果 detrend 是 False，则不进行去趋势。默认为 'constant'。
return_onesidedbool，可选: 如果 True，则为实数数据返回单边频谱。如果 False，则返回双边频谱。默认为 True，但是对于复数数据，始终返回双边频谱。
scaling{ ‘density’, ‘spectrum’ }，可选: 选择计算互谱密度（‘density’），其中如果 x 和 y 的单位为 V 并且 fs 的单位为 Hz，则 Pxy 的单位为 V**2/Hz，或选择计算互谱（‘spectrum’），其中 Pxy 的单位为 V**2。默认为 ‘density’。
axisint，可选: 计算两个输入的 CSD 的轴；默认值是沿最后一个轴（即 axis=-1）。
average{ ‘mean’, ‘median’ }，可选: 平均周期图时使用的方法。如果频谱是复数，则分别计算实部和虚部的平均值。默认为 ‘mean’。

1.2.0 版本新增。

返回:

fndarray: 采样频率的数组。
Pxyndarray: x, y 的互谱密度或互功率谱。

另请参见

periodogram: 简单，可选修改的周期图
lombscargle: 用于不均匀采样数据的 Lomb-Scargle 周期图
welch: Welch 方法的功率谱密度。[等效于 csd(x,x)]
coherence: Welch 方法的幅度平方相干性。

备注

按照惯例，Pxy 的计算方式是将 X 的共轭 FFT 乘以 Y 的 FFT。

如果输入序列的长度不同，则会用零填充较短的序列以匹配。

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

有关谱密度和（幅度）频谱的缩放的讨论，请参阅频谱分析部分的 SciPy 用户指南。

0.16.0 版本新增。

参考

[1]

P. Welch，“使用快速傅里叶变换来估计功率谱：一种基于对短的修改周期图进行时间平均的方法”，IEEE Trans。音频电声。卷 15，第 70-73 页，1967 年。

[2]

Rabiner, Lawrence R., 和 B. Gold。“数字信号处理的理论与应用” Prentice-Hall，第 414-419 页，1975 年

示例

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

生成两个具有一些共同特征的测试信号。

>>> fs = 10e3
>>> N = 1e5
>>> amp = 20
>>> freq = 1234.0
>>> 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)

计算并绘制互谱密度的幅度。

>>> f, Pxy = signal.csd(x, y, fs, nperseg=1024)
>>> plt.semilogy(f, np.abs(Pxy))
>>> plt.xlabel('frequency [Hz]')
>>> plt.ylabel('CSD [V**2/Hz]')
>>> plt.show()