Обработка данных электроэнцефалографии (ЭЭГ)

В примере представлен процесс цифровой обработки сигналов электроэнцефалографии (ЭЭГ), реализованный с применением специализированной библиотеки EngeeDSP. На примере реальных данных, записанных в стандартном формате EDF, последовательно иллюстрируются ключевые этапы анализ.

Функция engee.clear() выполняет очистку рабочего пространств:

engee.clear()

Для визуализации результатов доступны два режима отрисовки графиков:

Для статической, быстрой отрисовки - используйте gr().
Для интерактивной, динамической визуализации с возможностью масштабирования и просмотра значений - используйте plotlyjs().

Выберите одну из команд, закомментировав вторую. Оба бэкенда являются взаимоисключающими, и активирован будет последний вызванный.

#gr()
plotlyjs()

Plots.PlotlyJSBackend()

Чтение EDF файла

Подключим с помощью функции include файл "edfread.jl" для чтения EDF файлов:

include("$(@__DIR__)/edfread.jl")

edfread (generic function with 1 method)

Функция edfread предназначена для чтения данных в формате EDF (European Data Format) — стандартном формате хранения биомедицинских сигналов. Структура hdr, возвращаемая этой функцией, содержит полную метаинформацию о записи:

Общие параметры записи:

ver — версия формата EDF
patientID — идентификатор пациента
recordID — идентификатор записи
startdate и starttime — дата и время начала записи
bytes — размер заголовка в байтах
records — количество блоков данных в файле
duration — длительность одного блока в секундах
ns — количество каналов в записи

Параметры каждого канала:

labels — названия каналов
transducers — тип датчиков
physicalDims — физические единицы измерений
physicalMins и physicalMaxs — минимальные и максимальные физические значения
digitalMins и digitalMaxs — минимальные и максимальные цифровые значения
prefilters — применённые при записи фильтры
samples — количество отсчётов в одном блоке для каждого канала

В качестве примера обработки биомедицинских сигналов послужат данные электроэнцефалографии (ЭЭГ). Начнём с загрузки файла S004R02.edf для изучения его структуры и построения осциллограмм сигналов.

hdr, record = edfread("$(@__DIR__)/S004R02.edf")

