如何使用交叉谱密度计算两个相关信号的相移 [英] How to use the cross-spectral density to calculate the phase shift of two related signals

问题描述

我有两个信号，我希望一个信号响应另一个信号，但有一定的相移.

现在，我想计算相干或归一化的交叉频谱密度，以估计输入和输出之间是否存在因果关系，以找出该相干出现在哪些频率上.

例如参见这张图片(来自这里)，它似乎在频率10:

现在我知道我可以使用互相关来计算两个信号的相移，但是我如何使用相干性(频率为 10)来计算相移?

图片代码:

<代码>"""计算两个信号的相干性"""将numpy导入为np导入matplotlib.pyplot作为plt#在子图之间留出一些额外的空间plt.subplots_adjust(wspace = 0.5)nfft = 256dt = 0.01t = np.arange(0, 30, dt)nse1 = np.random.randn(len(t))#白噪声1nse2 = np.random.randn(len(t)) # 白噪声 2r = np.exp(-t/0.05)cnse1 = np.convolve(nse1, r, mode='same')*dt # 彩色噪声 1cnse2 = np.convolve(nse2，r，mode ='same')* dt#有色噪声2# 具有相干部分和随机部分的两个信号s1 = 0.01*np.sin(2*np.pi*10*t) + cnse1s2 = 0.01*np.sin(2*np.pi*10*t) + cnse2plt.subplot(211)plt.plot(t, s1, 'b-', t, s2, 'g-')plt.xlim(0,5)plt.xlabel('时间')plt.ylabel('s1和s2')plt.grid(真)plt.subplot(212)cxy, f = plt.cohere(s1, s2, nfft, 1./dt)plt.ylabel('一致性')plt.show()

.
.

对于它的价值，我已经添加了一个答案，也许是对的，也许是错误的.我不确定..

解决方案

让我尝试回答我自己的问题，也许有一天有可能对其他人有用，或者可能是(新)讨论的起点:

首先计算两个信号的功率谱密度，

 子图(121)psd(s1，nfft，1/dt)plt.title('signal1')子图(122)psd(s2, nfft, 1/dt)plt.title('signal2')plt.tight_layout()表演()

导致:

第二次计算互谱密度，即互相关函数的傅立叶变换:

  csdxy，fcsd = plt.csd(s1，s2，nfft，1./dt)plt.ylabel('CSD (db)')plt.title('信号1和2之间的交叉频谱密度')plt.tight_layout()表演()

哪个给:

通过使用互谱密度，我们可以计算相位，并且可以计算相干性(相干性会破坏相位).现在，我们可以将相干性和超过95％置信水平的峰结合起来

# 连贯性cxy，fcoh = cohere(s1，s2，nfft，1./dt)# 计算 95% 的置信水平edof = (len(s1)/(nfft/2)) * cxy.mean() # 等效自由度:(length(timeseries)/windowhalfwidth)*mean_coherencegamma95 = 1.-(0.05)**(1./(edof-1.))conf95 = np.where(cxy>gamma95)打印'gamma95'，gamma95，'edof'，edof# 绘制孪生图图，ax1 = plt.subplots()# 在 ax1 上绘制相干性ax1.plot(fcoh, cxy, 'b-')ax1.set_xlabel('频率(hr-1)')ax1.set_ylim([0,1])#使y轴标签和刻度标签与线条颜色匹配.ax1.set_ylabel('Coherence', color='b')对于 ax1.get_yticklabels() 中的 tl:tl.set_color('b')#在ax2上绘制相位ax2 = ax1.twinx()ax2.plot(fcoh [conf95]，phase [conf95]，'r.')ax2.set_ylabel('相位(度)'，color ='r')ax2.set_ylim([-200,200])ax2.set_yticklabels([-180,-135,-90,-45,0,45,90,135,180])对于 ax2.get_yticklabels() 中的 tl:tl.set_color('r')ax1.grid(真)#ax2.grid(真实)fig.suptitle('信号 1 和 2 之间的相干性和相位 (>95%)', fontsize='12')plt.show()

结果:

总而言之:在10分钟的时间段内，最一致的峰的相位约为1度(s1领先s2)(假设 dt 是一分钟的测量值)-> (10**-1)/dt

但是信号处理专家可能会纠正我，因为我60％的人确定我做对了

I've two signals, from which I expect that one is responding on the other, but with a certain phase shift.

Now I would like to calculate the coherence or the normalized cross spectral density to estimate if there is any causality between the input and output to find out on which frequencies this coherence appear.

See for example this image (from here) which seems to have high coherence at the frequency 10:

Now I know that I can calculate the phase shift of two signals using the cross-correlation, but how can I use the coherence (of frequency 10) to calculate the phase shift?

Code for image:

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

.
.
EDIT:

For what it's worth, I've add an answer, maybe it's right, maybe it's wrong. I'm not sure..

解决方案

Let me try to answer my own question and maybe one day it might be useful to others or function as a starting point for a (new) discussion:

Firstly calculate the power spectral densities of both the signals,

subplot(121)
psd(s1, nfft, 1/dt)
plt.title('signal1')

subplot(122)
psd(s2, nfft, 1/dt)
plt.title('signal2')

plt.tight_layout()
show()

resulting in:

Secondly calculate the cross-spectral density, which is Fourier transform of the cross-correlation function:

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

Which gives:

Than using the cross-spectral density we can calculate the phase and we can calculate the coherence (which will destroy the phase). Now we can combine the coherence and the peaks that rise above the 95% confidence level

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

result in:

To sum up: the phase of the most coherent peak is ~1 degrees (s1 leads s2) at a 10 min period (assuming dt is a minute measurement) -> (10**-1)/dt

But a specialist signal processing might correct me, because I'm like 60% sure if I've done it right

这篇关于如何使用交叉谱密度计算两个相关信号的相移的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！