[Из песочницы] Принцип анализа вариабельности сердечного ритма в MATLAB
Приветствую, Хабр! В этой публикации хочу представить свой опыт реализации алгоритма анализа ВСР человека в MATLAB. Теме анализа ВСР уделено достаточно внимания на хабре. (поиск по слову ЭКГ) однако, как мне показалось, некоторые моменты раскрыты слабо или вовсе не рассматриваются. В данной статье не уделяется много внимание объяснению явления ВСР и теории методов ее анализа. Подразумевается, что читатель подготовлен, а основной упор сделан на использование для целей анализа функций и процедур MATLAB.Анализ ВСР основан на определении последовательности RR-интервалов ЭКГ, т.е. между двумя последовательными сердечными сокращениями. Еще используют понятие NN-интервалов (normal-to-normal), то есть здесь учитываются промежутки только между нормальными сокращениями (в анализе не участвуют интервалы, записанные при нарушении сердечного ритма, а так же возникшие в результате внешних помех).
Верхний график рисунка иллюстрирует процесс детектирования длительности RR-интервалов по ЭКГ.Пациент после обследования получает кардиоинтервалограмму (КИГ), которая представляет собой совокупность RR-интервалов (нижний график рисунка, которые отображаются друг за другом. Здесь высота импульса равна длительности соответствующего RR-интервала в мс, а момент появления очередного импульса соответствует времени его детектирования.
В настоящее время существует несколько методов оценки ВСР. Среди них выделяют три группы:
методы временной области — опираются на статистические методы и направлены на исследование общей вариабельности; автокорреляционный анализ и корреляционная ритмография — интегральные показатели ВСР; методы частотной области — исследование периодических составляющих ВСР. В связи с этим различные методы анализа ВСР используют разные качественные и количественные критерии оценки. Иногда отмечается противоречие в интерпретации данных, полученных на основании разных методов оценки сердечного ритма. Поэтому актуальными являются методы суммарной оценки показателей ВСР.Р. М. Баевский предложил для комплексной оценки ритма сердца показатель активности регуляторных систем (ПАРС), который вычисляется в баллах на основании перечисленных методик. Т.е. качественный анализ ВСР должен быть проведен по всем трем методикам, а полученные данные используются для расчета показателя ПАРС.Исходными данными для анализа является запись последовательных моментов времени появления очередного R-зубца. Для анализа чаще используются 5-минутные записи, хотя есть и варианты анализа записей другой длительности.
Для построения КИГ необходимо знать не только моменты появления R-зубца на ЭКГ, но и длительность времени между двумя соседними R-зубцами. Для расчета этой длительности используется функция MATALAB diff, которая вычисляет конечные разности.
RR = diff (D_times); % Расчет RR-интервалов. Здесь D_times — моменты времени детектирования R-зубца, загруженные из файла данных times = D_times (2:(length (D_times))); % Итоговый вектор временных отметок. Далее следует определить возможные RR-интервалы, не являющиеся нормальными (получить NN-интервалы). Для этого сначала сформируем единичный вектор значений флагов для каждого RR-интервала: flags = ones (max (size (RR)),1); % Вектор флагов нормальных RR-интервалов. Теперь соответствующие «ненормальным» RR-интервалам элементы вектора flags (флаги) заменим нулями. Причем выявлять аномальные интервалы придаться вручную.Для этого построим КИГ исследуемого сигнала с помощью функции stem.
stem (RR, 'Marker','.', 'color', [0 0 0]),
title ('Уточнение синусового ритма', 'FontSize', 12);
xlabel ('Номер интервала','FontSize',12);
ylabel ('Длительность, сек.','FontSize',12);
На рисунке показан фрагмент построения с аномальным (не относящемуся к синусовому ритму) интервалом, учитывать значение которого в дальнейших расчетах нельзя.
Такие интервалы удаляются из КИГ, а участки их удаления будут интерполироваться линейным методом по величине соседних интервалов. Данный подход позволяет создать последовательность NN без нарушения общей динамики сигнала. Однако следует отметить, что результат анализа КИГ будет удовлетворительным при величине ложных интервалов не более 5–6.
Для восстановления удаленных интервалов воспользуемся функцией MATLAB для одномерной табличной интерполяции. Синтаксис: RRinterp = interp1(x, y, xi,»<метод>»).В данном случае принимаем x — время возникновения соседних ложному (или группе ложных) интервалов, y — их длительность, xi — время появления ложных интервалов, помеченных flags,»<метод>» — «linear». В результате переменная RRinterp возвращает длительность интерполированных интервалов.
На следующем шаге необходимо получить из неравномерной КИГ рис. 3 равномерно дискретизированный сигнал с помощью метода гладкой интерполяции. Для этого рекомендуется использовать интерполяцию кубическими сплайнами. Рекомендованной частотой отсчётов для равномерной дискретизации является 4 Гц.
fInt = 4; % Частота дискретизации интерполированного сигнала/
xSpline = min (times): 1/fInt: max (times); % Сетка интерполяции/
RRspline = interp1(times, RRint, xSpline, 'spline'); % Интерполяция кубическими сплайнами. times — время появления интервалов, RRint — длительность полученных NN-интервалов.
Для проведения частотного анализа обработанного сигнала КИГ RRspline необходимо избавиться от его постоянной составляющей, т.е. удалить линейный тренд. Эта достигается путем вычитания из отсчетов сигнала RRspline его среднего арифметического значения.RRdetrend = RRspline — mean (RRspline); Результат представлен на рисунке.
Для проведения частотного анализа (расчета спектральной плотности мощности) необходимо произвести процедуру дискретного преобразования Фурье (ДПФ). При этом нужно минимизировать утечку спектра и гребешковые искажения. Для этого применим метод взвешивания окном. Рекомендуется использовать окно Тьюки с 25% сглаживанием.
w = tukeywin (length (rrdetrend), 0.25); % Окно Тьюки с 25% косинусоидальными фронтами.
rrdetrend = (w'.*rrdetrend).*1000; % Умножение сигнала на оконную функцию для устранения растекания спектра и получение сигнала в мс.
А так же реализуем дополнение нулями.
Pspectr = fft (rrdetrend, 2048); % БПФ дополненное нулями до 2^11.
Результат расчета СПМ исследуемого сигнала КИГ, вычисленный на основании алгоритма БПФ, показан на рисунке.
Здесь описан периодограммный метод использования БПФ. Кроме него существуют метод Уэлча, авторегрессионный и др.По полученным данным СПМ вычисляем мощности компонент HF, LF, VLF сигнала, а так же общей мощности ТР. В соответствии с рекомендациями диапазоны будут ограничиваться следующими частотами.
BW = [0.003 0.04; … % Гц, VLF 0.04 0.15; … % Гц, LF 0.15 0.4]; % Гц, HF Для расчета мощности компонент используется функция trapz. Функция I = trapz (y) вычисляет интеграл, предполагая, что шаг интегрирования постоянен и равен единице; в случае, когда шаг h отличен от единицы, но постоянен, достаточно вычисленный интеграл умножить на h. В данном случае h = (fInt/N), где fInt = 4 Гц — частота дискретизации, N — число рассчитанных бинов БПФ.Полученных данных достаточно для расчета частотных показателей ПАРС по известным формулам.Определим геометрические параметры гистограммы. Для этого в соответствии с рекомендациями указываются диапазоны группировки кардиоинтервалов и используется команда hist.
dRR = 0.4: 0.05: 1.3; % Диапазоны группировки кардиоинтервалов
[ampRR, dRR] = hist (RRspline, dRR); % Расчет гистограммы
ampRR_proc = (ampRR.*100)./length (rrspline); % Получение данных в процентах
[yMo, yi] = max (ampRR);
Mo = dRR (yi); % Мода — наиболее часто встречающееся значение кардиоинтервала
AMo = (yMo*100)/length (RRspline); % Амплитуда моды
MxDMn = (max (RRspline) — min (RRspline)); % Вариационный размах
bar (dRR, ampRR_proc), grid on % Построение гистограммы.
title ('Гистограмма', 'FontSize', 12);
ylabel ('%','FontSize',12);
axis ([0.2 1.6 0 100])
Результат построения показан на рисунке.
Строим скатерограмму.
rr1 = RRspline (1: 2: end); % Абсцисса (х) RR (n)
rr2 = RRspline (2: 2: end); % Ордината (у) RR (n+1)
if length (rr1) ~= length (rr2)
if length (rr1) > length (rr2)
rr1 = rr1 (1: end — 1);
else rr2 = rr2 (1: end — 1);
end
end
plot (rr1, rr2,'Marker','.','LineStyle','none','Color',[0 0 0])
title ('Скатерограмма', 'FontSize', 12), grid on
axis ([0 1.6 0 1.6])
P.S.: Поскольку не являюсь медицинским специалистом в терминах допускаю возможные неточности.Описанный принцип анализа появился на основании найденных мною источников и уже представленных в открытом доступе наработок.
В силу особенностей жанра материал попытался изложить лаконично, опуская некоторые технические моменты MATLAB.Буду рад конструктивной критике и постараюсь ответить на возможные вопросы. Надеюсь, материал окажется кому-то полезен. Спасибо за внимание!
Источники Достаточно доступно и основное о методах исследованияБаевский Р.М. и др. Анализ вариабельности сердечного ритма при использовании различных электрокардиографических системСтандарты измерения и рекомендации по анализу ВСРБаза записей RR-интервалов от PhysioNetГотовые к использованию функции MATLAB: Software for Heart Rate Variability
