Методика оценки синхронизованности ритмов регуляторных систем для анализа сигналов в реальном времени

Ранее нашими коллегами был разработан метод изучения фазовой синхронизованности, ориентированный на анализ нестационарных сигналов, в частности низкочастотных (около 0,1 Гц) колебаний в сердечно-сосудистой системе [5–8], а также методы реконструкции систем с запаздывающей обратной связью по их сложным (в том числе хаотическим) реализациям [9–15]. Применение разработанных методов для исследования взаимодействия механизмов вегетативной регуляции кровообращения показало свою фундаментальную и клиническую значимость. Были получены новые фундаментальные сведения о принципах функционального взаимодействия вегетативной регуляции сердца и дистального сосудистого русла (оцениваемых по сигналу электрокардиограммы и фотоплетизмограммы — ЭКГ и ФПГ соответственно), имеющих важное клиническое значение [16–21].

В работе [22] была показана принципиальная возможность модификации уже зарекомендовавшего себя метода оценки синхронизованности 0,1 Гц- колебаний [16–21] на случай анализа сигналов в реальном времени. При этом метод предусматривал регистрацию только сигнала ФПГ. Однако, как будет показано далее, в предложенном в работе [22] виде реализация метода в виде конечного носимого устройства фактически не может быть выполнена, так как его ресурсоемкость требует использования мощного универсального микропроцессора и внешней оперативной памяти, делая длительное ношение устройства с питанием от встроенной батареи нецелесообразным. Поэтому целью данной работы является разработка специализированной методики, которую можно реализовать на современных микроконтроллерах.

Реализация методики оценки синхронизован-ности 0,1 Гц-колебаний в ритме сердца и кровенаполнении дистального сосудистого русла в виде аппаратно-программного комплекса, включающего носимое малогабаритное устройство с автономным питанием (по типу холтеровского монитора), а также распределенную базу данных для накопления и анализа истории записей, является перспективной задачей биоинженерии. Носимое устройство при этом должно обеспечивать:

• – регистрацию сигнала ФПГ в течение не менее суток за счет встроенных автономных источников питания;

• – мобильность обследуемого в ходе его повседневной активности;

  – расчет индекса синхронизованности 0,1 Гц-колебаний в реальном времени и сохранение результатов в энергонезависимой памяти для последующего сохранения на персональном компьютере и/ или в распределенной базе данных для дальнейшего анализа.

Удовлетворение совокупности этих технических требований требует разработки специализированного метода диагностики синхронизованности 0,1 Гц-колебаний, обеспечивающего анализ сигнала ФПГ в реальном времени. Специфическими особенностями такой методики должны являться: переход на анализ данных как можно меньшей разрядности с минимизацией частоты дискретизации; использование фильтров с более короткими наборами коэффициентов; использование методов обработки сигналов, сводящих к минимуму вычисление микропроцессором громоздких нелинейных математических формул, и др.

При этом будем полагать допустимым отклонение величины индекса, рассчитываемого подобной методикой, относительно предложенного ранее метода [5–8]. В данной работе будет оценена величина такого отклонения исходя.

Разработка метода, отличающегося перечисленными техническими особенностями, позволит использовать в носимом устройстве современные микроконтроллеры, которые не отличаются повышенным быстродействием и значительными объемами встроенной оперативной памяти для хранения коэффициентов фильтров и длинных экспериментальных выборок, однако обеспечивают низкое энергопотребление и имеют малые габариты.

Работоспособность метода продемонстрирована на примере анализа экспериментальных сигналов сердечно-сосудистой системы.

Материал и методы. В работе [22] предложена методика, подразумевающая анализ 24-битных выборок, сохраняемых в 32 разрядных ячейках памяти, получаемых при частоте дискретизации 100 Гц. Из

ФПГ выделялась информация о сердечном ритме и синтезировалась реализация эквидистантной карди-оинтервалограммы (КИГ) с частотой 100 Гц. Далее ряды ФПГ и КИГ фильтровались КИХ-фильтрами, содержащими 10001 коэффициент и с помощью свертки с импульсной характеристикой линейной системы — широкополосного фазовращателя, содержащей также 10001 коэффициент, получались сопряженные по Гильберту сигналы, которые использовались для выделения рядов мгновенных фаз. Эти ряды далее прореживались до 5 Гц, и разность фаз анализировалась на предмет наличия пологих участков — интервалов синхронизации.

