Coherence values do not agree when trying to calculate using different methods

Question

I'm trying to see if I can obtain identical coherence values for two signals when computing them two different ways:

Method 1: I first compute the cross-correlation and the two auto-correlations. Then I take the DFT of the results to obtain the cross-spectral density and the power spectral densities, from which the coherence can be computed.

Method 2: I simply apply the coherence function from scipy.signal.

import numpy as np
from matplotlib import pyplot as plt
from scipy.signal import correlate,coherence

t  = np.arange(0, 1, 0.001);
dt = t[1] - t[0]
fs=1/dt

f0=250
f1=400

x=np.cos(2*np.pi*f0*t)+np.cos(2*np.pi*f1*t)
y=np.cos(2*np.pi*f0*(t-.035))+np.cos(2*np.pi*f1*(t-.05))

fig, ax = plt.subplots(2,sharex=True)

# Coherence using method 1:
Rxx = correlate(x,x)
Pxx= abs(np.fft.rfft(Rxx))

Ryy = correlate(y,y)
Pyy = abs(np.fft.rfft(Ryy))

Rxy = correlate(x,y)
Pxy = abs(np.fft.rfft(Rxy))

f = np.fft.rfftfreq(len(Rxy))*fs
Coh1=np.divide(Pxy**2,np.multiply(Pxx,Pyy))

#Start at nonzero index, e.g. 50, since Pxx[0]=Pyy[0] where Coh[0] would be undefined
ax[0].plot(f[50:], Coh1[50:]) 
both equal zero
ax[0].set_xlabel('frequency [Hz]')
ax[0].set_ylabel('Coherence')

# Coherence using method 2:
f,Coh2=coherence(x,y)

ax[1].plot(f*fs,Coh2)
ax[1].set_xlabel('frequency [Hz]')
ax[1].set_ylabel('Coherence')

plt.show()

I obtain the graphs shown below.

The results don't make sense to me. Not only are they different, but neither produces coherence values equal to one at frequencies 250 and 400. The second method should yield an undefined coherence except at 250 and 400.

VMMF · Answer 1 · 2023-09-01T20:48:35.620

For your Method 1, instead of the Fourier transform of the correlation I would directly apply the cross power spectral density using scipy.signal.csd() which internally uses the Welch method to estimate the power spectral density, I would make sure I apply a Hanning Window to reduce spectral leakage, 50 % overlap to reduce magnitude square coherence estimator variance and I would try to make the time signals long enough so that I can have at least 32 averages (256 ideal), etc.

Not only are they different, but neither produces coherence values equal to one at frequencies 250 and 400. The second method should yield an undefined coherence except at 250 and 400.

I'm not sure I understand these assumptions.

Coherence or more specifically magnitude squared coherence (MSC) has multiple interpretations in different fields.

In the mechanical engineering world, a coherence value of 1 indicates that the waveforms’ phase difference is consistent for one or multiple samples, and a value of 0 means the difference in their phase has changed. Coherence indicates if the phase difference changes for the defined frequency

More generally coherence indicates if there is a linear relationship between 2 signals. If the MSC is 1 only at certain frequencies, then the system is linear time invariant at those frequencies. However be aware that numerical errors (be aware of the Matlab/Python eps value) and bias/variance of the MSC estimator may affect the MSC calculation

I would refer you to the following theoretical sources about coherence:

See below 2 alternatives (from pyplot and scipy) to your Method 2 and the plot I'm getting. Feel free to vary the duration, NFFT , window type and see how this varies the MSC estimator without changing the input signals. Be aware that the actual MSC doesn't change but the estimator we are getting of it does.

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

Fs = 1000
duration = 2 #seconds
time=np.arange(0,duration,1/Fs)

f0=250
f1=400

signal1 = np.cos(2*np.pi*f0*time)+np.cos(2*np.pi*f1*time)
 
plt.figure(0)
manager = plt.get_current_fig_manager()
manager.window.showMaximized()
plt.subplot(4, 1, 1)
plt.plot(signal1)
plt.title("Signal 1")
 
# signal 2
signal2 = np.cos(2*np.pi*f0*(time-.035))+np.cos(2*np.pi*f1*(time-.05))

 
plt.subplot(4, 1, 2) 
plt.plot(signal2)
plt.title("Signal 2")

NFFT = 128
avgWindow = np.kaiser(NFFT, 6) # 0 Rectangular, 5 Hamming, 6 Hanning, 8.6 Blackman

plt.subplot(4, 1, 3)

coh=plt.cohere(signal1,signal2,NFFT=NFFT,noverlap=NFFT//2,Fs=Fs,Fc=0,detrend='mean',window=avgWindow)

plt.subplot(4, 1, 4)
f, Cxy = sp_signal.coherence(signal1,signal2,fs=Fs,window=avgWindow,nfft=NFFT,noverlap=NFFT//2,detrend='constant') #Cxy = abs(Pxy)**2/(Pxx*Pyy)
plt.stem(f,Cxy)
plt.xticks(np.arange(0, f[-1]+1, f[-1]//25))
plt.xlabel('Frequency [Hz]')
plt.ylabel('Magnitude squared coherence')
 
# plot the coherence graph
plt.show()

Coherence values do not agree when trying to calculate using different methods

1 Answers1