Coherence Part 2 (Two noise signals, again)

# Load modules we'll need.

from scipy.io import loadmat
import matplotlib.pyplot as plt
import numpy as np
from scipy.signal import spectrogram

Make two noise signals, with multiple trials

N = 1000;
K = 100;
dt= 0.001;
T = N*dt;
x = np.random.randn(K,N)
y = np.random.randn(K,N)
t = np.arange(0,N)*dt

plt.plot(t,x[0,:])
plt.plot(t,y[0,:])
plt.xlabel('Time [s]');

Visualize the data across all trials

plt.imshow(x,                             # ... and show the image,
           extent=[min(t), max(t), K, 1],  # ... with meaningful axes,
           aspect='auto')                  # ... and a nice aspect ratio
plt.xlabel('Time [s]')
plt.ylabel('Trial #');
plt.title('All trials from E1');

Compute the cross-covariance, averaged across trials

cc_xy = "SOMETHING"                          # Compute cc for each trial, 
cc_xy = np.mean(cc_xy,0)                     # ... average over trials,
lags = np.arange(-N + 1, N)                  # ... create a lag axis,
plt.plot(lags * dt, cc_xy)                   # ... and plot the result.
plt.xlabel('Lag [s]')
plt.ylabel('Trial averaged cross-covariance');
plt.ylim([-0.1, 1]);

Compute the coherence

# Fourier transforms.
Xf = "SOMETHING"  # Compute Fourier transform of x for each trial
Yf = "SOMETHING"  # Compute Fourier transform of y for each trial

# Auto- and cross-spectra.
Sxx = "SOMETHING"  # Spectrum of x trials
Syy = "SOMETHING"  # ... and y trials
Sxy = "SOMETHING"  # ... and the cross spectrum

# Trial average.
Sxx = np.mean(Sxx,0)
Syy = np.mean(Syy,0)
Sxy = np.mean(Sxy,0)

# Calculate coherence.
cohr = "SOMETHING"

f = np.fft.fftfreq(N, dt)               # Define a frequency axis.
plt.plot(f, cohr.real)                  # Plot the coherence.
plt.ylim([0, 1.1])                      # ... with y-axis scaled,
plt.xlabel('Frequency [Hz]')            # ... and with axes labeled.
plt.ylabel('Coherence')
plt.title('Trial averaged coherence between two electrodes');