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

В данном примере продемонстрирован процесс загрузки, анализа и визуализации, записанных данных в формате EDF (European Data Format) - стандартном формате для хранения биомедицинских сигналов.

О формате EDF

EDF (European Data Format) - это открытый стандарт хранения и обмена многоканальными биосигналами, широко используемый в медицине и научных исследованиях. Его применяют для записи ЭЭГ, ЭКГ, ЭМГ, дыхательных сигналов, движений глаз, насыщения крови кислородом и других физиологических данных.

EDF разработан для обеспечения совместимости между оборудованием разных производителей и различным программным обеспечением. Благодаря фиксированной структуре этот формат легко обрабатывается многими инструментами анализа данных.

Использование EDF-формата

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

engee.clear()

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

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

edfread (generic function with 1 method)

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

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

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

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

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

Проверка на тестовых данных

Для проверки корректности работы алгоритмов чтения и обработки EDF-файлов используются стандартизированные тестовые данные с ресурса EDF/BDF Test Files.

Для демонстрации работы с EDF-форматом используется файл test_generator.edf. Этот файл содержит многоканальные данные для тестирования и верификации алгоритмов чтения.

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

структуру заголовка hdr с параметрами записи;
массив record, содержащий данные всех каналов.

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

(EDFHeader(0.0, "test file", "EDF generator", "02.10.08", "14.27.00", 4352, 900, 1.0, 16, ["F4", "F3", "X10", "FP2", "P4", "C4", "P3", "C3", "X9", "FP1", "F8", "F7", "DC01", "DC04", "DC03", "DC02"], ["AgAgCl", "AgAgCl", "AgAgCl", "AgAgCl", "AgAgCl", "AgAgCl", "AgAgCl", "AgAgCl", "AgAgCl", "AgAgCl", "AgAgCl", "AgAgCl", "Respiration", "SaO2", "BPM", ""], ["uV", "uV", "mV", "uV", "uV", "uV", "uV", "uV", "uV", "uV", "uV", "mV", "V", "%", "BPM", ""], [-3200.0, -3200.0, -1.6, -3200.0, -3200.0, -3200.0, -3200.0, -3200.0, -3200.0, -3200.0, -3200.0, -16.0, 0.0, -1200.0, -1200.0, -32768.0], [3200.0, 3200.0, 1.6, 3200.0, 3200.0, 3200.0, 3200.0, 3200.0, 3200.0, 3200.0, 3200.0, 16.0, 12.0, 1200.0, 1200.0, 32767.0], [-32768.0, -32768.0, -32768.0, -32768.0, -32768.0, -32768.0, -32768.0, -32768.0, -32768.0, -32768.0, -32768.0, -32768.0, 0.0, -32768.0, -32768.0, -32768.0], [32767.0, 32767.0, 32767.0, 32767.0, 32767.0, 32767.0, 32767.0, 32767.0, 32767.0, 32767.0, 32767.0, 32767.0, 32767.0, 32767.0, 32767.0, 32767.0], ["HP:0.015Hz", "HP:0.015Hz", "HP:0.015Hz", "HP:0.015Hz", "HP:0.015Hz", "HP:0.015Hz", "HP:0.015Hz", "HP:0.015Hz", "HP:0.015Hz", "HP:0.015Hz", "HP:0.015Hz", "HP:0.015Hz", "HP:0.015Hz", "HP:0.015Hz", "HP:0.015Hz", "HP:0.015Hz"], [200, 100, 200, 200, 50, 100, 200, 200, 200, 200, 200, 200, 200, 25, 25, 25]), [-799.9633783474479 -799.9633783474479 … 800.0610360875868 800.0610360875868; -800.7446402685588 -750.646219577325 … NaN NaN; … ; 60.004577706569066 60.004577706569066 … NaN NaN; 16384.0 16384.0 … NaN NaN])

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

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 пациента:           test file
Описание записи:       EDF generator
Дата/время начала:     02.10.08 14.27.00
Количество каналов:    16
Количество записей:    900
Общая длительность:    900.0 сек

Сформируем таблицу с характеристиками каждого канала: его название, физический диапазон значений, единицы измерения и частота дискретизации.

println(" № |  Канал   | Диапазон (мин/макс) | Ед.изм | Частота, Гц")
println("----------------------------------------------------------")

for ch in 1:hdr.ns
    label = hdr.labels[ch]

    # берём строку матрицы и убираем NaN (паддинг)
    row = record[ch, :]
    row = row[.!isnan.(row)]

    dataMin = round(minimum(row); digits = 2)
    dataMax = round(maximum(row); digits = 2)

    units = hdr.physicalDims[ch]
    units = units == "" ? "-" : units

    fs = round(hdr.samples[ch] / hdr.duration; digits = 2)

    println(
        lpad(ch, 2), " | ",
        rpad(label, 8), " | ",
        lpad(string(dataMin), 8), " / ",
        rpad(string(dataMax), 8), " | ",
        rpad(units, 6), " | ",
        fs
    )
end

 № |  Канал   | Диапазон (мин/макс) | Ед.изм | Частота, Гц