Легко подсчитать, что даже при реализации таких фильтров с 1-байтными коэффициентами для их хранения потребуется около 20 Кбайт энергонезависимой памяти, а при реализации операции свертки — будут постоянно заняты буфером еще порядка 160 Кбайт оперативной памяти (ОП): 4 буфера длиной по 1000132 разрядных значений. При этом расходы ОП на реализацию прочих этапов методики не превышают в каждый момент времени 3–4 Кбайт (2-3% от расхода ОП на реализацию фильтров и фазовращателя). При этом передовые микроконтроллеры с ядрами ARM, выполненные, например, по технологии picoPower technologies/lowpower/, обеспечивающей сниженное энергопотребление (на уровне около 90 мкА/МГц в активном режиме) и работу при напряжении питания от 1,6 В и удовлетворяющие требованиям к разрабатываемому носимому устройству, имеют не более 64 Кбайт встроенной ОП (семейство микроконтроллеров Atmel ATSAM4L с ядром ARM®CortexTM-M4). При этом объем ОП наиболее мощных (но уже не поддерживающих необходимые технологии сниженного энергопотребления) аналогов таких микроконтроллеров не превышает 160 Кбайт. Даже не проводя необходимое для количественного оценивания быстродействия моделирование, ясно, что и быстродействия таких процессоров (достигающего 48 MIPS) не хватит для реализации методики в реальном времени. Таким образом, оказывается, что технические особенности предложенной в работе [22] методики принципиально ограничивают реализацию этого подхода на базе существующих электронных компонент. Требуется разработка метода, учитывающего особенности архитектуры элементной базы.

В данной работе предложена методика диагностики синхронизованности 0,1 Гц-колебаний в регуляции ритма сердца и кровенаполнения дистального сосудистого русла в ходе анализа ФПГ, зарегистрированных разработанным нами ранее датчиком [23, 24]. Датчик обеспечивает регистрацию ФПГ с частотой дискретизации 120 Гц при разрядности квантования 24 бит (при указанной частоте квантования гарантируется линейность 19 старших бит выборок).

Методика включает последовательные этапы:

– предварительная обработка экспериментального сигнала ФПГ,

– получение эквидистантной 5 Гц КИГ,

– предварительная перевыборка ФПГ до 5 Гц,

– выделение исследуемых 0,1 Гц-колебаний,

– синтез сигналов мгновенных фаз и расчет в реальном времени количественной меры — суммарного процента фазовой синхронизации.

Для получения информации о вариабельности сердечного ритма по сигналу ФПГ идентифицировались моменты сердечных сокращений как вершины анакротических зубцов ФПГ (r-пики) и выделялась эк- видистантная 5 Гц КИГ. Детектирование r-пиков осуществлялось по быстро нарастающему переднему фронту пульсовой волны, который является типичным для всех основных форм ФПГ [25]. Для выделения начала переднего фронта пульсовой волны проводили фильтрацию ФПГ полосовым БИХ-фильтром Баттерворта [26] с узкой полосой пропускания от 0 до 1,5 Гц, которая обеспечила сглаживание высокочастотных составляющих пульсовых волн, связанных с артефактами и шумами радиофизической природы. Затем в скользящем окне, длительностью 1 секунда, методом скользящего среднего уменьшалось влияние тренда среднего, связанного с техническими особенностями процесса регистрации. Так как длительность переднего фронта пульсовой волны изменяется слабо, то в качестве аппроксимации расстояний между сокращениями сердца использовали последовательность расстояний между локальными минимумами фильтрованного сигнала, по которым строился неэквидистантный ряд КИГ. Эквидистантную КИГ получали с помощью кусочно-линейной аппроксимации и перевыборки до частоты 5 Гц.

Для понижения частоты дискретизации ФПГ до 5 Гц сигнал ФПГ предварительно фильтровался БИХ-фильтром нижних частот Баттерворта 1-го порядка с частотой среза 2 Гц. Такая фильтрация уменьшила влияние эффекта отражения частот при понижении частоты дискретизации сигнала ФПГ до 5 Гц, что, в свою очередь, обеспечило сжатие сигнала в 20 раз.

