import numpy as np
import matplotlib.pyplot as plt
from scipy.fft import fft, ifft
from scipy.io import loadmat
Chapter 11
Chapter 11
Analyzing Neural Time Series Data
Python code for Chapter 11 – 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
Figure 11.1
# Define parameters
= 1000 # sampling rate of 1 kHz
srate = np.arange(-1, 1+1/srate, 1/srate)
time = 10 # in Hz
freq = 2 # amplitude, or height of the sine wave
amp
# Create sine wave
= amp * np.sin(2 * np.pi * freq * time)
sine_wave
# Plot sine wave
plt.figure()
plt.plot(time, sine_wave)-1, 1])
plt.xlim([-3, 3]) # this adjusts the y-axis limits for visibility
plt.ylim(['My first sine wave!')
plt.title( plt.show()
Figure 11.2
# Define a sampling rate
= 500
srate
# List some frequencies
= [3, 10, 5, 15, 35]
frex
# List some random amplitudes
= [5, 15, 10, 5, 7]
amplit
# Phases
= [np.pi/7, np.pi/8, np.pi, np.pi/2, -np.pi/4]
phases
# Define time
= np.arange(-1, 1+1/srate, 1/srate)
time
# Loop through frequencies and create sine waves
= np.zeros((len(frex), len(time)))
sine_waves for fi in range(len(frex)):
= amplit[fi] * np.sin(2 * np.pi * frex[fi] * time + phases[fi])
sine_waves[fi, :]
# Plot each wave separately
plt.figure()for fi in range(len(frex)):
len(frex), 1, fi+1)
plt.subplot(=2)
plt.plot(sine_waves[fi, :], linewidth0, len(time), -max(amplit), max(amplit)])
plt.axis([
plt.tight_layout()
plt.show()
# Plot the result
plt.figure()sum(sine_waves, axis=0))
plt.plot(np.=True, axis='both', tight=True)
plt.gca().autoscale(enable'sum of sine waves')
plt.title( plt.show()
Figure 11.3
# Plot sum of sine waves plus random noise
plt.figure()sum(sine_waves + 5 * np.random.randn(*sine_waves.shape), axis=0))
plt.plot(np.0, 1020, -40, 50]) # this sets the x-axis and y-axis limits
plt.axis(['sum of sine waves plus white noise')
plt.title( plt.show()
Figure 11.4
= np.arange(-1, 1+1/srate, 1/srate)
time
# Create three sine waves
= np.sin(2 * np.pi * 3 * time)
s1 = 0.5 * np.sin(2 * np.pi * 8 * time)
s2 = s1 + s2
s3
# Plot the sine waves
plt.figure()for i, s in enumerate([s1, s2, s3], 1):
2, 3, i)
plt.subplot(
plt.plot(time, s)-1, 1])
plt.xlim([-1.6, 1.6])
plt.ylim([-1.5, 2, 0.5))
plt.yticks(np.arange(
# Plot power
2, 3, i+3)
plt.subplot(= fft(s) / len(time)
f = np.linspace(0, srate/2, int(np.floor(len(time)/2)+1))
hz abs(f[:len(hz)]) * 2)
plt.bar(hz, np.0, 11])
plt.xlim([0, 1.2])
plt.ylim([0, 11))
plt.xticks(np.arange(
plt.tight_layout() plt.show()
Figure 11.5
# Length of sequence
= 10
N # Random numbers
= np.random.randn(N)
data # Sampling rate in Hz
= 200
srate # Nyquist frequency
= srate / 2
nyquist # Initialize Fourier output matrix
= np.zeros(N, dtype=complex)
fourier # Frequencies in Hz
= np.linspace(0, nyquist, N//2+1)
frequencies # Time
= np.arange(N) / N
time
# Fourier transform is dot-product between sine wave and data at each frequency
for fi in range(N):
= np.exp(-1j * 2 * np.pi * (fi) * time)
sine_wave = np.sum(sine_wave * data)
fourier[fi] = fourier / N
fourier
# Plotting
= plt.figure(figsize=(12, 8))
fig
= fig.add_subplot(221)
ax1 '-o')
ax1.plot(data, -1, N])
ax1.set_xlim(['Time domain representation of the data')
ax1.set_title(
= fig.add_subplot(222, projection='3d')
ax2 //2+1]), np.abs(fourier[:N//2+1])**2, '-o', linewidth=3)
ax2.plot3D(frequencies, np.angle(fourier[:NTrue)
ax2.grid('Frequency (Hz)')
ax2.set_xlabel('Phase')
ax2.set_ylabel('Power')
ax2.set_zlabel('3-D representation of the Fourier transform')
ax2.set_title(
= fig.add_subplot(223)
ax3 abs(fourier[:N//2+1])**2)
ax3.bar(frequencies, np.-5, 105])
ax3.set_xlim(['Frequency (Hz)')
ax3.set_xlabel('Power')
ax3.set_ylabel('Power spectrum derived from discrete Fourier transform')
ax3.set_title(
= fig.add_subplot(224)
ax4 //2+1]))
ax4.bar(frequencies, np.angle(fourier[:N-5, 105])
ax4.set_xlim(['Frequency (Hz)')
ax4.set_xlabel('Phase angle')
ax4.set_ylabel(-np.pi, np.pi, np.pi/2))
ax4.set_yticks(np.arange('Phase spectrum derived from discrete Fourier transform')
ax4.set_title(
plt.tight_layout() plt.show()
Figure 11.6
# Compute sine waves and sum to recover the original time series
= np.zeros(N)
reconstructed_data for fi in range(N):
= fourier[fi] * np.exp(1j * 2 * np.pi * (fi) * time)
sine_wave += np.real(sine_wave)
reconstructed_data
# Plot original and reconstructed data
plt.figure()'-o', label='original data')
plt.plot(data, 'r-*', label='inverse Fourier transform data')
plt.plot(reconstructed_data, 0, N-1])
plt.xlim([
plt.legend() plt.show()
Figure 11.7
# Perform FFT
= fft(data) / N
fft_data
# Plotting
=(18, 6))
plt.figure(figsize
131)
plt.subplot(abs(fourier[:N//2+1])**2, '*-', label='time-domain Fourier')
plt.plot(frequencies, np.abs(fft_data[:N//2+1])**2, 'ro-', markersize=8, label='FFT')
plt.plot(frequencies, np.'Frequency (Hz)')
plt.xlabel('Power')
plt.ylabel('Power spectrum derived from discrete Fourier transform and from FFT')
plt.title(
plt.legend()
132)
plt.subplot(//2+1]), '*-', label='time-domain Fourier')
plt.plot(frequencies, np.angle(fourier[:N//2+1]), 'ro-', markersize=8, label='FFT')
plt.plot(frequencies, np.angle(fft_data[:N'Frequency (Hz)')
plt.xlabel('Phase')
plt.ylabel(-np.pi, np.pi+1, np.pi/2))
plt.yticks(np.arange('Phase spectrum derived from discrete Fourier transform and from FFT')
plt.title(
133)
plt.subplot('*-', label='Manual inverse Fourier transform')
plt.plot(reconstructed_data, 'ro-', markersize=8, label='ifft')
plt.plot(np.real(ifft(fft(data))), 'Time')
plt.xlabel('Amplitude')
plt.ylabel('Manual inverse Fourier transform and ifft')
plt.title(
plt.tight_layout() plt.show()
Figure 11.9
# List some frequencies
= [3, 10, 5, 7]
frex
# List some random amplitudes
= [5, 15, 10, 5]
amplit
# Phases
= [np.pi/7, np.pi/8, np.pi, np.pi/2]
phases
# Create a time series of sequenced sine waves
= 500
srate = np.arange(-1, 1+1/srate, 1/srate)
time = np.zeros(len(time) * len(frex))
stationary = np.zeros(len(time) * len(frex))
nonstationary
for fi in range(len(frex)):
= amplit[fi] * np.sin(2 * np.pi * frex[fi] * time + phases[fi])
temp_sine_wave += np.tile(temp_sine_wave, len(frex))
stationary *= (time + 1)
temp_sine_wave = fi * len(time)
start_idx = start_idx + len(time)
stop_idx = temp_sine_wave
nonstationary[start_idx:stop_idx]
# Plot stationary and non-stationary signals
=(10, 6))
plt.figure(figsize
221)
plt.subplot('r')
plt.plot(stationary, 1, len(stationary)])
plt.xlim(['stationary signal')
plt.title(
222)
plt.subplot(
plt.plot(nonstationary)1, len(nonstationary)])
plt.xlim(['non-stationary signal')
plt.title(
# Perform FFT and plot
= np.linspace(0, srate/2, len(nonstationary)//2+1)
frequencies = fft(nonstationary) / len(nonstationary)
fft_nonstationary = fft(stationary) / len(stationary)
fft_stationary 212)
plt.subplot(abs(fft_stationary[:len(frequencies)])*2, 'r', label='Power stationary')
plt.plot(frequencies, np.abs(fft_nonstationary[:len(frequencies)])*2, label='Power non-stationary')
plt.plot(frequencies, np.0, max(frex)*2])
plt.xlim([
plt.legend()
plt.tight_layout() plt.show()
Figure 11.10
# Load sample EEG data
= loadmat('../data/sampleEEGdata.mat')["EEG"][0, 0]
EEG = EEG["data"][46, :, 0]
eegdat4convol = EEG["srate"][0][0]
srate
# Create Gaussian
= np.arange(-1, 1+1/srate, 1/srate)
time = 5 / (2 * np.pi * 30)
s = np.exp(-time**2 / (2 * s**2)) / 30
gaussian
# Plot EEG data and Gaussian
=(10, 4))
plt.figure(figsize
211)
plt.subplot(
plt.plot(eegdat4convol)
212)
plt.subplot(
plt.plot(gaussian)
plt.show()
# Plot convolution results
=(10, 4))
plt.figure(figsize
211)
plt.subplot('same'))
plt.plot(np.convolve(eegdat4convol, gaussian,
212)
plt.subplot(abs(fft(np.convolve(eegdat4convol, gaussian, 'same'))))
plt.plot(np.
plt.show()
# Plot FFT results
=(10, 4))
plt.figure(figsize
211)
plt.subplot(abs(fft(eegdat4convol)))
plt.plot(np.
212)
plt.subplot(abs(fft(gaussian)))
plt.plot(np.
plt.show()
Figure 11.11
# Define parameters
= 1000
srate = np.arange(-0.5, 0.5, 1/srate)
time = 20
f = [15, 5]
fg = np.sin(2 * np.pi * f * time)
s
# Loop over Gaussian widths
for i in range(2):
= np.exp(-time**2 / (2 * (4 / (2 * np.pi * fg[i])**2))) / fg[i]
g
=(10, 8))
plt.figure(figsize
411)
plt.subplot(
plt.plot(time, s)'Sine wave (signal)')
plt.title(-0.5, 0.5])
plt.xlim([-1.1, 1.1])
plt.ylim([
412)
plt.subplot(
plt.plot(time, g)-0.5, 0.5])
plt.xlim(['Gaussian (kernel)')
plt.title(
413)
plt.subplot('same'))
plt.plot(time, np.convolve(s, g, -0.5, 0.5])
plt.xlim([-1.1, 1.1])
plt.ylim(['result of convolution')
plt.title(
427)
plt.subplot(= np.abs(fft(s))
fft_s = fft_s[:len(fft_s)//2+1] / np.max(fft_s[:len(fft_s)//2+1])
fft_s 501), fft_s, color='r')
plt.bar(np.arange(
= np.abs(fft(g))
fft_g = fft_g[:len(fft_g)//2+1] / np.max(fft_g[:len(fft_g)//2+1])
fft_g 501), fft_g)
plt.plot(np.arange(0, 40])
plt.xlim([0, 1.05])
plt.ylim(['individual power spectra')
plt.title(
428)
plt.subplot(501), fft_g * fft_s)
plt.bar(np.arange(0, 40])
plt.xlim([0, 0.035])
plt.ylim(['multiplied power spectra')
plt.title(
plt.tight_layout() plt.show()
Figure 11.12
= EEG["times"][0]
times = EEG["srate"][0][0]
srate
# Create Gaussian
= np.arange(-1, 1+1/srate, 1/srate)
time = 5 / (2 * np.pi * 30)
s = np.exp(-time**2 / (2 * s**2)) / 30
gaussian
# Plot EEG data, Gaussian, and result of convolution
=(10, 8))
plt.figure(figsize
411)
plt.subplot(
plt.plot(times, eegdat4convol)-1000, 1500])
plt.xlim([
412)
plt.subplot(
plt.plot(time, gaussian)-1, 1])
plt.xlim([
413)
plt.subplot('r')
plt.plot(times, eegdat4convol, 'same'))
plt.plot(times, np.convolve(eegdat4convol, gaussian, -1000, 1500])
plt.xlim([
= len(eegdat4convol)
nfft = np.abs(fft(eegdat4convol, nfft))
fft_s = fft_s[:nfft//2+1]
fft_s = np.linspace(0, srate/2, nfft//2+1)
f
427)
plt.subplot(/ np.max(fft_s), 'r')
plt.plot(f, fft_s
= np.abs(fft(gaussian, nfft))
fft_g = fft_g[:nfft//2+1]
fft_g / np.max(fft_g))
plt.plot(f, fft_g
0, 60])
plt.xlim([
428)
plt.subplot(* fft_g)
plt.plot(f, fft_s 0, 60])
plt.xlim([
plt.tight_layout() plt.show()