作者热门文章
- c - 在位数组中找到第一个零
- linux - Unix 显示有关匹配两种模式之一的文件的信息
- 正则表达式替换多个文件
- linux - 隐藏来自 xtrace 的命令
我有两个信号,我希望一个信号对另一个信号作出响应,但具有一定的相移。
现在我想计算相干性或归一化交叉谱密度,以估计输入和输出之间是否存在任何因果关系,以找出这种相干性出现在哪些频率上。
例如,这张图片(来自 here)似乎在频率 10 处具有高相干性:
现在我知道我可以使用互相关计算两个信号的相移,但我如何使用相干性(频率为 10)来计算相移?
图片代码:
"""
Compute the coherence of two signals
"""
import numpy as np
import matplotlib.pyplot as plt
# make a little extra space between the subplots
plt.subplots_adjust(wspace=0.5)
nfft = 256
dt = 0.01
t = np.arange(0, 30, dt)
nse1 = np.random.randn(len(t)) # white noise 1
nse2 = np.random.randn(len(t)) # white noise 2
r = np.exp(-t/0.05)
cnse1 = np.convolve(nse1, r, mode='same')*dt # colored noise 1
cnse2 = np.convolve(nse2, r, mode='same')*dt # colored noise 2
# two signals with a coherent part and a random part
s1 = 0.01*np.sin(2*np.pi*10*t) + cnse1
s2 = 0.01*np.sin(2*np.pi*10*t) + cnse2
plt.subplot(211)
plt.plot(t, s1, 'b-', t, s2, 'g-')
plt.xlim(0,5)
plt.xlabel('time')
plt.ylabel('s1 and s2')
plt.grid(True)
plt.subplot(212)
cxy, f = plt.cohere(s1, s2, nfft, 1./dt)
plt.ylabel('coherence')
plt.show()
.
.
编辑:
对于它的值(value),我添加了一个答案,也许是对的,也许是错的。我不确定..
最佳答案
让我尝试回答我自己的问题,也许有一天它可能对其他人有用或作为(新)讨论的起点:
首先计算两个信号的功率谱密度,
subplot(121)
psd(s1, nfft, 1/dt)
plt.title('signal1')
subplot(122)
psd(s2, nfft, 1/dt)
plt.title('signal2')
plt.tight_layout()
show()
导致:
其次计算互谱密度,即互相关函数的傅里叶变换:
csdxy, fcsd = plt.csd(s1, s2, nfft, 1./dt)
plt.ylabel('CSD (db)')
plt.title('cross spectral density between signal 1 and 2')
plt.tight_layout()
show()
给出:
我们可以使用交叉谱密度计算相位和相干性(这会破坏相位)。现在我们可以结合一致性和高于 95% 置信水平的峰值
# coherence
cxy, fcoh = cohere(s1, s2, nfft, 1./dt)
# calculate 95% confidence level
edof = (len(s1)/(nfft/2)) * cxy.mean() # equivalent degrees of freedom: (length(timeseries)/windowhalfwidth)*mean_coherence
gamma95 = 1.-(0.05)**(1./(edof-1.))
conf95 = np.where(cxy>gamma95)
print 'gamma95',gamma95, 'edof',edof
# Plot twin plot
fig, ax1 = plt.subplots()
# plot on ax1 the coherence
ax1.plot(fcoh, cxy, 'b-')
ax1.set_xlabel('Frequency (hr-1)')
ax1.set_ylim([0,1])
# Make the y-axis label and tick labels match the line color.
ax1.set_ylabel('Coherence', color='b')
for tl in ax1.get_yticklabels():
tl.set_color('b')
# plot on ax2 the phase
ax2 = ax1.twinx()
ax2.plot(fcoh[conf95], phase[conf95], 'r.')
ax2.set_ylabel('Phase (degrees)', color='r')
ax2.set_ylim([-200,200])
ax2.set_yticklabels([-180,-135,-90,-45,0,45,90,135,180])
for tl in ax2.get_yticklabels():
tl.set_color('r')
ax1.grid(True)
#ax2.grid(True)
fig.suptitle('Coherence and phase (>95%) between signal 1 and 2', fontsize='12')
plt.show()
结果:
总结:最相干峰的相位在 10 分钟内约为 1 度(s1 超前 s2)(假设 dt
是分钟测量值)-> (10 **-1)/dt
但是专业的信号处理可能会纠正我的错误,因为我有 60% 的把握做对了
关于python - 如何使用互谱密度计算两个相关信号的相移,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/21647120/
我是一名优秀的程序员,十分优秀!