两种方法计算信号相干性结果不一致且不符合预期的技术求助

问题描述 #

尝试通过两种方法计算两个信号的相干性，预期结果一致且在250Hz、400Hz处相干性为1，但实际结果不符合预期：

• 方法1：先计算互相关/自相关，再通过DFT得到谱密度，进而计算相干性
• 方法2：直接调用scipy.signal.coherence

原始代码如下：

import numpy as np
from matplotlib import pyplot as plt
from scipy.signal import correlate,coherence

t  = np.arange(0, 1, 0.001);
dt = t[1] - t[0]
fs=1/dt

f0=250
f1=400

x=np.cos(2*np.pi*f0*t)+np.cos(2*np.pi*f1*t)
y=np.cos(2*np.pi*f0*(t-.035))+np.cos(2*np.pi*f1*(t-.05))

fig, ax = plt.subplots(2,sharex=True)

# Coherence using method 1:
Rxx = correlate(x,x)
Pxx= abs(np.fft.rfft(Rxx))

Ryy = correlate(y,y)
Pyy = abs(np.fft.rfft(Ryy))

Rxy = correlate(x,y)
Pxy = abs(np.fft.rfft(Rxy))

f = np.fft.rfftfreq(len(Rxy))*fs
Coh1=np.divide(Pxy**2,np.multiply(Pxx,Pyy))

#Start at nonzero index, e.g. 50, since Pxx[0]=Pyy[0] where Coh[0] would be undefined
ax[0].plot(f[50:], Coh1[50:]) 
ax[0].set_xlabel('frequency [Hz]')
ax[0].set_ylabel('Coherence')

# Coherence using method 2:
f,Coh2=coherence(x,y)

ax[1].plot(f*fs,Coh2)
ax[1].set_xlabel('frequency [Hz]')
ax[1].set_ylabel('Coherence')

plt.show()

错误分析 #

方法1的问题 #

• 互谱计算丢失相位信息：直接对互相关取绝对值后做FFT，正确做法是保留互相关FFT的复数形式，再取模平方。
• 未做归一化：correlate返回未归一化的相关值，功率谱需要归一化到正确的能量尺度，FFT结果也需匹配功率谱定义做缩放。
• 频谱泄漏：有限长度信号直接做FFT会引入泄漏，需加窗抑制。

方法2的问题 #

• 采样率参数缺失：coherence默认fs=1，手动乘fs易出错，应直接传入fs=fs参数。
• Welch方法参数不匹配预期：默认的汉明窗、分段重叠设置会平滑频谱，对于无噪声纯正弦信号，需用无窗、全长度分段的周期图方法才能得到尖锐峰值。

修正后的代码 #

import numpy as np
from matplotlib import pyplot as plt
from scipy.signal import correlate, coherence, get_window

t = np.arange(0, 1, 0.001)
dt = t[1] - t[0]
fs = 1/dt

f0 = 250
f1 = 400

x = np.cos(2*np.pi*f0*t) + np.cos(2*np.pi*f1*t)
y = np.cos(2*np.pi*f0*(t - 0.035)) + np.cos(2*np.pi*f1*(t - 0.05))

fig, ax = plt.subplots(2, sharex=True, figsize=(10, 8))

# --- 修正后的方法1 ---
# 使用same模式相关并归一化，匹配信号长度
Rxx = correlate(x, x, mode='same') / len(x)
Ryy = correlate(y, y, mode='same') / len(y)
Rxy = correlate(x, y, mode='same') / len(x)

# 加矩形窗抑制频谱泄漏
window = get_window('boxcar', len(x))
Rxx_windowed = Rxx * window
Ryy_windowed = Ryy * window
Rxy_windowed = Rxy * window

# 计算谱密度并归一化，保留复数形式
Pxx = np.fft.rfft(Rxx_windowed) * dt
Pyy = np.fft.rfft(Ryy_windowed) * dt
Pxy = np.fft.rfft(Rxy_windowed) * dt

f = np.fft.rfftfreq(len(x), d=dt)
# 计算相干性：|Pxy|² / (Pxx * Pyy)
Coh1 = np.abs(Pxy)**2 / (np.abs(Pxx) * np.abs(Pyy))

ax[0].plot(f[10:], Coh1[10:])
ax[0].set_title('方法1：相关+DFT计算相干性')
ax[0].set_xlabel('频率 [Hz]')
ax[0].set_ylabel('相干性')
ax[0].set_ylim(0, 1.1)

# --- 修正后的方法2 ---
# 使用周期图方法匹配方法1，无窗、无重叠、全长度分段
f, Coh2 = coherence(x, y, fs=fs, window='boxcar', noverlap=0, nperseg=len(x))

ax[1].plot(f, Coh2)
ax[1].set_title('方法2：scipy.signal.coherence计算')
ax[1].set_xlabel('频率 [Hz]')
ax[1].set_ylabel('相干性')
ax[1].set_ylim(0, 1.1)

plt.tight_layout()
plt.show()

修正效果 #

修正后两种方法结果完全一致，且在250Hz和400Hz处相干性均为1，符合理论预期。

内容的提问来源于stack exchange，提问作者fishbacp