(EDFHeader(0.0, "X X X X", "Startdate 12-AUG-2009 X X BCI2000", "12.08.09", "16.15.00", 16896, 61, 1.0, 65, ["Fc5", "Fc3", "Fc1", "Fcz", "Fc2", "Fc4", "Fc6", "C5", "C3", "C1"  …  "Po7", "Po3", "Poz", "Po4", "Po8", "O1", "Oz", "O2", "Iz", "EDFAnnotations"], ["BCI2000", "BCI2000", "BCI2000", "BCI2000", "BCI2000", "BCI2000", "BCI2000", "BCI2000", "BCI2000", "BCI2000"  …  "BCI2000", "BCI2000", "BCI2000", "BCI2000", "BCI2000", "BCI2000", "BCI2000", "BCI2000", "BCI2000", ""], ["uV", "uV", "uV", "uV", "uV", "uV", "uV", "uV", "uV", "uV"  …  "uV", "uV", "uV", "uV", "uV", "uV", "uV", "uV", "uV", "-"], [-8092.0, -8092.0, -8092.0, -8092.0, -8092.0, -8092.0, -8092.0, -8092.0, -8092.0, -8092.0  …  -8092.0, -8092.0, -8092.0, -8092.0, -8092.0, -8092.0, -8092.0, -8092.0, -8092.0, -32768.0], [8092.0, 8092.0, 8092.0, 8092.0, 8092.0, 8092.0, 8092.0, 8092.0, 8092.0, 8092.0  …  8092.0, 8092.0, 8092.0, 8092.0, 8092.0, 8092.0, 8092.0, 8092.0, 8092.0, 32767.0], [-8092.0, -8092.0, -8092.0, -8092.0, -8092.0, -8092.0, -8092.0, -8092.0, -8092.0, -8092.0  …  -8092.0, -8092.0, -8092.0, -8092.0, -8092.0, -8092.0, -8092.0, -8092.0, -8092.0, -32768.0], [8092.0, 8092.0, 8092.0, 8092.0, 8092.0, 8092.0, 8092.0, 8092.0, 8092.0, 8092.0  …  8092.0, 8092.0, 8092.0, 8092.0, 8092.0, 8092.0, 8092.0, 8092.0, 8092.0, 32767.0], ["HP:0Hz LP:0Hz N:0Hz", "HP:0Hz LP:0Hz N:0Hz", "HP:0Hz LP:0Hz N:0Hz", "HP:0Hz LP:0Hz N:0Hz", "HP:0Hz LP:0Hz N:0Hz", "HP:0Hz LP:0Hz N:0Hz", "HP:0Hz LP:0Hz N:0Hz", "HP:0Hz LP:0Hz N:0Hz", "HP:0Hz LP:0Hz N:0Hz", "HP:0Hz LP:0Hz N:0Hz"  …  "HP:0Hz LP:0Hz N:0Hz", "HP:0Hz LP:0Hz N:0Hz", "HP:0Hz LP:0Hz N:0Hz", "HP:0Hz LP:0Hz N:0Hz", "HP:0Hz LP:0Hz N:0Hz", "HP:0Hz LP:0Hz N:0Hz", "HP:0Hz LP:0Hz N:0Hz", "HP:0Hz LP:0Hz N:0Hz", "HP:0Hz LP:0Hz N:0Hz", "HP:0Hz LP:0Hz N:0Hz"], [160, 160, 160, 160, 160, 160, 160, 160, 160, 160  …  160, 160, 160, 160, 160, 160, 160, 160, 160, 80]), [-16.0 -13.0 … 0.0 0.0; -11.0 -10.0 … 0.0 0.0; … ; 52.0 23.0 … 0.0 0.0; 12331.0 5140.0 … NaN NaN])

Следующим шагом является анализ метаданных EDF-файла. Ключевые параметры, которые мы извлечём для дальнейшей обработки, включают частоту дискретизации, количество каналов и длительность записи.

println("Версия формата:         ", hdr.ver)
println("ID пациента:            ", strip(hdr.patientID))
println("Описание записи:        ", strip(hdr.recordID))
println("Дата/время начала:      ", hdr.startdate, " ", hdr.starttime)
println("Количество каналов:     ", hdr.ns)
println("Количество записей:     ", hdr.records)
println("Общая длительность:     ", hdr.records * hdr.duration, " сек")

Версия формата:         0.0
ID пациента:            X X X X
Описание записи:        Startdate 12-AUG-2009 X X BCI2000
Дата/время начала:      12.08.09 16.15.00
Количество каналов:     65
Количество записей:     61
Общая длительность:     61.0 сек

Для наглядности также приведём ключевые характеристики первых трёх каналов.

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("    Частота дискрет.:   ", round(fs, digits=1), " Гц")
    println("    Длина сигнала:      ", length(channel_data), " отсчётов")
    println("    Длительность:       ", round(t[end], digits=2), " сек")
    println("    Амплитуда min/max:  ", round(minimum(channel_data), digits=2), " .. ", round(maximum(channel_data), digits=2))
end

  1. Fc5:
    Частота дискрет.:   160.0 Гц
    Длина сигнала:      9760 отсчётов
    Длительность:       61.0 сек
    Амплитуда min/max:  -93.0 .. 100.0
  2. Fc3:
    Частота дискрет.:   160.0 Гц
    Длина сигнала:      9760 отсчётов
    Длительность:       61.0 сек
    Амплитуда min/max:  -111.0 .. 197.0
  3. Fc1:
    Частота дискрет.:   160.0 Гц
    Длина сигнала:      9760 отсчётов
    Длительность:       61.0 сек
    Амплитуда min/max:  -165.0 .. 159.0

