CC BY
19УДК 004.04, 612.171.1
Адаптация фильтра нелокального усреднения для усиления звуков тонов сердца на фонокардиограммах плода и человека
Я. В. Костелей1, Д. С. Жданов, И. Г. Боровской
Представлено конструирование и обоснование работы фильтра pNLF, основанного на модификации фильтра нелокального усреднения, для решения задачи усиления звуков тонов сердца на фонокардиограммах плода и человека для последующей сегментации тонов и вычисления частоты сердечных сокращений, актуальных для фетального мониторинга. Приведены особенности подбора параметров фильтра и перспективы его оптимизации для конкретных задач. Предоставлены данные и результаты оценки фильтрации методом pNLF на фонокардиограммах плода и взрослого человека, зарегистрированных с использованием разных акустических датчиков.
Ключевые слова: фонокардиограмма, нелокальное усреднение, сегментация тонов сердца.
1. Введение
Мониторинг частоты сердечных сокращений человека применяется для диагностики патологических состояний и наблюдения за состоянием человека, контроля ведения спортивной деятельности, выявления эмоциональных всплесков для анализа поведения и интерактивных интерфейсов. В настоящее время для неинвазивной регистрации сердечных сокращений используются следующие технологии: электрокардиография, кардиотокография, плетизмография и фонокардиография.
Фонокардиография - это метод регистрации акустического сигнала сердца с его последующим анализом. Метод актуален ввиду того, что приборы, реализующие его, характеризуются простотой исполнения, отсутствием необходимости применения специализированных средств (проводящие гели, комплекты электродов и т.д.), применимостью для захвата сигналов артерий и сердцебиений плода, а также тем, что их применение не оказывает воздействия на организм пациента. Так же как и данные электрокардиографии и кардиотоко-графии, фонокардиограммы могут быть использованы для анализа состояния сердечнососудистой системы пациента по графику частоты сердечных сокращений, полученному при обработке и сегментации данных сигналов. Стоит отметить, что такое направление актуально для мониторинга состояния плода и матери.
В свою очередь, для задачи расчета частоты сердечных сокращений фонокардиограмма имеет ряд особенностей, затрудняющих детектирование отдельных тонов сердца. На ней может присутствовать множество фоновых сигналов: физиологические звуки (дыхание, икота, звуки пищеварения), звуки внешней среды (разговоры, шумы техники, музыка) и трение
1 Исследования проводились при финансовой поддержке Российского фонда фундаментальных исследований (РФФИ), грант 18-31-20012 мол_а_вед «Разработка математической модели и вычислительных алгоритмов для решения задачи адаптивной идентификации и подсчета сердечных сокращений на фонокардиограмме при их неопределенных характеристиках в условиях высокого уровня помех».
датчика по поверхности кожи. Амплитуда сердцебиений зависит как от постоянных параметров (толщина жировой прослойки и характеристики тканей организма), так и от изменяющихся во время проведения регистрации сигнала: выраженность удара сердца, степень прижатия датчика, удаление объекта регистрации (например, положение плода в утробе). Также формации сердечного биения (совокупность колебаний, входящих в один удар сердца, имеющая закономерности в количестве, выраженности амплитуды и последовательности колебаний) могут меняться не только для разных приборов, но и от фонокардиограммы к фо-нокардиограмме на одном приборе в зависимости от точки захвата сигнала и особенностей сердечно-сосудистой системы. На рис. 1 представлены различные формации тонов сердца для разных точек размещения акустического датчика относительно сердца.
а)
б)
в)
г)
Рис. 1. Примеры разных формаций одного сердечного биения на фонокардиограмме: (а)-(г) - участки фонокардиограмм, которые были захвачены акустическим датчиком с разных точек
захвата
Для решения задачи первичной сегментации фонокардиограммы на тоны сердца необходимо провести ее предварительную обработку, направленную на улучшение выраженности сердцебиений. При этом стоит отметить, что сердцебиения человека обычно регистрируются в области груди и более выражены по сравнению с сигналами сердцебиений плода, регистрируемыми с поверхности живота беременной. Наиболее часто используемым методом обработки фонокардиограмм является вейвлет-преобразование (материнские вейвлеты Морле, Койфлет 4, Добеши 4, Symlet 8), а целью настройки фильтра при использовании вейвлетов является выбор семейства вейвлетов и их порядка [1-3]. При этом, как было сказано ранее, формации комплекса первого тона могут меняться внутри даже одного непрерывного сигнала, что требует «аккуратного» подхода к подбору параметров фильтрации. Другим подходом в фильтрации акустического сигнала сердцебиений является разложение его на компоненты, например, разложение на эмпирические моды (преобразование Гилберта-Хуанга) [4-5] и анализ независимых компонент [6]. В данных методах необходимо подбирать количество уровней разложения и производить выбор уровня разложения для анализа. Также для решения задачи используют различные фильтры предварительной обработки, например, фильтр Калмана и FIR/ПR-фильтрация [7-8]. При этом стоит отметить, что задача фильтрации такого типа сигналов до сих пор остается актуальной, так как амплитудно-частотные параметры
реальных фонокардиограмм характеризуются изменениями даже внутри одной непрерывной регистрации, а уровень фонового шума и помех высок. Поэтому разработка новых подходов к адаптивной фильтрации фонокардиограмм является также актуальной.
В работе предлагается провести конструирование фильтра, настройка которого производится с помощью физиологических параметров сердцебиений и характеристик датчика, который позволяет повысить выраженность первого тона сердца на фонокардиограммах плода и человека с целью их дальнейшей сегментации. Конструируемый фильтр будет основан на методе нелокального усреднения [9].
2. Представление звуков сердцебиений на фонокардиограмме и их характеристики
Исходным сигналом для обработки выступает акустический сигнал сердцебиений плода или взрослого человека. Для построения фильтра необходимо учитывать физиологические параметры полезного сигнала и его особенности, которые зависят от характеристик датчика. Фонокардиограмма в простейшем представлении состоит из совокупности последовательно идущих тонов сердца: звука первого тона S1 и звука второго тона S2 (рис. 2).
Рис. 2. Представление сердцебиения на фонокардиограмме
Обычно на фонокардиограмме амплитуда S1 по сравнению с амплитудой S2 более выражена, но встречаются и обратные ситуации: S2 более выражен по амплитуде на протяжении всего сигнала, либо внутри одного непрерывного сигнала «лидирующая» амплитуда может меняться от S2 к S1 и обратно. Также один из тонов может проявляться и пропадать в рамках одной регистрации.
В терминологии электрокардиографии расстояние между двумя последовательно идущими первыми тонами сердца называется ЯЯ-интервалом, а пикообразное наиболее выраженное возрастание сигнала электрокардиограммы является частью комплекса S1 и называется сегментом Я. Обычно на фонокардиограммах в качестве сегмента Я принимается начало формации S1 [10], который называется сегментом S. Также при подсчете ЯЯ-интервалов допустимо фиктивно в качестве точки Я устанавливать выраженное колебание в формации S1 (рис. 2). ЯЯ-интервал - важный объект для выделения, так как является показателем, который позволяет рассчитать частоту сердечных сокращений.
Частота сердечных сокращений (ЧСС) в течение разных этапов взросления человека в норме колеблется в различных диапазонах [11]. Например, новорожденные в норме имеют ЧСС от 100 до 170 уд/мин. Дети старше 10 лет и взрослые - от 60 до 100 уд/мин. Спортсмены - от 40 до 60 уд/мин. Данные интервалы указаны для организма без патологий и в спокойном состоянии. Для сердцебиений плода с учетом патологических состояний можно вы-
делить интервал ЧСС от 100 до 200 уд/мин. Таким образом, для плода характерен разброс RR-интервала от 300 до 600 мс. Для среднестатистического взрослого - от 600 мс до 1000 мс. Это позволяет определить границы допустимого расстояния между двумя отдельными сердцебиениями.
Как было сказано выше, для построения графика ЧСС можно принять в качестве сегмента R определенное выраженное по амплитуде колебание комплекса S1, которое визуально детектируемо на последующем комплексе S1. Необходимость точного определения колебания, на котором находится сегмент R, связано с точностью определения ЧСС в точке, так как отклонение в 10-20 мс вносит значительную погрешность в расчете ЧСС. Для построения фильтра же допустимо более свободное представление RR-интервала, поэтому можно искусственно внести предположение о том, что сегмент R находится в середине комплекса S1.
Длина комплекса S1 составляет L мс, где L в соответствии с анализом фонокардиограмм, проведенным научным коллективом, и теории электрокардиографии составляет 120-150 мс [12]. В качестве амплитуды A первого тона S1 будет выступать амплитуда колебания комплекса S1, которая наиболее выражена (максимальна) по амплитуде (рис. 2).
В качестве другого важного параметра выступают частоты колебаний [tomin, tomax], которые соответствуют колебаниям, входящим в комплекс S1. Для построения фильтра диапазон частот укажет на период колебания сердцебиений T, который может быть использован в качестве окна фильтрации. Для фонокардиограмм значения диапазона частот [tomin, tomax] полезного сигнала сердцебиений будут зависеть от совокупности факторов: возраст, физические особенности и особенности аппаратной предварительной фильтрации, источник сигнала (сердце, артерия, аорта), место захвата сигнала относительно источника звука. Подавление колебаний, частоты которых находятся вне диапазона [tomin, tomax], позволяет исключить или заглушить часть «шумовых» составляющих, например, голос или часть физиологических звуков. В качестве метода предварительной обработки сигнала авторским коллективом использовался метод, описанный в его более ранних работах [13].
3. Модификация фильтра нелокального усреднения
Фильтр нелокального усреднения (non-local means, NLM) был создан для шумоподавления на фотографиях и отличался от «локального» подхода тем, что доля текущего анализируемого пикселя в усреднении значения целевого пикселя определяется мерой, которая указывает, насколько окрестность текущего пикселя подобна окрестности целевого пикселя [9] :
1
NL(u)( x) = —j exp(-L ( x)
(Ga * |u(X +.) -u(Y + .)|2)(0)
h
)u(Y )dY,
L(x) = }exp(- (Ga * lu(X + ) - u(Z + Л )(0))dZ,
h2
(Ga * \u(X +.) - u(Y + .)|2)(0) = { Ga (t) |u(x +1) - u(y +1)|2 dt,
Я2
(1)
где х - фильтруемый элемент массива и, х е Я , С(х) - нормализирующее значение, И - параметр фильтрации, Оа - ядро фильтра Гаусса со среднеквадратичным отклонением а. Фильтруемое значение в х устанавливается в (1) как среднее значение от всех пикселей изображения, чья гауссова окрестность подобна окрестности х.
Далее фильтр был адаптирован и стал использоваться для обработки одномерных сигналов. При выполнении работ особый интерес представляла обработка электрокардиограм-
мы [14-17] и фонокардиограммы [18]. Ниже представлен вид фильтра для одномерного вида из работы [14]:
Ф) = X ^г) • Лг)>
г ( геИ (
г (я) - X 4?, г),
(я )
X (у(? + 3) - У(г + 3))2
я, г) - ехр(-
3еА
г)
- ехР(- ,2),
(2)
2 Ьа^
где и(э) - новое значение после фильтрации в точке s, Z(s) - нормировочное значение для весовой функции, - весовая функция, определяемая через расстояние сР^^) окрестностей в ^ и I, определяет долю у(^), участвующую в усреднении текущего значения точки э, у(х) -значение сигнала в точке х, Щ(э) - диапазон в V, в котором будет происходить поиск подобных окрестностей I для окрестности э, определяется внутри диапазона [э-М; э+М], А - диапазон смещений, определяющий размах окрестностей [-Р; Р], Ьа - число элементов, входящих в А, равное 2Р+1, X - значение фильтра или значение полосы пропускания, в [14] равно среднеквадратичному отклонению шума на сигнале а и в (2) задается как 2а1, так как предполагается, что дисперсия равна d1~2o1. Визуальное представление получения одного весового значения представлено на рис. 3.
)
Рис. 3. Сравнение окрестности целевой точки э с окрестностью точки ^ из Щэ)
В [14-18] общая идея одномерной КЬМ-фильтрации заключалась в удалении шумовых составляющих из общей тенденции фонокардиограммы. Например, в работе [18] из фоно-кардиограммы исключалась составляющая фоновых дыхательных и других физиологических шумов. Визуально это соответствует тому, что были исключены низкоамплитудные высокочастотные колебания, а низкочастотные колебания, находящиеся между комплексами S1 и значительно уменьшающие выраженность S1, остаются на сигнале, что отвечает целям работы [18]. В работе интерес представляет использование фильтра нелокального усреднения для ослабления тех участков фонокардиограммы, которые не относятся к формации первого тона S1. По этим требованиям в области интереса пропускания фильтра будет попадать комплекс колебаний, который имеет выраженность по амплитуде и соответствует физиологическим и аппаратным характеристикам.
Как было уже рассмотрено выше, общей целью конструируемого фильтра является устранение фонового шума на фонокардиограмме с целью подготовки сигнала к сегментации тонов сердца. Поэтому можно обобщить, что целью фильтрации является подавление тех участков (формаций) сигнала, которые не соблюдают следующие условия:
- имеют выраженные по амплитуде колебания на протяжении 120-150 мс (длина комплекса S1);
- колебания формаций попадают в диапазон частот [штш, ^тах];
- между комплексами S1 RR-интервал по длине должен быть по продолжительности не менее чем допустимое значение минимальной частоты сердечных сокращений (Л/тт).
Таким образом, для достижения необходимого результата при анализе фильтруемой точки сигнала необходимо учитывать амплитуду колебания, в которую она входит. В (2) в качестве параметра полосы пропускания амплитуды сигнала выступало значение фильтра X. В [14] параметр фильтрации X определялся как среднеквадратичное отклонение шума относительно основной тенденции сигнала. В соответствии с условиями, определенными выше, в качестве полосы пропускания необходимо установить пороговый уровень амплитуды, который позволит выделить выраженные комплексы потенциального тона сердца. При этом нужно помнить о том, что амплитуда полезного сигнала может меняться с течением времени, то есть параметр фильтрации X должен динамически изменяться в соответствии с текущим уровнем амплитуды полезного сигнала. Поэтому он был определен как среднеквадратичное отклонение амплитуды сигнала относительно нулевого значения. Также если данный уровень должен изменяться динамически, то необходимо определить, внутри какого окна он должен быть рассчитан. Можно предположить, что если . - целевая точка фильтрации, то значение фильтра X должно быть рассчитано в диапазоне [л-////тш; л+////тт]. Таким образом, при нахождении точки . внутри комплекса S1 при расчете значения X будут участвовать колебания, полностью входящие в S1 и частично входящие в межтоновые участки справа и слева от S1. Если точка . будет находиться в середине межтонального участка, то в значение X будет внесена доля значений всего межтонального участка и части двух комплексов S1.
Из теории КЬМ-фильтрации известно [9], что для расчета степени подобия между двумя точками . и I необходимо учитывать не подобие конкретных значений в точках . и t, а подобие их окрестностей. Можно предположить, что для решаемой в работе задачи этой окрестностью могут являться колебания, в которые входит точка. Отдельные колебания характеризуются периодом Т, а также существует диапазон их частот, определяющийся свойствами полезного сигнала. Таким образом, достаточно установить Т в соответствии с минимальной частотой ^тш, которой будет соответствовать максимальный период. В выражении (2) значение Р будет устанавливаться как половина значения от Т. Для определения степени подобия между окрестностями . и t при расчете весового значения должны использоваться абсолютные значения (модуль) амплитуды в этих точках для того, чтобы производить сравнение отклонения от основной тенденции (нуля).
В выражении (2) коэффициент 2 при X2 внесем в масштабирующий коэффициент ¡и, так как это позволит регулировать степень сглаживания сигнала путем изменения доли сР(.^), необходимой для приближения к X2. Причем чем ближе отношение между и X2 к нулю,
тем большая доля текущего значения t будет вноситься в усреднение и будет ближе к ^(.,.)=1. Чем больше значение отношения между и X2, тем меньшая доля от t будет
вноситься в усреднение и тем большая доля значения в точке . будет использована для установки фильтруемого значения.
Также необходимо определить длину диапазона Щ(.) и правило, по которому будут выбираться элементы t из Щ(.). Для определения длины Щ(.) рассмотрим три ситуации: точка . находится в середине комплекса S1, точка . лежит в начале или в конце комплекса S1. Для сравнения точки . со всеми колебаниями, входящими в комплекс S1, в первом случае необходимо будет рассматривать диапазон [-0.5Ь, 0.5Ь], во втором случае - диапазон [0, Ь], в третьем - [-Ь, 0]. Таким образом, можно сделать вывод, что значение М в выражении (2) должно быть установлено равным Ь. Если точка . лежит внутри комплекса S1, ее можно будет сравнить со всем комплексом S1 вне зависимости от того, находится она в начале комплекса S1 или в конце. Так как для сравнения двух точек на сигнале важными являются амплитуды колебаний, в которые они входят, а не конкретные значения в точках, можно допу-
стить, что диапазон значений [я-М; я+М] можно разбить на узлы с шагом Р. Но стоит отметить, что при установке параметров М и Р нужно руководствоваться тем, что М должен быть кратен Р. Итоговый вид фильтра представлен в (3).
u(s) -
1
Z ]=-ш-1
m-1
2 w(s, tj ф) • v(t ; (s)),
M
tj (s) = s - M + jP, ш - —
Р
2 (ф + - - |v(t + -)|)2 2 w( s, t) - ехр(--^=—р---) - ехр(--^^),
(s)
цЛ2^)
2 v( s +1)2
Л2^) - 1 -
2 Штт +1
ш—1
z(s) - 2 w(s, tj (s)).
- ш+1
(3)
Параметры ЕЯтт, Р, М были установлены в соответствии с физиологическими параметрами сердцебиений плода и характеристиками датчика. Параметр ¡и методом подбора был установлен в значение 0.3. Результат работы фильтра представлен на рис. 4.
500 250 0
^ -250 -500
а
д
туд
и л п
27,7
27,72 27,4
Время, мин
27,76
Рис. 4. Результат работы полученного фильтра (3) на сигнале сердцебиений плода с полосой частот в [20 Гц, 40 Гц], красная кривая - исходный зашумленный сигнал, черная кривая - результат обработки
4. Обоснование и доработка фильтра
В конструируемом фильтре (3) в качестве параметров для определения расстояния (подобия) между двумя окрестностями МЛ2^^) и расчета значения фильтра X2 должны выступать амплитуды колебаний, которые находятся внутри окрестностей. Но весь механизм фильтра сравнивает и использует точечные значения амплитуд сигнала у(^), а не амплитуды колебаний. На каждом участке фонокардиограммы может присутствовать сумма гармоник с раз-
личными параметрами, поэтому при наложении двух участков друг на друга явно не прослеживается сравнение амплитуд колебаний.
Представим, что текущая гармоника описывается функцией:
где А - амплитуда гармоники. Тогда известно, что среднеквадратичное отклонение фонокар-диограммы равно:
Таким образом, среднеквадратичное отклонение гармоники примерно составляет 0.707 от амплитуды данной гармоники. Поэтому, учитывая разницу между значением А и X, погрешность, которую даст задание значения Т и присутствие нескольких гармоник на сигнале, можно предполагать зависимость ЩС2(э,^) неточной, но допустимой для задач фильтрации, так как гармоника с наибольшей амплитудой будет вносить основную долю в определение значения А , если параметр Т приближен к ее частоте.
Физически на реальном сигнале сердцебиений ожидается, что значение параметра X будет лежать на уровне или чуть выше, чем амплитуда колебаний между двумя биениями сердец, не содержащих выраженный комплекс S1. Параметр /л позволяет опускать или поднимать этот пороговый уровень. Динамическое определение параметра X позволит работать более деликатно на низкоамплитудных сигналах, когда разница между амплитудой межтоновых участков и амплитудной комплекса S1 будет менее выражена по сравнению с тем же соотношением на высокоамплитудных колебаниях (разное отношение сигнал-шум).
Другим вопросом является интерпретация параметра расстояния между двумя окрестностями ЩС2(э,^) как разницы между колебаниями внутри S1 и вне S1. Ранее было определено, что комплекс S1 имеет длину Ь, равную 120-150 мс, и что обычно комплекс S1 содержит несколько выраженных по амплитуде колебаний, а колебания на границах S1 имеют более низкую амплитуду по сравнению с «лидирующими» колебаниями (рис. 1, 2). Таким образом, можно предположить, что если в окрестность точки э входят все или часть «лидирующих» по амплитуде колебаний, а в окрестность точки I входят колебания, которые завершают или начинают комплекс S1, то существует разница в «усредненных» значениях амплитуд окрестностей э и Поэтому, если анализировать выражение (3), то можно понять, что фактически разница в «усредненных» значениях амплитуд будет заложена в параметр ЩС2(э,^), и она должна превышать значение /X в несколько раз. Для того чтобы добиться описанного выше эффекта, необходимо, чтобы Т был выбран таким образом, чтобы в него попадала как минимум половина одного «лидирующего» колебания. Представим, что длительность комплекса S1 составляет 150 мс, а Т равен 50 мс, следовательно, параметр Р равен 25 мс. Пусть точка э расположена в середине комплекса S1 (подобно рис. 5а), что соответствует точке R. Таким образом, в окрестность э [э-Р; э+Р] попали участки всех «лидирующих» колебаний комплекса S1, в две соседние окрестности [э-1Р; э] и [э; э+1Р] попала остальная часть «лидирующих» колебаний и часть комплекса S1, которая еще имеет выраженность по сравнению с межтоновыми участками. Общая длина окрестностей составит 4Р, что равняется 100 мс. Как видно на рис. 5а, на участке [э-Р; э+Р] амплитуды сигнала не меняются, на участке [э-1Р; э] сигнал практически не изменился, на участке [э; э+1Р] имеется «падение» амплитуды «последнего» колебания, так как далее следуют колебания подобного уровня.
у(?) - А $,т(а>г), со - —,
(4)
^ 1
Я(А, с) - — | (А 8т(сг))2ёг
V о
18т2(сгуг -V0.49А2 - 0.707А. (5)
о
0 2 4 6 8
Время, мс
а) б) в) г) д)
0 ж 4 2 ж 4 2 $ 4 4 s 6 6 s 8
Время, мс
Рис. 5. Изменение w(s,t) и н^э) на сигнале сердцебиений плода при изменении положения точки э. Верхний ряд: серая кривая - исходный участок сигнала, из которого будет браться точки э в нижнем ряду, черная кривая - результат его обработки. Нижний ряд: серые точки - значение пн(э^), черные точки - значение а) точка э находится в середине комплекса 81; б) точка э движется к концу
комплекса 81; в) точка э лежит в межтоновом участке, точки t лежат в области комплекса 81 и в межтоновом участке; г) точка э лежит в межтоновом участке, так же как и большинство точек V,
д) точка э находится в начале комплекса 81
Рассмотрим фактическое изменение н(э^) и при фильтрации фонокардиограммы
(рис. 5). При этом значимым для нас будет не только значение весовой функции н(э^), но и значение которое далее будет называться нормированным значением весовой
функции Нормированное значение весовой функции показывает, какая доля значения
v(t) фактически будет использоваться для усреднения. Так, на рис. 5а точка s находится в центре комплекса 81, значение пн(э,э) равно единице, значения стремятся к нулю, то
есть текущее значение v(s) вносится в фильтруемое значение и{() в полном объеме, а значения пн(£) не используются. На рис. 5б точка э находится в конце комплекса 81, пн(э,э) отдаляется от единицы, а пм(э^) вне комплекса 81 начинают вносить влияние в усреднение. На рис. 5в точка э лежит вне комплекса 81, примерно четыре значения пм(э^) лежат в комплексе 81 и стремятся к нулю, а остальные значения пн(э^) стремятся к уровню пн(э,э). На рис. 5г э находится в середине межтонового участка, а практически все значения пн(э^) и пн(э,э) находятся на одном уровне. На рисунке 4д э приближается к началу комплекса 81, а часть точек t лежат на следующем комплексе 81, пн(э^) в межтоновом участке находятся практически на уровне с пн(э,э), пн(э£) на комплексе 81 стремятся к нулю.
Из вышесказанного можно сделать вывод, что значимым показателем является пн(э,э), который позволяет оставлять без изменений значения «лидирующих» колебаний 81 и плавно понижать амплитуды колебаний комплекса первого тона ближе к межтоновым участкам. Если построить график пн(э,э) (рис. 6), то можно увидеть сигнал, который похож на огибающую обработанного сигнала с той разницей, что вне зависимости от текущего уровня амплитуды на сигнале в центрах комплексов 81 значение графика пн(э,э) будет стремиться к
единице. Это позволяет предположить, что этот график может быть полезен для сегментации тонов сердца на сигнале, на котором возможны резкие перепады амплитуды полезного сигнала.
Также стоит отметить, что если посмотреть на рис. 4, то можно увидеть, что на межтоновых участках фильтрованного сигнала отмечается падение амплитуды сигнала, но сам сигнал не соответствует форме исходного сигнала. Это искажение является следствием усреднения и погрешности, полученной в результате сравнения колебаний с разными параметрами фазы и частоты. Поэтому можно сделать предположение, что если не производить усреднение, а только изменять уровень амплитуды сигнала с помощью коэффициента nw(s,s), то это может позволить получить более низкое значение сигнал-шум и не допустить искажение сигнала. Тогда в конструируемый фильтр будут внесены следующие изменения:
= = 40. (6)
Z(s) Z (s) v '
Так как полученная модификация алгоритма фильтрации не соответствует понятию метода нелокального усреднения, так как из него исключено усреднение, а исходные параметры были подобраны в соответствии с физиологическими параметрами сердца, назовем его pNLF (phonocardiogram nonlocal filtration). Работа фильтра в виде (3) с изменениями (6) представлена на рис. 6.
0
27,7 27,72 27,4 27,76
Время, мин.
Рис. 6. Результаты работы фильтра рКЬБ и график нормированной весовой функции Верхний ряд: серая кривая - исходный зашумленный сигнал, черная кривая - результат обработки фильтром рКЬМ. Нижний ряд: значения нормированной весовой функции
5. Запуск фильтра на реальных данных
Оценка работы сконструированного фильтра (3) с изменениями (6) на реальных данных была осложнена тем, что интерес для работы представляла проверка работы фильтра на фо-нокардиограммах, полученных с различных приборов и с различным уровнем ЧСС. Но при этом существует ограниченное число источников данных фонокардиограмм, доступных для анализа. Также для оценки подавления межтоновых участков в результате работы фильтра необходимо было иметь разметку анализируемых сигналов, где были бы обозначены межтоновые участки, комплексы S1 и S2 и исключены участки, где тоны сердцебиений визуально не определимы экспертом. Первично были взяты акустические сигналы сердцебиений из соревнования «Classifying Heart Sounds Challenge» [19]. В этом источнике данных существует архив сигналов сердцебиений без патологий длиною пару секунд, а также csv-файл с разметкой комплексов S1 и S2. В соответствии с данной разметкой была произведена ручная сегментация комплексов S1 и S2 (данные и их описание представлены в [20]). В ходе сегментации были получены навыки разметки данного типа сигналов по следующим критериям (которые также основаны на физиологических критериях):
- комплексы S1 и S2 выражены по амплитуде, имеют повторяющиеся формации;
- расстояние между S1 и S2 одного биения обычно не превышает расстояния между текущим S2 и последующим S1, особенно на сигналах сердцебиений взрослых;
- S2 обычно по длительности меньше комплекса S1;
- соседние RR-интервалы не должны отличаться друг от друга на более чем примерно 3 уд/мин (гарантировано физиологически для нормального стабильного ритма), ЧСС имеет выраженную плавную тенденцию.
Предполагается, что конструированный фильтр должен оставлять без изменений потенциальные комплексы S1, визуально выраженные на сигнале, поэтому визуальная разметка выраженных формаций на сигнале по критериям, описанным выше, была принята достаточной для оценки отношения сигнал-шум на реальных сигналах. В данных [19] не указан датчик или датчики, с помощью которых были захвачены акустические сигналы, но определенный по спектру сигналов диапазон частот [wmin, ^max] охватывает полезный сигнал сердцебиений для всех сигналов данного репозитория. Также были выбраны только те сигналы, на которых присутствуют хотя бы 10 сердцебиений. Оценка такого небольшого числа полезных формаций является не полной, но позволяет сформировать примерное представление о работе фильтра. Все характеристики сигналов, которые будут описаны в работе, указаны в [20].
В качестве другого источника данных были взяты несколько сигналов из [21], которые используются в рамках исследования [22]. На сигналах представлены сердцебиения плода. Также для оценки были использованы деперсонализированные сигналы плода и беременной (в одном сигнале на разных полосах частот), которые были зарегистрированы авторским коллективом с помощью приложения и прибора FetalCare [20]. Все описанные выше сигналы были размечены и приведены к одинаковой частоте дискретизации.
Расчет характеристик сигналов и оценка эффективности фильтра были произведены с помощью следующих величин:
AmpSD -
N-1 A 2 у Ai
¿0 N '
SNRdB -20logio( dSNRdB - SNR
AmpSD
heart
LostAmp<>/0 -
pNLF AmpSD
AmpSDnoise
- SNRreal,
pNLF
AmpSD
REAL
100%,
где N - число колебаний, входящих в комплекс S1, S2 (hearttone) или межтонового шума (noise), At - амплитуда t-колебания, входящего в комплекс S1, S2 (hearttone) или межтонового шума (noise), AmpSD - среднеквадратичное отклонение амплитуд колебаний комплекса от нулевой оси, SNRdB - отношение сигнал-шум, dSNRdB - разница между отношением сигнал-шум до и после обработки фильтром, LostAmp% - потери амплитуды комплекса после проведения pNLF-фильтрации.
Основными параметрами, с помощью которых можно оценить выраженность полезного сигнала, являются отношения сигнал-шум SNRdB - необходимо сравнить отношение сигнал-шум до и после обработки (dSNR%). Но получение абсолютного числа кратного изменения затрудняется тем, что амплитуда полезного сигнала постоянно изменяется внутри сигнала. Поэтому для оценки работы фильтра на участках сигналов было рассчитано математическое ожидание и среднеквадратичное отклонение показателей, указанных в (7). Хотя это не является абсолютно верным подходом ввиду того, что это соединит не только погрешность физического изменения амплитуды полезного сигнала, но и погрешность в разметке и разных уровнях подавления участков сигнала сконструированным фильтром, но этот подход дает примерное представление об уровне «шумоподавления» и разбросе значений рассчитанных величин. Все рассчитанные параметры представлены в [20]. В табл. 1 указаны диапазоны математического ожидания показателей.
Таблица 1. Оценка эффективности сконструированного фильтра
База сигналов Разница SNR (dSNRdb) Потери (Ьо81Ашр%0) Потери шума (Ьо81Ашр8В%)
Classifying Heart Sounds Challenge [19] 12-14 дБ 0-8 % (новое значение: 92-100 % от амплитуды до фильтрации) 76-78 % (новое значение: 18-24 % от амплитуды до фильтрации)
Cesarelli и др. Physionet [21] 8-11 дБ 2-19 % (новое значение: 81-98 % от амплитуды до фильтрации) 66-72 % (новое значение: 28-34 % от амплитуды до фильтрации)
FetalCareDB (плод)[20] 9-11 дБ 1-8 % (новое значение: 92-99 % от амплитуды до фильтрации) 68-75 % (новое значение: 25-32 % от амплитуды до фильтрации)
FetalCareDB (мама) [20] 10-11 дБ 4-12 % (новое значение: 88-96 % от амплитуды до фильтрации) 73-75 % (новое значение: 25-27 % от амплитуды до фильтрации)
Стоит отметить, что параметры фильтрации были согласованы для применения на данных авторского коллектива, поэтому параметры оценки могут быть улучшены при «дона-стройке» параметров конструированного фильтра для сигналов [19] и [21]. В целом можно сделать выводы о том, что применение сконструированного фильтра позволяет улучшить
выраженность полезного сигнала, при этом изменение амплитуды полезного сигнала сердцебиений незначительно.
Также было определено, что реализация сконструированного фильтра, проведенная авторским коллективом, представляет хорошую скорость работы алгоритмов фильтрации: на ЭВМ с Intel Core Í5-4440 CPU 3.40GHz, 8Gb ОЗУ средняя скорость обработки 500 отчетов фонокардиограммы при частоте дискретизации в 500 Гц алгоритмом pNLF составила менее 1 мс [20].
6. Выводы и перспективы
Обработка акустических сигналов сердцебиений плода и человека является актуальным направлением, поэтому полученный в работе фильтр может применяться для первичной обработки с целью усиления полезного сигнала сердцебиений для его дальнейшей сегментации. В особенности это даст хорошие результаты на невыраженных сигналах сердцебиений, например, сигналах сердцебиений плода или сигналах сердцебиений взрослого, захваченных c артерии. Оценка результатов запуска полученного фильтра pNLF на реальных данных сердцебиений показала хороший уровень усиления полезного сигнала (тонов сердца). Также скорость работы алгоритма pNLF показала возможность его использования в потоковой обработке и внедрения в устройства мониторинга. Полученный алгоритм на данный момент внедряется в приложение FetalCare, а далее будет проводиться дальнейший анализ его применения на реальных данных сердцебиений плода с целью выявления критериев к улучшению и оптимизации работы алгоритма.
В качестве перспективы работы над фильтром pNLF можно предложить внесение адаптивных компонент. Например, изменение параметра RRmm на текущий уровень ЧСС. Также в зависимости от текущего уровня амплитуд полезного сигнала можно изменять динамически показатель Можно также составить механизмы более точного определения отношения между параметрами [«min, «max], T, P, m и M для применения алгоритма для фонокардио-грамм разных датчиков. Особую перспективу применения представляет график w(s,s)IZ(s) для сегментации сердцебиений на сигнале и выравнивание потери амплитуды комплекса S1. Также полученный комплекс алгоритмов может быть использован на акустических данных других ритмичных периодичных сигналов.
Литература
1. Sbrollini A., Strazza A., Caragiuli M., Mozzoni C., Tomassini S., Agostinelli A., Burattini L. Fetal Phonocardiogram Denoising by Wavelet Transformation: Robustness to Noise II Computing in Cardiology Conference. 2017. V. 44. DOI: 10.22489/cinc.2017.331-075.
2. Koutsiana E., Hadjileontiadis L. J., Chouvarda I., Khandoker A. H. Fetal Heart Sounds Detection Using Wavelet Transform and Fractal Dimension II Frontiers in Bioengineering and Biotechnology. 2017. V. 5, № 49. DOI: 10.3389/fbioe.2017.00049.
3. Chourasia V. S., Mittra A. K. Wavelet-based denoising of fetal phonocardiographic signals II International Journal of Medical Engineering and Informatics. 2010. V. 2, № 2. P. 139-150. DOI:10.15 0 4/ijmei.2 010.031516.
4. Samieinasab M., Sameni R Fetal phonocardiogram extraction using single channel blind source separation II 2015 23rd Iranian Conference on Electrical Engineering. 2015. P. 78-83. DOI: 10.110 9/iraniancee.2 015.714 618 6.
5. Taralunga D. D., Ungureanu M., Hurezeanu B., Strungaru R. Fetal Heart Rate Estimation from Phonocardiograms Using an EMD Based Method II Recent Advances in Computer Sci-
ence Proceedings of the 19th International Conference on Computers (part of CSCC'15). 2015. P. 414-417.
6. Warbhe A. D., Dharaskar R. V., Kalambhe B. Extraction of Fetal Heart Sound Signals using Independent Component Analysis // Computational Intelligence Applications. 2010. P. 66-68.
7. Chetlur A., Sankar R., Moreno W.A., Hart, S. Trends in fetal monitoring through phonocardiography: Challenges and future directions // Biomedical Signal Processing and Control. 2017. V. 33. P. 289-305. DOI: 10.1016/j.bspc.2016.11.007.
8. Ismail S., Stddtqt I., Akram U. Localization and classification of heart beats in phonocardiography signals - a comprehensive review // EURASIP Journal on Advances in Signal Processing. 2018. № 26. DOI: 10.1186/s13 63 4-018-05 45-9.
9. Buades A., Coll B., Morel J. M. A review of image denoising algorithms, with a new one // Multiscale Modeling and Simulation: A SIAM Interdisciplinary Journal, Society for Industrial and Applied Mathematics. 2005. V. 4, № 2, P. 490-530.
10. Аед А.М., Исаков Р. В., Сушкова Л. Т., Аль-Хайдри В. А. Алгоритм построения кардиоин-тервалограммы на основе фонокардиограммы // Радиотехнические и телекоммуникационные системы. 2016. № 2. С. 34-43.
11. Большее А. С., Сидоров Д. Г., Овчинников С. А. Частота сердечных сокращений. Физио-лого-педагогические аспекты. Н. Новгород: ННГАСУ, 2017. 76 с.
12. Altunkaya S., Kara S., Gormus N., Herdem S. Time Domain Features of Heart Sounds for Determining Mechanical Valve Thrombosis // Lecture Notes in Electrical Engineering. 2011. P. 173-181. DOI: 10.1007/978-94-007-1192-1_15.
13. Zhdanov D. S., Bureev A. S., Kostelet Y. V., Khokhlova L. A., Dtkman E. Y. A Mobile Device for Assessing Fetal Status Based on Monitoring Cardiovascular System Parameters. Biomedical Engineering. 2018. V. 52, № 2. P. 87-91. DOI: 10.1007/s10527-018-9789-9.
14. Tracey B. H., Miller E. L. Nonlocal Means Denoising of ECG Signals // IEEE Transactions on Biomedical Engineering. 2012. V. 59, № 9. P. 2383-2386. DOI: 10.1109/tbme.2012.2208 964.
15. Lee Y., HwangD. Periodicity-based nonlocal-means denoising method for electrocardiography in low SNR non-white noisy conditions // Biomedical Signal Processing and Control. 2018. V. 39. P. 284-293. DOI: 10.1016/j.bspc.2017.08.00 6.
16. Ttan X., Li Y., Zhou H., Li X., Chen L., ZhangX. Electrocardiogram Signal Denoising Using Extreme-Point Symmetric Mode Decomposition and Nonlocal Means // Sensors. 2016. V. 16, № 10. DOI: 10.3390/s16101584.
17. Singh P., Shahnawazuddin S., Pradhan G. An Efficient ECG Denoising Technique Based on Non-local Means Estimation and Modified Empirical Mode Decomposition // Circuits, Systems, and Signal Proc. 2018. V. 37. P. 4527-4547. DOI: 10.1007/s00034-018-0777-9.
18. Рудницкий А. Г. Использование метода нелокального усреднения для разделения звуков сердца и звуков дыхания // Акустический журн. 2014. Т. 60, № 6. С. 688-695.
19. Bentley P., Nordehn Gl., CoimbraM., Mannor S., Getz R Datasets for Classifying Heart Sounds Challenge [Электронный ресурс] // Classifying Heart Sounds Challenge: [Электронный ресурс]. [2012]. URL: http://www.peterjbentley.com/heartchallenge /#dates (дата обращения 07.06.2020 г.).
20. Kosteley Y., Bureev A.; Zhdanov D., Yuriev S. Fetal and human PCG Data with the first and second heart tone segmentation for pNLF algorithm testing [Электронный ресурс] // Mendeley Data: [сайт]. [2020]. URL: doi.org/10.17 632/wk6rkmdgrf.1 (дата обращения 07.06.2020 г.).
21. Cesarelli M, Ruffo M, Romano M, Bifulco P. Fetal PCG Database [Электронный ресурс] // PhysioNet: [сайт]. [2020]. URL: doi.org/10.13026/C2PD5B (дата обращения 07.06.2020 г.).
22. Cesarelli M, Ruffo M, Romano M, Bifulco P. Simulation of foetal phonocardiographic recordings for testing of FHR extraction algorithms // Computer Methods Programs Biomed. 2012. V. 107, № 3. P. 513-23. DOI: 10.1016/j.cmpb.2011.11.008.
Статья поступила в редакцию 25.08.2021; переработанный вариант — 23.09.2021.
Костелей Яна Валерьевна
инженер-программист ООО «Диагностика +», старший преподаватель кафедры экономической математики, информатики и статистики Томского государственного университета систем управления и радиоэлектроники, e-mail: kosteleyyv@gmail. com.
Жданов Дмитрий Сергеевич
к.т.н., старший научный сотрудник ООО «Диагностика +» (634029, Томск, пр. Фрунзе, 11б, пом. 407), старший научный сотрудник Томского государственного университета (634050, Томск. пр. Ленина, 36), e-mail: d_s_zhdanov@mail.ru.
Боровской Игорь Георгиевич
д.ф.-м.н., профессор, заведующий кафедры экономической математики, информатики и статистики Томского государственного университета систем управления и радиоэлектроники (634050, Томск, пр. Ленина, 40), e-mail: igor.g.borovskoi@tusur.ru.
Non-local averaging filter adaptation for heart sounds amplification on fetus and humans' phonocardiograms
Y. V. Kosteley, D. S. Zhdanov, I. G. Borovskoy
The design and substantiation of the pNLF filter based on a non-local averaging filter modification is presented. The filter is used for solving the problem of heart sounds amplification on fetal and human phonocardiograms for further tone segmentation and heart rate calculation relevant for fetal monitoring. The parameters for choosing the filter and the prospects for its optimization for specific tasks are presented. The data and results of filtration performed by the pNLF on fetus and human phonocardiograms recorded by using different acoustic sensors are provided.
Keywords: phonocardiogram, non-local means, heart rate segmentation.
