Analyzing Rhythms Part 1

Load the data and look at it.

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('Data/EEG-1.mat')    # Load the EEG data
# Data available here
# https://github.com/Mark-Kramer/BU-MA665-MA666/tree/master/Data
#
EEG  = data['EEG'].reshape(-1) # Extract the EEG variable
t    = data['t'][0]            # ... and the t variablea

Look at it.

plt.plot(t,EEG)
plt.xlabel('Time [s]')
plt.ylabel('EEG');

Preliminaries

Q: What is the sampling interval (dt)?

dt = "SOMETHING"
print(dt)

Q. What is the sampling frequency (f0)?

f0 = "SOMETHING"
print(f0)

Q. What is the number of points in the data (N)?

N = "SOMETHING"
print(N)

Q. What is the total time of the observation (T)?

T = "SOMETHING"
print(T)

---

Compute the spectrum “by hand”

# Q. What is the Nyquist frequency and frequency resolution?
fNQ = "SOMETHING"
df  = "SOMETHING"

# Q. What is the frequency axis?
fj  = np.arange(0,fNQ,df)

# Then, compute the Fourier transform "by hand".
x = EEG
X = np.ndarray(np.size(fj), complex);
for j in range( np.size(fj) ):
    X[j] = "SOMETHING"

# And the spectrum,
Sxx = "SOMETHING"

# Plot it,
plot(fj, Sxx.real)
plt.xlim([0, 100])                          # Select frequency range
plt.xlabel('Frequency [Hz]')                # Label the axes
plt.ylabel('Power [$\mu V^2$/Hz]');

Compute the power spectrum using the FFT function.

x   = EEG
X   = np.fft.fft(x)                        # Compute Fourier transform of x
Sxx = "SOMETHING"                          # Compute spectrum
Sxx = Sxx[0:int(N / 2)].real               # Ignore negative frequencies

Define the frequency axis

df  = "SOMETHING"                  # Determine frequency resolution
fNQ = "SOMETHING"                  # Determine Nyquist frequency
faxis = np.arange(0,fNQ,df)        # Construct frequency axis

Plot the spectrum versus frequency

Q. What do you see?

plt.plot(faxis, Sxx)
plt.xlim([0, 100])                          # Select frequency range
plt.xlabel('Frequency [Hz]')                # Label the axes
plt.ylabel('Power [$\mu V^2$/Hz]');

Plot the spectrum versus frequency on a decibel scale.

Q. Now what do you see?

plt.figure()
plt.plot(faxis, [???])   # Plot the spectrum in decibels.
plt.xlim([0, 100])                               # Select the frequency range.
plt.ylim([-60, 0])                               # Select the decibel range.
plt.xlabel('Frequency [Hz]')                     # Label the axes.
plt.ylabel('Power [dB]');

Plot the spectrum versus frequency on a logarithmic frequency axis.

Q. And now what do you see?

plt.figure()
plt.[???]                                            # Log-log scale
plt.xlim([df, 100])                                  # Select frequency range
plt.ylim([-60, 0])                                   # ... and the decibel range.
plt.xlabel('Frequency [Hz]')                         # Label the axes.
plt.ylabel('Power [dB]');

Apply Hanning taper

Apply the Hanning taper to the data.

x_tapered  = ???
plt.figure()
plt.plot(t,x)
plt.plot(t,x_tapered);

Apply the Hanning taper and look at the spectrum.

xf_tapered  = "SOMETHING"                        # Compute Fourier transform of x.
Sxx_tapered = "SOMETHING"                        # Compute the spectrum,
Sxx_tapered = np.real(Sxx_tapered[:int(N / 2)])  # ... and ignore negative frequencies.

plt.figure()
plt.semilogx(faxis,10*np.log10(Sxx))         # Plot spectrum of untapered signal.  
plt.semilogx(faxis,10*np.log10(Sxx_tapered)) # Plot spectrum vs tapered signal.
plt.xlim([faxis[1], 100])                    # Select frequency range,
plt.ylim([-70, 20])                          # ... and the power range.
plt.xlabel('Frequency [Hz]')                 # Label the axes
plt.ylabel('Power [$\mu V^2$/Hz]')

Spectrogram

# Plot the spectrogram.

Fs = 1 / dt               # Define the sampling frequency,
interval = int(Fs)        # ... the interval size,
overlap = int(Fs * 0.95)  # ... and the overlap intervals

                          # Compute the spectrogram
f0, t0, Sxx0 = spectrogram(
    EEG,                  # Provide the signal,
    fs=Fs,                # ... the sampling frequency,
    nperseg=interval,     # ... the length of a segment,
    noverlap=overlap)     # ... the number of samples to overlap,
plt.pcolormesh(t0, f0, 10 * np.log10(Sxx0),
               cmap='jet')# Plot the result
plt.colorbar()            # ... with a color bar,
plt.ylim([0, 70])             # ... set the frequency range,
plt.xlabel('Time [s]')       # ... and label the axes
plt.ylabel('Frequency [Hz]');