Далее выполняется визуализация сигналов выбранных каналов ЭЭГ. Построение осциллограмм позволяет наглядно оценить характер активности мозга, наличие артефактов.

# Выбор каналов
Nk1 = 1
Nk2 = 17
eegNk1 = record[Nk1, :]
eegNk2 = record[Nk2, :]

# Параметры сигнала
fs = hdr.samples[1] / hdr.duration    # Частота дискретизации
Tn = 60.0     # Время наблюдения
N1 = Int64(Tn * fs)    # Длина выборки
dims = hdr.physicalDims[Nk1]

# Извлечение сигналов
y1 = eegNk1[1:N1]
y2 = eegNk2[1:N1]
t = range(0, step=1/fs, length=N1)

0.0:0.00625:59.99375

p1 = plot(t, y1, 
          title="Осциллограмма ЭЭГ (канал $Nk1)", 
          xlabel="Время, с", 
          ylabel="Амплитуда, $dims", 
          grid=true)

p2 = plot(t, y2, 
          title="Осциллограмма ЭЭГ (канал $Nk2)", 
          xlabel="Время, с", 
          ylabel="Амплитуда, $dims", 
          grid=true)

Обработка ЭЭГ сигнала

Спектральный анализ ЭЭГ сигнала

Далее выполняется переход к частотному анализу данных с помощью библиотеки EngeeDSP. Для выбранных каналов применяется быстро преобразование Фурье, которое преобразует сигналы из временного представления в частотное. Полученные спектры центрируются для удобства последующей визуализации и анализа частотных составляющих.

nsamp = length(y1)
df = fs / nsamp

freq_vec = -fs/2 : df : fs/2 - df

-80.0:0.016666666666666666:79.98333333333333

fft_y1 = EngeeDSP.Functions.fft(y1)
fft_y2 = EngeeDSP.Functions.fft(y2)

# Центрированный спектр
fft_Y1 = EngeeDSP.Functions.fftshift(fft_y1)
fft_Y2 = EngeeDSP.Functions.fftshift(fft_y2)

9600-element Vector{ComplexF64}:
              -421.0 + 0.0im
 -167.45381690274462 - 4.689631736961019im
 -127.26103917783621 - 39.072676398338444im
  126.26479397906041 - 51.17359412739461im
 -173.70319686493804 - 61.30850404859848im
  -49.36668159247347 + 265.5510857894187im
 -178.81333356772302 - 54.330408165395056im
   -176.633080970405 - 68.62540635410005im
   90.29976763091145 + 176.30961270506486im
  -127.8082895730895 + 30.20710593311196im
 -124.23375485648694 + 116.3338724542482im
  13.046073715318926 + 108.89195101464429im
  -2.054308247030349 + 174.3480287311877im
                     ⋮
  -2.054308247030349 - 174.3480287311877im
  13.046073715318926 - 108.89195101464429im
 -124.23375485648694 - 116.3338724542482im
  -127.8082895730895 - 30.20710593311196im
   90.29976763091145 - 176.30961270506486im
   -176.633080970405 + 68.62540635410005im
 -178.81333356772302 + 54.330408165395056im
  -49.36668159247347 - 265.5510857894187im
 -173.70319686493804 + 61.30850404859848im
  126.26479397906041 + 51.17359412739461im
 -127.26103917783621 + 39.072676398338444im
 -167.45381690274462 + 4.689631736961019im

После получения частотного представления сигнала выполняется визуализация его амплитудного спектра.

Amp1 = abs.(fft_Y1)./ nsamp

p3 = plot(
    freq_vec, Amp1,
    xlabel="Частота, Гц",
    ylabel="Амплитуда, $dims",
    title="Амплитудный спектр (канал $Nk1)",
    grid=true,
    legend=false
)

Для другого канала выполняется дополнительное преобразование спектра в децибелы.

