0
$\begingroup$

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: Coherence Result with artificial data

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()
$\endgroup$
1
$\begingroup$

Your colleague is correct.

In

Carter, G. Clifford. "Coherence and time delay estimation." Proceedings of the IEEE 75.2 (1987): 236-255.

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]$

enter image description here

$\endgroup$
  • 1
    $\begingroup$ 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. $\endgroup$ – Bax Menker May 9 '18 at 19:10
  • $\begingroup$ So, I guess the coherence calculated in my example is correct, right? $\endgroup$ – Bax Menker May 9 '18 at 19:11
  • $\begingroup$ I would say No, taking the noiseless coherence of something with respect to itself results in one isn't what coherence is about. It is a rather useless result, enough for what, circular reasoning? If you know that they are the same, why use the coherence? N=1 is explicitly ruled out. Period. $\endgroup$ – Stanley Pawlukiewicz May 9 '18 at 19:53
  • $\begingroup$ I think I did not explain myself clearly enough. Ok, I understand N=1 is ruled out. But what does N refer to? Reading the source carfully it seems to me N is the number of partitions. So if I have one input signal and one noisy output signal, as in my example code, I cut those in N partitions and compute the N input/output pairs' cross spectral densities. Those are then used to calculate my overall coherence of the original single input and noisy output signals or am I still wrong? Thanks for your patience :-) $\endgroup$ – Bax Menker May 9 '18 at 20:46
  • $\begingroup$ N independent. There are papers for overlap as well and 50% is recommended in a number of papers. The DFT should also be “large” $\endgroup$ – Stanley Pawlukiewicz May 9 '18 at 21:21

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.