Chapter 2

Analyzing Neural Time Series Data

Python code for Chapter 2 – converted from original Matlab by AE Studio (and ChatGPT)
Original Matlab code by Mike X Cohen

This code accompanies the book, titled “Analyzing Neural Time Series Data” (MIT Press).
Using the code without following the book may lead to confusion, incorrect data analyses, and misinterpretations of results.
Mike X Cohen and AE Studio assume no responsibility for inappropriate or incorrect use of this code.

Import necessary libraries

Type %pip install my_package in a cell to install any missing packages

import numpy as np
import matplotlib.pyplot as plt
from scipy.io import loadmat
from scipy.fft import fft, ifft
from scipy.signal import firls, filtfilt

Figure 2.1

# Load the sample EEG data
EEG = loadmat('../data/sampleEEGdata.mat')['EEG'][0, 0]

# Set the number of trials to analyze
nTrials = 6 # modifiable

# Initialize data array
data = np.zeros((nTrials, EEG['pnts'][0][0]))

# Define wavelet parameters
wavetime = np.arange(-1, 1 + 1/EEG['srate'][0][0], 1/EEG['srate'][0][0])
n_conv = len(wavetime) + EEG['pnts'][0][0] - 1
waveletfft = fft(np.exp(2 * 1j * np.pi * 10 * wavetime) * np.exp(-wavetime**2 / (2 * (5 / (2 * np.pi * 10))**2)) / 10, n_conv)
data10hz = np.zeros((nTrials, EEG['pnts'][0][0]))

# Plotting setup
fig, axs = plt.subplots(nTrials, 3, figsize=(15, 10))

# Loop over trials
for triali in range(nTrials):
    
    # Create single trial "ERP"
    data[triali, :] = 0.5 * np.sin(2 * np.pi * 6 * EEG['times'][0] / 1000 + 2 * np.pi * triali / nTrials - np.pi) + np.random.randn(EEG['pnts'][0][0]) / 6
    # Add non-phase-locked stimulus potential
    data[triali, 259:360] = data[triali, 259:360] + np.sin(2 * np.pi * 10 * EEG['times'][0][259:360] / 1000 + 2 * np.pi * triali / nTrials - np.pi) + np.random.randn(101) / 5
    
    # Plot data from this trial
    axs[triali, 0].plot(EEG['times'][0], data[triali, :])
    axs[triali, 0].set_xlim([-250, 850])
    axs[triali, 0].set_ylim([-2.2, 2.2])
    
    # Plot ERP from trial 1 to current
    axs[triali, 1].plot(EEG['times'][0], np.mean(data[:triali+1, :], axis=0))
    axs[triali, 1].set_xlim([-250, 850])
    axs[triali, 1].set_ylim([-2.2, 2.2])
    
    # Convolve with 10 Hz wavelet
    convolution_result_fft = ifft(waveletfft * fft(data[triali, :], n_conv)) * np.sqrt(5 / (2 * np.pi * 10))
    convolution_result_fft = convolution_result_fft[int(np.floor(len(wavetime)/2)) : -int(np.floor(len(wavetime)/2))]
    data10hz[triali, :] = np.abs(convolution_result_fft)**2
    
    # Plot 10 Hz power
    axs[triali, 2].plot(EEG['times'][0], np.mean(data10hz[:triali+1, :], axis=0))
    axs[triali, 2].set_xlim([-250, 850])
    axs[triali, 2].set_ylim([-0.1, 0.8])

plt.tight_layout()
plt.show()

Figure 2.2

This code involves performing convolution with a complex Morlet wavelet, which you will learn about in Chapters 10-13.

# Define parameters for the time-frequency analysis
srate = 1000
time = np.arange(0, 10 + 1/srate, 1/srate)
DCoffset = -0.5

# Create multi-frequency signal
a = np.sin(2 * np.pi * 10 * time)  # High frequency part
b = 0.1 * np.sin(2 * np.pi * 0.3 * time) + DCoffset  # Low frequency part

# Combine signals
data = a * b
data += (2 * np.sin(2 * np.pi * 3 * time) * np.sin(2 * np.pi * 0.07 * time) * 0.1 + DCoffset)

# Morlet wavelet convolution parameters
num_frex = 40
min_freq = 2
max_freq = 20

Ldata = len(data)
Ltapr = len(data)
Lconv1 = Ldata + Ltapr - 1
Lconv = 2**int(np.ceil(np.log2(Lconv1)))