Amp2 = abs.(fft_Y2)./nsamp
# амплитуда в dB (относительно 1 мкВ)
Amp2_dB = 20 .* log10.(Amp2 .+ eps())   # добавим eps чтобы избежать -Inf

p4 = plot(
    freq_vec, Amp2_dB,
    xlabel="Частота, Гц",
    ylabel="Амплитуда, дБ (отн. 1 $dims)",
    title="Амплитудный спектр в дБ (канал $Nk2)",
    grid=true,
    legend=false
)

Рассчитываем периодограмму для первого канала с помощью функции periodogram, чтобы оценить спектральную плотность мощности (СПМ). График показывает распределение энергии сигнала по частотам.

power,freq = EngeePhased.Functions.periodogram(y1, fs=fs,out=:data)

p5 = plot(vec(freq), power, title="Периодограмма ЭЭГ (канала $Nk1)", xlabel="Частота, Гц", ylabel="СПМ мкВ^2/Гц", grid=true)

На спектрах наблюдается выраженная сетевая помеха (60 Гц). Перед анализом ритмов ЭЭГ требуется её устранить путём фильтрации.

Фильтрация ЭЭГ сигнала

Результаты спектрального анализа выявили выраженный пик на частоте 60 Гц, соответствующий наводке от сети переменного тока. Для её устранения с помощью библиотеки EngeeDSP применяется режекторный фильтр Баттерворта 10-го порядка. Фильтр настроен на подавление полосы частот от 57 Гц до 63 Гц с резонансной частотой 60 Гц, что позволяет устранить сетевую помеху.

f_center = 60.0            # резонансная частота, Гц
bw = 6.0                   # полоса пропускания, Гц
f1 = f_center - bw/2       # нижняя граничная частота,Гц
f2 = f_center + bw/2       # верхняя граничная частота, Гц

wn1 = f1 / (fs/2)
wn2 = f2 / (fs/2)

order = 10  # порядок фильтра

b, a = EngeeDSP.Functions.butter(order, [wn1, wn2], "stop")

2-element Vector{Any}:
 [0.4698292983238282 6.6907669963375955 … 6.6907669963375955 0.4698292983238282]
 [1.0 13.168543560078367 … 3.3802180890456617 0.22073956956346077]

Сигналы обоих каналов обрабатываются режекторным фильтром через функцию filter библиотеки EngeeDSP. Это позволяет подавить сетевую помеху на 60 Гц.

filt_y1 = EngeeDSP.Functions.filter(b, a, y1)
filt_y2 = EngeeDSP.Functions.filter(b, a, y2)

9600-element Vector{Float64}:
   1.8793171932953128
   2.01519766223684
 -11.161084001659852
 -22.223250803947558
 -17.419966985844184
  -5.378208460928653
 -15.976140468269232
 -17.489645349799662
 -12.366007559676508
 -19.02299470607572
 -20.260811820031957
   1.8435917754153994
  20.34792173591789
   ⋮
 -37.158101916029956
 -26.422698608435752
  -3.218062926172324
  -3.2807346339035037
 -10.572428911745424
  -8.136393100657186
  -9.502412037543898
 -22.962737096877152
 -23.034456665437208
 -29.53225578379856
 -45.0151022958248
 -50.55831627614732

После подавления сетевой помехи выполняется построение осциллограммы ЭЭГ-сигнала.

p6 = plot(
    t, filt_y1,
    xlabel = "Время, c",
    ylabel = "Амплитуда, мкВ",
    title  = "Осциллограмма ЭЭГ сигнала после фильтрации (канал $Nk1)",
    grid = true,
    legend = false
)

С помощью функции EngeeDSP.Functions.fft выполняется быстрое преобразование Фурье, затем спектр центрируется. Построенный амплитудный спектр позволяет визуально оценить отсутствие пика на частоте 60 Гц

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 = "Частота, Гц)",
    ylabel = "Амплитуда, $dims",
    title  = "Спектр ЭЭГ сигнала после фильтрации (канал $Nk1)",
    grid = true,
    legend = false
)