----------------------------------------------------------
 1 | F4       |  -799.96 / 800.06   | uV     | 200.0
 2 | F3       |  -800.74 / 800.84   | uV     | 100.0
 3 | X10      |     -0.8 / 0.8      | mV     | 200.0
 4 | FP2      |  -3200.0 / 3200.0   | uV     | 200.0
 5 | P4       |   -798.3 / 798.4    | uV     | 50.0
 6 | C4       |   -798.3 / 798.4    | uV     | 100.0
 7 | P3       |  -798.99 / 799.08   | uV     | 200.0
 8 | C3       |   -798.3 / 798.4    | uV     | 200.0
 9 | X9       |   -798.3 / 798.4    | uV     | 200.0
10 | FP1      |  -799.96 / 800.06   | uV     | 200.0
11 | F8       |  -692.74 / 692.83   | uV     | 200.0
12 | F7       |     -4.0 / 4.0      | mV     | 200.0
13 | DC01     |      0.0 / 6.0      | V      | 200.0
14 | DC04     |    100.0 / 100.0    | %      | 25.0
15 | DC03     |     60.0 / 60.0     | BPM    | 25.0
16 | DC02     |  16384.0 / 16384.0  | -      | 25.0

Для проверки корректности чтения и интерпретации данных выполним сравнение загруженных метаданных с эталонной информацией, приведённой на странице EDF/BDF Test Files.

 signal label  waveform       physical range         f         sf
 --------------------------------------------------------------------
    1    F4     block          +800uV/-800uV          1Hz       200Hz
    2    F3     triangle       +800uV/-800uV          3Hz       100Hz
    3    X10    impulse        +0.8mV/-0.8mV          5Hz       200Hz
    4    FP2    noise          +3200uV/-3200uV        -Hz       200Hz
    5    P4     sine           +800uV/-800uV          1Hz        50Hz
    6    C4     sine           +800uV/-800uV          2Hz       100Hz
    7    P3     sine           +800uV/-800uV          3Hz       200Hz
    8    C3     sine           +800uV/-800uV          4Hz       200Hz
    9    X9     sine           +800uV/-800uV          8Hz       200Hz
   10    FP1    sine           +800uV/-800uV         16Hz       200Hz
   11    F8     sine           +800uV/-800uV         32Hz       200Hz
   12    F7     triangle       +4mV/-4mV              5Hz       200Hz
   13    DC01   sine square    +6V/-0V                5Hz       200Hz
   14    DC04   DC             +100%                  -Hz        25Hz
   15    DC03   DC             +60BPM                 -Hz        25Hz
   16    DC02   DC             +16384                 -Hz        25Hz

Таким образом, загруженные данные соответствуют описанию тестового файла, что подтверждает корректность работы функции edfread.

Построим осциллограммы первых 5 секунд многоканальной записи, чтобы сравнить полученные графики с тестовым изображением.

t_max = 5.0     # ограничение по времени, с
nchan = size(record, 1) # кол-во каналов

plt = plot(
    layout = (nchan, 1),
    size   = (1000, 200*nchan),
    margin = 20*Plots.px
)

for ch in 1:nchan
    # Частота дискретизации канала
    fs = hdr.samples[ch] / hdr.duration

    # Максимальное число отсчетов для канала
    n_time = min(Int(round(t_max * fs)), hdr.samples[ch] * hdr.records)

    # Время и сигнал
    t = (0:n_time-1) ./ fs
    y = record[ch, 1:n_time]

    # Единицы измерения
    units = hdr.physicalDims[ch]
    units = units == "" ? "-" : units

    plot!(
        plt[ch],
        t, y,
        label  = "$(hdr.labels[ch])",
        xlabel = "Время, с",
        ylabel = units,
        legend = :topright
    )
end

display(plt)

Для тестирования работы функции edfread с расширенным форматом EDF+ загрузим файл test_generator_2.edf.

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

(EDFHeader(0.0, "X X X X", "Startdate 10-DEC-2009 X X test_generator", "10.12.09", "12.44.02", 3328, 600, 1.0, 12, ["squarewave", "ramp", "pulse", "ECG", "noise", "sine1Hz", "sine8Hz", "sine85Hz", "sine15Hz", "sine17Hz", "sine50Hz", "EDFAnnotations"], ["", "", "", "", "", "", "", "", "", "", "", ""], ["uV", "uV", "uV", "uV", "uV", "uV", "uV", "uV", "uV", "uV", "uV", ""], [-1000.0, -1000.0, -1000.0, -1000.0, -1000.0, -1000.0, -1000.0, -1000.0, -1000.0, -1000.0, -1000.0, -1.0], [1000.0, 1000.0, 1000.0, 1000.0, 1000.0, 1000.0, 1000.0, 1000.0, 1000.0, 1000.0, 1000.0, 1.0], [-32768.0, -32768.0, -32768.0, -32768.0, -32768.0, -32768.0, -32768.0, -32768.0, -32768.0, -32768.0, -32768.0, -32768.0], [32767.0, 32767.0, 32767.0, 32767.0, 32767.0, 32767.0, 32767.0, 32767.0, 32767.0, 32767.0, 32767.0, 32767.0], ["", "", "", "", "", "", "", "", "", "", "", ""], [200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 51]), [99.99237048905174 99.99237048905174 … -99.96185244525817 -99.96185244525817; -99.96185244525817 -98.95475700007621 … 98.0086976424812 98.98527504386978; … ; 99.99237048905174 0.015259021896781633 … -99.96185244525817 0.015259021896781633; 0.37633325703822385 0.1568780041199359 … NaN NaN])

