Python SciPy signal.coherence用法及代碼示例

本文簡要介紹 python 語言中 scipy.signal.coherence 的用法。

用法:

scipy.signal.coherence(x, y, fs=1.0, window='hann', nperseg=None, noverlap=None, nfft=None, detrend='constant', axis=-1)#

使用 Welch 方法估計 discrete-time 信號 X 和 Y 的幅度平方相幹估計值 Cxy。

Cxy = abs(Pxy)**2/(Pxx*Pyy)，其中像素和派伊是 X 和 Y 的功率譜密度估計，以及像素是 X 和 Y 的交叉譜密度估計。

參數 ：：

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

axis： 整數，可選

為兩個輸入計算相幹性的軸；默認值在最後一個軸上(即 axis=-1 )。

返回 ：：

f： ndarray

采樣頻率數組。

Cxy： ndarray

x 和 y 的幅度平方相幹性。

注意：

適當的重疊量取決於窗口的選擇和您的要求。對於默認的 Hann 窗口，50% 的重疊是在準確估計信號功率和不過度計算任何數據之間的合理權衡。較窄的窗口可能需要較大的重疊。

參考：

[1]

P. Welch，“使用快速傅裏葉變換來估計功率譜：一種基於對短的修改周期圖進行時間平均的方法”，IEEE Trans。音頻電聲。卷。 15，第 70-73 頁，1967 年。

[2]

Stoica、Petre 和 Randolph Moses，“信號的頻譜分析”Prentice Hall，2005 年

例子：

>>> 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, Cxy = signal.coherence(x, y, fs, nperseg=1024)
>>> plt.semilogy(f, Cxy)
>>> plt.xlabel('frequency [Hz]')
>>> plt.ylabel('Coherence')
>>> plt.show()