Для наглядной оценки работы режекторного фильтра выполняется сравнение амплитудных спектров сигнала до и после подавления помехи. На одном графике отображаются исходный спектр и спектр после фильтрации.

p8 = plot(freq_vec, Amp1, 
    label = "Без фильтрации"
)

plot!(freq_vec, Amp_y1_filt,
    label = "После фильтрации"
)

xlabel!("Частота, Гц")
ylabel!("Амплитуда, $dims")
title!("Амплитудный спектр (канала $Nk1)")

Выделение ритмов головного мозга

Для извлечения альфа-ритма (8-12 Гц) из очищенного ЭЭГ-сигнала применяется полосовой КИХ-фильтр (конечной импульсной характеристики). С помощью функции EngeeDSP.Functions.fir1 проектируется фильтр 200-го порядка с заданной полосой пропускания. Затем функция EngeeDSP.Functions.filter применяет этот фильтр к сигналу filt_y1, выделяя компоненты в диапазоне альфа-ритма и формируя новый сигнал alpha_y1.

f_low  = 8.0    # нижняя граничная частота, Гц
f_high = 12.0   # верхняя граничная частота, Гц

wn = [f_low, f_high] ./ (fs/2)

fir_order = 200 # порядок фильтра

b_alpha = EngeeDSP.Functions.fir1(fir_order, wn)

alpha_y1 = EngeeDSP.Functions.filter(b_alpha, [1.0], filt_y1)

9600-element Vector{Float64}:
  -1.2642257727728393e-17
  -0.0014772799882990105
  -0.005539546070600877
  -0.015279098406798688
  -0.028659168288808372
  -0.041483379140291855
  -0.045984402374753056
  -0.0444267402821369
  -0.03603825284940148
  -0.022191236447088324
  -0.007060365828623284
   0.004409769212829775
   0.015129422954688276
   ⋮
  -4.904243363704513
  -8.618264246193597
 -10.850745206304234
 -11.218257659774965
  -9.665142972402778
  -6.472970349632717
  -2.2091004739029922
   2.372356741708738
   6.4585955829921815
   9.322794813175799
  10.464216600544383
   9.691176482531793

Построена осциллограмма выделенного альфа-ритма - alpha_y1.

# Осциллограмма альфа-ритма
p9 = plot(
    t, alpha_y1,
    xlabel = "Время, с",
    ylabel = "Амплитуда, мкВ",
    title  = "Альфа-ритм (8–12 Гц), (канал $Nk1)",
    legend = false,
    grid   = true
)

display(p9)

Для дальнейшего анализа заданы характерные частотные диапазоны основных ритмов мозговой активности:

Дельта-ритм: 0.5 – 4 Гц
Тета-ритм: 4 – 8 Гц
Альфа-ритм: 8 – 12 Гц
Бета-ритм: 13 – 30 Гц

Эти диапазоны будут использоваться для настройки полосовых фильтров с целью выделения соответствующих ритмических компонент из ЭЭГ-сигнала.

delta_band = (0.5, 4.0)
theta_band = (4.0, 8.0)
alpha_band = (8.0, 12.0)
beta_band  = (13.0, 30.0)

(13.0, 30.0)

Для удобства выделения различных ритмов ЭЭГ создана функция band_fir_filter, которая реализует типовую операцию полосовой КИХ-фильтрации с использованием библиотеки EngeeDSP

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

band_fir_filter (generic function with 1 method)

С помощью ранее созданной функции band_fir_filter из ЭЭГ-сигнала последовательно выделяются три основные ритмические компоненты: дельта-ритм (0.5–4 Гц), тета-ритм (4–8 Гц) и бета-ритм (13–30 Гц).

# Фильтрация первого сигнала
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)

