Python SciPy signal.csd用法及代码示例

本文简要介绍 python 语言中 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。

参数 ：：

x： array_like

测量值的时间序列

y： array_like

测量值的时间序列

fs： 浮点数，可选

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

window： str 或 tuple 或 数组，可选

想要使用的窗口。如果窗户是一个字符串或元组，它被传递给scipy.signal.get_window生成窗口值，默认为DFT-even。看scipy.signal.get_window获取窗口列表和所需参数。如果窗户是数组，它将直接用作窗口，其长度必须为nperseg。默认为 Hann 窗口。

nperseg： 整数，可选

每个段的长度。默认为 None，但如果 window 是 str 或 tuple，则设置为 256，如果 window 是 数组，则设置为窗口的长度。

noverlap: int, optional：

段之间重叠的点数。如果None,noverlap = nperseg // 2.默认为None.

nfft： 整数，可选

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

detrend： str 或函数 或False， 可选的

指定如何去除每个段的趋势。如果scipy.signal.detrend是一个字符串，它作为类型参数scipy.signal.detrend函数。如果它是一个函数，它接受一个段并返回一个去趋势的段。如果scipy.signal.detrend是False，没有去趋势。默认为‘constant’。

return_onesided： 布尔型，可选

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

scaling： { ‘density’, ‘spectrum’ }，可选

选择计算交叉谱密度 (‘density’)，其中 Pxy 的单位为 V**2/Hz，计算交叉谱 (‘spectrum’)，其中 Pxy 的单位为 V**2，如果测量 x 和 y以 V 为单位，fs 以 Hz 为单位。默认为‘density’

axis： 整数，可选

为两个输入计算 CSD 的轴；默认值在最后一个轴上(即 axis=-1 )。

average： { ‘mean’, ‘median’ }，可选

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

返回 ：：

f： ndarray

采样频率数组。

Pxy： ndarray

x,y 的交叉谱密度或交叉功率谱。

注意：

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

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

适当的重叠量取决于窗口的选择和您的要求。对于默认的 Hann 窗口，50% 的重叠是在准确估计信号功率和不过度计算任何数据之间的合理权衡。较窄的窗口可能需要较大的重叠。

参考：

[1]

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

[2]

拉宾纳、劳伦斯 R. 和 B. Gold。 “数字信号处理的理论与应用”Prentice-Hall, pp. 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()