# Should scipy.signal.coherence be 1 for single input and output signals?

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?

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


a closed form expression for the pdf of the estimate $\mid \hat{\gamma^2}\mid$, the $p(\mid \hat{\gamma^2}\mid |\; N,\mid \gamma^2 \mid)$, is given and the confidence pdf for $N=1$ is flat on $[0,1]$
• Thank you very much Stanley for your help! The source you pointed me at is indeed answering most of my questions! The source also states that "[...] we first partition each time-limited realization [time series] into N segments [...]" (p.237). Therefore, it seems to me taking one single time series for both, input and output, is actually enough. Later the source states that Welch's method (used by scipy.signal.coherence) fits into this approach (p.238). I'd get a coherence = 1 if I didn't slice my time series into N>1 parts. – Bax Menker May 9 '18 at 19:10