python - クロススペクトル密度を使用して、関連する 2 つの信号の位相シフトを計算する方法

Question

2 つの信号があり、一方が他方に応答していると予想されますが、特定の位相シフトがあります。

ここで、コヒーレンスまたは正規化されたクロス スペクトル密度を計算して、入力と出力の間に因果関係があるかどうかを推定し、このコヒーレンスがどの周波数で現れるかを調べたいと思います。

たとえば、周波数 10 でコヒーレンスが高いと思われる この画像 (ここから) を参照してください。

これで、相互相関を使用して 2 つの信号の位相シフトを計算できることがわかりましたが、(周波数 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()

.
.
編集：

価値があるので、答えを追加しました。おそらく正しいか、間違っているかもしれません。わからない..

Answer 1

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

相関する 2 つの信号のパワー スペクトル密度:

2 つの信号のコヒーレンスと位相 (10 Hz に拡大):

そして、クロス スペクトル密度の実部と虚部 (!) は次のとおりです。

Answer 2

不確実性を含むクロススペクトル分析を説明するJupyter Notebookを用意しました。

スクリーンショット: