I am trying to calculate the coherence between input and output signals. I thought I could work with a single input and a single output time series and calculate the coherence $\gamma^2$ between them. I did so with the scipy.signal.coherence function (See code and results below).
Now, a colleague challanged my procedure, saying that the coherence needs to be calculated by averaging several input and output signals, since otherwise the coherence will simply be $\gamma^2=1$. I wanted to prove otherwise, after all I used a scipy function and I regard those as trustworthy.
So I started out with the general definition of the coherence as [1][2] [3]:
$$
\gamma^2 = \frac{|\mathbb{E}[X(-j\omega)Y(j\omega)]|^2}{\mathbb{E}[|X(j\omega)|^2]~~\mathbb{E}[|Y(j\omega)|^2]}, \quad \hat{\gamma}^2 = \frac{|\overline{S}_{XY}|^2}{\overline{S}_{XX}\overline{S}_{YY}},
$$
where $\overline{\cdot}$ indicates a mean value, $|\cdot|$ the magnitude of a complex number, ${S}_{XX}$ and ${S}_{YY}$ are the energy spectral density for $X$ and $Y$ respectively and ${S}_{XY}$ is the cross spectral denisty between $X$ and $Y$. And $X(j\omega)$ being the Fourier Transform of a time series $X(t)$ and $X(-j\omega)$ being its complex conjugate. Let's define those as:
$$ X(j\omega) = x_1+ix_2, \quad X(-j\omega) = x_1 -ix_2, \quad Y(j\omega) = y_1+iy_2. $$
Now, assuming I do not have several signals, but only a single input signal $X$ and single output signal $Y$, I drop the mean $\overline{\cdot}$ and end up with the following: $$ \gamma^2 = \frac{|{S}_{XY}|^2}{{S}_{XX}{S}_{YY}} = \frac{|X(-j\omega)Y(j\omega)|^2}{|X(j\omega)|^2|Y(j\omega)|^2} = \frac{\left(|X(-j\omega)||Y(j\omega)|\right)^2}{|X(j\omega)|^2|Y(j\omega)|^2}= \frac{\left(\sqrt{x^2_1+x^2_2}\sqrt{y^2_1+y^2_2}\right)^2}{\left(x^2_1+x^2_2\right)\left(y^2_1+y^2_2\right)}=1. $$
In case I did not make a mistake, my scipy.signal.coherence function should give me a coherence = 1. But as shown below, it does not. Instead the coherence is close to zero everywhere, except at the signals' natural frequencies. Please note that I get $\gamma^2=1$ only if I chose $X=Y$.
I must be missing something here but I just don't see it. Any ideas?
Many thanks in advance! Max
Result:
Code:
from scipy import signal
import matplotlib.pyplot as plt
import numpy as np
from scipy.fftpack import fft
fs = 5e3 # sampling frequency
N = 10000 # number of observations
T = N/fs # max time T
amp = 20 # sine wave amplitude input
amp2 = 2*amp # sine wave amplitude output
freq = 500.0 # sine wave frequency input
freq2 = 500 # sine wave frequency output
time = np.arange(N) / fs # time vector t_0 up to T
# create data
x = amp*np.sin(2*np.pi*freq*time) # input
y = amp2*np.sin(2*np.pi*freq2*time)+np.random.normal(scale = 10, size = len(x)) + 100 # output
fig, axes = plt.subplots(1,3, figsize=(15,5))
# plot time signals
axes[0].plot(x[:400], label = 'Input', alpha = 0.7)
axes[0].plot(y[:400], label = 'Output', alpha = 0.7)
axes[0].legend()
axes[0].grid(linestyle='--')
axes[0].set_xlabel('Time')
# plot fft input
xfft = fft(x)
N = len(x) # number of sample points
T = 1 / fs # sample spacing
xfreq = np.linspace(0.0, 1.0/(2.0*T), N//2)
axes[1].plot(xfreq, 2.0/N * np.abs(xfft[0:N//2]), label = 'Input', alpha = 0.7)
# plot fft output
yfft = fft(y)
N = len(y) # number of sample points
T = 1 / fs # sample spacing
yfreq = np.linspace(0.0, 1.0/(2.0*T), N//2)
axes[1].plot(xfreq, 2.0/N * np.abs(yfft[0:N//2]), label ='Output', alpha = 0.7)
axes[1].legend()
axes[1].set_xlim(1,1000)
axes[1].grid(linestyle='--')
axes[1].set_xlabel('Frequency [Hz]')
# plot coherence:
f, Cxy = signal.coherence(x, y, fs)
axes[2].semilogy(f, Cxy, color='black')
axes[2].set_xlabel('Frequency [Hz]')
axes[2].set_ylabel('Coherence')
axes[2].grid(linestyle='--')
plt.show()
