2
$\begingroup$

My question is somewhat similar to one I posted previously

I understand how coherence is defined in terms of spectral densities, but I would be elated if I could explain the basic concept of coherence, in terms of orthogonality, to my linear algebra students who have no signal processing background.

I have two signals, x and y, sampled at 1000 Hz for one second each. Below I enter the signals and plot the coherence $Coh_{xy}$.

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

dt = .001
t = np.arange(0, 1, dt)
N=len(t)
fs=int(1/dt)

noise=np.random.uniform(low=-.0000125, high=.0000125, size=(N,))

f0=300
f1=450

x = np.cos(2 * np.pi * f0 * t) + np.cos(2 * np.pi * f1 * t)
y = np.cos(2 * np.pi * f0 * (t**1.005))+np.cos(2*np.pi*f1*(t+5*noise))

fig, ax = plt.subplots()

f,Coh_xy=coherence(x,y,fs=fs,nfft=len(x))

ax.plot(f,Coh_xy)
ax.set_xlabel('frequency')
ax.set_ylabel('coherence Cxy')
ax.set_ylim(0, 1.2)
ax.grid(True)
plt.show()

enter image description here

print(abs(Coh_xy[f0]))
print(abs(Coh_xy[f1]))

0.7405039912715564
0.9999417256272132

Given how the signals are constructed, the peak coherence values make sense to me.

Denoting the signals as $\vec{x}$ and $\vec{y}$ in $\mathbb{R}^{1000}$, I can compute their DFTs using the 1000-by-1000 DFT Matrix, $W$. (I use the convention that columns are normalized.)

The frequency resolution is 1 Hz, so the frequencies $f_0=300$ and $f_1=450$ correspond to bins 300 and 450, respectively. To determine the spectral component of a signal at $f_0=300$, I could project the signal vector onto columns 300 and 700 (=1000-300) of $W$ and add the two results. For example, if $\vec{w}_0$ and $\vec{w}_1$ denote columns 300 and 700 of $W$, respectively, then the spectral component of $\vec{x}$ at frequency bin 300 is

\begin{equation*} \vec{x}_f=(\vec{x} \bullet \vec{w}_0)\vec{w}_0+(\vec{x} \bullet \vec{w}_1)\vec{w}_1. \end{equation*}

In a similar manner, I could compute the spectral components of $\vec{y}$. Below, I create and plot the spectral components of each signal at $f_0=300$ and $f_1=450$. In addition, I create scatter plots of $\vec{x}_f$ versus $\vec{y}_f$ at $f_0=300$ Hz and at $f_1=450$. For each scatter plot, I also compute the correlation coefficient.

def DFT_matrix(N):
    i, j = np.meshgrid(np.arange(N), np.arange(N))
    omega = np.exp( - 2 * np.pi * 1J / N )
    W = np.power( omega, i * j )  # / np.sqrt(N)
    return W

W=DFT_matrix(N)
F=f0
w0=W[:,F]/np.linalg.norm(W[:,F])
w1=W[:,fs-F]/np.linalg.norm(W[:,F])
x_f=np.real(np.dot(x,w0)*w0+np.dot(x,w1)*w1)
y_f=np.real(np.dot(y,w0)*w0+np.dot(y,w1)*w1)
print('correlation coefficient of x_f and y_f at f0 = 300: '+str(np.corrcoef(x_f,y_f)
[0,1]))

fig, ax = plt.subplots(3)
n_min=10
n_max=50
ax[0].plot(t[n_min:n_max],x_f[n_min:n_max],color='b')
ax[0].grid(color='grey', linestyle='-', linewidth=2)
ax[1].plot(t[n_min:n_max],y_f[n_min:n_max],color='r')
ax[1].grid(color='grey', linestyle='-', linewidth=2)
ax[2].scatter(x_f,y_f,color='b')
ax[2].set_xlabel('x_f')
ax[2].set_ylabel('y_f',rotation=0)
plt.show()

enter image description here

F=f1
w0=W[:,F]/np.linalg.norm(W[:,F])
w1=W[:,fs-F]/np.linalg.norm(W[:,F])
x_f=np.real(np.dot(x,w0)*w0+np.dot(x,w1)*w1)
y_f=np.real(np.dot(y,w0)*w0+np.dot(y,w1)*w1)
print('correlation coefficient of x_f and y_f at f1 = 450: '+str(np.corrcoef(x_f,y_f)
[0,1]))

fig, ax = plt.subplots(3)
n_min=10
n_max=50
ax[0].plot(t[n_min:n_max],x_f[n_min:n_max],color='b')
ax[0].grid(color='grey', linestyle='-', linewidth=2)
ax[1].plot(t[n_min:n_max],y_f[n_min:n_max],color='r')
ax[1].grid(color='grey', linestyle='-', linewidth=2)
ax[2].scatter(x_f,y_f,color='b')
ax[2].set_xlabel('x_f')
ax[2].set_ylabel('y_f',rotation=0)
plt.show()

enter image description here

This leads me to my question.

Is it reasonable to say the following: When the square of the correlation coefficient between two spectral components, plotted as ordered pairs in $\mathbb{R}^2$, is close to one, the spectral components are close to being linear functions of one another (constant multiples even). When the spectral components are close to being linear functions of one another, the coherence between the two original signals at the corresponding frequency will be close to one.

My question is based upon a point raised in an editorial found in International Journal of Psychophysiology, 57 (2005) 83–85:

Coherence of a signal pair at a given frequency can be conceptualized as 1) passing each signal through a narrow band filter tuned to the given frequency, 2) computing the Pearson-product moment correlation coefficient of the filtered signals, and 3) squaring the resulting correlation coefficient. The amplitudes of coherence functions range from a minimum of zero to a maximum of +1.0. In formal terms, the coherence function is defined as the squared magnitude of the cross-spectral density of the signal pair divided by the product of the power spectra of the two signals.

$\endgroup$

0

Your Answer

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

Browse other questions tagged or ask your own question.