СПОСОБ АДАПТИВНОГО ПОДАВЛЕНИЯ ПОМЕХ В ЭЛЕКТРОКАРДИОСИГНАЛЕ
METHOD OF ADAPTIVE SUPPRESSION OF NOISE IN ELECTROCARDIOSIGNAL

Изобретение относится к медицине, в частности к кардиологии, и может быть использовано для адаптивного подавления помех в электрокардиосигнале (ЭКС).
Эффективность автоматического анализа ЭКС определяется точностью измерения его амплитудных и временных параметров, достоверностью обнаружения и распознавания его отдельных элементов. В свою очередь, основной причиной погрешностей измерений и ошибок обнаружения (распознавания) элементов ЭКС являются помехи, различные по своему происхождению, интенсивности, спектральным и статистическим характеристикам, взаимодействию с полезным сигналом.
В алгоритмах и технических средствах автоматического анализа ЭКС устанавливается допустимый уровень помех - заранее заданный количественный показатель, характеризующий уровень остаточных помех в электрокардиосигнале, при котором последующие процедуры автоматического анализа могут быть выполнены с заявленным качеством.
Традиционно для подавления помех на ЭКС в основном применяются способы, основанные на частотной линейной фильтрации [1], что связано с наличием подходящего математического аппарата, простотой интерпретации и расчета линейных фильтров. В то же время использование методов линейной фильтрации не позволяет получить допустимый уровень остаточных помех без искажений ЭКС. В этих случаях оптимальными способами подавления помех, позволяющими получить сигнал с допустимым уровнем помех, являются нелинейные фильтры [2]. Среди группы нелинейных фильтров получили распространение в настоящее время фильтры, основанные на частотно-временном представлении сигнала. Особенностью этих способов является разложение сигнала на элементарные частотные компоненты, обработка (фильтрация) отдельных компонент и последующее восстановление сигнала.
Предполагается, что такой подход обеспечивает максимально возможное удаление помех при минимальном искажении ЭКС.
Известен способ подавления помех на основе непрерывного вейвлет-преобразования [3], заключающийся в регистрации сигнала, его непрерывном вейвлет-преобразовании, пороговой обработке (thresholding) вейвлет-коэффициентов и восстановлении сигнала без помех.
К недостатку известного способа подавления помех на основе непрерывного вейвлет-преобразования относится недостаточно высокая эффективность помехоподавления, обусловленная невозможностью адаптации базового вейвлета к структуре и локальным особенностям зарегистрированного ЭКС.
Действительно, в известном способе подавления помех на основе непрерывного вейвлет-преобразования применяется фиксированный базовый вейвлет (Gaussian) для разложения ЭКС и поэтому данный подход не позволяет учесть все локальные особенности конкретного исследуемого ЭКС и сопровождающих его помех.
Наиболее близким к предлагаемому изобретению является способ подавления помех и удаления артефактов в электрокардиосигнале [4], заключающийся в том, что осуществляют регистрацию электрокардиосигнала, разложение его на частотные составляющие, восстановление и вывод электрокардиосигнала.
Как следует из формулы изобретения известного способа подавления помех и удаления артефактов в электрокардиосигнале сначала осуществляется разложение электрокардиосигнала на частотные составляющие, затем подавление помех и удаление артефактов посредством фильтрации частотных составляющих и восстановление электрокардиосигнала посредством суммирования. Как следует из описания, известный способ подавления помех и удаления артефактов в электрокардиосигнале основан на разложении ЭКС на частотные составляющие с помощью декомпозиции на эмпирические моды (ДЭМ) [5]. Эмпирические моды - это монокомпонентные составляющие сигнала, модулированные по амплитуде и частоте, т.е. их амплитуда и частота меняются во времени. Моды не имеют строгого аналитического описания, но должны удовлетворять двум условиям [6, 7]:
- общее число экстремумов и число пересечений нуля должны отличаться не более чем на единицу;
- среднее значение двух огибающих: верхней, интерполирующей локальные максимумы, и нижней, интерполирующей локальные минимумы, должно быть приближенно равно нулю.
Декомпозиции на эмпирические моды является адаптивным методом анализа нестационарных сигналов. Базис, используемый для разложения ЭКС (набор эмпирических мод), конструируется непосредственно из самого регистрируемого ЭКС. Это позволяет учесть все его локальные особенности, форму информативных участков, присутствие помех различного вида и интенсивности. Кроме адаптивности, ДЭМ обладает и другими важными для практических приложений свойствами: ортогональностью и полнотой [5, 6].
На фигуре 1 приведена схема известного способа подавления помех и удаления артефактов в электрокардиосигнале.
На фигуре 2 приведен идеальный ЭКС с обозначением элементов (зубцов и интервалов).
На фигуре 3 (а, б) приведены кардиоцикл зарегистрированного (а) и зашумленного (б) ЭКС.
На фигуре 4 приведена схема разложения ЭКС на эмпирические моды.
На фигуре 5 (а, б, в, г, д, е, ж, з) приведены результаты разложения зашумленного ЭКС на эмпирические моды (а, б, в, г, д, е, ж) и остаток - глобальный тренд сигнала (з).
На фигуре 6 (а, б) приведены результаты подавления помех в ЭКС известным способом подавления помех и удаления артефактов в электрокардиосигнале. На фигуре 6,а представлен восстановленный ЭКС при удалении первой высокочастотной моды, на фигуре 6,б - восстановленный ЭКС при удалении первых двух высокочастотных мод.
Согласно фигуре 1 в известном способе подавления помех и удаления артефактов в ЭКС выполняются следующие этапы.
1. Регистрация ЭКС. На этом этапе осуществляется регистрация ЭКС пациента. При регистрации ЭКС решаются вопросы:
- электрического взаимодействия между электродами и поверхностью кожи в месте расположения электродов;
- усиления ЭКС;
- аналого-цифрового преобразования ЭКС;
- электрической безопасности пациента;
- сопряжения устройства регистрации ЭКС с компьютером.
После этапа регистрации на вход компьютера поступает ЭКС пациента (см. фигуру 3,а).
2. Разложение сигнала на частотные составляющие. На этом этапе осуществляется разложение ЭКС на частотные составляющие с помощью декомпозиции на эмпирические моды (см. фигуру 5).
Схема разложения ЭКС на частотные составляющие приведена на фигуре 4 и включает в себя следующие этапы.
j
i
1. Определение локальных экстремумов (максимумов и минимумов) ЭКС f(t):
j
i
j
i-1
j
i
j
i+1
- значение i-го отсчета f(t) является локальным максимумом, если выполняется условие f(t)< f(t)≥ f(t);
j
i
j
i-1
j
i
j
i+1
значение i-го отсчета f(t)) является локальным минимумом, если выполняется условие f(t)>f(t)≤f(t).
j
i
j
i
j
i
2. Определение верхней e(t) и нижней g(t) огибающих ЭКС с помощью кубической сплайн-интерполяции [8] по найденным локальным экстремумам f(t):
в
в
в
в
где α, b, c, d - коэффициенты для каждого значения i-го отсчета верхней огибающей ЭКС;
н
н
н
н
α, b, c, d - коэффициенты для каждого значения i-го отсчета нижних огибающих ЭКС.
3. Вычисление среднего значения огибающих ЭКС в соответствии с выражением:
j
i
где h(t) - среднее значение огибающих ЭКС;
j
i
j
i
e(t) и g(t) - верхняя и нижняя огибающая исходных ЭКС соответственно.
4. Вычисление остатка ЭКС по формуле
<mrow><msub><mi>s</mi><mi>j</mi></msub><mo stretchy="false">(</mo><msub><mi>t</mi><mi>i</mi></msub><mo stretchy="false">)</mo><mo>=</mo><msub><mi>f</mi><mi>j</mi></msub><mo stretchy="false">(</mo><msub><mi>t</mi><mi>i</mi></msub><mo stretchy="false">)</mo><mo>−</mo><msub><mi>h</mi><mi>j</mi></msub><mo stretchy="false">(</mo><msub><mi>t</mi><mi>i</mi></msub><mo stretchy="false">)</mo><mo>,</mo><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mo stretchy="false">(</mo><mn>4</mn><mo stretchy="false">)</mo></mrow>
j
i
где s(t) - остаток частотных составляющих ЭКС.
5. Вычисление значения критерия останова. В качестве критерия останова декомпозиции используется значение нормализованной квадратичной разности, определяемое как [7]:
<mrow><mi>S</mi><mi>D</mi><mo>=</mo><mstyle displaystyle="true"><munderover><mo>∑</mo><mn>1</mn><mi>n</mi></munderover><mrow><mrow><mo>[</mo><mrow><mfrac><mrow><msup><mrow><mo stretchy="false">(</mo><mo>|</mo><msub><mi>f</mi><mi>j</mi></msub><mo stretchy="false">(</mo><msub><mi>t</mi><mi>i</mi></msub><mo stretchy="false">)</mo><mo>−</mo><msub><mi>S</mi><mi>j</mi></msub><mo stretchy="false">(</mo><msub><mi>t</mi><mi>i</mi></msub><mo stretchy="false">)</mo><mo>|</mo><mo stretchy="false">)</mo></mrow><mn>2</mn></msup></mrow><mrow><msub><mi>f</mi><mi>j</mi></msub><msup><mrow><mo stretchy="false">(</mo><msub><mi>t</mi><mi>i</mi></msub><mo stretchy="false">)</mo></mrow><mn>2</mn></msup></mrow></mfrac></mrow><mo>]</mo></mrow></mrow></mstyle><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mo stretchy="false">(</mo><mn>5</mn><mo stretchy="false">)</mo></mrow>
6. Проверка условия останова. На этом этапе осуществляется сравнение значения остатка ЭКС со значением нормализованной квадратичной разности. Если:
j
i
- SD>s(t), то переходят к выполнению действия 1;
j
i
j
i
j
i
SD<s(t) и h(t)>s(t), то переходят к выполнению следующего действия.
k
i
j
i
7. Вывод частотных составляющих ЭКС. На этом этапе осуществляется вывод эмпирических мод m(t) и остатка S(t) ЭКС. На фигуре 5 (а, б, в, г, д, е, ж) приведены моды зашумленного ЭКС. Частотная составляющая, показанная на фигуре 5,з, является глобальным трендом ЭКС и дальнейшему разложению не подлежит. Глобальный тренд представляет собой функцию монотонную (постоянную) функцию, которая не может быть разложена далее на моды [6].
3. Фильтрация частотных составляющих ЭКС. На этом этапе осуществляется фильтрация частотных составляющих ЭКС путем удаления высокочастотных мод ЭКС.
По мнению авторов известного способа подавления помех и удаления артефактов в электрокардиосигнале удалением высокочастотных мод, полученных с помощью декомпозиции ЭКС, обеспечивается подавление помех в ЭКС.
4. Восстановление и вывод ЭКС. Результатом этого этапа является восстановление ЭКС, который представляет собой сумму неудаленных мод (оставшихся) ЭКС и вывод ЭКС (см. фигуру 6).
Для оценки возможностей известного способа подавления помех и удаления артефактов в электрокардиосигнале необходимо оценить его помехоподавление. Для этого авторами предлагаемого изобретения осуществляется
выделение фрагмента зарегистрированного ЭКС с минимальным уровнем помех. Этот фрагмент считается эталонным ЭКС (см. фигуру 3,а);
- зашумление эталонного ЭКС (см. фигуру 3,б);
- подавление помех в зашумленном ЭКС.  При этом в соответствии с алгоритмом известного способа подавления помех и удаления артефактов в электрокардиосигнале получены два сигнала (см. фигуру 6): на фигуре 6,а показан восстановленный ЭКС при удалении первой высокочастотной моды, на фигуре 6,б представлен восстановленный ЭКС при удалении первых двух высокочастотных мод.
i
i
Сравнение эталонного ЭКС (см. фигуру 3,а) с восстановленными ЭКС (см. фигуры 6,а и 6,б) показывает их различие. Действительно, из анализа восстановленных ЭКС следует, что применение известного способа подавления помех и удаления артефактов в электрокардиосигнале искажает форму эталонного ЭКС. Для количественной оценки качества помехоподавления используется среднеквадратическая ошибка (СКО) отклонения восстановленного ЭКС y, от эталонного ЭКС x, [9]:
<mrow><mi>С</mi><mi>К</mi><mi>О</mi><mo>=</mo><mrow><mo>(</mo><mrow><msqrt><mrow><mfrac><mrow><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>i</mi><mo>=</mo><mn>1</mn></mrow><mi>n</mi></munderover><mrow><msup><mrow><mo stretchy="false">(</mo><msub><mi>x</mi><mi>i</mi></msub><mo>−</mo><msub><mi>y</mi><mi>i</mi></msub><mo stretchy="false">)</mo></mrow><mn>2</mn></msup></mrow></mstyle></mrow><mrow><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>i</mi><mo>=</mo><mn>1</mn></mrow><mi>n</mi></munderover><mrow><msubsup><mi>x</mi><mi>i</mi><mn>2</mn></msubsup></mrow></mstyle></mrow></mfrac></mrow></msqrt></mrow><mo>)</mo></mrow><mo>×</mo><mn>1</mn><mn>0</mn><mn>0</mn><mi>%</mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mo stretchy="false">(</mo><mn>6</mn><mo stretchy="false">)</mo></mrow>
i
где х - значения отсчетов исходного (эталонного) ЭКС (см. фигуру 3,а);
j
y - значения отсчетов восстановленного ЭКС (см. фигуру 6);
n - объем выборки.
Рассчитанная по выражению (6) среднеквадратическая ошибка для сигнала, приведенного на фигуре 6,а, составляет 37,8%; а для сигнала, приведенного на фигуре 6,б, - 28,1%. Полученные значения СКО указывают на недостаточное помехоподавление известного способа подавления помех и удаления артефактов в электрокардиосигнале. При этом изображения восстановленных ЭКС (см. фигуру 6) свидетельствуют об искажении и невысокой точности восстановления ЭКС.
Таким образом, проведенный анализ известного способа подавления помех и удаления артефактов в электрокардиосигнале доказал невозможность эффективного подавления помех в электрокардиосигнале посредством удаления высокочастотных мод частотных составляющих ЭКС. Выявлено, что при удалении высокочастотных мод в ЭКС и последующей реконструкции ЭКС происходит его искажение (см. фигуру 6).
Следовательно, известный способ подавления помех и удаления артефактов в электрокардиосигнале не обеспечивает качественного подавления помех в ЭКС, т.е. подавления помех с минимальным искажением полезного сигнала.
Недостатком известного способа подавления помех и удаления артефактов в электрокардиосигнале является низкая эффективность подавления помех, заключающаяся в высоком уровне остаточных помех и искажении ЭКС.
Изобретение направлено на повышение точности восстановления электрокардиосигнала и уменьшение уровня остаточных помех в электрокардиосигнале.
Это достигается тем, что в способе подавления помех и удаления артефактов в электрокардиосигнале, заключающемся в регистрации электрокардиосигнала, разложении его на частотные составляющие, восстановлении и выводе электрокардиосигнала, дополнительно осуществляют выделение кардиоцикла электрокардиосигнала, выделение участка изолинии на сегменте TP кардиоцикла электрокардиосигнала, построение поверхностей энергетической плотности выделенного кардиоцикла электрокардиосигнала и участка изолинии на основе преобразования Гильберта-Хуанга в системе координат энергия-частота-время, формирование пороговой поверхности из поверхности энергетической плотности участка изолинии, сравнение поверхности энергетической плотности выделенного кардиоцикла электрокардиосигнала со сформированной пороговой поверхностью.
Сущность предлагаемого способа адаптивного подавления помех в электрокардиосигнале заключается в том, что подавление помех в электрокардиосигнале осуществляется за счет выделения кардиоцикла (КЦ) электрокардиосигнала и участка изолинии кардиоцикла электрокардиосигнала на сегменте TP, построения поверхностей энергетической плотности (ПЭП) кардиоцикла электрокардиосигнала и выделенного участка изолинии, формирования и сравнения ПЭП выделенного кардиоцикла с пороговой поверхностью.
Главной отличительной особенностью предлагаемого способа адаптивного подавления помех в электрокардиосигнале от известного способа является выделение кардиоцикла ЭКС и участка изолинии на сегменте TP, сравнение ПЭП выделенного кардиоцикла со сформированной поровой поверхностью.
В предлагаемом изобретении выделение КЦ и участка изолинии на сегменте TP осуществляется исходя из следующих предположений.
Во-первых, сегмент TP представляет собой участок изолинии, в течение периода времени которого электрическая активность сердца отсутствует (за исключением зубца U) [10]. Авторами предлагаемого изобретения выделяется на сегменте TP участок, не содержащий зубец U. Авторы считают, что так как сердце на выделенном участке изолинии сегмента TP электрически не активно, то все регистрируемые сигналы а этом участке являются помехами. Выделение участка изолинии позволит адаптивно сформировать помеху кардиоцикла электрокардиосигнала.
Во-вторых, авторы считают маловероятным кардинальное изменение помехи в течение одного кардиоцикла. Поэтому выделенный участок изолинии является реальной помеховой составляющей кардиоцикла электрокардиосигнала.
Далее осуществляется сравнение пороговой поверхности с ПЭП кардиоцикла ЭКС; выделение области помехи под пороговой поверхностью; удаление области помехи. Эти операции, по мнению авторов предлагаемого изобретения, позволяют уменьшить уровень остаточных помех и искажений в ЭКС.
На фигуре 7 приведена схема предлагаемого способа адаптивного подавления помех в электрокардиосигнале.
На фигуре 8 (а, б) приведены выделенный участок изолинии кардиоцикла ЭКС на сегменте TP (а) и результаты разложения участок изолинии кардиоцикла ЭКС на эмпирические моды (б, в, г) и остаток - глобальный тренд сигнала (д).
На фигуре 9 приведена схема алгоритма построения поверхности энергетической плотности ЭКС.
На фигуре 10 (а, б) приведены поверхности энергетической плотности выделенного кардиоцикла ЭКС (а) и участка изолинии (б).
На фигуре 11 приведены сечение пороговой поверхности энергетической плотности выделенного участка изолинии в системе координат энергия-частота (а) и пороговая поверхность в системе координат энергия-время-частота (б).
На фигуре 12 приведен алгоритм сравнения поверхности энергетической плотности выделенного кардиоцикла электрокардиосигнала с пороговой поверхностью.
На фигуре 13 приведен результат подавления помех в ЭКС предлагаемым способом адаптивного подавления помех в электрокардиосигнале.
На фигуре 14 (а, б) приведены поверхности энергетической плотности в системе координат энергия-время-частота с обозначением элементов ЭКС (зубцов Р, Т, кардиокомплекса QRS и помех).
Анализ схем алгоритмов, приведенных на фигурах 1, 7, 9, 12, показывает, что предлагаемый способ адаптивного подавления помех в электрокардиосигнале основывается совершенно на другом подходе к подавлению помех и позволяет устранить указанный недостаток известного способа [4].
Авторы предлагают вместо удаления отдельных частотных составляющих ЭКС осуществлять выделение КЦ ЭКС и участка изолинии на сегменте TP, построение поверхностей энергетической плотности (ПЭП) кардиоцикла ЭКС и участка изолинии, формирование пороговой поверхности и сравнение ее с ПЭП кардиоцикла ЭКС.
Кроме того, визуальное представление ПЭП кардиоцикла ЭКС в системе координат энергия-время-частота, по мнению авторов, открывает возможность выявления новых диагностических признаков в электрокардиосигнале.
Рассмотрим особенности реализации предлагаемого способа адаптивного подавления помех в электрокардиосигнале (см. фигуру 7). Первый этап («Регистрация ЭКС») аналогичен этапу известного способа [4] (см. стр.4).
Для выполнения этапа «Выделение кардиоцикла ЭКС» используется алгоритм, приведенный в [11]. Согласно способу выделения начала кардиоцикла начало кардиоцикла ЭКС находится на расстоянии 1/4 интервала R-R от вершины зубца R. Известный алгоритм обеспечивает надежное выделение кардиоцикла ЭКС. На фигуре 3,а приведен выделенный кардиоцикл ЭКС известным способом [11].
«Выделение участка изолинии на сегменте ТР». На этом этапе осуществляется выделение участка изолинии на сегменте TP кардиоцикла ЭКС (см. фигуру 8а).
Сегмент TP представляет собой отрезок времени от конца зубца Т до начала зубца Р следующего кардиоцикла. В этот отрезок времени электрическая активность сердца отсутствует [12]. Следовательно, на участке изолинии [TP/2, Р] от середины сегмента TP до начала зубца Р следующего кардиоцикла (см. фигуру 26) отсутствуют зубец U и любые другие возможные информативные составляющие сигнала, а присутствуют только помехи. Поэтому участок изолинии [TP/2, Р] может быть использован в качестве источника помехи для одного кардиоцикла.
Таким образом, на этом этапе осуществляется адаптивное формирование помехи кардиоцикла электрокардиосигнала за счет выделения реальной помеховой составляющей кардиоцикла ЭКС.
Выделенный участок изолинии [TP/2, Р] на сегменте TP используется в предлагаемом изобретении для последующего адаптивного формирования пороговой поверхности. Под адаптивностью в данном случае понимается формирование собственной пороговой поверхности для каждого выделенного кардиоцикла ЭКС.
Следующие два этапа («Разложение кардиоцикла ЭКС на частотные составляющие» и «Разложение участка изолинии ЭКС на частотные составляющие») аналогичны этапу известного способа [4] «Разложение сигнала на частотные составляющие» (см. стр.4-6). На этих этапах осуществляется разложение кардиоцикла ЭКС и участка изолинии на частотные составляющие с помощью [6, 7] декомпозиции на эмпирические моды.
Этапы «Построение ПЭП кардиоцикла ЭКС» и «Построение ПЭП участка изолинии ЭКС» осуществляют на основе преобразования Гильберта-Хуанга в системе координат энергия-частота-время. Поверхность энергетической плотности (ПЭП) частотных составляющих электрокардиосигнала представляет собой распределение мгновенной энергии частотных составляющих ЭКС в каждой точке частотно-временной плоскости.
Преобразование Гильберта-Хуанга включает в себя последовательное выполнение двух действий: декомпозицию сигнала на эмпирические моды (ДЭМ) и преобразование Гильберта для каждой моды. Значения мгновенных частот представляют собой спектр Гильберта, который позволяет [13]
- определить по его виду наличие и характер амплитудной и частотной модуляций;
- идентифицировать временные и частотные диапазоны, где концентрируется энергия исходного сигнала;
- построить трехмерную поверхность энергетической плотности в системе координат энергия-частота-время.
После того как получены частотные составляющие кардиоцикла ЭКС (см. фигуру 5) и участка изолинии (см. фигуру 8) с помощью ДЭМ, выполняется преобразование Гильберта и построение трехмерной поверхности.
Построение трехмерной поверхности ЭКС в системе координат энергия-частота-время осуществляется посредством преобразования Гильберта частотных составляющих (эмпирических мод) сигнала [7, 13]. Схема построения ПЭП приведена на фигуре 9 и включает в себя следующие этапы:
1. Получение сопряженного по Гильберту ЭКС:
<mrow><msub><mi>M</mi><mi>k</mi></msub><mo stretchy="false">(</mo><msub><mi>t</mi><mi>i</mi></msub><mo stretchy="false">)</mo><mo>=</mo><mfrac><mn>1</mn><mi>π</mi></mfrac><mstyle displaystyle="true"><mrow><mo>∫</mo><mrow><mfrac><mrow><msub><mi>m</mi><mi>k</mi></msub><mo stretchy="false">(</mo><mi>τ</mi><mtext></mtext><mo stretchy="false">)</mo></mrow><mrow><msub><mi>t</mi><mi>i</mi></msub><mo>−</mo><mi>τ</mi></mrow></mfrac></mrow></mrow></mstyle><mi>d</mi><mi>τ</mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mo stretchy="false">(</mo><mn>7</mn><mo stretchy="false">)</mo></mrow>
k
i
где m(t) - эмпирическая мода, подвергнутая преобразованию Гильберта;
k
i
k
i
M(t) - сопряженный по Гильберту ЭКС, соответствующий моде m(t);
k - номер моды;
τ - независимая переменная.
2. Определение аналитического (комплексного) значения для каждой моды ЭКС:
<mrow><mi>Z</mi><mo>=</mo><msub><mi>m</mi><mi>k</mi></msub><mo stretchy="false">(</mo><msub><mi>t</mi><mi>i</mi></msub><mo stretchy="false">)</mo><mo>+</mo><mi>q</mi><msub><mi>M</mi><mi>k</mi></msub><mo stretchy="false">(</mo><msub><mi>t</mi><mi>i</mi></msub><mo stretchy="false">)</mo><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mo stretchy="false">(</mo><mn>8</mn><mo stretchy="false">)</mo></mrow>
<mrow><mi>q</mi><mo>=</mo><msqrt><mrow><mo>−</mo><mn>1</mn></mrow></msqrt></mrow>
где
- мнимая единица.
3. Определение значения мгновенной амплитуды для каждой моды ЭКС:
<mrow><msub><mi>α</mi><mi>k</mi></msub><mo stretchy="false">(</mo><mi>t</mi><mo stretchy="false">)</mo><mo>=</mo><msqrt><mrow><msub><mi>m</mi><mi>k</mi></msub><msup><mrow><mo stretchy="false">(</mo><msub><mi>t</mi><mi>i</mi></msub><mo stretchy="false">)</mo></mrow><mn>2</mn></msup><mo>+</mo><msub><mi>M</mi><mi>k</mi></msub><msup><mrow><mo stretchy="false">(</mo><msub><mi>t</mi><mi>i</mi></msub><mo stretchy="false">)</mo></mrow><mn>2</mn></msup></mrow></msqrt><mtext>                             </mtext><mo stretchy="false">(</mo><mn>9</mn><mo stretchy="false">)</mo></mrow>
4. Определение значения мгновенной частоты для каждой моды ЭКС:
<mrow />
5. Представление ЭКС в виде поверхности в системе координат амплитуда-частота-время. После выполнения преобразования Гильберта в каждой эмпирической моде сигнал в виде поверхности в системе координат амплитуда-частота-время может быть выражен следующим образом [6, 13]:
k
6. Преобразование поверхности в системе координат амплитуда-частота-время к виду поверхности в системе координат энергия-частота-время. Чтобы представить энергетическую плотность сигнала и построить поверхность в системе координат энергия-частота-время мгновенной амплитуде α{t) в выражении (9) придают квадратичную форму:
Трехмерная поверхность сигнала в координатах энергия-частота-время представляет собой поверхность энергетической плотности (ПЭП). Переход от представления сигнала в системе координат амплитуда-частота-время (см. выражение 11) к представлению сигнала в системе координат энергия-частота-время (см. выражение 12) обеспечивает более точное формирование пороговой поверхности. Визуальное представление ПЭП, построенной на основе преобразования Гильберта-Хуанга, характеризует распределение мгновенной энергии ЭКС и участка изолинии в каждой точки частотно-временной плоскости и, по мнению авторов, открывает возможность выделения новых диагностических признаков в ЭКС.
На фигуре 10 (а, б) приведены ПЭП кардиоцикла ЭКС (а) и участка изолинии (б), рассчитанные по выражению (12) согласно алгоритму построения ПЭП (см. фигуру 9).
<mrow><mover accent="true"><mrow><msub><msup><mi>E</mi><mo>′</mo></msup><mrow><msub><mi>f</mi><mi>j</mi></msub></mrow></msub></mrow><mo stretchy="true">¯</mo></mover><mrow><mo>(</mo><mi>t</mi><mo>)</mo></mrow></mrow>
j
k
На следующем этапе «Формирование пороговой поверхности» осуществляется вычисление среднего значения энергии
ПЭП участка изолинии для каждой частоты f в течение заданного промежутка времени [1, t] согласно следующему выражению:
<mrow><mover accent="true"><mrow><msub><msup><mi>E</mi><mo>′</mo></msup><mrow><msub><mi>f</mi><mi>j</mi></msub></mrow></msub></mrow><mo stretchy="true">¯</mo></mover><mrow><mo>(</mo><mi>t</mi><mo>)</mo></mrow><mo>=</mo><mfrac><mn>1</mn><mrow><msub><mi>t</mi><mi>k</mi></msub></mrow></mfrac><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>i</mi><mo>=</mo><mn>1</mn></mrow><mrow><msub><mi>t</mi><mi>k</mi></msub></mrow></munderover><mrow><msub><msup><mi>E</mi><mo>′</mo></msup><mrow><msub><mi>f</mi><mi>j</mi></msub></mrow></msub><mo stretchy="false">(</mo><msub><mi>t</mi><mi>i</mi></msub></mrow></mstyle><mo stretchy="false">)</mo><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mo stretchy="false">(</mo><mn>1</mn><mn>3</mn><mo stretchy="false">)</mo></mrow>
1
где t - отсчет, принятый за середину сегмента TP;
k
t - отсчет, принятый за начало зубца Р следующего кардиоцикла.
<mrow><mover accent="true"><mrow><msub><msup><mi>E</mi><mo>′</mo></msup><mrow><msub><mi>f</mi><mi>j</mi></msub></mrow></msub></mrow><mo stretchy="true">¯</mo></mover><mrow><mo>(</mo><mrow><msub><mi>t</mi><mi>i</mi></msub></mrow><mo>)</mo></mrow></mrow>
<mrow><mover accent="true"><msup><mi>E</mi><mo>′</mo></msup><mo stretchy="true">¯</mo></mover><mrow><mo>(</mo><mrow><msub><mi>f</mi><mi>j</mi></msub></mrow><mo>)</mo></mrow></mrow>
Рассчитанные средние значения энергии
ПЭП участка изолинии по формуле (13) позволяют получить функцию распределения
среднего значения энергии по частоте:
<mrow><mover accent="true"><msup><mi>E</mi><mo>′</mo></msup><mo stretchy="true">¯</mo></mover><mrow><mo>(</mo><mrow><msub><mi>f</mi><mi>j</mi></msub></mrow><mo>)</mo></mrow><mo>=</mo><mi>ψ</mi><mo stretchy="false">(</mo><mover accent="true"><mrow><msub><msup><mi>E</mi><mo>′</mo></msup><mrow><mtext> </mtext><msub><mi>f</mi><mi>j</mi></msub></mrow></msub></mrow><mo stretchy="true">¯</mo></mover><mo stretchy="false">(</mo><mi>t</mi><mo stretchy="false">)</mo><mo stretchy="false">)</mo><mi></mi><mtext>                                   </mtext><mo stretchy="false">(</mo><mn>14</mn><mo stretchy="false">)</mo></mrow>
где t - const.
<mrow><mover accent="true"><mrow><msub><msup><mi>E</mi><mo>′</mo></msup><mrow /></msub></mrow><mo stretchy="true">¯</mo></mover><mrow><mo>(</mo><mrow><msub><mi>f</mi><mi>j</mi></msub></mrow><mo>)</mo></mrow></mrow>
Результатом расчета выражения (14) является функция распределения
характеризующая усредненное сечения ПЭП участка изолинии (см. фигуру 11а).
<mrow><mover accent="true"><msup><mi>E</mi><mo>′</mo></msup><mo stretchy="true">¯</mo></mover><mrow><mo>(</mo><mrow><msub><mi>f</mi><mi>j</mi></msub></mrow><mo>)</mo></mrow></mrow>
После вычисления функции распределения среднего значения энергии по частоте
, формирование пороговой поверхности сводится к сумме средних значений энергий:
<mrow><msub><mi>V</mi><mrow><msub><mi>f</mi><mi>j</mi></msub></mrow></msub><mo stretchy="false">(</mo><msub><mi>t</mi><mi>i</mi></msub><mo stretchy="false">)</mo><mo>=</mo><mover accent="true"><msup><mi>E</mi><mo>′</mo></msup><mo stretchy="true">¯</mo></mover><mo stretchy="false">(</mo><msub><mi>f</mi><mi>j</mi></msub><mo stretchy="false">)</mo><mo>*</mo><msub><mi>t</mi><mi>k</mi></msub><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mi></mi><mo stretchy="false">(</mo><mn>1</mn><mn>5</mn><mo stretchy="false">)</mo></mrow>
k
1
k
где t - значение величины временного интервала выделенного кардиоцикла ЭКС [t, t].
<mrow><msub><mi>V</mi><mrow><msub><mi>f</mi><mi>j</mi></msub></mrow></msub><mo stretchy="false">(</mo><msub><mi>t</mi><mi>i</mi></msub><mo stretchy="false">)</mo></mrow>
Результатом вычисления формулы (15) является пороговая поверхность
выделенного участка изолинии кардиоцикла электрокардиосигнала в системе координат энергия-частота-время, представленная на фигуре 11,б.
Таким образом, сформированная пороговая поверхность находится в одной системе координат и одном масштабе времени с ПЭП кардиоцикла ЭКС, что позволяет сравнить две поверхности без дополнительных их преобразований.
Следующий этап «Сравнение ПЭП выделенного кардиоцикла электрокардиосигнала со сформированной пороговой поверхностью» заключается в совмещении сформированной пороговой поверхности помехи с ПЭП кардиоцикла ЭКС, выделении области помехи под пороговой поверхностью и ее удаление.
fj
i
Выделение и удаление области помехи осуществляется путем применения правила обработки сигналов на основе пороговых процедур (пороговая обработка или thresholding). Пороговая обработка применяется к каждому i-му элементу ПЭП Е(t) и заключается в принятии решения о сохранении значения энергии рассматриваемого элемента ПЭП или ее пересчете [14].
Среди известных простейшим в реализации и наиболее эффективным решающим правилом для подавления помех на ЭКС является правило мягкого трешолдинга (soft thresholding) [15]. Для выполнения этапа сравнения ПЭП кардиоцикла ЭКС и пороговой поверхности необходимо сформулировать правило объемного мягкого трешолдинга, которое может быть выражено следующим образом:
<mrow><msub><mi>E</mi><mrow><msub><mi>f</mi><mi>j</mi></msub></mrow></msub><mo stretchy="false">(</mo><msub><mi>t</mi><mi>i</mi></msub><mo stretchy="false">)</mo><mrow><mo>{</mo><mrow><mtable><mtr><mtd><mrow><msub><mi>E</mi><mrow><msub><mi>f</mi><mi>j</mi></msub></mrow></msub><msub><mrow /><mrow><msub><mrow /><mtext></mtext></msub></mrow></msub><mo stretchy="false">(</mo><msub><mi>t</mi><mi>i</mi></msub><mo stretchy="false">)</mo><mo>−</mo><msub><msup><mi>E</mi><mo>′</mo></msup><mrow><mtext></mtext><msub><mi>f</mi><mi>j</mi></msub></mrow></msub><mo stretchy="false">(</mo><msub><mi>t</mi><mi>i</mi></msub><mo stretchy="false">)</mo><mo>,</mo><mi></mi><mi>е</mi><mi>с</mi><mi>л</mi><mi>и</mi><mi></mi><msub><mi>E</mi><mrow><msub><mi>f</mi><mi>j</mi></msub></mrow></msub><mo stretchy="false">(</mo><msub><mi>t</mi><mi>i</mi></msub><mo stretchy="false">)</mo><mo>></mo><msub><msup><mi>E</mi><mo>′</mo></msup><mrow><msub><mi>f</mi><mi>j</mi></msub><msub><mrow /><mtext></mtext></msub></mrow></msub><mo stretchy="false">(</mo><msub><mi>t</mi><mi>i</mi></msub><mo stretchy="false">)</mo></mrow></mtd></mtr><mtr><mtd><mrow><mi></mi><mn>0,</mn><mi></mi><mtext> </mtext><mi>е</mi><mi>с</mi><mi>л</mi><mi>и</mi><mi></mi><mi></mi><msub><mi>E</mi><mrow><msub><mi>f</mi><mi>j</mi></msub></mrow></msub><mo stretchy="false">(</mo><msub><mi>t</mi><mi>i</mi></msub><mo stretchy="false">)</mo><mo><</mo><msub><msup><mi>E</mi><mo>′</mo></msup><mrow><msub><mi>f</mi><mi>j</mi></msub><msub><mrow /><mtext></mtext></msub></mrow></msub><mo stretchy="false">(</mo><msub><mi>t</mi><mi>i</mi></msub><mo stretchy="false">)</mo><mi></mi></mrow></mtd></mtr></mtable></mrow></mrow><mtext>                  </mtext><mo stretchy="false">(</mo><mn>16</mn><mo stretchy="false">)</mo></mrow>
<mrow><msub><mi>E</mi><mrow><msub><mi>f</mi><mi>j</mi></msub></mrow></msub><mo stretchy="false">(</mo><msub><mi>t</mi><mi>i</mi></msub><mo stretchy="false">)</mo></mrow>
где
- ПЭП выделенного кардиоцикла ЭКС.
Под объемным мягким трешолдингом понимается уменьшение значения энергии ПЭП кардиоцикла ЭКС на величину значения энергии пороговой поверхности помехи (значения энергии, содержащую только помеховую составляющую сигнала).
Таким образом согласно выражению (16), по мнению авторов, путем выполнения правила объемного мягкого трешолдинга можно удалить помеховую составляющую кардиоцикла ЭКС. Выполнение данного этапа предлагаемого способа позволит уменьшить уровень помех в ЭКС (см. фигуру 13).
Заключительный этап «Восстановление и вывод электрокардиосигнала» аналогичен этапу известного способа [4]. Выполнение данного этапа необходимо для восстановления ЭКС с целью последующего измерения уровня помех и проверки его соответствия допустимому уровню.
Согласно выражению (6), для количественной оценки помехоподавления ЭКС предлагаемым способом осуществляется вычисление значения среднеквадратической ошибки отклонения восстановленного ЭКС от эталонного ЭКС. У восстановленного ЭКС, показанного на фигуре 13, значение среднеквадратической ошибки (СКО) равно 6,2%, что существенно меньше СКО известного способа подавления помех и удаления артефактов в электрокардиосигнале.
Анализ результатов работы известного и предлагаемого способа показывает, что СКО отклонения восстановленного ЭКС от эталонного в предлагаемом способе адаптивного подавления помех в электрокардиосигнале в более чем четыре раза (!) меньше СКО известного способа подавления шумов и удаления артефактов в электрокардиосигнале.
Следовательно, за счет формирования пороговой поверхности участка изолинии кардиоцикла ЭКС и его сравнения с ПЭП кардиоцикла ЭКС в предлагаемом способе адаптивного подавления помех в электрокардиосигнале обеспечивается более эффективное помехоподавление и более точное восстановление ЭКС. Восстановленный таким образом ЭКС (см. фигуру 13) является более приближенным к исходному ЭКС (см. фигуру 3,а), по сравнению с восстановленным ЭКС известного способа (см. фигуру 6).
Приведенное описание способа адаптивного подавления помех в электрокардиосигнале, по мнению авторов предлагаемого изобретения, показывает, что предлагаемый способ адаптивного подавления помех в электрокардиосигнале позволяет устранить недостаток известного способа подавления помех и удаления артефактов в электрокардиосигнале: повысить точность восстановления электрокардиосигнала и уменьшить уровень остаточных помех в электрокардиосигнале.
Другим преимуществом предлагаемого способа адаптивного подавления помех в электрокардиосигнале является визуальное представление ПЭП ЭКС в системе координат энергия-время-частота, что позволяет выявить новые диагностические признаки в электрокардиосигнале. На фигуре 14 (а, б) показаны ПЭП кардиоцикла ЭКС, согласно которым видно, что на поверхности в системе координат энергия-частота-время отчетливо выделяются зубцы Р, Т, кардиокомплекс QRS и помехи, присутствующие при регистрации электрокардиосигнала.
Эффективное подавление помех при регистрации и обработке электрокардиосигнала ведет к увеличению достоверных заключений и, следовательно, к повышению эффективности диагностики и лечения сердечно-сосудистых заболеваний.
Источники информации
1. Рангайян P.M. Анализ биомедицинских сигналов. Практический подход / Пер. с англ. под ред. А.П. Немирко. М.: ФИЗМАТЛИТ, 2007 - 440 с.
2. Щербаков М.А. Нелинейная фильтрация сигналов и изображений: Учебн. пособие. Пенза: ПТУ, 1999 - 166 с.
3. Method for signal denoising using continuous wavelet: патент США №20100010780 / Hailong Zhu/14.10.2010 - 8 p.
4. Denoising and artifact rejection for cardiac signal in a senses system: патент США №20080269628 / Detlef W., Honghuan Z., Bryon P. / 30.10.2008 - 9 p.
5. Huang N.E., Shen Z., Long S.R. The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis // Proc. R. Soc. Lond. A., Vol.454, 1998 - pp.903-995.
6. Huang N.E., Shen S.S. The Hilbert-Huang transform and its applications - World Scientific Publication, 2005 - 323 p.
7. Huang N.E., Attoh-Okine Nii O. The Hilbert-Huang transform in engineering - Taylor and Francis, 2005 - 95 p.
8. Корнейчук Н.П. Сплайны в теории приближения. - М.: Наука, Главная редакция физико-математической литературы, 1984 - 352 с.
9. Cohen A., Zigel Y. Compression of Multichannel ECG Through Multichannel Long Term Prediction. IEEE BME magazine, Vol.17, 1998 - pp.109-115.
10. Орлов B.H. Руководство по электрокардиографии. 3-е изд. - M.: Медицинское информационное агентство, 2003 - 528 с.
11. Способ выделения начала кардиоцикла и устройство для его осуществления: патент РФ №2294139 / А.Н.Митрошин, О.Н.Бодин, Д.С.Логинов, И.О.Жулев, В.В.Прошкин // Изобретения. Полезные модели: официальный бюллетень - №6, от 27.02.2007.
12. Зудбинов Ю.И., Азбука ЭКГ. Изд-во - М.: «Медицина для вас», 2003 - 160 с.
13. Sharpley R.C., Vatchev V.D. Analysis of the intrinsic mode functions. Department of Mathematics University of South Carolina, 2004 - 31 p.
14. Donoho D.L., Johnstone I.M. Neo-classical minimax problems, thresholding, and adaptation. Bernoulli, 1996 - pp.39-62.
15. Donoho D.L. De-noising by soft-thresholding // IEEE Trans, on Inform. Theory, 1995 - pp.613-627.