9600-element Vector{Float64}:
   0.0032654379858672727
   0.007422172824834053
   0.014801690250850997
   0.01165754148359467
  -0.0031911525904436255
  -0.02693768353113933
  -0.025993363024740108
  -0.014631477180530721
  -0.0074227076199889115
  -0.00622183076701969
  -0.0002976026605465516
   0.003940445906677193
   0.018343208019203788
   ⋮
   2.4802134675647785
   5.3437746345359844
   6.485972782545103
   6.051318593077311
   2.7836681642821213
  -4.1124061256879605
 -12.035547405119221
 -15.21028471340803
  -9.187943810884414
   4.28918838990696
  17.187421337929987
  20.690486728872056

Аналогично, для второго сигнала производится выделение основных ритмов мозга с использованием той же функции полосовой фильтрации. Из сигнала filt_y2 получают соответствующие компоненты: дельта-, тета-, альфа- и бета-ритмы.

# Фильтрация второго сигнала
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)

9600-element Vector{Float64}:
  -0.0008163594964668182
  -0.0011922511153292235
   0.00497267421639784
   0.012843450207834648
   0.009921533738226282
  -0.0044660220059472245
  -0.008916031800116834
  -0.004127870318159943
  -0.002099704456340565
  -0.005045801669564937
  -0.01255520255552034
  -0.02411127441779878
  -0.02587110552048713
   ⋮
  -4.483366022584483
  -1.6915932100886193
   3.3538114998521458
   8.008266102529634
   7.965540927375576
   1.1699021121243658
  -8.904363025305269
 -14.764620414594361
 -10.857934387212628
   1.3189661692662014
  13.757997160850445
  17.945613846215153

Визуализация ритмов головного мозга для первого сигнала.

p10 = plot(
    plot(t, delta_y1,
         xlabel = "Время, с",
         ylabel = "Амлитуда, мкВ",
         title  = "Дельта-ритм (0.5–4 Гц)",
         legend = false,
         grid   = true),
    plot(t, theta_y1,
         xlabel = "Время, с",
         ylabel = "Амлитуда, мкВ",
         title  = "Тета-ритм (4–8 Гц)",
         legend = false,
         grid   = true),
    plot(t, alpha_y1,
         xlabel = "Время, с",
         ylabel = "Амлитуда, мкВ",
         title  = "Альфа-ритм (8–12 Гц)",
         legend = false,
         grid   = true),
    plot(t, beta_y1,
         xlabel = "Время, с",
         ylabel = "Амлитуда, мкВ",
         title  = "Бета-ритм (13–30 Гц)",
         legend = false,
         grid   = true),
    layout = (4, 1),
    size   = (900, 900),
    margin = 40*Plots.px
)

display(p10)

При анализе электроэнцефалограммы сигнал разделяется на стандартные частотные диапазоны, каждый из которых связан с определёнными состояниями мозга:

Дельта-ритм (0.5–4 Гц) – самый медленный ритм. Преобладает во время глубокого сна.
Тета-ритм (4–8 Гц) – связан с состоянием дремоты, релаксации и процессами памяти.
Альфа-ритм (8–12 Гц) – ритм спокойного бодрствования, наиболее выражен при закрытых глазах.
Бета-ритм (13–30 Гц) – ассоциирован с активным бодрствованием, сосредоточением и умственной деятельностью.

Корреляционный анализ

Переходим к разделу корреляционной обработки данных. Для количественной оценки внутренней структуры и периодичности выделенных ритмов используется корреляционный анализ.

Автокорреляционная функция

С помощью функции EngeeDSP.Functions.xcorr вычисляется автокорреляция сигнала alpha_y1. Функция возвращает два массива:

R_aa — значения коэффициентов корреляции для каждого сдвига.
lag_aa — соответствующие величины временного сдвига в отсчетах.

maxlag = length(alpha_y1) - 1

R_aa, lag_aa = EngeeDSP.Functions.xcorr(alpha_y1, alpha_y1, maxlag)

