CC BY
109
Алгоритмы фильтрации сигналов биоэлектрической природы
12 2 И.А. Тарасова , А.В. Леонова , С.А. Синютин
волгоградский государственный технический университет, г. Волгоград
2
Южный федеральный университет, факультет электроники и
приборостроения, г. Таганрог
Для выявления событий в электрофизиологических сигналах применяют различные формы их представления. Так, например, в [1] для исследования вызванных потенциалов предлагалась аппроксимация ЭЭГ синусоидальной функцией, умноженной на экспоненциально затухающий множитель. Однако данный подход приводит к расхождению ряда при отрицательных значениях аргумента (времени). Кроме того, необходимо точно знать момент начала события.
В настоящее время специалистами [2, 3] был предложен более адекватный метод представления сигналов, называемый представлением в
Метод представления в базисе Г абора рассматривается как
электроэнцефалографии. Кроме того, базис Габора активно используется при анализе изображений.
Рассмотрим синтез эффективных фильтров, основанных на дискретномвейвлет-преобразовании (ДВП). На примере удаления низкочастотной (НЧ) и высокочастотных (ВЧ) составляющих для тестовых сигналов ЭКГ и ЭЭГ синтезированы соответствующие фильтры.
Рассмотрим синтез полосового фильтра для ЭЭГ на основе ДВП. Эталонный сигнал ^Ее;с!( ^) аддитивно смешаем с помехами низкой и высокой
базисе функций Габора: gY (t) = K (у) exp
перспективный, для задач электрофизиологии, особенно
частоты, соответственно
ideal
EEGnoise
Параметры фильтра при экстремальном значении функционала:
N-1 , ч 2
Я vldeal (к Т )- vldeal (к Т ))
/ ,' wave v, EEGnoise \ sample/ EEG \ sample// к=0
Yn = argmin n-1 •
Y,er Я sideal (к -T n
/ ‘У wave/,- EEGnoise \ sample) I к=0
В результате можно получить параметры полосовых фильтров для вейвлетовДобеши, симплетов и койфлетов:
Y = {Nwawe = 4, Nzero = 2 NwaveDec = 6,WaveName = ' Db' }
Y = {Nwawe = 4 Nzero = 2 NwaveDec = 6,WaveName = ' Sym ' }
Y = {Nwawe = 2, Nzero = 2, NwaveDec = 6,WaveName =' Coif' }
Аналогично для ФНЧ:
Y = {Nwawe = 4 Nzero = 2, NwaveDec = 2,WaveName = 'Db' }
Y = {Nwawe = 5, Nzero = 2, NwaveDec = 2,WaveName = ' Sym ' }
Y = {Nwawe = 2, Nzero = 2, NwaveDec = 2,WaveName = 'Coif' }
Для полосового фильтра, применяемого при фильтрации ЭКГ набор параметров:
Y> = {Nwawe = 4 Nzero = 4 NwaveDec = 7,WaveName = 'Db' }
Y> = {Nwawe = 4 Nzero = 4, NwaveDec = 7,WaveName = 'Sym ' }
Yi = {Nwawe = 2, Nzero = 3, NwaveDec = 7,WaveName = 'Coif' }
Оптимальный набор y для сигнала ЭКГ ФНЧ для вейвлетаДобеши выглядит так:
Y> = {Nwawe = 7 Nzero = 4, NwaveDec = 4,WaveName = 'Db' }
Аналогично для симплетов и койфлетов:
Y> = {Nwawe = 7 Nzero = 4 NwaveDec = 4,WaveName = 'Sym ' }
Yi = {Nwawe = 3 Nzero = 3 NwaveDec = 3,WaveName = 'Coif' }
Для полосовых фильтров выигрыш от использования фильтров ДВП по сравнению с КИХ-фильтрами, основанных на вейвлетахДобеши и симплетов:
imprInv
для фильтров, основанных на койфлетах:
FIRFiltLen
imprInv
waveDec
л
(2)
wave
l=0 2 l = N
wave
V l=0 ^ l=N z
V z
у
Для разложения по функциям Габора в рамках быстрого алгоритма воспользуемся следующими соображениями.
В рамках реального времени полный перебор приводит к существенным вычислительным затратам, потому необходима редукция вычислительной сложности алгоритма разложения.
Предлагаемые методы уменьшения времени работы алгоритма:
1) свертка сигнала с функцией Габора методом быстрой корреляции;
2) разбиение интервала на независимые сегменты, определяемые спадом гауссианы;
3) учет ограниченности спектра и потому отсутствие необходимости полного рассмотрения на отрезке [0,® 1є] (как на всем отрезке, так и в
пределах отдельных гауссиан).
Для разложения по функциям Габора применяют модифицированный алгоритм «Поиска совпадений» MatchingPursuit (MP), который сводится к шагам, показанным на блок-схеме (Рисунок 1). Прекращаем разложение в момент достижения энергетической «выборки» сигнала на некоторую заданную величину ЗЖ, т.е. проверяем условие (/п,/п )<ЗЖ. Если условие
не выполнено, то возвращаемся на шаг поиска максимума свертки. В противном случае вычисления прекращаем.
n := 0
инициализация
переменных
fo (t)^ f (t)
n := n +1
Vy єГ: y, := argmax(gy,f)
вычисляем скалярное произведение для всех функций Г из словаря и ищем такие параметры у,, которые обеспечивают максимум свертки
а :=
( gy,,f ) ( gy ’ gy)
вычисляем нормировочный коэффициент ап по формуле, следующей из принципа наименьших квадратов
L+1 (t)^ fn (t)-angyXt)
определяем новое значение раскладываемой функции f, путем вычитания найденной элементарной д
Рисунок 1 - Применение алгоритма MatchingPursuit для разложения по функциям Г абора
Ng
В результате получаем решение в виде f (t) = ^ algy (t), гдеNGab = n.
i=0
при прямом переборе всех точек:
Имеем следующее допустимое пространство определения функций Габора, которое является сеточным и определяет вычислительные затраты
t е [0,t ], At = т , ;
ч_ sample’
U е [0,tmax], AU = Tsample;
^ е [0,2tmax1, ^ = Tsample; е [0,fösample
], Afö= 112tmax"
Естественная предельная алгоритмическая сложность разложения сигнала методом MP (поиск одной функции), определяемая параметрами АЦП по операциям умножения для всего временного отрезка Totr получается
порядка 0( 2 N0tr4nNotr), эквивалентно 0( No4tr).
При применении алгоритмов быстрой корреляции для поиска
V^. еГ;^л := argmax(gY ,f) алгоритмическая сложность снижается до
0( ^о1г1°8 ^о1г), что может быть использовано при оптимизации
вычислительных операций и зачастую позволяет уложить время обработки сигналов в рамки реального.
Обработка электрокардиосигнала, полученного в результате холтеровскогомониторирования, заключается в фильтрации, прореживании и вычисления ЧСС (вычисление ЧСС основано на алгоритме выделения R-зубца). Поскольку в общем случае на снимаемый с пациента электрокардиосигнал могут влиять различные электрические факторы (сетевая помеха 50 Гц, импульсные помехи), то существует необходимость снизить их влияние на сигнал, используемый для анализа алгоритмом выделения Я-зубца. Общий алгоритм выделения Я-зубца представлен на рисунке 1, в него не входят цифровые фильтры, т.к. их реализация зависит от разрядности и вычислительной мощности используемой аппаратной платформы.
Фильтрация НЧ и ВЧ компонент производится достаточно традиционно, а именно, с помощью цифровых рекурсивных или нерекурсивных фильтров. Такие фильтры можно синтезировать с помощью ЕШегВе81^Тоо1Вох пакета Ма1;ЬаЬ [4]. Интерес представляет нелинейный фильтр, используемый для выделения зоны QRS комплекса и подчеркивания Я зубца. Этот фильтр одновременно уменьшает амплитуду сетевой помехи 50Гц. Выражения для фильтра имеют следующий вид:
(^ - X-2 )2 (Х,._2д- X,_2а_2 )2
28 28
(1)
U, = U,-1 + (р 2U-1), (2)
edge, = р - U. (3)
где X -входные отсчеты кардио-сигнала, индекс w - является константой, w -коэффициент указывающий на смещение отсчета кардио-сигнала
относительно 0-го элемента, нулевым считается центральный отсчет во временном окне из последних (2w+3) отсчетов, данный коэффициент зависит
от частоты дискретизации сигнала, например для частоты съема ЭКГ 500 Гц, коэффициент w равен 10.
Адаптация граничных (пороговых) значений алгоритма в зависимости от изменения уровня сигнала
Рисунок 1 - Общий алгоритм выделения Я-зубца
Вычисление переменной Р[ производится по следующей формуле:
р = р -
Гі і-1
{іпВи//егс-ю-2 — тВи//егс—ю )2 (тВи£Тегс+ю-2 — тВи//егс+ю )
2
2
; (4)
где тБи//ег - массив отсчетов электрокардиосигнала, индексы с, w -являются константами, с - центральный элемент массива, w - коэффициент,
зависящий от частоты дискретизации сигнала, для частоты съема ЭКГ 500 Гц, коэффициент w равен 10 (это необходимо для эффективной режекции частоты 50 Гц).
Блок фиксирования R-зубца содержит пользовательскую процедуру, в которой может выполняться измерение RR-интервалов, расчет ЧСС (частота сердечных сокращений) и статистических параметров ЭКГ для анализа ВСР (вариабельности сердечного ритма).
Подстройка граничных параметров заключается в динамической корректировке порогового значения threshold отражающего порог фиксации R-зубца. Такая подстройка необходима для адаптации алгоритма к различным амплитудам электрокардиосигнала, так как сложно и невозможно установить среднее значение размаха сигнала ЭКГ для различных пациентов. Кроме того амплитуда сигнала может изменяться и у одного и того же пациента с течением времени. На основании проведенных расчетов, разработан алгоритм подстройки амплитуды порога для выделения R-зубца. Суть алгоритма заключается в определении величины фронта (оценки производной) входного сигнала, и установки граничного значения порога для фиксации R-зубца в процентах от этого значения, в данном случае порог установлен на 20% от максимальной величины фронта [5]. Описанные алгоритмы позволяют устройствам для проведения
холтеровскогомониторирования удовлетворять требованиям к встраиваемому программному обеспечению, к которому в свою очередь предъявляются жесткие требования ко времени выполнения, а их нарушение может привести к нарушению функционирования устройства в целом.
Работа выполнена при финансовой поддержке Министерства образования и науки РФ в рамках реализации ФЦП «Научные и научнопедагогические кадры инновационной России» на 2009-2013 гг., мероприятие 1.4, соглашение от 14.11.2012 г. № 14.A18.21.2081.
Литература
1. Жадин М. Н. Биофизические механизмы формирования электроэнцефалограммы / М. Н. Жадин. - М.: Наука, 1984. - 197 с.
2. Durka P. Matching Pursuit and Unification in EEG Analysis/ P. Durka -Norwood:Artech House, 2007. - 204 p.
3. Акулов Л.Г. Алгоритмические особенности представления электрофизиологических временных рядов в базисе функций Габора / Л.Г. Акулов, И.А. Тарасова, Ю.П. Муха // Биомедицинская радиоэлектроника. - 2010. - № 6. - C. 31-37.
4. А.Б.Сергиенко. Цифровая обработка сигналов [Текст] // СПб, Изд-во: Питер, 2002. — 608 с.
5. Синютин С.А. Анализ RR интервального ряда водителя в условиях сильных помех с помощью Wavelet преобразования [Текст] // Инженерный вестник Дона. - 2012. - №3 (16.08.12).
