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