Аналогично предыдущему файлу, извлекаем и анализируем ключевые параметры записи.

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 10-DEC-2009 X X test_generator
Дата/время начала:      10.12.09 12.44.02
Количество каналов:     12
Количество записей:     600
Общая длительность:     600.0 сек

Сформируем таблицу с характеристиками каждого канала: его название, физический диапазон значений, единицы измерения и частота дискретизации.

println(" № |   Канал    | Диапазон (мин/макс) | Ед.изм | Частота, Гц")
println("------------------------------------------------------------")

for ch in 1:hdr.ns-1
    label = hdr.labels[ch]

    # берём строку матрицы и убираем NaN (паддинг)
    row = record[ch, :]
    row = row[.!isnan.(row)]

    dataMin = round(minimum(row); digits = 2)
    dataMax = round(maximum(row); digits = 2)

    units = hdr.physicalDims[ch]
    units = units == "" ? "-" : units

    fs = round(hdr.samples[ch] / hdr.duration; digits = 2)

    println(
        lpad(ch, 2), " | ",
        rpad(label, 10), " | ",
        lpad(string(dataMin), 8), " / ",
        rpad(string(dataMax), 8), " | ",
        rpad(units, 6), " | ",
        fs
    )
end

 № |   Канал    | Диапазон (мин/макс) | Ед.изм | Частота, Гц
------------------------------------------------------------
 1 | squarewave |   -99.96 / 99.99    | uV     | 200.0
 2 | ramp       |   -99.96 / 98.99    | uV     | 200.0
 3 | pulse      |     0.02 / 99.99    | uV     | 200.0
 4 | ECG        |   -17.32 / 61.48    | uV     | 200.0
 5 | noise      |     0.02 / 98.99    | uV     | 200.0
 6 | sine1Hz    |   -99.96 / 99.99    | uV     | 200.0
 7 | sine8Hz    |   -99.78 / 99.81    | uV     | 200.0
 8 | sine85Hz   |   -99.96 / 99.99    | uV     | 200.0
 9 | sine15Hz   |   -99.96 / 99.99    | uV     | 200.0
10 | sine17Hz   |   -99.96 / 99.99    | uV     | 200.0
11 | sine50Hz   |   -99.96 / 99.99    | uV     | 200.0

Для проверки корректности чтения и интерпретации данных выполним сравнение загруженных метаданных с эталонной информацией, приведённой на странице EDF/BDF Test Files.

 signal label/waveform  amplitude    f       sf
---------------------------------------------------
   1    squarewave        100 uV    0.1Hz   200 Hz
   2    ramp              100 uV    1 Hz    200 Hz
   3    pulse             100 uV    1 Hz    200 Hz
   4    ECG               100 uV    1 Hz    200 Hz
   5    noise             100 uV    - Hz    200 Hz
   6    sine 1 Hz         100 uV    1 Hz    200 Hz
   7    sine 8 Hz         100 uV    8 Hz    200 Hz
   8    sine 8.5 Hz       100 uV    8.5Hz   200 Hz
   9    sine 15 Hz        100 uV   15 Hz    200 Hz
  10    sine 17 Hz        100 uV   17 Hz    200 Hz
  11    sine 50 Hz        100 uV   50 Hz    200 Hz

Таким образом, загруженные данные соответствуют описанию тестового файла, что подтверждает корректность работы функции edfread.

Для проверки корректности функционирования загрузки и интерпретации данных EDF+ построим графики первых 10 секунд записи.

t_max = 10.0     # ограничение по времени, с
nchan = size(record, 1)-1

plt = plot(
    layout = (nchan, 1),
    size   = (1000, 200*nchan),
    margin = 30*Plots.px
)

for ch in 1:(nchan)
    # Частота дискретизации канала
    fs = hdr.samples[ch] / hdr.duration

    # Максимальное число отсчетов для канала
    n_time = min(Int(round(t_max * fs)), hdr.samples[ch] * hdr.records)

    # Время и сигнал
    t = (0:n_time-1) ./ fs
    y = record[ch, 1:n_time]

    # Единицы измерения
    units = hdr.physicalDims[ch]
    units = units == "" ? "-" : units

    plot!(
        plt[ch],
        t, y,
        label  = "$(hdr.labels[ch])",
        xlabel = "Время, с",
        ylabel = units,
        legend = :topright
    )
end

display(plt)

Заключение

В данном примере был рассмотрен принцип работы с данными в форматах EDF и EDF+. На примере тестовых файлов (test_generator.edf и test_generator_2.edf), взятых с EDF/BDF Test Files.

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