2
$\begingroup$

Suppose I have 2 seconds of data sampled at 1000 Hz., and two sinusoids differing only in amplitude and phase shift:

import numpy as np
from matplotlib import pyplot as plt
from scipy.signal import csd,periodogram

t = np.arange(0, 2, 0.001);
x = np.sin(2*np.pi*200*t)
y = 3*np.sin(2*np.pi*200*(t-np.pi/2))
fig,ax = plt.subplots()
ax.cohere(x, y)

The resulting plot of coherence is shown below.

enter image description here

Q1: Am I correct that the coherence is identically equal to one throughout and that the frequency axis is really going from 0 to 500 Hz, the Nyquist frequency?

Q2: For simple sinusoids, does a coherence identically equal to one mean they essentially differ in amplitude and have identical frequencies and a fixed phase difference?

Now I tried to compute the coherence using the cross-spectral density of x and y, along with their respective power spectral densities:

f,Pxx=periodogram(x)
f,Pyy=periodogram(y)
f,Pxy=csd(x,y)

Coh=np.zeros(len(f))
for i in range(len(f)):
    Coh[i]=abs(Pxy[i])**2/(Pxx[i]*Pyy[i])

fig,ax = plt.subplots()
ax.plot(f,Coh)

This resulted in the graph below.

Q3: Shouldn't the graphs be identical? What am I missing?

enter image description here

$\endgroup$

3 Answers 3

3
$\begingroup$

I suspect it's something different between csd and periodogram. If, instead, I use csd for all calculations:

f,Pxx=csd(x , x )
f,Pyy=csd(y , y)
f,Pxy=csd(x , y)

Coh=np.zeros(len(f))
for i in range(len(f)):
    Coh[i]=abs(Pxy[i])**2/(abs(Pxx[i]*Pyy[i]))

fig,ax = plt.subplots()
ax.plot(f,Coh)

then I get the figure below, which seems much closer to what you're expecting to see.

Coherence calculated using csd only

periodogram probably uses some averaging which doesn't help in the coherence calculation. By default, it looks like csd uses no averaging, but when called to calculate periodogram it uses averaging.

$\endgroup$
1
  • 1
    $\begingroup$ Ah! That's very helpful. Thanks! $\endgroup$
    – fishbacp
    Apr 1, 2022 at 12:31
2
$\begingroup$

What you are trying to do here is ill defined. At most frequencies $P_{xx}$, $P_{yy}$ and $P_{xy}$ are essentially zero, so the coherence is mostly $\frac{0}{0}$ which is undefined.

What happens in the actual code is highly dependent on how exactly the estimators are implemented (spectral leakage, windowing, hop size, scaling, etc) and good old fashioned numerical noise.

In your second graph, your Y-axis scale appears to be $10^{56}$ which is a strong indicator of some serious numerical problem.

$\endgroup$
0
$\begingroup$

As Peter indicated there is a difference between csd and periodogram. However, from his answer it was not clear when and why are averages used. So I will expand on that.

Note that periodogram uses Welch with a hardcoded noverlap=0 and a default window='boxcar'. As can be seen in the last portion of the function source code

return welch(x, fs=fs, window=window, nperseg=nperseg, noverlap=0,
             nfft=nfft, detrend=detrend, return_onesided=return_onesided,
             scaling=scaling, axis=axis)

In turn, Welch function internally uses csd of the input with itself, window='hann' and noverlap=None (which translates to a 50% segment overlap) as can be seen in the last portion of the function source code

freqs, Pxx = csd(x, x, fs=fs, window=window, nperseg=nperseg,
                 noverlap=noverlap, nfft=nfft, detrend=detrend,
                 return_onesided=return_onesided, scaling=scaling,
                 axis=axis, average=average)

Finally CSD uses _spectral_helper internal function.

It is easy to see that both CSD and periodogram end up calling this same function. The difference is that they call it with different input arguments.

Besides the differences already explained above, they key difference is that CSD calls _spectral_helper with nperseg=None, which means _spectral_helper defaults to nperseg=256 samples, computing a few averages because your input has 2000 samples. On the other hand periodogram calls _spectral_helper with nperseg=2000 and no averages are computed.

In other words if you compute the code below your are doing exactly the same as if you were calling periodogram

import numpy as np
from matplotlib import pyplot as plt
from scipy.signal import csd,periodogram

t = np.arange(0, 2, 0.001)
x = np.sin(2*np.pi*200*t)
y = 3*np.sin(2*np.pi*200*(t-np.pi/2))
fig,ax = plt.subplots()
ax.cohere(x, y)

f,Pxx=csd(x , x, window='boxcar',noverlap=0,nperseg=len(x))
f,Pyy=csd(y , y, window='boxcar',noverlap=0,nperseg=len(y))
# f,Pxx=periodogram(x)
# f,Pyy=periodogram(y)
f,Pxy=csd(x , y, window='boxcar',noverlap=0)

Coh=np.zeros(len(f))
for i in range(len(f)):
    Coh[i]=abs(Pxy[i])**2/(abs(Pxx[i]*Pyy[i]))

fig,ax = plt.subplots()
ax.plot(f,Coh)
plt.show()

But there is no way to call periodogram to obtain the same results as in the default CSD because the inputs required to do this were hardcoded and you don't have the option to change them

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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