CC BY
149
X УПРАВЛЕНИЕ В МЕДИЦИНЕ И БИОЛОГИИ
УДК 615.471:617.7
О ТОЧНОСТИ СПЕКТРАЛЬНЫХ МЕТОДОВ РАСЧЕТА ПОКАЗАТЕЛЕЙ ВАРИАБЕЛЬНОСТИ СЕРДЕЧНОГО РИТМА
А. Н. Калиниченко,
научный сотрудник, доцент
Санкт-Петербургский государственный электротехнический университет «ЛЭТИ»
Представлены результаты экспериментального исследования ошибки определения оценок частотных показателей вариабельности сердечного ритма при использовании наиболее распространенных методов спектрального анализа: на основе быстрого преобразования Фурье и на основе модели авторегрессии. Показано, что как точность, так и статистическая устойчивость общепринятых спектральных показателей вариабельности, рассчитанных по отдельной реализации сигнала, относительно невелики, что необходимо учитывать при практическом использовании данной методики.
The results of experimental estimation of the heart rate variability spectral parameters accuracy and validity are presented. Two groups of the spectral analysis methods were considered: FFT based and parametric (based on autoregressive model). It is shown that all the common methods provide comparatively high level of the spectral parameters estimation errors.
Введение
Методика анализа вариабельности сердечного ритма (ВСР) предполагает исследование статистических и частотных свойств сигнала, образуемого последовательностью интервалов времени между смежными сокращениями сердца, относящимися к фоновому ритму. Установлено [1, 2], что анализ ВСР позволяет получить количественные индикаторы активности различных отделов вегетативной нервной системы человека, что оказывается полезным для решения многих медицинских диагностических задач. Интервалы времени между последовательными сокращениями сердца (КК-ин-тервалы) измеряются по ЭКГ (обычно — как расстояния между вершинами К-зубцов QКS-комп-лексов). Если за период регистрации ЭКГ в сигнале наблюдались нарушения сердечного ритма или помехи, вызвавшие ошибки в определении КК-ин-тервалов, то участки сигнала, соответствующие этим событиям, должны исключаться из процедуры анализа. В связи с этим при анализе ВСР в качестве входных данных принято использовать не непосредственно последовательность КК-интерва-лов, а так называемую последовательность NN-интервалов (т. е. интервалов между смежными QRS-комплексами фонового ритма или комплексами класса «норма») [1, 2]. Так как анализируемая последовательность интервалов времени не
является равномерно дискретизованным сигналом, традиционно используемым при частотном анализе, то к ней нельзя непосредственно применить общепринятые методы цифровой обработки. Поэтому на первом этапе анализа из последовательности NN-интервалов с использованием гладкой интерполяции получают равномерно дискретизованный сигнал (частота дискретизации обычно составляет 4 Гц), который рассматривается как описание восстановленной функции управления сердечным ритмом [1, 2]. Именно исследование этой функции, порождающей последовательность регистрируемых КК-интервалов, является предметом анализа ВСР. Международными стандартами по анализу ВСР [1] в качестве основных спектральных показателей рекомендованы три параметра, представляющих собой суммарные мощности в трех диапазонах частот: VLF (0,003-0,04 Гц), LF (0,04-0,15 Гц) и HF (0,15-0,4 Гц), а также ряд показателей, производных от перечисленных параметров.
Спектральные показатели ВСР рассчитываются с использованием как непараметрических (классических) методов, основанных на дискретном преобразовании Фурье (или на алгоритме ускоренного вычисления этого преобразования — быстром преобразовании Фурье (БПФ)), так и параметрических методов, в основе которых лежит модель авторегрессии — скользящего среднего [1]. При-
менение указанных методов связано с рядом ограничений и допущений, каждое из которых оказывает влияние на точность и статистическую устойчивость получаемых оценок [3].
В частности, использование методов на основе БПФ подразумевает следующие условия.
• Входной сигнал рассматривается как бесконечный и периодический с периодом повторения, равным длительности анализируемого фрагмента. Это отражается в дискретности получаемого спектра и его ограниченном разрешении по частоте, которое тем хуже, чем короче взятый для анализа фрагмент.
• Взятие фрагмента сигнала конечной длительности в неявном виде предполагает использование прямоугольного окна, что ведет к появлению эффекта спектральной утечки, т. е. растекания энергии периодических компонент сигнала по частотному диапазону.
• Представление анализируемого сигнала в дискретной форме приводит к тому, что получаемый спектр является бесконечным и периодическим.
• Достоверная оценка какой-либо частотной составляющей может быть получена лишь при условии, что анализируемый фрагмент сигнала включает, по меньшей мере, 4-5 периодов колебаний с соответствующей частотой.
Для параметрических методов, в свою очередь, характерны следующие особенности.
• Рассчитанный частотный спектр не является обратимым, т. е. по нему не может быть точно восстановлен исходный сигнал.
• Получаемая спектральная оценка существенно зависит от выбора порядка модели. Кроме того, необходимость подбора оптимального порядка исключает возможность полностью автоматического расчета спектральных оценок.
• Для моделирования низкочастотных компонент анализируемого сигнала требуется использовать достаточно большие значения порядка модели, что приводит к высокой вычислительной сложности процедур анализа.
Помимо этого для обеих групп методов характерны следующие проблемы.
• Использование методов спектрального анализа предполагает выполнение условия стационарности исследуемого процесса, что для сигнала сердечного ритма может считаться справедливым лишь условно и на ограниченных отрезках времени.
• Получаемые спектральные оценки носят случайный характер, что связано со случайностью анализируемого процесса [3].
Из перечисленных условий видно, что спектральные параметры, получаемые по фрагментам сигнала относительно небольшой длительности (один из наиболее распространенных подходов при анализе ВСР), могут служить лишь приблизительными оценками истинных частотных свойств изучаемого сигнала. С другой стороны, метод получения показателей VLF, LF и HF как суммарных мощностей в определенных частотных диапазонах
соответствует использованию усреднения по частоте, которое само по себе должно повышать статистическую устойчивость результатов [3]. Целью данной работы являлось исследование наиболее часто используемых методик расчета спектральных показателей ВСР с точки зрения оценки точности и статистической устойчивости получаемых частотных параметров.
Методы исследования
Представленные в настоящей работе исследования проводились путем моделирования сигнала сердечного ритма и алгоритмов его анализа в программной среде МАТЪАВ.
В целях исследования точности вычисления спектральных оценок была создана модель сигнала, обладающего частотными свойствами, характерными для сигнала ВСР в норме (концентрация мощности в диапазонах частот VLF, LF и HF) [1, 4]. Для этого в трех указанных диапазонах частот были сформированы три огибающие функции спектральной плотности мощности (СПМ) в форме гауссовых кривых, центр каждой из которых расположен в середине соответствующего диапазона частот, а значения на границах диапазона составляют приблизительно 1% от максимума для данной кривой. Масштаб каждой из трех кривых был выбран таким, чтобы величина площади под ней равнялась 500, что соответствует порядку наблюдаемых в норме значений показателей ВСР VLF, LF и HF, измеренных в миллисекундах в квадрате.
Выражение для отсчетов полученной функции имеет вид
„-1 3 - Лм-к.
500 1 г 2/г 1 ,1Ч
Р =-/2=ХХ—1е 1 1 ^[мс/Гц], (1)
л/2п ^0 1=1 —1
где N = 512 — выбранное число отсчетов по частоте в диапазоне частот от 0 до /тах = 0,5 Гц; А/ = /тах/N = 0,5/512 — шаг по частоте; I = 1, 2, 3 — номер диапазона частот (соответственно для диапазонов VLF, LF и HF); — £ = (2 - /11 )/6 — среднеквадратическое отклонение (СКО) для соответствующей гауссовой кривой (здесь /;1 и /2 — границы £-го диапазона).
По отсчетам модельной СПМ были восстановлены значения модулей амплитудного спектра
А =\12А/Р] [мс]. (2)
Рассчитанные значения компонент амплитудного спектра были использованы для создания модельного сигнала функции управления сердечным ритмом:
К N -1
= ЕЕ А ^[2п((А/+ф 1)+с [мс]
к=0 /=0
где К = 102 400 — выбранное число отсчетов сигнала; Т = 10 мс — интервал дискретизации (соот-
ветствует частоте дискретизации 100 Гц); С = = 1000 мс — произвольно выбранное среднее значение RR-интервала, а = rand (j) — фазовый
сдвиг для j -й синусоиды, полученный как случайное число с равномерным распределением в диапазоне от 0 до 1 (случайный фазовый сдвиг использован для предотвращения эффектов интерференции синусоид с частотами, кратными А/). В результате была получена реализация модельного сигнала, имеющая длительность 1024 с (или 17 мин и 4 с), с частотой дискретизации 100 Гц.
Последовательность RR-интервалов, соответствующая рассчитанному по формуле (3) сигналу, была сгенерирована с использованием концепции модели порождения сигнала сердечного ритма, известной как IPFM (Integral Pulse Frequency Modulation — интегральная импульсно-частотная модуляция) [5], в соответствии с описанным ниже алгоритмом.
1. n := 1 (начальное значение номера RR-интервала).
2. k := 1 (начальное значение номера отсчета сигнала).
3. r := 0 (начальное значение первого RR-интервала).
4. r := r + T (накопление значения очередного интервала).
5. Если r < sk, то переход к п. 8.
6. Rn := sk-1 (запись значения n-го RR-интервала).
7. r := sk - r (начальное значение следующего RR-интервала).
8. k := k +1 (переход к очередному отсчету сигнала).
9. Если k < K, то переход к п. 4, в противном случае — завершение процедуры.
Полученная последовательность RR-интервалов содержит 1024 значения. Далее с использованием интерполяции кубическими сплайнами была восстановлена исходная функция управления в виде равномерно дискретизованного сигнала с частотой /д = 4 Гц (частота, рекомендуемая международным стандартом по ВСР [1]). Число отсчетов полученной реализации составило 4096 = 212, что позволяет непосредственно использовать алгоритм БПФ для расчета СПМ сигнала по всей реализации целиком.
На рис. 1 показаны графики СПМ модельного сигнала и сигнала, восстановленного по последовательности RR-интервалов, а также приведены значения спектральных показателей ВСР, рассчитанные для обеих кривых СПМ. Как можно видеть, СПМ, рассчитанная по всей реализации модельного сигнала, почти совпадает с исходной СПМ модели, а значения параметров близки к заданным для модели. Некоторая потеря мощности для показателя HF (около 6%) объясняется сглаживающим эффектом интерполяции с использованием кубических сплайнов.
На рис. 2, а показан фрагмент сгенерированной последовательности RR-интервалов, а на
Р
3.5
3!
2.5 2!
1.5 1
0,5 0'
0,1 0,2 0,3 0,4 /, Гц
■ Рис. 1. Графики СПМ модельного сигнала и сигнала, восстановленного по модельной последовательности ВВ-интервалов
рис. 2, б — график равномерно дискретизованной функции, восстановленной по данному фрагменту.
Достаточно точное совпадение как СПМ, так и спектральных параметров с теоретически заданными значениями объясняется в данном случае тем, что при таком способе получения модельного сигнала все синусоиды, составляющие сигнал согласно формуле (3), целое число раз укладывались в генерируемый фрагмент, что исключает появление эффекта спектральной утечки. Таким образом, было показано, что модельный сигнал обладает заданными для модели частотными свойствами и может рассматриваться в качестве образца реализации сигнала сердечного ритма, стационарной на всем протяжении.
При практическом использовании методик оценки частотных параметров ВСР принято выполнять анализ либо по коротким фрагментам сигнала (обычно — 5 мин), получаемым в ходе процедур функциональных исследований, либо по суточным записям, получаемым в ходе холтеров-ского мониторного наблюдения [1, 2]. В данном исследовании рассматривался только первый подход.
Исследовались методы, наиболее часто применяемые при спектральном анализе, рекомендованные международным и российским стандартами по анализу ВСР [1-3].
1. Непараметрические (классические) методы на основе дискретного преобразования Фурье:
— периодограммный метод с использованием четырех типов окон — прямоугольного, Бартлетта (треугольного), фон Ханна и Тьюки (10 %-80 %-10 %);
— метод Уэлча с разбиением пятиминутного фрагмента на три перекрывающихся сегмента длительностью по 3 мин с использованием окна фон Ханна [2].
• 104, мс2/Гц
а) RR(n), мс
1200 ------------
1000
800
600
400
200
0
0
50
100
150
200
250
300 n
б) S(t), мс
Рис. 2. Фрагмент ритмограммы для модельной последовательности RR-интервалов (а) и график равномерно дискретизованной функции, восстановленной по данному фрагменту (б)
2. Параметрические методы на основе авторегрессионной модели со значениями порядка 20, 30 и 50:
— метод Берга;
— ковариационный метод;
— модифицированный ковариационный метод.
Графики использованных оконных функций
приведены на рис. 3. Выбор окон определялся следующими соображениями. Использование прямоугольного окна соответствует операции взятия конечного фрагмента без какого-либо сглаживания сигнала по краям, что ведет к существенной спектральной утечке, но обеспечивает наилучшее возможное разрешение по частоте. Окно Бартлетта (треугольное) является часто используемым простейшим сглаживающим окном. Окно фон Ханна обеспечивает высокую степень снижения спектральной утечки (правда, за счет заметного ухудшения разрешения спектра) [2]. Окно Тьюки с 10%-ми косинусными склонами представляет собой часто используемое эмпирическое окно, обеспечивающее достижение практичного компромисса между спектральной утечкой и потерей спектрального разрешения.
Для получения оценок статистической устойчивости результатов из модельной последовательности RR-интервалов случайным образом выбирались 20 фрагментов продолжительностью по 5 мин каждый. Эти фрагменты были использованы для вычисления СПМ и спектральных параметров с помощью каждого из перечисленных методов. К каждому из анализируемых фрагментов применяли следующие процедуры предварительной обработки:
— восстановление функции управления сердечным ритмом методом гладкой интерполяции с использованием кубических сплайнов;
— получение равномерно дискретизованного сигнала (1200 отсчетов с частотой / = 4 Гц);
— удаление линейного тренда с использованием метода наименьших квадратов.
Дальнейшие процедуры для методов на основе БПФ включали следующие этапы (с учетом различий между периодограммным методом и методом Уэлча):
— умножение фрагмента на оконную функцию;
— дополнение нулями до N = 2048 отсчетов;
— расчет дискретного преобразования Фурье с помощью алгоритма БПФ;
01230100480102002302002302000023010023010001230100230200002301002301000123010023010000230100230100002301002301000023010123010000230200230100012302002301000123020023010201230100230100002301002301
■ Рис. 3. Окна, использованные при исследовании непараметрических методов
а)4 Р • 104, мс2/Гц
3.5
3
2.5 2
1.5 1
0,5
0
в)
4
3.5 З
2.5 2
1.5 1
0,5
0
0
б)
Р • 104, мс2/Гц
0 0,1 0,2 0,3
Р • 104, мс2/Гц
0,4
4
3.5 З
2.5 2
1.5 1
0,5
0,5 f, Гц 0 0"
г)
0,5 f, Гц
4
3.5 З
2.5 2
1.5 1
0,5
0
Р • 104, мс2/Гц
------СПМ модели
----— средняя СПМ
для 20 фрагментов:
VLF = 509 мс2, LF = 484 мс2, HF = 472 мс2
0,5 f, Гц
0
0,1
0,2
0,З
0,4
0,5 f, Гц
Рис. 4. Графики оценок СПМ, полученных для одного и того же фрагмента модельного сигнала непараметрическими методами: а и б — периодограммный метод с прямоугольным окном и окном фон Ханна; в — метод Уэлча с окном фон Ханна; г — усредненная оценка СПМ для периодограммного метода с использованием прямоугольного окна
— расчет СПМ в диапазоне частот от 0 до 0,5 Гц:
Р/ ) = Т/ X { [X{)] + 1т2 [X/)]},
к =0
где к = 4 Гц, /к = кА/ (А/ = /д/1), к = 0, 1, 2, ..., 1/2; X (/к) и р (/к) — соответственно значения компонент разложения по Фурье и СПМ для частот /к;
— домножение СПМ на корректирующий коэффициент, учитывающий потери мощности при использовании окна:
к = 1 ’ тт 5>2<п>
1Уи1 п=1
где Ыт — размер окна в отсчетах; ю(п) — отсчеты оконной функции;
— домножение СПМ на корректирующий коэффициент, учитывающий потери мощности при дополнении нулями:
к0 = 1/ Т0 ,
где 10 — число отсчетов фрагмента сигнала до дополнения нулями.
Для параметрических методов СПМ рассчитывалась непосредственно с помощью процедур МАТЪАВ, реализующих соответствующие методы.
Показатели VLF, LF и HF вычислялись интегрированием СПМ по соответствующим диапазонам частот с использованием формулы прямоугольников [3]:
VLF = А/ХР(/к), 0,003</к <0,04Гц;
к
LF = А/Xр(/к), 0,04</к < 0,15Гц;
к
ОТ = А/Xр(/к), 0,15 </к <0,4Гц.
к
На рис. 4 и 5 приведены примеры графиков оценок СПМ, полученных для одного и того же фраг-
а)
4
3.5 3
2.5 2
1.5 1
0,5
0
Р • 104, мс2/Гц
------СПМ модели
----— СПМ фрагмента 5 мин:
VLF = 270 мс2, LF = 597 мс2, HF = 426 мс2
б) Р • 104, мс2/Гц
4
0
0,1
0,2
0,3
0,4
0,5 /, Гц
3.5 3
2.5 2
1.5 1
0,5
0
------СПМ модели
----— СПМ фрагмента 5 мин:
VLF = 322 мс2, LF = 558 мс2, ОТ = 421 мс2
в) Р • 104, мс2/Гц 4
г)
Р • 104, мс2/Гц
0,5 /, Гц
3.5 3
2.5 2
1.5 1
0,5
0
1 1 1 1 4 1 1 1 1
СПМ модели — СПМ фрагмента 5 мин: ' І' VLF = 268 мс2, LF = 598 мс2, ' ОТ = 436 мс2 1 1 1 1 ' 1 1 ‘ 1 1 -Iі • 1 1 3.5 3 2.5 2 1.5 — Iі Iі Iі 1 1 1 1 1 1 ' 1 СПМ модели средняя СПМ для 20 фрагментов: VLF = 451 мс2, LF = 510 мс2, ОТ = 475 мс2
"К > Р\ 1 1 1 Ч
! \У У . 0,5 1 \ 1 ! \У і
0
0,1 0,2 0,3 0,4 0,5 /, Гц
0,1
0,2
0,3
0,4
0,5 /, Гц
0
0
■ Рис. 5. Графики оценок СПМ, полученных для одного и того же фрагмента модельного сигнала параметрическими методами: а и б — метод Берга со значениями порядка 20 и 50; в — модифицированный ковариационный метод с порядком, равным 20; г — усредненная оценка СПМ для метода Берга с порядком модели, равным 50
мента модельного сигнала с использованием трех различных вариантов непараметрических и параметрических методов соответственно. На графиках также указаны рассчитанные по этим СПМ значения параметров VLF, ЬЕ и НЕ. Для контроля на всех графиках пунктирной линией показана истинная СПМ модельного сигнала. Кроме того, для периодограммного метода с использованием прямоугольного окна и для метода Берга с порядком модели, равным 50, были рассчитаны оценки СПМ, усредненные по 20 случайно выбранным фрагментам модельного сигнала (см. рис. 4, г, 5, г).
Можно видеть, что оценки СПМ, полученные таким способом, оказываются существенно ближе к истинной СМП модели. Однако этот результат получен благодаря стационарности процесса на всем протяжении модельного сигнала, что в реальных условиях анализа ВСР практически недостижимо.
Результаты
Анализ приведенных на рис. 4 и 5 оценок СПМ, рассчитанных для отдельного фрагмента модельного сигнала, показывает, что как погрешность
■ Таблица 1
Параметр ВСР Показатель точности Периодограммный метод Метод Уэлча с окном фон Ханна
Окно
прямоугольное Бартлетта фон Ханна Тьюки (10 % -- 80 % - 10%)
VLF Среднее значение, мс2 509 551 555 526 537
Максимальная ошибка, % 57 113 125 71 96
СКО, % 31 47 49 34 41
ЬЕ Среднее значение, мс2 484 469 468 482 468
Максимальная ошибка, % 18 39 45 25 36
СКО, % 9 25 27 12 21
НЕ Среднее значение, мс2 472 470 471 470 469
Максимальная ошибка, % 11 25 23 13 12
СКО, % 6 12 13 6 8
■ Таблица 2
Параметр ВСР Показатель точности Метод Берга Ковариационный метод Модифицированный ковариационный метод
Порядок модели Порядок модели Порядок модели
20 30 50 20 30 50 20 30 50
VLF Среднее значение, мс2 387 382 451 393 394 462 386 384 455
Максимальная ошибка, % -66 -64 -63 -66 -63 -62 -66 -63 -63
СКО, % 32 30 35 32 29 36 32 30 36
ЬЕ Среднее значение, мс2 553 539 510 554 538 510 554 539 510
Максимальная ошибка, % 28 24 21 31 29 30 27 23 23
СКО, % 9 9 9 10 11 12 9 9 10
НЕ Среднее значение, мс2 478 476 475 473 474 474 478 475 474
Максимальная ошибка, % 11 11 11 14 14 14 12 11 12
СКО, % 6 5 5 6 5 6 5 5 6
расчета параметров VLF, LF и HF, так и разброс получаемых значений достаточно велики. В целях получения статистических оценок этих погрешностей для параметров, полученных для 20 случайно выбранных фрагментов, рассчитывались следующие показатели:
— среднее значение (характеризует наличие смещения величин оценок);
— максимальная относительная ошибка (показывает наихудший случай);
— СКО (характеризует разброс значений).
В табл. 1 и 2 сведены результаты полученных оценок точности и статистической устойчивости соответственно для непараметрических и параметрических методов.
Выводы
На основании имеющихся результатов могут быть сделаны следующие выводы.
1. Все исследованные методы дают достаточно большую ошибку, особенно для значения показателя VLF, соответствующего наиболее низкочастотным компонентам сигнала сердечного ритма.
Литература
1. Heart rate variability. Standards of Measurement, Physiological interpretation and clinical use // Circulation. 1996. Vol. 93. P. 1043-1065.
2. Анализ вариабельности сердечного ритма при использовании различных электрокардиографических систем: Методические рекомендации // Вестник аритмологии. 2001. № 24. С. 65-87.
3. С. Л. Марпл-мл. Цифровой спектральный анализ и его приложения. М.: Мир, 1990.
2. Из непараметрических методов наименьшую ошибку обеспечивает применение простейшего периодограммного метода с использованием прямоугольного окна.
3. Параметрические методы позволяют достоверно оценить компоненту VLF только при относительно высоких значениях порядка модели (не менее 30).
4. Относительно большие величины разброса значений и максимальной ошибки, полученные для большинства методов, означают, что спектральные показатели ВСР, рассчитанные по отдельно взятой пятиминутной реализации сигнала, обладают невысокой статистической устойчивостью.
Основным практическим выводом из приведенных результатов исследований может служить рекомендация рассматривать рассчитанные значения спектральных показателей ВСР лишь как приблизительные оценки, дающие ориентировочное представление об интенсивности периодических процессов в соответствующих частотных диапазонах сердечного ритма.
4. McSharry P. E., Clifford G., Tarassenko L., Smith L. A.
A dynamical model for generating synthetic electrocardiogram signals // IEEE Transactions on Biomedical Engineering. 2003. Vol. 50. N 3. P. 289-294.
5. Hyndman B.W., Mohn R. K. A Pulse Modulator Model of Pacemaker Activity: Digest of the 10th Int. Conf. on Med. & Biol. Eng. Dresden, 1973. P. 223.
