Algorithmic Support for Spectral Processing of Cardiograms

Full Text

Компьютерное моделирование сердечнососудистой системы человека является актуальной задачей современной науки. Построение максимально полной модели позволило бы интенсифицировать научные исследования, например, за счет возможности проводить виртуальные эксперименты. В настоящее время известно множество подходов к моделированию сердечно-сосудистой системы [1- 3], так или иначе основанных на декомпозиции сложной биологической системы сердца. Все эти модели достаточно сложны и, как правило, не могут представить сердце во всем его многообразии. Поэтому заслуживает внимания модель сердца как «черного ящика», о состоянии которого судят по ограниченному набору доступных для наблюдения параметров. Этот подход широко применяется в медицинской практике, примером является электрокардиография. Существует также множество систем искусственного интеллекта [4], которые, как в представленной статье, необходимо оснастить проектом с открытым исходным кодом. Лучшие практики разработки программного обеспечения и ИТ-архитектуры были рас- смотрены на примерах инновационных стартапов [5]. Приведенные методы и приемы будут аналогичным образом представлены в следующем исследовании. Современной тенденцией функциональной диагностики является получение максимальной информации при минимальном воздействии на организм пациента. Таким методом неинвазивного исследования является метод регистрации электрокардиограмм, который широко применяется в современной медицине. Развиваются спектральные методы диагностики [6-9], основанные на преобразовании Фурье и вейв- лет-преобразовании. Как указано в [10], оценка изменчивости сердечного ритма, называемая также анализом вариабельности ритма сердца (ВРС), как клиническая практика, получила развитие с на- чала 60-х гг. ХХ в. Этому процессу способствовало применение для исследования методов математической статистики, алгоритмов обработки биологических сигналов и развитие физиологической интерпретации полученных данных. В дальнейшем ВРС выделился в самостоятельный неинвазивный метод в кардиологии. Метод активно развивается в настоящее время [11-15]. Методы анализа кардиоинтервалов разделяются на визуальные и математические методы. Математические методы делятся на три больших класса: 1) исследование общей вариабельности (ста- тистические методы или временной анализ); 2) исследование периодических составляющих ВРС (частотный анализ); 3) исследование внутренней организации динамического ряда кардиоинтервалов (автокорреляционный анализ). Для целей идентификации нарушений сер- дечного ритма методом исследования выбран спектральный (частотный) анализ кратковременных записей ЭКГ, вплоть до одного периода сердечных сокращений. Цель исследования - поиск ответов на следующие вопросы: 1) как изменяется спектр ЭКГ при нарушении сердечного ритма; 2) какова динамика изменения спектра ЭКГ (продолжительность переходного процесса); 3) какое число гармоник спектра необходимо для надежной фиксации факта нарушения сердечного ритма; 4) возможно ли использование спектра сердечных сокращений как многомерной функции состояния сердца. Спектральный анализ представляет собой чувствительный инструмент исследований, основанный на преобразовании Фурье. Это преобразование по известной для некоторого сигнала функции времени f(t) позволяет построить функцию частоты F(ω), описывающую этот же сигнал. Преобразование Фурье имеет некоторые математические ограничения для исходной функции времени, но его можно выполнить над любым физическим сигналом[1] [16; 17], в частности над кардиосигналом [9; 8]. Как известно, электрокардиографический сигнал является периодическим с периодом Т, поэтому этот сигнал можно разложить в ряд Фурье вида f t( )= a0 + (an·cos +nt bn·sin nt); n=1 2 =n · ,n (1) T где T T 1 2 a0 = T 0 f t dt a( ) ; n = T 0 f t( )·cos nt dt· ; T 2 bn = T 0 f t( )·sin nt dt· . (2) В связи с заменой аналоговых приборов цифровыми фактическим стандартом записи электрокардиографических сигналов является EDF-формат, где сигнал f t( ) представлен сеточной функцией f t( )k - набором значений напряжения в некоторые фиксированные моменты времени tk . Таким образом, в настоящей работе решается задача спектрального анализа такого сигнала. Алгоритм решения достаточно очевиден: 1) извлечение информации из EDF-файла; 2) определение периода электрокардиографического сигнала; 3) разложение этого сигнала в ряд Фурье; 4) определение необходимого количества членов ряда для представления сигнала с заданной точностью. Европейский формат данных (EDF) был введен в 1992 г. как стандарт для записей ЭЭГ и PSG (сна) [19]. В процессе применения было установлено, что формат имеет ограничения для применения в других областях (миография, вызванные потенциалы, кардиология). Основное ограничение - это то, что EDF требует непрерывности записи. Поэтому был разработан формат EDF+, в котором это ограничение снято и сохранены все остальные спецификации формата EDF без изменений. EDF+ позволяет хранить несколько нес- межных записей в одном файле. Это единственная несовместимость с EDF. Используя EDF+, все сигналы, аннотации и события, записанные за один сеанс с использованием одной системы записи, можно безопасно хранить вместе в одном файле. EDF+ также может хранить только события и аннотации без какихлибо сигналов. Такая гибкость позволяет выбирать оптимальное сочетание. Для обработки был выбран пакет MATLAB. За основу была взята функция edfread2, исходный скрипт доступен по [15]. Формат вызова функции: [hdr, dataMat] = edfread2(fname). Здесь fname - строковая константа, полное имя файла с расширением, например, ‘any_name.edf', файл должен быть в текущем каталоге или должен быть указан полный путь. Функция читает все записи из файла. Заголовок содержит общую информацию (идентификатор пациента, начало записи, конец записи и т.д.), эту информацию функция возвращает структуре hdr. Данные всех сигналов возвращаются в матрицу данных dataMat, каждый столбец этой матрицы представляет собой дискретный набор всех значений одного из сигналов. Первые 256 байт EDF-файла содержат общую информацию о самом формате, данные о пациенте, данные о записях сигналов, в том числе количестве сигналов (ns). Эти данные дополняются блоками по 256 байт на каждый сигнал, где указывается тип сигнала по характеру содержащейся информации (например, температура тела, кардиограмма и т.д.), амплитуда аналоговая и цифровая, продолжительность (duration) сигнала в секундах, количество дискретных значений. Таким образом, заголовочный блок содержит 256 + (ns * 256) байт. За заголовочным блоком следует массив дискретных значений сигналов, на каждое значение выделено 2 байта, значение дискретное, представлено целым числом со знаком. Положение отдельных записей в заголовочном блоке следующее (каждый байт представляет собой символ в коде ASCII): 8 байт: версия формата данных (по умолчанию 0); 80 байт: идентификационные данные пациента; 80 байт: идентификационные данные записи; 8 байт: дата начала записи (dd.mm.yy); 8 байт: время начала записи (hh.mm.ss); 8 байт: количество байт в заголовочном блоке; 44 байта зарезервировано; 8 байт: количество записей данных (-1, если неизвестно, каждая запись может содержать несколько сигналов); 8 байт: продолжительность записи сигнала в секундах; 4 байта: количество сигналов в записи (ns); далее множитель (ns ×) означает, что параметр записан для каждого сигнала: ns × 16 байт: ns меток сигналов (имеются в виду названия сигналов, например 'ECGV2Ref' или ‘Body temp’); ns × 80 байт: тип преобразователя (например, электрод Ag/AgCl); ns × 8 байт: единица измерения физической величины (mV - милливольты, degree - градусы Цельсия, и т.д.); ns × 8 байт: физический минимум сигнала; ns × 8 байт: физический максимум сигнала; ns × 8 байт: цифровой минимум сигнала; ns × 8 байт: цифровой максимум сигнала; ns × 80 байт: параметры фильтрации при записи сигнала (например, полоса пропускания фильтра); ns × 8 байт: количество цифровых значений (nr) при записи сигнала (произведение времени записи на частоту дискретизации; например, при записи продолжительностью 3330 с и частоте дискретизации 200 Гц здесь будет записано 666000); ns × 32 байта зарезервировано. Массив дискретных значений сигналов представляет собой матрицу nr×ns (nr строк и ns столбцов), каждый столбец - это сигнал. Таким образом, описанная функция полностью извлекает информацию из EDF-файла. Первый этап обработки сигнала - визуальный анализ кардиограммы в целом, выявление участков кардиограммы нормального ритма для данного пациента и характерных аритмий. На рис. 1 показана осциллограмма сигнала, восстановленного из EDF-файла 01_ГУСА.edf тестовой базы критических состояний 2019 г.[2] Кардиограмма на первый взгляд подобна белому шуму, однако это впечатление создается из-за большой длительности записи. Тем не менее на кардиограмме даже в этом виде можно выделить, по меньшей мере, три различных участка, обозначенные на рис. 1 как 1, 2, 3. Подробнее рассмотреть кардиограмму можно изменением масштаба времени. Нап- ример, выбором 400 точек на участке 1 получается часть сигнала, показанная на рис. 2. Здесь уже наблюдается характерный вид кардиограммы, позволяющий выделить один период сердечного ритма (точки данных отмечены символом «*»). Этот период представлен на рис. 3. Период ограничен точками с индексами 80233 и 80350. Полученные из EDF-файлов кардиограммы можно разложить в ряд Фурье с достаточно большим числом гармоник (спектр). Существует несколько алгоритмов разложения в ряд Фурье[3] [16; 17] - быстрое, дискретное и т. д. Эти алгоритмы обычно требуют фиксированного и кратного 2N числа точек на период, при выполнении этого и ряда других условий такие алгоритмы обеспечивают выигрыш в скорости вычислений. В случае кардиосигналов наблюдается значительный разброс числа точек на период из-за вариации частоты сердечных сокращений. Дискретизация кардиосигнала с 2N числа точек на период в этом случае потребует аппроксимации сигнала многочленами и вычисления значений аппроксимирующей функции на сетке 2N значений аргумента. Такие преобразования искажают исходный сигнал и нивелируют выигрыш в быстродействии. Для оценки совокупности положительных и отрицательных сторон необходимо достаточно сложное дополнительное исследование. В свете изложенного, ввиду небольшого числа точек за период и предположительно небольшого числа вычисляемых гармоник, было отдано предпочтение прямому вычислению интегралов по формуле (2) методом трапеций, который по точности и объему вычислений является компромиссом между другими методами (прямоугольников, парабол и т.д.). Расчетные формулы принимают вид 2 1 m-1 =n T · ;n a0 = T ( fk + fk+1)(· tk+1 -tk ); 2 k=1 1 m-1 an = T k=1( fk·cos +n kt fk+1·cos n kt +1)(· tk+1 -tk ); 1 m-1 bn = T k=1( fk·sin +n kt fk+1·sin n kt +1)(· tk+1 -tk ), где m - число точек на период кардиосигнала. Контроль корректности преобразования производился вычислением частичной суммы ряда Фурье по формуле (1). Разложение кардиосигнала в ряд Фурье с восстановлением функции вычислением частичной суммы ряда, построением графиков и спектра было реализовано в скрипте (все переменные в Workspace должны быть сохранены). Результаты экспериментального определения оптимального числа гармоник приведены ниже. При n = 5 гармониках на 3-м участке (см. рис. 1) показаны на рис. 4 и рис. 5. Как видно из рис. 4, пяти гармоник явно недостаточно для хорошей аппроксимации кардиосигнала. При n = 10 получается картина, представленная на графиках рис. 6, при n = 20 (рис. 7), при n = 40 (рис. 8). Сравнение рисунков показывает, что наи- лучший результат получается при числе гармоник, примерно равном половине числа значений кардиосигнала за период, в рассматриваемом случае это соответствует n = 26. Экспериментальное подтверждение данной гипотезы представлено на рис. 9. В этом случае наблюдается совпадение частичной суммы ряда Фурье и значений кардиосигнала. Значительно меньшее число гармоник, равно как и значительно большее, вызывает искажение формы восстановленного кардиосигнала. Таким образом, число гармоник при разложении в ряд Фурье должно равняться половине числа значений кардиосигнала за период. Следующий этап - сравнение спектров на различных участках кардиограммы. Кардиограммы, использованные в данной работе, были предоставлены РОХМИНЭ. Спектры представлены на рис. 10-14. Визуальный сравнительный анализ спектров кардиограмм позволяет сделать предположение, что при нормальной форме кардиосигнала в спектре присутствует значительное количество высших гармоник, при этом их амплитуда примерно равна амплитуде первой гармоники и даже превосходит ее. По мере развития аритмии происходит уменьшение амплитуды высших гармоник по сравнению с первой гармоникой, число гармоник значительной величины уменьшается. В процессе исследования рассмотрена спецификация цифровых кардиограмм в EDFформате и применение существующей функ- ции системы MATLAB для извлечения этой информации. Проведено разложение в классический ряд Фурье на одном периоде кардиосигнала. Определено, что максимальная точность описания кардиосигнала достигается при числе гармоник, равном половине числа точек дискретизации кардиосигнала в течение периода. Корректность работы разработанного для спектрального анализа скрипта проверялась восстановле- нием кардиосигнала по его спектру и сравнением с исходным сигналом. Подтверждена явная зависимость спектра от формы кардиосигнала, что позволяет сделать вывод о применимости метода спектрального анализа для идентификации нарушений сердечного ритма. При этом выявление законо- мерностей требует большого объема экспериментов и интерпретации их результатов профильными специалистами с целью диагностического использования. Тем не менее спектр электрических сигналов сердечных сокращений можно рассматривать как многомерную функцию состояния сердца. Авторы отмечают значительную трудоемкость анализа кардиосигналов «вручную». Для развития работы по спектральным исследованиям требуется существенная автоматизация процедуры обработки исходной кардиограммы и проектирование понятного конечному пользователю интерфейса. Без развития функционала программного обеспечения исследование динамики изменения спектра ЭКГ представляется затруднительным.

About the authors

Denis A. Andrikov

RUDN University

Email: andrikovdenis@mail.ru
ORCID iD: 0000-0003-0359-0897
SPIN-code: 8247-7310

Ph.D. of Technical Sciences, Associate Professor of the Department of Mechanics and Control Processes, Academy of Engineering

Moscow, Russia

Sinan V. Kurbanov

RUDN University

Author for correspondence.
Email: ya.sinan@yandex.ru
ORCID iD: 0009-0005-6632-9102

Postgraduate student of the Department of Mechanics and Control Processes, Academy of Engineering

Moscow, Russia

References

Supplementary files