Processing of electroencephalography (EEG) data
The example shows the process of digital signal processing of electroencephalography (EEG), implemented using the specialized EngeeDSP library. Using the example of real data recorded in the standard EDF format, the key stages of the analysis are consistently illustrated.
Function engee.clear() cleans up workspaces:
engee.clear()
Two graph rendering modes are available to visualize the results.:
-
For static, fast rendering, use
gr(). -
For interactive, dynamic visualization with the ability to scale and view values, use
plotlyjs().
Select one of the commands by commenting out the second one. Both backends are mutually exclusive, and the last one called will be activated.
# gr()
plotlyjs()
Reading an EDF file
We will connect using the function include the file "edfread.jl" for reading EDF files:
include("$(@__DIR__)/edfread.jl")
Function edfread It is designed to read data in the EDF (European Data Format) format, the standard format for storing biomedical signals. Structure hdr The information returned by this function contains the full meta information about the record:
General recording parameters:
-
ver— EDF format version -
patientID— patient ID -
recordID— record ID -
startdateandstarttime— date and time of the start of recording -
bytes— header size in bytes -
records— the number of data blocks in the file -
duration— the duration of one block in seconds -
ns— the number of channels in the recording
Parameters of each channel:
-
labels— channel names -
transducers— type of sensors -
physicalDims— physical units of measurement -
physicalMinsandphysicalMaxs— minimum and maximum physical values -
digitalMinsanddigitalMaxs— minimum and maximum digital values -
prefilters— filters applied during recording -
samples— the number of samples in one block for each channel
Electroencephalography (EEG) data will serve as an example of biomedical signal processing. Let's start by uploading the file S004R02.edf to study its structure and build waveforms of signals.
hdr, record = edfread("$(@__DIR__)/S004R02.edf")
The next step is to analyze the metadata of the EDF file. The key parameters that we will extract for further processing include the sampling rate, the number of channels, and the recording duration.
println("Format version: ", hdr.ver)
println("Patient ID: ", strip(hdr.patientID))
println("Description of the record: ", strip(hdr.recordID))
println("Start date/time: ", hdr.startdate, " ", hdr.starttime)
println("Number of channels: ", hdr.ns)
println("Number of entries: ", hdr.records)
println("Total duration: ", hdr.records * hdr.duration, " sec")
For clarity, we also present the key characteristics of the first three channels.
for i in 1:min(3, size(record, 1))
channel_data = record[i, :]
fs = hdr.samples[i] / hdr.duration
t = (1:length(channel_data)) ./ fs
println(" $(i). $(hdr.labels[i]):")
println(" The frequency is discretionary.: ", round(fs, digits=1), " Hz")
println(" Signal length: ", length(channel_data), " counts")
println(" Duration: ", round(t[end], digits=2), " sec")
println(" Min/max amplitude: ", round(minimum(channel_data), digits=2), " .. ", round(maximum(channel_data), digits=2))
end
Next, the signals of the selected EEG channels are visualized. The construction of waveforms allows you to visually assess the nature of brain activity and the presence of artifacts.
# Channel selection
Nk1 = 1
Nk2 = 17
eegNk1 = record[Nk1, :]
eegNk2 = record[Nk2, :]
# Signal Parameters
fs = hdr.samples[1] / hdr.duration # Sampling rate
Tn = 60.0 # Observation time
N1 = Int64(Tn * fs) # Sample Length
dims = hdr.physicalDims[Nk1]
# Signal extraction
y1 = eegNk1[1:N1]
y2 = eegNk2[1:N1]
t = range(0, step=1/fs, length=N1)
p1 = plot(t, y1,
title="EEG waveform (channel $Nk1)",
xlabel="Time, from",
ylabel="Amplitude, $dims",
grid=true)
p2 = plot(t, y2,
title="EEG Waveform ($Nk2 channel)",
xlabel="Time, from",
ylabel="Amplitude, $dims",
grid=true)
EEG signal processing
Spectral analysis of the EEG signal
Next, we switch to frequency data analysis using the EngeeDSP library. For the selected channels, a fast Fourier transform is applied, which converts the signals from a time representation to a frequency representation. The obtained spectra are centered for the convenience of subsequent visualization and analysis of the frequency components.
nsamp = length(y1)
df = fs / nsamp
freq_vec = -fs/2 : df : fs/2 - df
fft_y1 = EngeeDSP.Functions.fft(y1)
fft_y2 = EngeeDSP.Functions.fft(y2)
# Centered spectrum
fft_Y1 = EngeeDSP.Functions.fftshift(fft_y1)
fft_Y2 = EngeeDSP.Functions.fftshift(fft_y2)
After receiving the frequency representation of the signal, its amplitude spectrum is visualized.
Amp1 = abs.(fft_Y1)./ nsamp
p3 = plot(
freq_vec, Amp1,
xlabel="Frequency, Hz",
ylabel="Amplitude, $dims",
title="Amplitude spectrum (channel $Nk1)",
grid=true,
legend=false
)
For the other channel, an additional conversion of the spectrum to decibels is performed.
Amp2 = abs.(fft_Y2)./nsamp
# amplitude in dB (relative to 1 mv)
Amp2_dB = 20 .* log10.(Amp2 .+ eps()) # add eps to avoid -Inf
p4 = plot(
freq_vec, Amp2_dB,
xlabel="Frequency, Hz",
ylabel="Amplitude, dB (rel. 1 $dims)",
title="Amplitude spectrum in dB (channel $Nk2)",
grid=true,
legend=false
)
We calculate the periodogram for the first channel using the function periodogram to estimate the spectral power density (SPM). The graph shows the frequency distribution of the signal energy.
power,freq = EngeePhased.Functions.periodogram(y1, fs=fs,out=:data)
p5 = plot(vec(freq), power, title="EEG periodogram ($Nk1 channel)", xlabel="Frequency, Hz", ylabel="SPM MV^2/Hz", grid=true)
Pronounced network interference (60 Hz) is observed on the spectra. Before analyzing EEG rhythms, it is required to eliminate it by filtering.
Filtering the EEG signal
The results of the spectral analysis revealed a pronounced peak at a frequency of 60 Hz, corresponding to a tip from the AC network. To eliminate it using the EngeeDSP library, a 10th-order Butterworth notch filter is used. The filter is configured to suppress the frequency band from 57 Hz to 63 Hz with a resonant frequency of 60 Hz, which eliminates network interference.
f_center = 60.0 # resonant frequency, Hz
bw = 6.0 # bandwidth, Hz
f1 = f_center - bw/2 # lower limit frequency, Hz
f2 = f_center + bw/2 # upper limit frequency, Hz
wn1 = f1 / (fs/2)
wn2 = f2 / (fs/2)
order = 10 # filter order
b, a = EngeeDSP.Functions.butter(order, [wn1, wn2], "stop")
The signals of both channels are processed by the notch filter through the function filter EngeeDSP libraries. This makes it possible to suppress network interference at 60 Hz.
filt_y1 = EngeeDSP.Functions.filter(b, a, y1)
filt_y2 = EngeeDSP.Functions.filter(b, a, y2)
After network interference is suppressed, an EEG waveform is constructed.
p6 = plot(
t, filt_y1,
xlabel = "Time, c",
ylabel = "Amplitude, MV",
title = "EEG waveform after filtering (channel $Nk1)",
grid = true,
legend = false
)
Using the function EngeeDSP.Functions.fft A fast Fourier transform is performed, then the spectrum is centered. The constructed amplitude spectrum allows you to visually assess the absence of a peak at a frequency of 60 Hz.
fft_filt_y1 = EngeeDSP.Functions.fft(filt_y1)
fft_filt_y1_shift = EngeeDSP.Functions.fftshift(fft_filt_y1)
Amp_y1_filt = abs.(fft_filt_y1_shift)./nsamp
p7 = plot(
freq_vec, Amp_y1_filt,
xlabel = "Frequency, Hz)",
ylabel = "Amplitude, $dims",
title = "The spectrum of the EEG signal after filtering (channel $Nk1)",
grid = true,
legend = false
)
To visually evaluate the operation of the notch filter, the amplitude spectra of the signal are compared before and after interference suppression. One graph shows the original spectrum and the spectrum after filtering.
p8 = plot(freq_vec, Amp1,
label = "Without filtering"
)
plot!(freq_vec, Amp_y1_filt,
label = "After filtering"
)
xlabel!("Frequency, Hz")
ylabel!("Amplitude, $dims")
title!("Amplitude spectrum ($Nk1 channel)")
Isolation of brain rhythms
A band-pass FIR filter (finite impulse response) is used to extract the alpha rhythm (8-12 Hz) from the purified EEG signal. Using the function EngeeDSP.Functions.fir1 A 200th order filter with a given bandwidth is being designed. Then the function EngeeDSP.Functions.filter applies this filter to the signal filt_y1 by isolating components in the alpha rhythm range and forming a new signal alpha_y1.
f_low = 8.0 # lower limit frequency, Hz
f_high = 12.0 # upper limit frequency, Hz
wn = [f_low, f_high] ./ (fs/2)
fir_order = 200 # filter order
b_alpha = EngeeDSP.Functions.fir1(fir_order, wn)
alpha_y1 = EngeeDSP.Functions.filter(b_alpha, [1.0], filt_y1)
An oscillogram of the selected alpha rhythm is constructed - alpha_y1.
# The alpha rhythm waveform
p9 = plot(
t, alpha_y1,
xlabel = "Time, from",
ylabel = "Amplitude, MV",
title = "Alpha rhythm (8-12 Hz), (channel $Nk1)",
legend = false,
grid = true
)
display(p9)
For further analysis, the characteristic frequency ranges of the main rhythms of brain activity are set.:
-
Delta rhythm: 0.5 – 4 Hz
-
Theta rhythm: 4 – 8 Hz
-
Alpha Rhythm: 8 – 12 Hz
-
Beta Rhythm: 13 – 30 Hz
These ranges will be used to adjust bandpass filters in order to isolate the corresponding rhythmic components from the EEG signal.
delta_band = (0.5, 4.0)
theta_band = (4.0, 8.0)
alpha_band = (8.0, 12.0)
beta_band = (13.0, 30.0)
To make it easier to distinguish different EEG rhythms, a function has been created band_fir_filter, which implements a typical bandpass FIR filtering operation using the EngeeDSP library
function band_fir_filter(sig, fs, band, order)
f1, f2 = band
wn = [f1, f2] ./ (fs/2)
b = EngeeDSP.Functions.fir1(order, wn)
y = EngeeDSP.Functions.filter(b, [1.0], sig)
return y
end
Using a previously created function band_fir_filter Three main rhythmic components are sequentially distinguished from the EEG signal: delta rhythm (0.5–4 Hz), theta rhythm (4-8 Hz) and beta rhythm (13-30 Hz).
# Filtering the first signal
delta_y1 = band_fir_filter(filt_y1, fs, delta_band, fir_order)
theta_y1 = band_fir_filter(filt_y1, fs, theta_band, fir_order)
beta_y1 = band_fir_filter(filt_y1, fs, beta_band, fir_order)
Similarly, for the second signal, the main rhythms of the brain are isolated using the same bandpass filtering function. From the signal filt_y2 The corresponding components are obtained: delta, theta, alpha and beta rhythms.
# Filtering the second signal
delta_y2 = band_fir_filter(filt_y2, fs, delta_band, fir_order)
theta_y2 = band_fir_filter(filt_y2, fs, theta_band, fir_order)
alpha_y2 = band_fir_filter(filt_y2, fs, theta_band, fir_order)
beta_y1 = band_fir_filter(filt_y2, fs, beta_band, fir_order)
Visualization of brain rhythms for the first signal.
p10 = plot(
plot(t, delta_y1,
xlabel = "Time, from",
ylabel = "Amlituda, mkV",
title = "Delta rhythm (0.5–4 Hz)",
legend = false,
grid = true),
plot(t, theta_y1,
xlabel = "Time, from",
ylabel = "Amlituda, mkV",
title = "Theta rhythm (4-8 Hz)",
legend = false,
grid = true),
plot(t, alpha_y1,
xlabel = "Time, from",
ylabel = "Amlituda, mkV",
title = "Alpha rhythm (8-12 Hz)",
legend = false,
grid = true),
plot(t, beta_y1,
xlabel = "Time, from",
ylabel = "Amlituda, mkV",
title = "Beta rhythm (13-30 Hz)",
legend = false,
grid = true),
layout = (4, 1),
size = (900, 900),
margin = 40*Plots.px
)
display(p10)
When analyzing an electroencephalogram, the signal is divided into standard frequency ranges, each of which is associated with specific brain conditions.:
-
Delta rhythm (0.5–4 Hz) – the slowest rhythm. Predominates during deep sleep.
-
Theta rhythm (4-8 Hz) – is associated with the state of drowsiness, relaxation and memory processes.
-
Alpha rhythm (8-12 Hz) – the rhythm of calm wakefulness is most pronounced when the eyes are closed.
-
Beta rhythm (13-30 Hz) – is associated with active wakefulness, concentration and mental activity.
Correlation analysis
Let's move on to the correlation data processing section. Correlation analysis is used to quantify the internal structure and frequency of the selected rhythms.
The autocorrelation function
Using the function EngeeDSP.Functions.xcorr the autocorrelation of the signal is calculated alpha_y1. The function returns two arrays:
-
R_aa— the values of the correlation coefficients for each shift. -
lag_aa— the corresponding values of the time shift in the counts.
maxlag = length(alpha_y1) - 1
R_aa, lag_aa = EngeeDSP.Functions.xcorr(alpha_y1, alpha_y1, maxlag)
tau_aa = vec(lag_aa) ./ fs
Visualization of the autocorrelation function of the alpha rhythm
R_aa_norm = R_aa ./ maximum(R_aa)
plot(
tau_aa, R_aa_norm,
xlabel = "The shift, from",
ylabel = "ACF, mkV^2",
title = "Alpha Rhythm ACF (Channel $Nk1)",
grid = true,
legend = false
)
The cross-correlation function
To analyze the relationship of alpha activity in different brain regions, the cross-correlation function (VCF) between the alpha rhythm signals of the first and second channels is calculated.
R_ab, lag_ab = EngeeDSP.Functions.xcorr(alpha_y1, alpha_y2, maxlag)
tau_ab = vec(lag_ab) ./ fs
Visualization of the mutual correlation of alpha rhythms
R_ab_norm = R_ab ./ maximum(R_ab)
plot(
tau_ab, R_ab_norm,
xlabel = "The shift, from",
ylabel = "VKF, mkV^2",
title = "The mutual correlation function of alpha rhythms",
grid = true,
legend = false
)
Conclusion
In this example, the processing of EEG data recorded in the standard EDF format using the EngeeDSP library was demonstrated. All the work is based on its functions and executed in a logical sequence.:
-
Data download: The data was read from the EDF file.
-
Spectral analysis: Using the functions
EngeeDSP.Functions.fftandEngeeDSP.Functions.fftshiftThe signals were converted to the frequency domain, which revealed the presence of network interference. -
Filtering: To suppress 60 Hz interference, a Butterworth notch filter was designed and applied using the following functions
EngeeDSP.Functions.butterandEngeeDSP.Functions.filter. -
Highlighting rhythms: The main rhythms of the brain (delta, theta, alpha, beta) were isolated using bandpass FIR filters created by the function
EngeeDSP.Functions.fir1and applied viaEngeeDSP.Functions.filter. -
Correlation analysis: To study the temporal structure and relationships of the signals, the autocorrelation (ACF) and cross-correlation (VCF) functions were calculated using
EngeeDSP.Functions.xcorr.