Далее с помощью цифрового фильтра с конечной импульсной характеристикой Хэмминга 101 значение и полосой пропускания [0,06; 0,14] Гц выделялись спектральные составляющие сигналов ФПГ и КИГ с частотой 0,1 Гц, отражающие активность исследуемых систем.

С помощью специализированного подхода строился аналитический сигнал. Для этого реализовано преобразование Гильберта в виде линейной системы, импульсная характеристика которой такова, что система представляет собой широкополосный фазовращатель, на -π/2 восстанавливались мгновенные фазы составляющих сигналов ФПГ и КИГ с частотой 0,1 Гц.

Далее с помощью модификации для анализа в реальном времени автоматической процедуры, предложенной в [5–8], по мгновенной разности фаз φ (t) исследуемых ритмов осуществлялась детекция областей фазовой синхронизации. В окне, имеющем ширину b, с помощью метода наименьших квадратов проводилась линейная аппроксимация зависимости φ (t). В результате для момента времени ti, соответствующего середине окна, получаем угловой коэффициент наклона аппроксимирующей прямой αi. Сдвигая окно на одну точку вдоль временной реализации φ (t), вычисляем угол наклона αi+1 для момента времени ti+1 и т.д. Областям фазовой синхронизации, на которых зависимость φ (t) имеет почти горизонтальный вид, должны соответствовать участки с малым значением α. Введем пороговое значение коэффициента наклона α, равное ±a, полагая, что синхронизация колебаний возможна лишь при |α|≤|a|. Будем также считать, что вторым необходимым условием синхронизации является достаточно большая протяженность области малых значений α. Ее длина должна быть не менее величины l. При описанном подходе конечная ширина скользящего окна не позволяет исследовать синхронизацию сигналов на их границах — начальном участке протяженности b/2 и конечном участке той же длины. Выбраны сле- дующие значения параметров: b=20 с, |a|=0.01, l=20 с. Суммарная длительность всех областей синхронизации выражается в процентах от длительности всей записи — суммарный процент фазовой синхронизации (S).

Анализ описанной новой методики показывает ее принципиальную реализуемость на базе современных микроконтроллеров для мобильных приложений, поддерживающих технологию picoPower. Действительно, так как после предварительной полосовой фильтрации высокоамплитудные тренды будут в основном удалены, то разрядность квантования может быть снижена до 16 бит без потери полезной информации. Тогда, с учетом децимации сигнала и соответствующего усечения импульсной характеристики фильтра и фазовращателя, расход ОП при реализации этих элементов не превысит 0.8 Кбайт (при длине ИХ 101). При этом децимация сигналов на начальных этапах обработки позволяет получить, как минимум, 20-кратный выигрыш в эффективной скорости выполнения операций (5 Гц против 100 Гц). Таким образом, разработанный алгоритм допускает реализацию на микроконтроллерах, представляющих семейство Atmel ATSAM4L, что дает принципиальную возможность использования носимого устройства с автономным питанием, обеспечивающего суточный анализ синхронизованности исследуемых ритмов в реальном времени без ограничения мобильности обследуемого.

Эксперимент. Работоспособность предложенной методики обработки сигналов в реальном времени продемонстрирована в натурном эксперименте при количественной оценке фазовой синхронизованности между подсистемами нервной регуляции деятельности ССС у 9 добровольцев. Проводилась регистрация 120 Гц 24-битного ФПГ разработанным носимым устройством для длительного мониторирования кровенаполнения конечностей [24]. Одновременно проводилась регистрация 250 Гц 16-битных временных реализаций ЭКГ и ФПГ с помощью эталонного прибора «Медиком Энцефалан-10». На рис. 1а представлены участки временных реализаций ФПГ, зарегистрированные разработанным носимым устройством (рис. 1а, тонкая линия) и эталонным прибором «Медиком Энцефалан-10» (рис. 1а, толстая линия), на рис. 1б показаны их спектры мощности.

