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()
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()
We will connect using the function include
include("$(@__DIR__)/edfread.jl")
Function edfreadhdr
General recording parameters:
-
ver -
patientID -
recordID -
startdatestarttime -
bytes -
records -
duration -
ns
Parameters of each channel:
-
labels -
transducers -
physicalDims -
physicalMinsphysicalMaxs -
digitalMinsdigitalMaxs -
prefilters -
samples
Electroencephalography (EEG) data will serve as an example of biomedical signal processing. Let's start by uploading the file S004R02.edf
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("Дата/время начала: ", 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: ", 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="Осциллограмма ЭЭГ (канал $Nk1)",
xlabel="Time, from",
ylabel="Амплитуда, $dims",
grid=true)
p2 = plot(t, y2,
title="Осциллограмма ЭЭГ (канал $Nk2)",
xlabel="Time, from",
ylabel="Амплитуда, $dims",
grid=true)
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="Амплитуда, $dims",
title="Амплитудный спектр (канал $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="Амплитуда, дБ (отн. 1 $dims)",
title="Амплитудный спектр в дБ (канал $Nk2)",
grid=true,
legend=false
)
We calculate the periodogram for the first channel using the function periodogram
power,freq = EngeePhased.Functions.periodogram(y1, fs=fs,out=:data)
p5 = plot(vec(freq), power, title="Периодограмма ЭЭГ (канала $Nk1)", xlabel="Frequency, Hz", ylabel="СПМ мкВ^2/Гц", grid=true)
Pronounced network interference (60 Hz) is observed on the spectra. Before analyzing EEG rhythms, it is required to eliminate it by filtering.
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
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 = "Осциллограмма ЭЭГ сигнала после фильтрации (канал $Nk1)",
grid = true,
legend = false
)
Using the function EngeeDSP.Functions.fft
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 = "Амплитуда, $dims",
title = "Спектр ЭЭГ сигнала после фильтрации (канал $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!("Амплитуда, $dims")
title!("Амплитудный спектр (канала $Nk1)")
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.fir1EngeeDSP.Functions.filterfilt_y1alpha_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 = "Альфа-ритм (8–12 Гц), (канал $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
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
# 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
# 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.
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.
Using the function EngeeDSP.Functions.xcorralpha_y1
-
R_aa -
lag_aa
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 = "Shift, from",
ylabel = "ACF, mkV^2",
title = "АКФ альфа-ритма (канал $Nk1)",
grid = true,
legend = false
)
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 = "Shift, from",
ylabel = "VKF, mkV^2",
title = "The mutual correlation function of alpha rhythms",
grid = true,
legend = false
)
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.fftEngeeDSP.Functions.fftshift -
Filtering: To suppress 60 Hz interference, a Butterworth notch filter was designed and applied using the following functions
EngeeDSP.Functions.butterEngeeDSP.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.fir1EngeeDSP.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