Module 1: Signal Processing Workbook
Welcome to Module 1! This workbook covers essential signal processing concepts for EEG analysis. You'll learn about sampling theory, frequency analysis, filtering techniques, and practical noise reduction methods. Each exercise builds your understanding of how to process and analyze EEG signals effectively.
Learning Objectives:
• Understand sampling theory and the Nyquist theorem
• Apply Fourier transforms to analyze EEG frequencies
• Design and implement digital filters for EEG data
• Perform time-frequency analysis
• Remove common noise sources from EEG recordings
Topic 1: Digital Signal Processing
Complete the following code to demonstrate sampling theory and aliasing:
import numpy as np
import matplotlib.pyplot as plt
# Original continuous signal parameters
true_frequency = 40 # Hz (beta band)
duration = 1 # second
# Create "continuous" signal (very high sampling rate)
continuous_fs = 10000 # Hz
t_continuous = np.linspace(0, duration, int(continuous_fs * duration))
continuous_signal = np.sin(2 * np.pi * * t_continuous)
# Sample at different rates
fs_low = 50 # Hz (below Nyquist)
fs_good = 100 # Hz (above Nyquist)
# Nyquist frequency for each sampling rate
nyquist_low = / 2
nyquist_good = / 2
# Create sampled versions
t_low = np.arange(0, duration, )
t_good = np.arange(0, duration, 1/ )
# Sample the continuous signal
sampled_low = np.sin(2 * np.pi * true_frequency * )
sampled_good = np.sin(2 * np.pi * true_frequency * )
Questions:
- What is the Nyquist frequency for a 256 Hz sampling rate?
- Which EEG frequency bands would be affected by aliasing with a 50 Hz sampling rate?
Calculate the apparent frequency when a signal is undersampled:
f_apparent = |f_true - n × f_sampling|
Given a 45 Hz beta wave sampled at 60 Hz:
| Parameter |
Value |
Your Calculation |
| True frequency (f_true) |
45 Hz |
- |
| Sampling frequency (f_s) |
60 Hz |
- |
| Nyquist frequency |
? |
Hz |
| Apparent frequency |
? |
Hz |
This aliased frequency would appear in which EEG band?
Topic 2: Frequency Domain Analysis
Complete the code to perform FFT on EEG data and identify frequency components:
import numpy as np
from scipy import signal
def analyze_eeg_spectrum(eeg_signal, sampling_rate):
"""
Compute frequency spectrum of EEG signal
"""
# Length of signal
n = len(eeg_signal)
# Compute FFT
fft_values = np.fft. (eeg_signal)
# Get frequency bins
freqs = np.fft. (n, sampling_rate)
# Get magnitude (only positive frequencies)
magnitude = np.abs(fft_values[:n//2]) *
freqs = freqs[: ]
# Convert to power spectral density
psd = magnitude **
return freqs, psd
# Find dominant frequency
def find_peak_frequency(freqs, psd, band_range):
"""Find peak frequency within a specific band"""
# Extract band of interest
band_mask = (freqs >= band_range[0]) & (freqs <= )
band_freqs = freqs[ ]
band_psd = psd[band_mask]
# Find peak
peak_idx = np. (band_psd)
peak_freq = band_freqs[peak_idx]
return peak_freq
Match the frequency bands with their typical EEG characteristics:
1. Delta (0.5-4 Hz)
2. Theta (4-8 Hz)
3. Alpha (8-13 Hz)
4. Beta (13-30 Hz)
5. Gamma (30-100 Hz)
A. Eyes closed, relaxation
B. Deep sleep
C. Active thinking, focus
D. Drowsiness, meditation
E. Cognitive processing
Match: 1-___, 2-___, 3-___, 4-___, 5-___
Calculate power in specific frequency bands:
def calculate_band_power(psd, freqs, band):
"""Calculate total power in a frequency band"""
# Define standard EEG bands
bands = {
'delta': (0.5, 4),
'theta': (4, 8),
'alpha': (8, 13),
'beta': (13, 30),
'gamma': (30, 50)
}
# Get band limits
low_freq, high_freq = bands[ ]
# Find indices for band
idx = np.where((freqs >= ) &
(freqs <= ))[0]
# Calculate power (area under PSD curve)
# Using trapezoidal integration
band_power = np. (freqs[idx], psd[idx])
# Calculate relative power (percentage)
total_power = np.trapz( , )
relative_power = (band_power / total_power) *
return band_power, relative_power
For a typical relaxed EEG recording with eyes closed, rank these bands by expected power (highest to lowest):
1. ________________
2. ________________
3. ________________
4. ________________
5. ________________
Topic 3: Digital Filtering for EEG
Design different types of filters for EEG processing:
from scipy import signal
def design_eeg_filters(fs):
"""Design common EEG filters"""
# High-pass filter (remove DC drift)
hp_cutoff = 0.5 # Hz
hp_order = # Filter order (even number)
# Low-pass filter (anti-aliasing)
lp_cutoff = 50 # Hz
lp_order = 6
# Notch filter (remove power line interference)
notch_freq = # Hz (US: 60, Europe: 50)
notch_q = 30 # Quality factor
# Design Butterworth filters
nyquist = fs / 2
# High-pass
b_hp, a_hp = signal.butter(hp_order, hp_cutoff/ ,
btype=' ')
# Low-pass
b_lp, a_lp = signal.butter(lp_order, /nyquist,
btype='low')
# Band-pass for alpha band
alpha_low = 8
alpha_high = 13
b_bp, a_bp = signal.butter(4, [alpha_low/nyquist, /nyquist],
btype=' ')
# Notch filter
b_notch, a_notch = signal. (notch_freq, notch_q, fs)
return (b_hp, a_hp), (b_lp, a_lp), (b_bp, a_bp), (b_notch, a_notch)
Select all correct statements about digital filters:
Apply filters correctly to avoid phase distortion:
def apply_filter(data, b, a, method='filtfilt'):
"""
Apply digital filter to EEG data
Parameters:
- data: EEG signal
- b, a: filter coefficients
- method: 'filter' or 'filtfilt'
"""
if method == 'filter':
# Forward filtering only (causes phase shift)
filtered = signal. (b, a, data)
elif method == 'filtfilt':
# Forward-backward filtering (zero phase shift)
filtered = signal. (b, a, data)
return filtered
# Compare filter effects
def compare_filter_methods(signal_data, b, a):
# Apply both methods
forward_only = apply_filter(signal_data, b, a, ' ')
zero_phase = apply_filter(signal_data, b, a, ' ')
# Calculate phase shift
phase_shift = np.angle(np.fft.fft(forward_only)) - np.angle(np.fft.fft(signal_data))
return forward_only, zero_phase, phase_shift
When should you use each filtering method?
| Scenario |
Best Method |
| Real-time EEG processing |
|
| Offline ERP analysis |
|
| Phase-sensitive connectivity analysis |
|
Topic 4: Time-Frequency Analysis
Implement STFT for time-frequency analysis of EEG:
def compute_stft(signal, fs, window_size, overlap):
"""
Compute Short-Time Fourier Transform
Parameters:
- signal: EEG data
- fs: sampling rate
- window_size: window length in samples
- overlap: overlap between windows (0-1)
"""
# Calculate hop size
hop_size = int(window_size * (1 - ))
# Create window function
window = signal.windows. (window_size)
# Calculate STFT parameters
f, t, Zxx = signal.stft(signal, fs=fs,
window= ,
nperseg= ,
noverlap=int(window_size * ))
# Convert to power
power = np.abs(Zxx) **
return f, t, power
# Choose appropriate window size
def select_window_size(fs, freq_resolution_needed):
"""
Calculate window size for desired frequency resolution
Frequency resolution = fs / window_size
"""
window_size = int( / freq_resolution_needed)
# Make it even
if window_size % 2 != 0:
window_size +=
return window_size
For a 250 Hz sampling rate, what window size gives 1 Hz frequency resolution?
Frequency resolution = Sampling rate / Window size
Compare wavelet transform with STFT for EEG analysis:
import pywt
def wavelet_analysis(signal, fs, wavelet='morl'):
"""
Perform continuous wavelet transform
Parameters:
- signal: EEG data
- fs: sampling rate
- wavelet: wavelet type ('morl', 'cmor', 'mexh')
"""
# Define scales for frequency range 1-50 Hz
frequencies = np.arange(1, 51, 1)
# Convert frequencies to scales
# scale = central_frequency * fs / frequency
central_freq = pywt.central_frequency( )
scales = (central_freq * fs) /
# Compute CWT
coefficients, freqs = pywt. (signal,
scales,
wavelet,
sampling_period=1/ )
# Calculate power
power = np.abs(coefficients) ** 2
return power, freqs
# Time-frequency resolution trade-off
def resolution_comparison():
"""
Compare time-frequency resolution of different methods
"""
# STFT: Fixed window → constant resolution
stft_time_res = # Proportional to window size
stft_freq_res = # Inversely proportional to window size
# Wavelet: Variable window → adaptive resolution
# High frequency: time resolution, frequency resolution
# Low frequency: time resolution, frequency resolution
Match the analysis method with its best use case:
| Use Case |
Best Method |
| Detecting brief gamma bursts |
____________________ |
| Analyzing steady-state responses |
____________________ |
| Tracking slow delta oscillations |
____________________ |
| Event-related spectral perturbations |
____________________ |
Topic 5: Practical Noise Reduction
Identify and remove different types of noise from EEG:
def identify_powerline_interference(signal, fs, powerline_freq=60):
"""
Detect powerline interference in EEG signal
"""
# Compute spectrum
freqs, psd = signal.periodogram(signal, fs=fs)
# Find power at powerline frequency
powerline_idx = np.argmin(np.abs(freqs - ))
powerline_power = psd[ ]
# Compare to surrounding frequencies
surrounding_idx = np.arange(powerline_idx-5, powerline_idx+6)
surrounding_idx = surrounding_idx[surrounding_idx != powerline_idx]
surrounding_power = np. (psd[surrounding_idx])
# Calculate SNR at powerline frequency
snr_powerline = powerline_power /
# Threshold for detection
if snr_powerline > : # Typical threshold
return True, snr_powerline
else:
return False, snr_powerline
def remove_baseline_drift(signal, fs, cutoff=0.5):
"""
Remove slow baseline drift from EEG
"""
# Design high-pass filter
nyquist = fs / 2
b, a = signal.butter( , cutoff/nyquist, btype=' ')
# Apply zero-phase filter
filtered = signal. (b, a, signal)
return filtered
def detect_movement_artifacts(signal, threshold_factor=5):
"""
Detect large amplitude movement artifacts
"""
# Calculate signal statistics
signal_std = np. (signal)
signal_mean = np.mean(signal)
# Find outliers
artifact_mask = np.abs(signal - signal_mean) > (threshold_factor * )
return artifact_mask
Match each artifact type with its characteristic and removal method:
| Artifact |
Frequency Range |
Removal Method |
| Eye blinks |
|
____________________ |
| Muscle (EMG) |
|
____________________ |
| Powerline |
|
____________________ |
Build a comprehensive preprocessing pipeline for EEG data:
class EEGPreprocessor:
def __init__(self, fs, powerline_freq=60):
self.fs = fs
self.nyquist = fs / 2
self.powerline_freq = powerline_freq
def preprocess_pipeline(self, raw_eeg):
"""
Complete preprocessing pipeline
"""
processed = raw_eeg.copy()
# Step 1: Remove baseline drift
hp_freq = # Hz
b_hp, a_hp = signal.butter(4, hp_freq/self.nyquist, btype='high')
processed = signal.filtfilt(b_hp, a_hp, processed)
# Step 2: Remove powerline interference
b_notch, a_notch = signal.iirnotch( ,
Q= ,
fs=self.fs)
processed = signal.filtfilt( , , processed)
# Step 3: Low-pass filter (anti-aliasing)
lp_freq = # Hz (for EEG)
b_lp, a_lp = signal.butter(6, lp_freq/self.nyquist, btype='low')
processed = signal.filtfilt(b_lp, a_lp, processed)
# Step 4: Detect and mark bad segments
artifact_mask = self.detect_artifacts(processed)
return processed, artifact_mask
def detect_artifacts(self, signal_data, threshold=5):
"""
Detect various artifacts
"""
# Z-score based detection
z_scores = np.abs((signal_data - np.mean(signal_data)) / np.std(signal_data))
# Mark artifacts
artifacts = z_scores >
# Extend artifact regions
for i in range(len(artifacts)):
if artifacts[i]:
# Mark surrounding samples too
start = max(0, i - int(0.1 * self.fs)) # 100ms before
end = min(len(artifacts), i + int(0.1 * self.fs)) # 100ms after
artifacts[start:end] =
return artifacts
Order these preprocessing steps correctly:
1. ____________________
2. ____________________
3. ____________________
4. ____________________
5. ____________________
Options: Low-pass filter, Artifact detection, Baseline correction,
Notch filter, Segmentation/Epoching
Project: Build a Frequency Analyzer
Create a comprehensive frequency analyzer for EEG signals:
class EEGFrequencyAnalyzer:
def __init__(self, fs=256):
self.fs = fs
self.bands = {
'delta': (0.5, 4),
'theta': (4, 8),
'alpha': (8, 13),
'beta': (13, 30),
'gamma': (30, 50)
}
def analyze_signal(self, eeg_signal, method='welch'):
"""
Perform complete frequency analysis
"""
results = {}
# Compute power spectrum
if method == 'welch':
freqs, psd = signal. (eeg_signal,
fs=self.fs,
nperseg=self.fs*2) # 2-second windows
else:
freqs, psd = signal.periodogram(eeg_signal, fs=self.fs)
results['frequencies'] = freqs
results['psd'] = psd
# Calculate band powers
band_powers = {}
total_power = np. (freqs, psd)
for band_name, (low, high) in self.bands.items():
# Find band indices
band_idx = np.where((freqs >= low) & (freqs <= high))[0]
# Calculate absolute and relative power
abs_power = np.trapz(psd[band_idx], freqs[band_idx])
rel_power = (abs_power / total_power) * 100
band_powers[band_name] = {
'absolute': abs_power,
'relative': rel_power,
'peak_freq': freqs[band_idx[np. (psd[band_idx])]]
}
results['band_powers'] = band_powers
# Calculate spectral features
results['spectral_edge_95'] = self.spectral_edge(freqs, psd, 0.95)
results['mean_frequency'] = self.mean_frequency(freqs, psd)
results['spectral_entropy'] = self.spectral_entropy(psd)
return results
def spectral_edge(self, freqs, psd, percentage):
"""Find frequency below which percentage of power is contained"""
cumsum_power = np. (psd)
total_power = cumsum_power[-1]
edge_idx = np.where(cumsum_power >= total_power * percentage)[0][0]
return freqs[edge_idx]
def mean_frequency(self, freqs, psd):
"""Calculate mean frequency weighted by power"""
return np.sum(freqs * psd) / np. (psd)
def spectral_entropy(self, psd):
"""Calculate spectral entropy (measure of signal complexity)"""
# Normalize PSD to probability distribution
psd_norm = psd / np.sum(psd)
# Calculate entropy
entropy = -np.sum(psd_norm * np. (psd_norm + 1e-15))
return entropy
Design a visualization for your frequency analyzer. Draw a sketch showing:
Include: PSD plot, band power bars, and spectral features display
Answer these questions about applying your frequency analyzer:
- What window length (in seconds) would you use for Welch's method to get good frequency resolution for delta waves?
- If alpha power increases from 20% to 40% when eyes close, what is this phenomenon called?
- Calculate the frequency resolution for a 5-second EEG segment sampled at 250 Hz:
Show your calculation:
- Which spectral feature would best distinguish between sleep stages?
Answer Key
Check your answers after completing all exercises.
Exercise 1: true_frequency, fs_low, fs_good, 1/fs_low, fs_good, t_low, t_good
Q1: B (128 Hz), Q2: Beta and Gamma bands would be affected
Exercise 2: Nyquist = 30 Hz, Apparent frequency = 15 Hz (|45-60| = 15)
This appears in the Beta band (13-30 Hz)
Exercise 3: fft, fftfreq, 2/n, n//2, 2, band_range[1], band_mask, argmax
Matching: 1-B, 2-D, 3-A, 4-C, 5-E
Exercise 4: band, low_freq, high_freq, trapz, psd, freqs, 100
Ranking: Alpha, Beta, Theta, Delta, Gamma
Exercise 5: 4, 60 (or 50), nyquist, high, lp_cutoff, alpha_high/nyquist, band, iirnotch
All statements are correct
Exercise 6: lfilter, filtfilt, filter, filtfilt
Real-time: filter, Offline ERP: filtfilt, Phase-sensitive: filtfilt
Exercise 7: overlap, hann, window, window_size, overlap, 2, fs, 1
250 samples gives 1 Hz resolution
Exercise 8: wavelet, frequencies, cwt, fs, constant, constant, good, poor, poor, good
Gamma bursts: Wavelet, Steady-state: STFT, Delta: STFT, ERSP: Wavelet
Exercise 9: powerline_freq, powerline_idx, mean, surrounding_power, 3, 4, high, filtfilt, std, signal_std
Eye blinks: 1-4 Hz/ICA, EMG: >20 Hz/Low-pass, Powerline: 50/60 Hz/Notch
Exercise 10: 0.5, self.powerline_freq, 30, b_notch, a_notch, 50, threshold, True
Order: 1. Baseline correction, 2. Notch filter, 3. Low-pass filter, 4. Artifact detection, 5. Segmentation
Exercise 11: welch, trapz, argmax, cumsum, sum, log2
Window length: 4 seconds for delta, Phenomenon: Alpha blocking/enhancement,
Resolution: 0.2 Hz, Best feature: Delta/theta power ratio
Congratulations! 🎉
You've completed Module 1 on Signal Processing Fundamentals! You now understand sampling theory, frequency analysis, digital filtering, and noise reduction techniques essential for EEG analysis.
Next step: Move on to Module 2 - Neuroscience for EEG