frex = np.logspace(np.log10(min_freq), np.log10(max_freq), num_frex)

# Initialize time-frequency representation matrix
tf = np.zeros((num_frex, len(data)))
datspctra = fft(data, Lconv)

s = 4 / (2 * np.pi * frex)
t = np.arange(-((len(data) - 1) / 2) / srate, ((len(data) - 2) / 2) / srate + 1/srate, 1/srate)

# Perform wavelet convolution
for fi in range(len(frex)):
    wavelet = np.exp(2 * 1j * np.pi * frex[fi] * t) * np.exp(-t**2 / (2 * s[fi]**2))
    m = ifft(datspctra * fft(wavelet, Lconv), Lconv)
    m = m[:Lconv1]
    m = m[int(np.floor((Ltapr - 1) / 2)) : -int(np.ceil((Ltapr - 1) / 2))]
    tf[fi, :] = np.abs(m)**2

# Plot the signals and time-frequency representation
fig, axs = plt.subplots(2, 2, figsize=(12, 8))

axs[0, 0].plot(a)
axs[0, 0].set_xlim([1 * 1000, 8 * 1000])
axs[0, 0].set_ylim([-1, 1])
axs[0, 0].set_xticks(np.arange(1 * 1000, 9 * 1000, 1000))
axs[0, 0].set_xticklabels(np.arange(1, 9))
axs[0, 0].set_title('10 Hz signal, DC=0')

axs[0, 1].plot(b)
axs[0, 1].set_xlim([1 * 1000, 8 * 1000])
axs[0, 1].set_ylim([-1, 1])
axs[0, 1].set_xticks(np.arange(1 * 1000, 9 * 1000, 1000))
axs[0, 1].set_xticklabels(np.arange(1, 9))
axs[0, 1].set_title(f'.3 Hz signal, DC={DCoffset}')

axs[1, 0].plot(data)
axs[1, 0].set_xlim([1 * 1000, 8 * 1000])
axs[1, 0].set_ylim([-1, 1])
axs[1, 0].set_xticks(np.arange(1 * 1000, 9 * 1000, 1000))
axs[1, 0].set_xticklabels(np.arange(1, 9))
axs[1, 0].set_title('Time-domain signal')

im = axs[1, 1].imshow(tf, aspect='auto', extent=[1, len(data), 0, num_frex], origin='lower')
axs[1, 1].set_xlim([1 * 1000, 8 * 1000])
axs[1, 1].set_yticks(np.arange(1, num_frex, 8))
axs[1, 1].set_yticklabels(np.round(frex[::8]))
axs[1, 1].set_xticks(np.arange(1 * 1000, 9 * 1000, 1000))
axs[1, 1].set_xticklabels(np.arange(1, 9))
axs[1, 1].set_title('Time-frequency representation')

plt.tight_layout()
plt.show()

Figure 2.3

# Select the channel to plot
chan2plot = 'Pz' # you can pick any electrode (type `EEG['chanlocs'][0]['labels']` in a new cell for all electrodes)

# Compute ERP (time-domain trial average from selected electrode)
erp = np.squeeze(np.mean(EEG['data'][EEG['chanlocs'][0]['labels']==chan2plot, :, :], axis=2))

# Low-pass filter parameters
nyquist = EEG['srate'][0][0] / 2
filter_cutoff = 40  # Hz
trans_width = 0.1  # Transition width, in fraction of 1

# Filter design
ffrequencies = [0, filter_cutoff, filter_cutoff * (1 + trans_width), nyquist] / nyquist
idealresponse = [1, 1, 0, 0]
filterweights = firls(101, ffrequencies, idealresponse)
filtered_erp = filtfilt(filterweights, 1, erp)

# Plotting the ERP
plt.figure(figsize=(10, 5))
plt.plot(EEG['times'][0], filtered_erp, 'k.-', label='256 Hz')

# Down-sample and plot
times2plot = np.arange(-200, 1001, 40)
downsampled_erp = filtered_erp[::5]
downsampled_times = EEG['times'][0][::5]
plt.plot(downsampled_times, downsampled_erp, 'mo-', label='50 Hz')

plt.xlim([-200, 1000])
plt.ylim([-4, 10])
plt.gca().invert_yaxis()
plt.xlabel('Time (ms)')
plt.ylabel('Voltage (µV)')
plt.title('ERP from electrode {}'.format(chan2plot))
plt.legend()

plt.show()