如何使用互谱密度计算两个相关信号的相移答案

【问题标题】：How to use the cross-spectral density to calculate the phase shift of two related signals如何使用互谱密度计算两个相关信号的相移
【发布时间】：2014-03-06 00:12:00
【问题描述】：

我有两个信号，我预计其中一个会响应另一个信号，但会有一定的相移。

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

例如看这张图片（来自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()

.
.
编辑：

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

【问题讨论】：

问题是您不知道如何根据相干性计算相移，还是您的代码不起作用？前者不是这里的主题问题。
我的问题是我不知道如何从相干计算相移。如果您能给我一些提示或建议如何做到这一点（编程或数学），我们将不胜感激
我不是 Python 程序员，但我很确定你不会。根据matplotlib.org/api/mlab_api.html 的文档，plt.cohere 计算由信号模量归一化的交叉谱密度的模平方，这就是为什么它是纯实值的。如果两个相同频率的正弦信号之间存在相移，则信号之间的互相关将是振荡的，并具有与之相关的相移，并且该相移在经过傅里叶变换后将保留，但随后被破坏通过取模数。
所以假设它只是在 en.wikipedia.org/wiki/Spectral_density#Cross-spectral_density 处对公式的绝对值进行离散化版本，你不会找到相移。有什么理由不能只对两个单独的信号进行 FFT 并使用从中获取的相位差吗？
为了澄清起见，当您说“一个响应另一个，但有一定的相移”时，您的意思是您期望两个信号相同（噪声除外），但两者之间有时间延迟？

标签： python matplotlib signal-processing cross-correlation spectral-density

【解决方案1】：

我准备了一份Jupyter Notebook，其中解释了交叉光谱分析，包括其不确定性。

截图：

【讨论】：

@KylerBrown 它一定是 github 服务器上的问题。对我来说，链接再次起作用
亲爱的 Mattjin，我真的很喜欢您解决问题的方法。不过，我对您模型的一个方面表示怀疑：如何在您的模型中考虑采样频率单位？我无法理解您将理论上以秒为单位的信号（从您的示例合成）转换为以小时为单位的相干图的方式。例如，如果我要向您的模型添加一些每月采样的数据，我应该如何解释您的绘图的标签和图例？
嗨@Mattijn。感谢您提供此解决方案。但是，您能否就如何估计相干性而非相位的不确定性提供一些见解？

【解决方案2】：

我不确定，在@Mattijn 的答案中计算相位变量的位置。

您可以从实数和实数之间的角度计算相移交叉谱密度的虚部。

from matplotlib import mlab

# First create power sectral densities for normalization
(ps1, f) = mlab.psd(s1, Fs=1./dt, scale_by_freq=False)
(ps2, f) = mlab.psd(s2, Fs=1./dt, scale_by_freq=False)
plt.plot(f, ps1)
plt.plot(f, ps2)

# Then calculate cross spectral density
(csd, f) = mlab.csd(s1, s2, NFFT=256, Fs=1./dt,sides='default', scale_by_freq=False)
fig = plt.figure()
ax1 = fig.add_subplot(1, 2, 1)
# Normalize cross spectral absolute values by auto power spectral density
ax1.plot(f, np.absolute(csd)**2 / (ps1 * ps2))
ax2 = fig.add_subplot(1, 2, 2)
angle = np.angle(csd, deg=True)
angle[angle<-90] += 360
ax2.plot(f, angle)

# zoom in on frequency with maximum coherence
ax1.set_xlim(9, 11)
ax1.set_ylim(0, 1e-0)
ax1.set_title("Cross spectral density: Coherence")
ax2.set_xlim(9, 11)
ax2.set_ylim(0, 90)
ax2.set_title("Cross spectral density: Phase angle")

plt.show()

fig = plt.figure()
ax = plt.subplot(111)

ax.plot(f, np.real(csd), label='real')
ax.plot(f, np.imag(csd), label='imag')

ax.legend()
plt.show()

要关联的两个信号的功率谱密度：

两个信号的相干性和相位（放大到 10 Hz）：

这里是交叉谱密度的实部和虚部（！）：

【讨论】：

np.angle(csd) 返回一个全零向量，因为交叉谱是实值，我猜。
@fccoelho：交叉谱密度具有复数值。只需执行示例代码并打印 csd 的值： array([ 1.90723361e-04 +0.00000000e+00j, 1.98803902e-03 +1.75715712e-03j, ...
你是对的。我错误地查看了相干谱。

【解决方案3】：

让我尝试回答我自己的问题，也许有一天它可能对其他人有用或作为（新）讨论的起点：

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

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% 的把握我是否做得对

【讨论】：