How to use the cross-spectral density to calculate the phase shift of two related signals

Question

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..

Answer 1

I've prepared a Jupyter Notebook which explains the cross-spectral analysis including its uncertainty.

screenshot:

Answer 2

I am not sure, where the phase variable was calculated in the answer of @Mattijn.

You can calculate the phase shift from the angle between the real and the imaginary part of the cross-spectral density.

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

The power spectral density of the two signals to be correlated:

The coherence and the phase of the two signals (zoomed in to 10 Hz):

And here the real and imaginary(!) part of the cross spectral density: