scipy 尝试使用不同方法计算时，相干值不一致

我试着看看，当用两种不同的方法计算两个信号时，我是否可以获得相同的相干值：

**方法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:]) 
both equal zero
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()

我得到了下面所示的图表。
结果对我来说毫无意义。它们不仅不同，而且在频率250和400处都不产生等于1的相干值。第二种方法应产生除250和400以外的未定义相干性。

对于您的方法1，我将使用scipy.signal.csd()直接应用互功率谱密度，而不是相关性的傅里叶变换，该方法在内部使用Welch方法来估计功率谱密度，我将确保应用汉宁窗口来减少频谱泄漏，50%的重叠，以减少幅度平方相干估计方差，我会尝试使时间信号足够长，以便我可以有至少32个平均值（256个理想值），等等。
它们不仅不同，而且在频率250和400处都不产生等于1的相干值。第二种方法应产生除250和400以外的未定义相干性。
我不确定我是否理解这些假设。
相干性或更具体地说幅度平方相干性（MSC）在不同领域有多种解释。
在机械工程领域，相干值为1表示波形的相位差对于一个或多个样本是一致的，值为0表示它们的相位差已经改变。相干性指示相位差是否针对所定义的频率而改变
更一般地，相干性指示2个信号之间是否存在线性关系。如果MSC仅在某些频率下为1，则系统在这些频率下是线性时不变的。但是，请注意，MSC估计器的数值误差（注意Matlab/Python eps值）和偏差/方差可能会影响MSC计算
我想请你参考以下关于连贯性的理论来源：

1987 Coherence and Time Delay Estimation, G. Clifford Carter
Coherence Graph
Coherence Mathematics
Coherence
2010 Random Data Analysis and Measurement Procedures 4th Ed, Julius S. Bendat, Allan G. Piersol
请参阅下面的2个替代方案（从pyplot和scipy）到你的方法2和我得到的图。您可以随意改变持续时间、NFFT、窗口类型，看看这如何在不改变输入信号的情况下改变MSC估计器。请注意，实际的MSC不会改变，但我们从中获得的估计值会改变。

import numpy as np
import matplotlib.pyplot as plt
from scipy import signal as sp_signal

Fs = 1000
duration = 2 #seconds
time=np.arange(0,duration,1/Fs)

f0=250
f1=400

signal1 = np.cos(2*np.pi*f0*time)+np.cos(2*np.pi*f1*time)
 
plt.figure(0)
manager = plt.get_current_fig_manager()
manager.window.showMaximized()
plt.subplot(4, 1, 1)
plt.plot(signal1)
plt.title("Signal 1")
 
# signal 2
signal2 = np.cos(2*np.pi*f0*(time-.035))+np.cos(2*np.pi*f1*(time-.05))

 
plt.subplot(4, 1, 2) 
plt.plot(signal2)
plt.title("Signal 2")

NFFT = 128
avgWindow = np.kaiser(NFFT, 6) # 0 Rectangular, 5 Hamming, 6 Hanning, 8.6 Blackman

plt.subplot(4, 1, 3)

coh=plt.cohere(signal1,signal2,NFFT=NFFT,noverlap=NFFT//2,Fs=Fs,Fc=0,detrend='mean',window=avgWindow)

plt.subplot(4, 1, 4)
f, Cxy = sp_signal.coherence(signal1,signal2,fs=Fs,window=avgWindow,nfft=NFFT,noverlap=NFFT//2,detrend='constant') #Cxy = abs(Pxy)**2/(Pxx*Pyy)
plt.stem(f,Cxy)
plt.xticks(np.arange(0, f[-1]+1, f[-1]//25))
plt.xlabel('Frequency [Hz]')
plt.ylabel('Magnitude squared coherence')
 
# plot the coherence graph
plt.show()

scipy 尝试使用不同方法计算时，相干值不一致

1条答案

相关问题

热门标签

最新问答