Далее на всех этапах проводилось сопоставление результатов количественной оценки фазовой синхронизации, полученных с помощью разработанного метода по сигналу ФПГ, зарегистрированному разработанным носимым устройством, и результатов, полученных разработанным ранее методом [5–8].

На рис. 2 показаны 0,1 Гц составляющие ФПГ и КИГ, выделенные разработанным методом и прямоугольным фильтром, использующимся в эталонном методе.

На рис. 3 показаны мгновенные фазы 0,1 Гц составляющих ФПГ (рис. 3б) и КИГ (рис. 3а), выделенные разработанным методом (тонкая линия) и с помощью преобразования Гильберта (толстая линия), использующимся в эталонном методе, и разности мгновенных фаз. На рис. 3 (в, г) представлена разность мгновенных фаз КИГ (рис. 3в) и мгновенных фаз ФПГ (рис. 3г), полученных двумя методами. На рис. 3в и 3г фазы колеблются около постоянного значения и характеризуются скачками на 2pi.

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

Рис. 1. Участки временных реализаций ФПГ (а), зарегистрированные разработанным носимым устройством (тонкая линия) и эталонным прибором «Медиком Энцефалан-10» (толстая линия), их нормированные спектры мощности соответственно (б)

Рис. 2. Участки 0,1 Гц составляющих КИГ (а) и ФПГ (б), выделенные КИХ фильтром Хэмминга длиной 101 значение (тонкие линии) и с помощью прямоугольного цифрового фильтра (толстые линии) добровольца A

Рис. 3. Мгновенные фазы 0,1 Гц составляющих КИГ (а) и ФПГ (б) полученные КИХ фильтром Гильберта 101 значение (тонкие линии) и с помощью цифрового преобразования Гильберта (толстые линии) добровольца A. Разность мгновенных фаз КИГ (в) и ФПГ (г), полученных двумя описанными методами

Рис. 4. Разности мгновенных фаз, полученные с помощью предложенного ранее метода по сигналам ФПГ и ЭКГ (толстая линия) и адаптированного для работы в реальном времени метода по унивариантному сигналу ФПГ (толстая линия) для добровольцев A и B. Штриховкой слева-снизу направо-наверх отмечены участки фазовой синхронизации, детектированные по разности фаз, полученной в реальном времени. Штриховка слева-сверху направо-вниз — участки фазовой синхронизации, детектированные ранее разработанным методом. Интервалы времени, в которых пересекаются области фазовой синхрони-зованности, детектированные обоими методами, заштрихованы наклонными клетками чественная оценка фазовой синхронизации (S) для добровольца A эталонным методом составила 39%, разработанным методом 45%, для добровольца B — 65 и 62% соответственно.

В областях фазового захвата разности фаз, полученные обоими методами, ведут себя очень похоже и характеризуются колебаниями в пределах pi. Области фазовой синхронизованности, одновременно детектированные обоими методами, на рис. 4 отмечены двойной штриховкой. Видно, что длительные участки фазовой синхронизации колебаний, детектированные по разностям фаз, полученным двумя способами, перекрываются. В целях большего перекрытия участков фазовой синхронизации для разработанной методики необходимо провести коррекцию параметров методики детекции областей фазовой синхронизации.

С помощью ранее разработанного метода и предложенного метода проводилась количественная оценка показателя фазовой синхронизованности изучаемых 0,1 Гц-колебаний. Средняя ошибка составила примерно 11±4%.

Заключение. Разработана новая методика количественной оценки фазовой синхронизованности 0,1 Гц-колебаний в ритме сердца и дистальном сосудистом русле, ориентированная на обработку 120 Гц 24-битного сигнала ФПГ в реальном времени на базе носимого устройства по типу холтеровского монитора. В частности, анализ технических особенностей реализации метода показывает возможность его реализации на базе 32-битных RISC-микроконтроллерах Atmel семейства SAM4L.

Разработанная методика включает специализированные методики обработки сигналов в реальном времени:

– методику получения эквидистантной КИГ по сигналу ФПГ,

– методику сжатия ФПГ по частоте дискретизации,

– методику фильтрации сигналов,

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

– методику автоматического поиска участков фазовой синхронизации по разности фаз.

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