([-2.0790449801918666e-12, -0.01431658111363992, -0.0691432963450996, -0.21981182733943047, -0.49881005826080144, -0.883644715207206, -1.2554774059748002, -1.498024805586091, -1.4948728866986274, -1.1694355430628611  …  -1.1694355430941927, -1.4948728866654877, -1.4980248056160164, -1.2554774059928664, -0.8836447151822767, -0.4988100582848648, -0.2198118273108758, -0.06914329639196168, -0.014316581095402736, -1.1522072915149645e-11], [-9599 -9598 … 9598 9599])

tau_aa = vec(lag_aa) ./ fs

19199-element Vector{Float64}:
 -59.99375
 -59.9875
 -59.98125
 -59.975
 -59.96875
 -59.9625
 -59.95625
 -59.95
 -59.94375
 -59.9375
 -59.93125
 -59.925
 -59.91875
   ⋮
  59.925
  59.93125
  59.9375
  59.94375
  59.95
  59.95625
  59.9625
  59.96875
  59.975
  59.98125
  59.9875
  59.99375

Визуализация автокорреляционной функции альфа-ритма

R_aa_norm = R_aa ./ maximum(R_aa)

plot(
    tau_aa, R_aa_norm,
    xlabel = "Сдвиг, с",
    ylabel = "АКФ, мкВ^2",
    title  = "АКФ альфа-ритма (канал $Nk1)",
    grid = true,
    legend = false
)

Взаимнокорреляционная функция

Для анализа связи альфа-активности в разных областях мозга вычисляется взаимная корреляционная функция (ВКФ) между сигналами альфа-ритма первого и второго каналов.

R_ab, lag_ab = EngeeDSP.Functions.xcorr(alpha_y1, alpha_y2, maxlag)

tau_ab = vec(lag_ab) ./ fs

19199-element Vector{Float64}:
 -59.99375
 -59.9875
 -59.98125
 -59.975
 -59.96875
 -59.9625
 -59.95625
 -59.95
 -59.94375
 -59.9375
 -59.93125
 -59.925
 -59.91875
   ⋮
  59.925
  59.93125
  59.9375
  59.94375
  59.95
  59.95625
  59.9625
  59.96875
  59.975
  59.98125
  59.9875
  59.99375

Визуализация взаимной корреляции альфа-ритмов

R_ab_norm = R_ab ./ maximum(R_ab)

plot(
    tau_ab, R_ab_norm,
    xlabel = "Сдвиг, с",
    ylabel = "ВКФ, мкВ^2",
    title  = "Взаимнокорреляционная функция альфа-ритмов",
    grid   = true,
    legend = false
)

Заключение

В данном примере была продемонстрирована обработка данных ЭЭГ, записанных в стандартном формате EDF, с использованием библиотеки EngeeDSP. Вся работа построена на её функциях и выполнена в логической последовательности:

Загрузка данных: Данные были считаны из EDF-файла.
Спектральный анализ: С помощью функций EngeeDSP.Functions.fft и EngeeDSP.Functions.fftshift сигналы были переведены в частотную область, что позволило выявить наличие сетевой помехи.
Фильтрация: Для подавления помехи 60 Гц был спроектирован и применён режекторный фильтр Баттерворта с использованием функций EngeeDSP.Functions.butter и EngeeDSP.Functions.filter.
Выделение ритмов: Основные ритмы мозга (дельта, тета, альфа, бета) были выделены с помощью полосовых КИХ-фильтров, созданных функцией EngeeDSP.Functions.fir1 и применённых через EngeeDSP.Functions.filter.
Корреляционный анализ: Для изучения временной структуры и связей сигналов была рассчитана автокорреляционная (АКФ) и взаимнокорреляционная (ВКФ) функции с помощью EngeeDSP.Functions.xcorr.

Обработка данных электроэнцефалографии (ЭЭГ)