Taking a simple example of the coherence between two signals, each composed of two sinusoids of frequency 100 and 200 Hz, I expect the coherence to show two peaks in the spectrum located at precisely these two frequencies. The script below can be used to study the effect of adding noise An =0.1 or no noise An = 0. before taking the coherence. This produces in the case of added noise With Noise

and in the case of no noise enter image description here

It is clear that with no noise there are many false peaks in coherence.

Question: Can someone show me mathematically speaking why this is the case?

import numpy as np
import matplotlib.pyplot as plt
from scipy import signal

t0, t1 = 0, 2.
fs = 2000
dt = 1./fs
t = np.arange(t0,t1-dt,dt)

#f1, f2 = 50, 50
#x1 = np.sin(2*np.pi*f1*t)
#x2 = np.sin(2*np.pi*f2*t+2*np.pi*0.)

An = 0.1            # Noise amplitude
f1, f2 = 100, 200   # Signal frequencies
p1, p2 = 0.25, 0.5  # Phase shifts (times pi)

x1 = np.cos(2*np.pi*f1*t) + np.sin(2*np.pi*f2*t) + An*np.random.randn(len(t));
x2 = 0.5*np.cos(2*np.pi*f1*t-np.pi*p1) + 0.35*np.sin(2*np.pi*f2*t-np.pi*p2) + An*np.random.randn(len(t));

fcoh, coh = signal.coherence(x1, x2, fs, window='hamming', nperseg=400, noverlap=0)

fig2, ax2 = plt.subplots(1,1)
ax2.plot(fcoh, coh**2,"k-")
ax2.set_xlabel(r"$f$ [Hz]")

without noise, the spectrum levels outside your 2 tones are very low, just the quantization noise of the tones.

Adding noise introduces a flat spectral level across the full spectrum. The normalization is well behaved.

Coherence is usually used in moderate to low SNR situations where you can average a lot.

If you look in :

Carter GC. Coherence and time delay estimation. Proceedings of the IEEE. 1987 Feb;75(2):236-55.

there are equations for the pdf of the sample coherence as a function of independent samples. The number of independent samples required for reasonable confidence intervals can be high, on the order of several hundred independent samples at each frequency bin.

  • $\begingroup$ Thank you for your answer! I am not sure I completely understand your first sentence - you say the spectrum levels are low apart from the two tones, but that is clearly not the case as there are several peaks between 0.8-1.0. So I have probably misunderstood something. Could you please clarify? $\endgroup$ – Dipole Jul 23 '18 at 22:41
  • $\begingroup$ Those peaks are extraneous. You didn’t place any tones at those frequencies. There is nothing there to cohere. The peaks are meaningless. Get the Carter paper. $\endgroup$ – Stanley Pawlukiewicz Jul 24 '18 at 4:38

To elaborate on Stanley's answer:

Even with a minuscule amount of noise say An = 0.0000001 you will get a clean graph like your first graph. Python's implementation of the Welch method certainly uses finite-precision floating point numbers and quantization errors in the input and from the calculation steps such as windowing and Fast Fourier Transform (FFT) propagate to the results. From SciPy docs:

Cxy = abs(Pxy)**2/(Pxx*Pyy), where Pxx and Pyy are power spectral density estimates of X and Y, and Pxy is the cross spectral density estimate of X and Y.

The normalizing division in the equation means that the values of Pxx and Pyy can be arbitrarily close to 0 while still giving large Cxy.

  • $\begingroup$ Thank you for your answer! Yes I found that limiting behavior of An interesting! However I still dont understand the role the noise plays. Are you saying that the reason for the false peaks in the coherence spectrum is due to quantization errors, and that introducing noise, however small, will now act to smoothen out the effect of quantization errors? $\endgroup$ – Dipole Jul 23 '18 at 22:45
  • $\begingroup$ Yes. The added noise acts as dither, albeit unnecessarily large, which decorrelates the quantization noise in the output of each operation from its input. The noise does not repeat itself between windows of the Welch method. Windowed white noise in time domain is correlated noise in frequency domain. Gaussian white noise is Gaussian white noise in both domains. Multiplication (windowing) in time domain is convolution in the frequency domain. Floating point number quantization is difficult to analyze but you can find by testing how much dither is needed. It depends on the signals. $\endgroup$ – Olli Niemitalo Jul 25 '18 at 5:06

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.