Analyzing Rhythms Part 2b (Autocovariance)

# Load modules we'll need.
from scipy.io import loadmat
import matplotlib.pyplot as plt
import numpy as np
from scipy.signal import spectrogram

# Load the data.
data = loadmat('AC_Example.mat')  # Load the data,
# Data available here
# https://github.com/Mark-Kramer/BU-MA665-MA666/tree/master/Data
#
d = data['d'][0]                # ... from the first electrode.
t = data['t'][0]              # Load the time axis
N = np.size(d,0)              # Store number of observations.
dt = t[1]-t[0]                # Store sampling interval.

Load the data and look at it.

Q. Do you see rhythms?

Conclusions

# Code to compute the spectrum.

Compute the spectrum.

Q. What do you find? What rhythms are present in the data?

Conclusions

# Code to compute the autocovariance.

Compute the autocovariance.

Q. What do you find? Is it consistent with the spectrum?

Conclusions

Example: Spectrum = FT{Autocovariance}

Make a simple signal.

dt = 0.001
N  = 1000
t  = np.arange(0,N)*dt
d  = np.sin(2*np.pi*10*t) + np.random.randn(N)

Compute the autocovariance with a modifications: circular-shift the data

rxx = [];
lags = np.arange(-int(N/2),int(N/2));
for idx,L in enumerate(lags):
    rxx = np.append(rxx, 1/N*np.sum(np.roll(d,L) * d))
plt.plot(lags, rxx)
plt.xlabel('Lags [indices]');
plt.ylabel('Autocovariance rxx');

Compute the spectrum via FT and directly from the data.

# Compute the spectrum from the FT{rxx}
Sxx_via_rxx = 2*dt*np.fft.fft(rxx)

# Compute the spectrum from the data.
T   = t[-1];
xf  = np.fft.fft(d);
Sxx = 2 * dt ** 2 / T * (xf * xf.conj()) 

plt.plot(10*np.log10(Sxx))
plt.plot(10*np.log10(Sxx_via_rxx), 'o')
plt.xlabel('Freq Index')
plt.legend({'Direct Sxx', 'Sxx via rxx'})
plt.xlim([0,150]);