ЭЭГ фильтры и методы фильтрации сигнала

Цифровые фильтры играют важную роль в анализе ЭЭГ. Функция фильтров состоит в подавлении данной частотной компоненты в сигнале и пропускании первоначального сигнала за исключением этих подавленных компонент. В отличие от аналоговых фильтров работа цифровых выполняется математическими действиями на сигнал.

  • Режекторный фильтр очень точный фильтр, который уменьшает сигнал определенной частоты. В ЭЭГ используется режекторный фильтр 50 (60) Гц для устранения шума электрической сети комнаты.
  • Фильтр высоких частот подавляет низкие частоты ЭЭГ-сигнала, пропуская высокие. Фильтр характеризуется нижней граничной частотой, которую принято обозначать в виде эквивалентной постоянной времени (в секундах).
  • Фильтр низких частот — аналоговый или цифровой фильтр, который подавляет высокие частоты ЭЭГ-сигнала и пропускает низкие.
  • Фильтр подавления синфазного сигнала (CMR)

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

КИХ и БИХ-фильтры

В зависимости от алгоритмов расчета существуют два типа фильтров: фильтры с конечной импульсной характеристикой (КИХ) (finite impulse response — FIRи фильтры с бесконечной импульсной характеристикой (БИХ) (infinite impulse response — IIR). Последние также называются рекурсивными фильтрами. Импульсная характеристика цифрового фильтра — пропускание фильтра, когда на его входе есть единичный сигнал (единичная величина 1 во время t = 0). Ответ КИХ-фильтра имеет конечную продолжительность. Ответ БИХ-фильтра теоретически (но не фактически) продолжается бесконечно из-за рекурсивной природы фильтра.

Для фильтров КИХ текущий отсчет (уn) рассчитывается исключительно от текущих и предыдущих входных значений (хn, хn-1, хn-2). Напротив, в рекурсивных фильтрах в дополнение к значению входного сигнала также используют значения предыдущих результирующих сигналов. Слово «рекурсивное» буквально означает «бегущий назад» и относится к факту, что предварительно расчетные значения выходного сигнала возвращаются при вычислении последующих выходных сигналов. Поэтому выражение рекурсивного фильтра содержит не только входные значения (хn, хn-1, хn_2), но и выходные (уn_р уn2). Для достижения требуемой амплитудно-частотной характеристики рекурсивные фильтры имеют более низкий порядок по сравнению с аналогичными нерекурсивными, и, таким образом, используется меньше элементарных вычислительных операций. Чтобы отфильтровать 50 (60) Гц шум сетевой наводки, применяется режекторный фильтр.

Цифровая фильтрация сигналов

Основные понятия и определения

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

y(n)=\sum_{i=1}^{I} b_{i} x(n-i)+\sum_{j=1}^{J} a_{j} y(n-j),

где х(n) и у(n) — соответственно отсчеты входной и выходной последовательностей, a bi, aj — коэффициенты фильтра. Приведенная формула называется обобщенным разностным уравнением цифрового фильтра и полностью определяет его.

Если aj = 0 для любых значений j, то фильтр называется нерекурсивным, или КИХ-фильтром (фильтром с конечной импульсной характеристикой). В противном случае фильтр называется рекурсивным, или БИХ-фильтром (фильтром с бесконечной импульсной характеристикой).

Фильтры, все коэффициенты которых являются константами, называются фильтрами с постоянными коэффициентами. Фильтры, коэффициенты которых могут изменяться во времени, называются адаптивными.

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

Частотные свойства ЦФ определяются его амплитудно-частотной и фазово-частотной характеристиками (АЧХ и ФЧХ). АЧХ показывает зависимость коэффициента передачи (усиления) фильтра от частоты, а ФЧХ — зависимость вносимого фильтром фазового сдвига (задержки) от частоты.

Классификация цифровых фильтров по виду амплитудно-частотных характеристик

Важнейшей характеристикой цифрового фильтра (ЦФ) является его АЧХ, являющаяся модулем комплексной передаточной функции и показывающая зависимость коэффициента передачи фильтра от частоты. Передаточная функция рассчитывается с использованием специализированного математического аппарата, называемого z-преобразованием. ЦФ принято классифицировать по виду их АЧХ.

Диапазон частот АЧХ, где ее значения близки к единице, принято называть полосой пропускания фильтра. Диапазон частот, где АЧХ приближается к нулю, называют полосой задержки, а область частот в промежутке между полосами пропускания и задержки называют переходной полосой. Фильтр может иметь по нескольку полос каждого из перечисленных видов. Точки, где АЧХ пересекает √2/2 уровень (примерно 0,7, что соответствует передаче половины мощности сигнала), называются частотами среза фильтра. На рис. 5 приведен пример амплитудно-частотной характеристики фильтра, на которой отмечены перечисленные характерные элементы АЧХ.

Пример АЧХ цифрового фильтра
Рис. 5. Пример АЧХ цифрового фильтра (фильтра нижних частот).

На рисунке отмечены характерные полосы частот, а также частота среза Fср

Примеры амплитудно-частотных характеристик наиболее распространенных типов фильтров:
Рис. 6. Примеры амплитудно-частотных характеристик наиболее распространенных типов фильтров

а — фильтр нижних частот, ФНЧ; б — фильтр верхних частот, ФВЧ; в — полосовой пропускающий фильтр; г — полосовой заграждающий фильтр; д — фильтр-резонатор; е — режекторный фильтр

Примеры амплитудно-частотных характеристик дифференцирующего (слева) и интегрирующего (справа) ЦФ.
Рис. 7. Примеры амплитудно-частотных характеристик дифференцирующего (слева) и интегрирующего (справа) ЦФ.

Пунктирными линиями показаны идеальные характеристики, сплошными — реальные. Фильтры рассчитаны на частоту дискретизации fn = 500 Гц

На рис. 6 приведены примеры амплитудно-частотных характеристик нескольких наиболее распространенных типов фильтров.

Фильтры, АЧХ которых приведены на рис. 6, решают задачу пропускания одних частотных составляющих сигнала и подавления других. Помимо таких ЦФ широкое применение в цифровой обработке сигналов находят фильтры, моделирующие операции дифференцирования и интегрирования (см. ниже). В идеальном случае амплитудно-частотная характеристика дифференциатора должна иметь вид линейно возрастающей функции, а интегратора гиперболы4. Характеристики реализуемых на практике цифровых дифференциаторов и интеграторов отличаются от идеальных и приближаются к ним только в определенном диапазоне частот, который называют областью дифференцирования или областью интегрирования (рис. 7).

Методы анализа цифровых фильтров

Ранее линейный цифровой фильтр был определен его обобщенным разностным уравнением:

y(n)=\sum_{i=1}^{I} b_{i} x(n-i)+\sum_{j=1}^{J} a_{j} y(n-j) (2.1)

где х(n) и у(n) — соответственно отсчеты входной и выходной последовательностей, bi, aj — коэффициенты фильтра. Как можно видеть, процедура цифровой фильтрации представляет собой комбинацию из грех элементарных операций: сложения, умножения и задержки (сдвига). Очень удобной и наглядной формой описания ЦФ является структурная схема, основу которой составляют графические обозначения перечисленных операций (см. рис. 8, где а — сложение двух отсчетов; б — умножение отсчета на некоторую величину k; в — сдвиг (задержка) отсчета на m тактов).

Рассмотрим для примера два простых фильтра, заданные уравнениями:

у(n) = 1/(8x(n)) + 1/(4x(n – I)) + 1/(4x(n – 2)) + + 1/(4x(n – 3)) + 1/(8x(n – 4)),

у(n) = у(n – 1) + 1/(2x(n)) + 1/(2x(n – 1)).

Графические обозначения элементарных операций, используемые в структурных схемах ЦФ
Рис. 8. Графические обозначения элементарных операций, используемые в структурных схемах ЦФ
Структурные схемы нерекурсивного (сверху) и рекурсивного (снизу) цифровых фильтров
Рис. 9. Структурные схемы нерекурсивного (сверху) и рекурсивного (снизу) цифровых фильтров

Первый из них является КИХ-фильтром нижних частот и известен как фильтр Хсмминга, а второй БИХ-фильтром и представляет собой формулу численного интегрирования по правилу трапеций5. На рис. 9 показаны структурные схемы обоих фильтров. Как видим, главной отличительной чертой БИХ-фильтра является наличие обратной связи, т. с. рекурсии, с чем и связано его название «рекурсивный». Структурные схемы ЦФ нс только являются наглядной формой их описания, но могут также использоваться для эквивалентных преобразований фильтров.

О воздействии ЦФ на сигнал можно судить по его импульсной и переходной характеристикам, которые являются откликом соответственно на единичный импульс и на единичное ступенчатое воздействие. На рис. 10 показаны примеры импульсных и переходных характеристик рассмотренных ЦФ. Обратите внимание на то, что импульсная характеристика КИХ-фильтра повторяет набор его коэффициентов (рис. 10, б), а импульсная характеристика БИХ-фильтра бесконечна (рис. 10, в), с чем и связано название этого фильтра. Ряд графиков на рис. 10, а показывает соответственно единичное импульсное и ступенчатое воздействия, подаваемые на входы фильтров.

Пример импульсных (слева) и переходных (справа) характеристик рассмотренных в настоящем разделе фильтров
Рис. 10. Пример импульсных (слева) и переходных (справа) характеристик рассмотренных в настоящем разделе фильтров

Частотные свойства ЦФ описываются его передаточной функцией, которая позволяет определить АЧХ и ФЧХ. АЧХ показывает зависимость коэффициента передачи (усиления) от частоты, а ФЧХ — зависимость вносимого фильтром фазового сдвига (задержки) от частоты.

Передаточная функция ЦФ может быть рассчитана по его разностному уравнению с использованием математического аппарата z-преобразования, которое широко используется при анализе импульсных и дискретных систем678. Смысл z-преобразования можно кратко описать следующим образом.

Пусть задана дискретная последовательность х(n), n = 0, 1,2, … …. N. Тогда ^-преобразованием этой последовательности называется выражение

X(z)=\sum_{n=0}^{N} x(n)z^{-n}

Для ЦФ, заданного разностным уравнением вида (2.1), z-преобразование его передаточной функции имеет вид

H(z)=\sum_{i=0}^{I} b_{i} z^{-i} /\left(1-\sum_{j=1}^{J} a_{j} z^{-j}\right)

Аналитическое выражение для зависимости передаточной функции от частоты H(w) может быть получено путем подстановки z = ejwT, где w — угловая частота в радианах, а Т = 1/fд — интервал дискретизации в секундах (здесь fд — частота дискретизации). Функция H(w)— комплексная. Ее модуль | H(\omega)=\sqrt{\operatorname{Re}^{2}[H(\omega)]+\operatorname{Im}[H(\omega)]} и аргумент = \angle H(\omega)=\operatorname{arctg}\{\operatorname{Im}[H(\omega)] / \operatorname{Re}|H(\omega)|\} дают выражения для АЧХ и ФЧХ фильтра соответственно. Для перехода к более удобной в практических задачах циклической частоте f следует выполнить подстановку w = 2πf.

Рассчитаем передаточные функции для рассмотренных выше ЦФ. Выражение для z-преобразования для КИХ-фильтра имеет вид

H(z) = 1/8z0 + 1/4z-1 + 1/4z-2 + 1/4z-3 + 1/8z-4.

Подставим z = ejwT:

H(z) = 1/8 + 1/4e-jwT   + 1/4e-2jwT + 1/4e-3jwT + 1/8e-4jwT = 1/4e-2jwT(1/2e2jwT +ejwT + 1 + 4e-4jwT + 1/2e-2jwT) = 1/4e-2jwT[cos(2wT) + 2cos(wT) + 1].

Последнее преобразование выполнено на основании одной из формул Эйлера, связывающих тригонометрические функции синус и косинус с комплексной экспонентой:

cos x = (ejx + e-jx) / 2, sin x = (ejx – e-jx) / 2j.

Тогда выражения для модуля и аргумента передаточной функции (АЧХ и ФЧХ фильтра) соответственно будут иметь вид

|H(w)| = l/4[cos(2wT) + 2cos(wT) +1];

∠H(w) = -2wТ.

С целью упрощения математических преобразований их обычно выполняют для нормализованной частоты w = 2π радиан, или f = 1 Гц, что соответствует Т = 1 с. По окончании преобразований переходят обратно к реальным частотам, учитывая частоту дискретизации f д. В этом случае выражения для ЛЧХ и ФЧХ для циклической частоты принимают вид

|Н(f )| =  1/4|сos(4πf /f д) + 2соs(2πf /f д) + 1];

∠H(f) = 4πf /f д.

Выполним аналогичные преобразования для рассматриваемого БИХ-фильтра:

H(z)=0,5 \frac{1+z^{-1}}{1-z^{-2}}

H(\omega)=0,5 \frac{1+e^{-j \omega T}}{1-e^{-j \omega T}}=0,5 \frac{e^{-j \omega T / 2}\left(e^{j \omega T / 2}+e^{-j \omega T / 2}\right)}{e^{-j \omega T / 2}\left(e^{j \omega T / 2}-e^{-j \omega T / 2}\right)}=

=0,5 \frac{\cos (\omega T / 2)}{j \sin \left(\omega T / 2\right)}=-0,5 j \operatorname{ctg}(\omega T / 2)

Здесь в преобразованиях использована вторая формула Эйлера (для синуса). Тогда окончательные выражения для АЧХ и ФЧХ фильтра примут вид

|H(f)|=0,5ctg(πf / f д);

∠H(f) = -π2.

На рис. 11 приведены графики АЧХ и ФЧХ обоих фильтров для случая f n = 200 Гц (рис. 11, а и б соответственно). Как видно из приведенных характеристик, первый фильтр является фильтром нижних частот (сглаживающим фильтром), в то время как АЧХ второго на начальном участке частотного диапазона близка к характеристике идеального интегратора. Важно отметить также, что в соответствии с теоремой отсчетов преобразования с дискретным сигналом корректны только в пределах диапазона частот от 0 Гц до поэтому в данном примере оси частот ограничены полосой от 0 до 100 Гц.

Графики ЛЧХ (слева) и ФЧХ (справа) КИХ- и БИХ-фильтров
Рис. 11. Графики ЛЧХ (слева) и ФЧХ (справа) КИХ- и БИХ-фильтров

Нетрудно заметить, что числитель и знаменатель выражения для передаточной функции ЦФ являются полиномами переменной z. Для каждого из них могут быть найдены корни, которые для числителя носят название нулей, а для знаменателя — полюсов. Нули и полюса передаточной функции представляют собой комплексные числа и равняются таким значениям переменной z, при которых передаточная функция обращается соответственно в нуль и в бесконечность.

Передаточная функция для КИХ-фильтра содержит только выражение для числителя, поэтому для фильтров этого типа существуют только нули — корни полинома для числителя. Формально для КИХ-фильтра также определен один полюс со значением z = 0. но его обычно не рассматривают. В приведенном примере КИХ-фильтра полином в числителе передаточной функции имеет четвертую степень. Решение уравнения дает следующие значения нулей:

z10 = j; z20 = -j; z3,40 = -1

Для БИХ-фильтров существуют и нули, и полюса. В рассмотренном примере БИХ-фильтр имеет один ноль и один полюс, которые соответственно принимают следующие значения:

z10 = 1; z10 = -1..

Широкое применение при анализе цифровых систем находят графики, показывающие расположение нулей и полюсов на комплексной плоскости (карты нулей и полюсов). Нули на графиках принято обозначать кружочками, а полюса звездочками. На графиках всегда отображают единичную окружность, очерчиваемую проведенным из начала координат вектором с длиной, равной единице (единичным вектором). Единичная окружность определяется уравнением z = ejwT. Угол поворота вектора против часовой стрелки начиная от положительного направления вещественной оси можно интерпретировать как нормализованную угловую частоту Если, учитывая фактически используемую частоту дискретизации сигнала fд, перейти к реальным частотам (как показано выше для ЛЧХ и ФЧХ), то полному обороту

вектора соответствует изменение частоты от 0 Гц до fn. Нетрудно продемонстрировать9, что при приближении конца вектора к какому-нибудь из нулей модуль передаточной функции уменьшается, а при приближении к полюсам — возрастает. Благодаря отмеченным обстоятельствам визуальный анализ карты нулей и полюсов позволяет, не выполняя сложных аналитических преобразований, не только составить представление о виде ЛЧХ и ФЧХ ЦФ, но также приблизительно рассчитать их графики, пользуясь циркулем, транспортиром и линейкой.

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

На рис. 12 приведены карты нулей и полюсов: а — для КИХ-фильтра; б — для БИХ-фильтра. Сопоставляя их с ЛЧХ фильтров (см. рис. 11) и учитывая описанную интерпретацию вращения единичного вектора, можно легко проследить причины возникновения максимумов и минимумов ЛЧХ на определенных частотах. Анализ устойчивости с помощью приведенного выше правила показывает, что оба фильтра устойчивые. Отметим, что КИХ-фильтры устойчивы всегда, что относится к важнейшим их достоинствам.

Карты нулей и полюсов КИХ- (а) и БИХ-фильтров (б)
Рис. 12. Карты нулей и полюсов КИХ- (а) и БИХ-фильтров (б)

В рассмотренных примерах ЦФ аналитические выражения для ЛЧХ и ФЧХ удалось вывести из передаточной функции относительно просто, используя формулы Эйлера. Однако на практике часто используются фильтры с таким набором коэффициентов, что вывод аналитических формул оказывается чрезвычайно громоздкой процедурой. В настоящее время эта проблема без труда преодолевается благодаря наличию специальных программ расчета АЧХ и ФЧХ цифровых систем.

Одним из основных достоинств КИХ-фильтров является относительная простота обеспечения линейности фазовой характеристики (т.е. постоянства временной задержки для составляющих всех частот), что, в частности, может быть очень важным при анализе электрофизиологических сигналов, так как нелинейность фазовой характеристики может приводить к изменению взаимного расположения элементов сигнала или даже к появлению «лишних» колебаний.

Линейность фазовых характеристик КИХ-фильтров обеспечивается при соблюдении одного из четырех условий симметрии импульсной характеристики1011:

  • симметричность при нечетном числе коэффициентов;
  • симметричность при четном числе коэффициентов;
  • антисимметричность при нечетном числе коэффициентов;
  • антисимметричность при четном числе коэффициентов.

В табл. 1 сравниваются характеристики КИХ- и ВИХ-фильтров. В практических задачах для выбора какого-либо из этих двух видов ЦФ необходимо исходить из конкретных требований и условий их применения.

Таблица 1

ХарактеристикиТип цифровых фильтров
КИХБИХ
ДостоинстваВсегда устойчивы.

Невозможно накопление ошибок округления.

Простота обеспечения линейности ФЧХ

Относительно небольшое число коэффициентов для обеспечения заданного вида передаточной функции
НедостаткиБольшое число коэффициентов для обеспечения заданного вида передаточной функцииВозможна неустойчивость.

Возможно накопление ошибок округления.

Сложность или иногда невозможность обеспечения линейности ФЧХ

Фильтры Баттерворта

Основным достоинством рекурсивных фильтров является возможность получения крутого спада ЛЧХ при использовании относительно небольшого числа коэффициентов фильтра. Поэтому, несмотря на наличие серьезных недостатков (возможная неустойчивость, накопление ошибок округления, нелинейность фазово-частотной характеристики), рекурсивные фильтры применяются очень часто. Среди них одними из самых распространенных являются так называемые фильтры Баттерворта. Их отличительной особенностью является монотонное изменение АЧХ в переходной полосе частот и практически полное отсутствие пульсаций в полосах пропускания и задержки. В некоторых применениях это оказывается очень важным, так как позволяет гарантировать постоянство коэффициента передачи фильтра во всей рабочей полосе частот и обеспечить надежное подавление устраняемых фильтром частотных составляющих.

Вид амплитудно-частотной характеристики фильтра нижних частот Баттерворта определяется, исходя из условия:

|H(\omega)|^{2}=\frac{1}{1+\left(\omega / \omega_{\mathrm{c}}\right)^{2}}

где wc— частота среза фильтра в радианах (уровень √2/2 на АЧХ), а N порядок фильтра, определяющий число коэффициентов в рекурсивной и нерекурсивной частях разностного уравнения, которое в этом случае будет иметь вид:

y(n)=\sum_{i=1}^{I} b_{i} x(n-i)+\sum_{j=1}^{J} a_{j} y(n-j)

Чем выше порядок фильтра Баттерворта, тем более крутой спад имеет его АЧХ в переходной полосе. В рекурсивных фильтрах достаточно крутой спад может быть достигнут при относительно небольшом числе коэффициентов. Однако крутой спад АЧХ не всегда можно рассматривать как положительное свойство фильтра. Рассмотрим пример фильтрации фрагмента ЭКГ с использованием фильтров Баттерворта. На рис. 13 приведены ЛЧХ двух фильтров нижних частот, имеющих одинаковую частоту среза 30 Гц: с порядком 2 и с порядком 8. Последний фильтр имеет АЧХ, быстро спадающую в переходной полосе.

На графике «а» на рис. 14 показан фрагмент ЭКГ длительностью 2 с, частота дискретизации сигнала здесь составляет 250 Гц. Графики «б» и «в» представляют тот же сигнал, прошедший соответственно через ФНЧ 2-го и 8-го порядков, ЛЧХ которых приведены на рис. 13. На нижнем из графиков отчетливо наблюдается так называемый «звон» — высокочастотные затухающие колебания, следующие сразу за участком сигнала, содержащим резкий всплеск (в данном случае — QRS-комплекс ЭКГ). При последующем анализе ЭКГ этот «звон» может быть ошибочно интерпретирован как S-волна ЭКГ, что может привести к неверному диагностическому заключению.

Амплитудно-частотные характеристики фильтров нижних частот Баттерворта
Рис. 13. Амплитудно-частотные характеристики фильтров нижних частот Баттерворта

2-го порядка (1) и 8-го порядка (2), имеющих одну и ту же частоту среза 30 Гц. Частота дискретизации равняется 250 Гц

Иллюстрация возникновения «звона» в результате фильтрации
Рис. 14. Иллюстрация возникновения «звона» в результате фильтрации

а) исходный фрагмент ЭКГ; 6) результат обработки БИХ-фильтром 2-го порядка; в) результат обработки БИХ-фильтром 8-го порядка (в сигнале наблюдается «звон»)

Характеристики фильтра верхних частот Баттерворта
Рис. 15. Характеристики фильтра верхних частот Баттерворта

4-го порядка с частотой среза 1 Гц: а) амплитудно-частотная характеристика; б) фазовочастотная характеристика (нелинейная)

Иллюстрация возникновения фазовых искажений в результате фильтрации

Рис. 16. Иллюстрация возникновения фазовых искажений в результате фильтрации

Одним из наиболее существенных недостатков БИХ-фильтров является почти всегда имеющая место нелинейность фазово-частотной характеристики, что означает непостоянство вносимой фильтром задержки для разных частот. Это может вызывать фазовые искажения сигналов, что в некоторых приложениях оказывается неприемлемым. Одним из вариантов выхода из данного затруднения является использование КИХ-фильтров, для которых нетрудно обеспечить линейность ФЧХ за счет выбора симметричного набора коэффициентов. Однако существует решение, позволяющее избежать фазовых искажений и при использовании БИХ-фильтров.

Рассмотрим задачу устранения низкочастотной помехи (дрейфа изоэлектрической линии) из ЭКГ. Так как все полезные составляющие электрокардиосигнала лежат в диапазоне частот выше I Гц, низкочастотную помеху можно существенно ослабить фильтром верхних частот с частотой среза около 1 Гц, нс затрагивая при этом саму ЭКГ. Пример АЧХ и ФЧХ такого фильтра показан на рис. 15. В данном случае предлагается использовать фильтр Баттерворта 4-го порядка. Как видно из графика «б» на рис. 15, фазово-частотная характеристика фильтра является нелинейной.

График «а» на рис. 16 показывает фрагмент ЭКГ, содержащий значительную низкочастотную помеху. Результат использования предложенного выше ФВЧ показан на графике «б». Дрейф изолинии был успешно устранен фильтром, но в то же время на графике хорошо видны искажения формы кардиоцикла, причиной которых являются фазовые искажения, вносимые фильтром. Данный недостаток может быть устранен за счет применения так называемой фильтрации с нулевым фазовым сдвигом. Суть этого приема заключается в том, что фрагмент сигнала фильтруется одним и тем же фильтром сперва в прямом, а затем в обратном направлении. Возникающие при этом фазовые искажения взаимно компенсируются и результирующий фазовый сдвиг равен нулю для всех частот. На графике «в» представлен результат такой фильтрации. Как можно видеть, заметных искажений сигнала при этом не наблюдается.

При использовании фильтрации с нулевым фазовым сдвигом следует учитывать два важных обстоятельства:

  1. так как один и тот же фильтр применяется дважды, эквивалентная амплитудно-частотная характеристика всей процедуры в целом определяется как квадрат АЧХ исходного фильтра. Частота среза при этом смещается в сторону верхних частот, а спад в переходной полосе становится круче (см. рис. 17). Для получения желаемого результата фильтрации необходимо так выбирать АЧХ фильтра, чтобы результирующая АЧХ соответствовала требуемой;
  2. фильтрация с нулевым фазовым сдвигом может быть выполнена только в тех случаях, когда весь сигнал известен заранее. Это означает, что данный прием невозможно использовать при обработке сигналов в реальном масштабе времени.

Практически такой же результат может быть получен при использовании КИХ-фильтра с линейной фазово-частотной характеристикой. На рис. 18 приведен пример амплитудно-частотной характеристики и импульсной характеристики (набора коэффициентов) КИХ-фильтра 250-го порядка, имеющего примерно такую же АЧХ, что и рассмотренный ранее ФВЧ Баттерворта. Результат использования этого фильтра показан на графике «г» рис. 16. Также, как и в случае использования фильтрации с нулевой фазой, низкочастотная помеха успешно устраняется, а существенные искажения сигнала отсутствуют. Однако использование КИХ-фильтров часто ограничивается вычислительной мощностью применяемого компьютера (данный фильтр имеет 250 коэффициентов, в то время как примерно эквивалентный БИХ-фильтр — всего 10 коэффициентов). Кроме того, КИХ-фильтр вносит существенную задержку в обрабатываемый сигнал (в данном примере она составляет 125 интервалов дискретизации), что не всегда приемлемо в системах анализа сигналов в реальном масштабе времени.

Иллюстрация изменения вида АЧХ фильтра при использовании фильтрации с нулевым фазовым сдвигом
Рис. 17. Иллюстрация изменения вида АЧХ фильтра при использовании фильтрации с нулевым фазовым сдвигом

исходная АЧХ (/) и эквивалентная АЧХ при двукратной фильтрации (2)

Характеристики КИХ-фильтра верхних частот 250-го порядка с частотой среза 1 Гц
Рис. 18. Характеристики КИХ-фильтра верхних частот 250-го порядка с частотой среза 1 Гц

а) амплитудно-частотная характеристика; б) импульсная характеристика (набор коэффициентов); здесь верхняя часть импульса, отображающего центральный коэффициент, равный 0,9936, нс поместилась на графике в выбранном масштабе

Цифровое дифференцирование и интегрирование

Дифференцирование и интегрирование являются важнейшими математическими операциями как для решения теоретических проблем, так и при практическом использовании математических методов. Для случаев, когда входные данные представлены последовательностью чисел (то есть являются дискретными), были разработаны методы численной математики, позволяющие, в частности, реализовать операции дифференцирования и интегрирования с любой заданной точностью. Многие из этих методов можно интерпретировать как варианты цифровой фильтрации 121314, в связи с чем представляет интерес их рассмотрение с точки зрения частотного подхода. В настоящем разделе рассмотрены некоторые наиболее часто используемые практические формулы, приведены их частотные характеристики, а также представлены примеры их применения в задачах цифровой обработки биомедицинских сигналов.

Для рассмотрения операции дифференцирования с позиций частотного подхода необходимо исследовать ее воздействие на гармонический (синусоидальный) сигнал, который в наиболее общем виде описывается комплексной экспонентой:

x(t) =  ejwT = cos wt + j sin wt,

где w — угловая частота в радианах.

Производная от гармонического сигнала определяется выражением: x’(t) = ejwT/dt = jwejwT.

Рассмотрим некоторую функцию Н(о) = x’(t)/x(t) = jw, которая показывает связь между входом и выходом операции дифференцирования и может интерпретироваться как передаточная функция дифференциатора по отношению к гармоническому сигналу. Модуль и аргумент данной функции являются амплитудно-частотной и фазово-частотной характеристиками и определяются выражениями:

|H(w)| = w

∠H(w) = π/2.

Это означает, что усиление дифференциатора по отношению к непрерывному гармоническому сигналу прямо пропорционально частоте, а вносимый фазовый сдвиг всегда равен π/2.

Как известно, операция вычисления производной для непрерывной функции x(t) определяется выражением:

x^{\prime}(t)=\lim _{\Delta t \rightarrow 0} \frac{x(t+\Delta t)-x(t)}{\Delta t}

где Δt — некоторый интервал времени. Если входные данные представлены последовательностью равноотстоящих отсчетов, следующих с интервалом Т, то можно предположить, что аналогом вышеприведенной формулы для случая дискретного сигнала является выражение:

y(n)= \frac{x(n)-x(n-1)}{T}= \frac{1}{T} x (n-1)

где х(n) и у(n) — соответственно отсчеты входной и выходной последовательностей, а n — целочисленный порядковый номер отсчета. Нетрудно заметить, что данная формула представляет собой нерекурсивный цифровой фильтр с двумя коэффициентами, равными 1/Т и -1/Т. Используя z-преобразование, получим выражение для передаточной функции данного фильтра:

Н( z) = (1- z-1)/Т;

Н(w) = (1 – e-jw)/Т = e-jw/2(ejw/2e-jw/2)/Т = 2je-jw/2sin(w/2)/T.

Соответственно выражения для ЛЧХ и ФЧХ будут иметь вид:

|Н(w)| = 2sin(w/2)/T;

∠H(w) = π/2 – w/2.

Кривая 2 на рис. 19 показывает вид AЧХ этого фильтра в диапазоне нормализованной частоты Fн от О до 0,5, а пунктирной линией (кривая 1) показана АЧХ идеального дифференциатора. Под нормализованной частотой понимается частота, деленная на используемую частоту дискретизации:

Fн = F/Fд = FT.

Использование нормализованной частоты является общепринятым способом упрощения математических выкладок на этапе анализа свойств цифровых фильтров и эквивалентно принятию предположения, что интервал дискретизации Т = 1 с. Тогда выражение для модуля передаточной функции принимает вид:

|H(w)| = 2sin(w/2).

При этом значению Fн= 0.5 соответствует частота свертывания (частота Найквиста) Fc = Fд/2. Графики обрываются при значениях Fн = 0,5, т. к. в соответствии с теоремой отсчетов, цифровая обработка сигнала на более высоких частотах некорректна из-за эффекта наложения частот.

Можно видеть, что формула первой разности отсчетов обладает свойствами дифференциатора только в диапазоне нормализованной частоты от 0 до примерно 0,1. Этот диапазон называют областью дифференцирования. На более высоких частотах кривая существенно отклоняется от идеальной. Фактически это приводит к усилению высокочастотных составляющих сигнала, в связи с чем часто говорят, что за пределами области дифференцирования данная формула усиливает шумы. Так как данное свойство обычно считается недостатком, формулу первой разности отчетов на практике применяют редко.

Амплитудно-частотные характеристики
Рис. 19. Амплитудно-частотные характеристики идеального дифференциатора (1). формулы первой разности отсчетов (2) и формулы центральной первой разности отсчетов (3)
Характеристики дифференцирующих фильтров
Рис. 20. Характеристики дифференцирующих фильтров

а) амплитудно-частотные характеристики идеального дифференциатора (1) и дифференцирующего фильтра 50-го порядка (2); б) импульсная характеристика (набор коэффициентов) дифференцирующего фильтра 50-го порядка

Значительно большее распространение получила формула центральной первой разности отсчетов:

y(n)= \frac{x(n)-x(n-2)}{2T}= \frac{1}{2T} x (n) - \frac{1}{2T} x (n-2)

Выражения для АЧХ и ФЧХ в этом случае имеют вид:

|Н(w)| = sin(w)/T;

∠H(w) = π/2 – w.

Соответствующий график АЧХ представлен на рис. 19 кривой 3, которая показывает, что данный фильтр имеет примерно такую же область дифференцирования, что и формула первой разности отсчетов, но усиление шумов проявляется в гораздо меньшей степени. В связи с этим, центральная первая разность часто используется для цифрового дифференцирования в тех случаях, когда нет жестких требований к математической точности данной операции.

Более точно цифровое дифференцирование может быть реализовано с использованием специально рассчитанных КИХ-фильтров. При этом число коэффициентов может составлять несколько десятков. На рис. 2.20 приведены амплитудно-частотная и импульсная характеристики цифрового дифференцирующею фильтра 50-го порядка. Отметим. что импульсная характеристика КИХ-фильтра совпадает с набором его коэффициентов. Данный фильтр имеет АЧХ близкую к характеристике идеального дифференциатора в диапазоне нормализованной частоты от 0 до приблизительно 0,14, а далее кривая резко спадает, что исключает эффект усиления помех.

Одним из наиболее широко известных применении цифрового дифференцирования в обработке биомедицинских сигналов является решение задачи выделения R-зубца ЭКГ на фоне остальных компонентов кардиоцикла и помех. График на рис. 21 демонстрирует пример фрагмента ЭКГ продолжительностью 2 с, содержащего типичные для ЭКГ помехи (сетевая наводка 50 Гц. широкополосная миографическая помеха, дрейф изоэлектрической линии). Кроме того, присутствующие в ЭКГ Т-зубцы также затрудняют выделение R-зубцов. Частота дискретизации сигнала в данном случае составляет 250 Гц, что является типичным значением для многих приборов и систем, предназначенных для контроля сердечного ритма. Справа от графика ЭКГ приведен график спектральной плотности мощности (график «б»), рассчитанной для данного фрагмента сигнала. Дрейф изолинии проявляется в спектре в виде высоких значений на частотах близких к нулю, сетевая наводка дает в спектре пик на частоте 50 Гц, а миографический шум — примерно равномерный фоновый уровень во всем диапазоне частот. Собственно компоненты кардиоцикла дают основной вклад в спектр в диапазоне от 0.5 Гц до 50 Гц. Составляющие спектра ЭКГ в области частот выше 7-10 Гц обусловлены в основном вкладом QRS-комплекса. Требуется подобрать такую процедуру фильтрации, которая позволила бы в максимальной степени подавить все компоненты сигнала, за исключением QRS-комплекса, что позволило бы обеспечить его надежное обнаружение на следующем этапе обработки сигнала. Так как QRS-комплекс является наиболее высокочастотным компонентом кардиоцикла, для выделения QRS-комплекса часто используют дифференцирующие фильтры.

Графики «л» и «г» демонстрируют результат применения формулы первой разности отсчетов. При этом оказались подавленными низкочастотные помехи и Т-волны ЭКГ. Однако отчетливо проявился свойственный данному фильтру эффект усиления шумов. С точки зрения решения задачи выделения R-зубца ЭКГ данный фильтр не дает удовлетворительною результата. Использование центральной первой разности (графики «д» и «е») лишь несущественно улучшает результат за счет меньшего влияния эффекта усиления помех.

Применение описанного выше дифференциатора 50-го порядка (графики «ж» и «з») дает наилучший результат: существенно подавлены все виды помех, а QRS-комплекс отчетливо выделен в сигнале. Отметим существенную задержку, вносимую этим фильтром в сигнал, что связано с большим числом коэффициентов фильтра. Применив несложные пороговые правила к полученному на выходе фильтра сигналу, можно получить вполне работоспособный обнаружитель R-зубца.

Недостатком последнего из описанных фильтров является его относительно высокая вычислительная сложность, что может сделать невозможным его применение в портативных устройствах контроля ЭКГ в реальном масштабе времени. Практически такой же результат с точки зрения решения задачи обнаружения QRS-комплекса ЭКГ может быть получен с использованием значительно более простого в вычислительном отношении фильтра. Примером такого простого фильтра может служить сглаживающий дифференциатор, определяемый разностным уравнением вида:

y(n)= 1/8[х(n – 8) + х(n – 7) + х(n – 6) + х(n – 5) – х(n – 3) – х(n -2) – х(n – 1) — х(n )].

На рис. 22 приведена амплитудно-частотная характеристика этого фильтра (кривая 2). Пунктиром на графике показана АЧХ идеального дифференциатора (кривая 1). Видно, что данный фильтр в диапазоне от О Гц до 15 Гц обладает свойствами дифференциатора. Далее характеристика спадает, т. е. фильтр пропускает преимущественно ту часть спектральных составляющих ЭКГ, которые относятся к QRS-комплексу. Еще одним полезным свойством фильтра является нулевой коэффициент передачи на частоте 50 Гц (это справедливо только для частоты дискретизации сигнала 250 Гц!), что означает одновременное решение как задачи выделения QRS-комплекса, так и задачи подавления сетевой наводки. Графики «и» и ок» рис. 21 демонстрируют результат применения данного фильтра. Как можно видеть, задача выделения QRS-комплекса ЭКГ на фоне остальных волн кардиоцикла и помех решается в данном случае достаточно успешно.

Фрагмент ЭКГ, содержащий помехи (о) и график спектральной плотности мощности для этого фрагмента (б).
Рис. 21. Фрагмент ЭКГ, содержащий помехи (о) и график спектральной плотности мощности для этого фрагмента (б).

Фрагмент ЭКГ и соответствующие графики спектральной плотности мощности после применения различных процедур цифровой фильтрации: формулы первой разности отсчетов (в, г); формулы центральной первой разности отсчетов (д, е) дифференциатора 50-го порядка (ж, з); сглаживающего дифференциатора (и, к)

Амплитудно-частотные характеристики идеального дифференциатора (1) и сглаживающего дифференциатора (2) для частоты дискретизации равной 250 Гц
Рис. 22. Амплитудно-частотные характеристики идеального дифференциатора (1) и сглаживающего дифференциатора (2) для частоты дискретизации равной 250 Гц
 Амплитудно-частотные характеристики идеального дифференциатора 2-го порядка (1) и формулы второй разности отсчетов (2) для нормализованной частоты
Рис. 23. Амплитудно-частотные характеристики идеального дифференциатора 2-го порядка (1) и формулы второй разности отсчетов (2) для нормализованной частоты

Взятие второй производной от гармонического сигнала дает:

x”(t) = – w2e-jvt

Тогда амплитудно-частотная и фазово-частотная характеристики соответственно определяются выражениями:

|Н(w)| = w2;

∠H(w) = π.

Кривая / на рис. ХХХ+5 показывает вид АЧХ операции дифференцирования 2-го порядка для нормализованной частоты.

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

z(n)= \frac{x(m)-x(n-1)}{T}

\begin{aligned} y(n)=\frac{z(n)-z(n-1)}{T}=& \frac{[x(n)-x(n-1)]-[x(n-1)-x(n-2)]}{T^{2}}=\\ &=\frac{1}{T^{2}} x(n)-\frac{2}{T^{2}} x(n-1)-\frac{1}{T^{2}} x(n-2) \end{aligned}

Используя z-преобразование, получим выражение для передаточной функции данного фильтра:

Н( z) = (1- 2z-1+ z-2)/Т2;

Н(w) = (1 – 2e-jw + e-j2w)/Т2 = e-jw(ejw – 2 – e-jw)/Т2 = 2e-jw[cos(w) – 1]/T2.

Соответственно выражения для АЧХ и ФЧХ будут иметь вид:

|Н(w)| = 2[cos(w) – 1]/T2;

∠H(w) = –w

Так как в случае использования нормализованной частоты Т = 1, то получим:

|Н(w)| = 2|cos(w) – 1|.

Кривая 2 на рис. 23 показывает вид АЧХ для формулы второй разности отсчетов. Как видно из графика, эта формула дает результат близкий к дифференцированию 2-го порядка в диапазоне нормализованной частоты от 0 до примерно 0,15.

С другой стороны, такая кривая соответствует АЧХ фильтра верхних частот. Это последнее свойство можно, например, использовать для оценки уровня широкополосного (миографического) шума в ЭКГ.

Как известно, основная доля мощности ЭКГ сосредоточена в области частот не выше 50 Гц15. Если выбрать частоту дискретизации равной 250 Гц, то формула второй разности практически не будет пропускать полезные составляющие электрокардиосигнала, а будет усиливать частоты, лежащие в более высокочастотном диапазоне. Так как спектр широкополосной помехи примерно равномерен во всем диапазоне частот, то по мощности высокочастотной части диапазона (заведомо не имеющей отношения к полезной ЭКГ) можно судить о содержании широкополосных помех в сигнале. На рис. 24 приведен пример использования второй разности отсчетов для решения такой задачи. На графике «о» показан фрагмент ЭКГ. в котором на участке от 0,65 с до 0.9 с наблюдается всплеск миографической помехи. Графики «б» и «в» демонстрируют соответственно результат применения формулы второй разности отсчетов и операции взятия модуля полученного сигнала. Как можно видеть, формула второй разности позволила выделить широкополосную помеху. Для большей устойчивости такого выделителя помех желательно далее применить фильтрацию нижних частот, для чего, например, можно использовать сглаживание по определенному числу последовательных точек:

y(n)=\frac{1 }{N}\sum_{i=0}^{N-1} x(n-i)

На графике «г» рис. 24 показан результат применения этого фильтра для случая N = 10. Пунктирная горизонтальная линия на этом же графике показывает возможное положение порогового значения, по превышению которого можно определять участки сигнала, где уровень миографической помехи недопустимо высок.

Пример использования второй разности отсчетов для оценки интенсивности миографической помехи в ЭКГ:
Рис. 24. Пример использования второй разности отсчетов для оценки интенсивности миографической помехи в ЭКГ

а) фрагмент ЭКГ, в котором наблюдается всплеск миографической помехи; б, в) результат применения формулы второй разности отсчетов и операции взятия модуля полученного сигнала; г) результат применения сглаживающего фильтра и возможный уровень порога (горизонтальная пунктирная линия) для определения зашумленных участков сигнала

Амплитудно-частотные характеристики идеального интегратора (1) и цифровых интегрирующих фильтров:
Рис. 25. Амплитудно-частотные характеристики идеального интегратора (1) и цифровых интегрирующих фильтров

формулы прямоугольников (2); формулы трапеций (3) и формулы Симпсона (4)

Иллюстрация принципа вычисления интеграла по формулам
Рис. 26. Иллюстрация принципа вычисления интеграла по формулам

прямоугольников (графики «а» и «г»); трапеций (графики «б» и «д») и Симпсона (графики «в» и «е»)

Аналогично операции дифференцирования для рассмотрения операции интегрирования с позиций частотного подхода необходимо исследовать ее воздействие на гармонический сигнал16.

Интеграл от гармонического сигнала определяется выражением:

{\int x(t) d t}={\int e^{jwt} d t}=\frac{1}{j \omega}e^{jwt}

Тогда передаточную функцию операции интегрирования можно выразить как:

H(\omega)=\frac{\int x(t) d t}{x(t)}=\frac{1}{j \omega}

Модуль и аргумент данной функции (АЧХ и ФЧХ) определяются выражениями:

|Н(w)| = 1/w;

∠H(w) = –π/2

Это означает, что усиление дифференциатора по отношению к непрерывному гармоническому сигналу обратно пропорционально частоте, а вносимый фазовый сдвиг всегда равен –π/2. Кривая 1 на рис. 25 показывает амплитудно-частотную характеристику идеального (аналогового) интегратора от нормализованной частоты.

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

Формула прямоугольников. Разностное уравнение, соответствующее формуле, имеет вид:

y(n) = y(n — 1) + Тх(n — 1),

где х(n) и y(n) — соответственно отсчеты входной и выходной последовательностей, Т — интервал дискретизации, а n — целочисленный порядковый номер отсчета.

Интеграл по данной формуле рассчитывается как сумма прямоугольников, образованных вертикальными прямыми, проходящими через точки взятия отсчетов и имеющих высоту, соответствующую значению отсчета х(n – 1), предшествующего текущему х(n). Графики «л» и «г» на рис. 26 иллюстрируют принцип расчета интеграла по формуле прямоугольников. На верхнем из этих графиков пунктиров показана исходная непрерывная функция, на которой точками отмечены моменты взятия дискретных отсчетов. В данном случае интервал дискретизации равен Т = 0,05 с, что соответствует частоте дискретизации Гд — 20 Гц. Заливкой отмечен прямоугольник, площадь которого рассчитывается для отсчета х(n), соответствующего здесь моменту времени t = 0,2 с. Как можно видеть, расхождение между площадью прямоугольника и площадью под соответствующим участком исходной кривой оказывается достаточно существенным. Этим объясняется заметная погрешность оценки интеграла, что иллюстрирует график «г» на рис. 2.26. Здесь сплошная линия показывает истинное значение интеграла, а точки, соединенные пунктирной линией, соответствуют результату оценки интеграла с использованием формулы прямоугольников.

Так как формула прямоугольников представляет собой рекурсивный цифровой фильтр, то с использованием z-преобразования можно рассчитать его передаточную функцию:

H(z)={T} \frac{z^{-1}}{1-z^{-1}}

H(\omega)=\frac{T}{1-e^{-j \omega T}}= \frac{Te^{-j \omega T / 2}}{e^{j \omega T / 2}-e^{-j \omega T / 2}}= \frac{Te^{-j \omega T / 2}}{2 j \sin \left(\omega T / 2\right)}

Соответственно АЧХ и ФЧХ фильтра описываются выражениями:

|Н(w)| = T/2 sin (wT/2) ;

∠H(w) = –π/2 – wT/2

Кривая 2 на рис. 25 показывает вид АЧХ формулы прямоугольников для нормализованной частоты. Из графика видно, что область совпадения данной кривой с характеристикой идеального интегратора (кривая 1) охватывает диапазон нормализованной частоты от 0 до примерно 0,1. Это означает, что удовлетворительные результаты интегрирования могут быть получены только для сигналов, частотный состав которых нс выходит за пределы 10 % от частоты дискретизации сигнала.

Более качественный результат может быть получен с помощью формулы трапеций, определяемой выражением:

y(n) = y(n – 1) + Т|х(n -1) + х(n – I )|.

Принцип вычисления интеграла по этой формуле и получаемый результат иллюстрируются на графиках «б» и «д» рис. 26. Как хорошо видно из графиков, погрешность вычисления в этом случае значительно ниже, чем для формулы прямоугольников. Найдем передаточную функцию формулы трапеций:

H(z)=\frac{T}{2} \frac{1+z^{-1}}{1-z^{-1}}

H(\omega)=\frac{T}{2} \frac{1+e^{-j \omega T}}{1-e^{-j \omega T}}=

=\frac{T}{2} \frac{e^{j \omega T / 2}+e^{-j \omega T / 2}}{e^{j \omega T / 2}-e^{-j \omega T / 2}}=\frac{T}{2} \frac{2 \cos (\omega T / 2)}{2 j \sin \left(\omega T / 2\right)}=\frac{T}{2 j} \operatorname{ctg}(\omega T / 2)

Соответственно АЧХ и ФЧХ фильтра описываются выражениями:

|Н(w)| = T/2 ctg (wT/2) ;

∠H(w) = –π/2 

Кривая 3 на рис. 25 показывает вид АЧХ формулы трапеций для нормализованной частоты.

Одной из наиболее популярных формул численного интегрирования является формула Симпсона:

y(n) = y(n – 2) + T/3 |х(n -1) + 4х(n -1) + х(n – 2)|.

Формула основана на интерполяции сигнала полиномом второй степени по трем последовательным точкам. Принцип вычисления интеграла по этой формуле и получаемый результат иллюстрируются на графиках «в» и «е» рис. 26. В данном случае видно, что достигается очень хорошее приближение к истинному результату. Однако результат получается с частотой отсчетов вдвое ниже исходной, так как для вычисления каждого следующего значения используются три подряд идущих отсчета. Найдем передаточную функцию формулы Симпсона:

H(z)=\frac{T}{3} \frac{1+4 z^{-1}+z^{-2}}{1-z^{-2}}
H(\omega)=\frac{T}{3} \frac{1+4 e^{-j \omega T}+e^{-2 j \omega T}}{1-e^{-2 j \mathcal{H} T}}=\frac{T}{2} \frac{e^{j \omega T}+4+e^{-j \omega T}}{e^{j \omega T}-e^{-j \omega T}}=
=\frac{T}{3} \frac{4+2 \cos (\omega T)}{2 j \sin (\omega T)}=\frac{T}{3} \frac{2+\cos (\omega T)}{j \sin (\omega T)}

Соответственно АЧХ и ФЧХ для формулы Симпсона описываются выражениями:

|H(w)|=\frac{3}{} \frac{2+\cos (\omega T)}{\sin (\omega T)}

∠H(w) = –π/2 

Кривая 4 на рис. 25 показывает вид АЧХ данного фильтра для нормализованной частоты. Как видно из графика, АЧХ формулы Симпсона демонстрирует очень хорошее совпадение с характеристикой идеального интегратора для значений нормализованной частоты от 0 до 0,2, но начиная с частоты примерно 0,3, она резко уходит вверх. На практике это означает, что формулу Симпсона нельзя использовать для расчета интеграла от сигналов, содержащих компоненты с частотой выше 25 % от частоты дискретизации.

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

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

Применение цифровой фильтрации в медицинских приборах и системах

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

  • устранение помех;
  • выделение информативных составляющих;
  • анализ сигналов по частотным диапазонам.

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

Рассмотрим пример использования ЦФ при анализе электроэнцефалограммы, т. е. записи электрических потенциалов мозга, снимаемых с поверхности головы человека. Выделяют четыре основных ритма ЭЭГ 17:

  • δ (дельта)     с частотой 0,5…3 Гц и амплитудой 40…300 мкВ;
  • θ (тета)       с частотой 4…6 Гц и амплитудой 40…300 мкВ;
  • α (альфа)     с частотой 8…13 Гц и амплитудой до 100 мкВ;
  • β (бета)       с частотой 14…40 Гц и амплитудой до 15 мкВ.

Использование полосовой цифровой фильтрации позволяет выделить каждый из ритмов отдельно, после чего по полученным сигналам вычисляются некоторые численные показатели, характеризующие дан ный ритм. На рис. 27 слева показаны исходные сигналы ЭЭГ по 19 каналам, а справа — эти же сигналы, пропущенные через полосовой ЦФ. рассчитанный на выделение полосы частот, соответствующей a-ритму ЭЭГ. На рис. 28 показана АЧХ данного фильтра. Можно видеть, что в результате фильтрации в сигнале отчетливо проявились участки, где присутствует о-ритм.

Слева — исходные сигналы ЭЭГ по 19 каналам, а справа — эти же сигналы, пропущенные через фильтр, выделяющий а-ритм
Рис. 27. Слева — исходные сигналы ЭЭГ по 19 каналам, а справа — эти же сигналы, пропущенные через фильтр, выделяющий а-ритм
АЧХ полосового ЦФ, рассчитанного на выделение полосы частот, соответствующей a-ритму ЭЭГ
Рис. 28. АЧХ полосового ЦФ, рассчитанного на выделение полосы частот, соответствующей a-ритму ЭЭГ

Footnotes

  1. Аифичер Э. С.. Джервис Б. У. Цифровая обработка сигналов: практический подход / Пер. с англ. — М.: ИД «Вильямс», 2004. — 992 с.
  2. Рабинер Л., Гоулд Б. Теория и применение цифровой обработки сигналов. — М.: Мир, 1978. — 848 с.
  3. Хемминг Р. В. Цифровые фильтры. — М.: Сов. радио, 1980. — 224 с.
  4. Микрокомпьютерные медицинские системы: Проектирование и применение / Под ред. У.Томпкинса, Дж. Уэбстера: Пер. с англ. — М.: Мир, 1983. – 544 с.
  5. Микрокомпьютерные медицинские системы: Проектирование и применение / Под ред. У.Томпкинса, Дж. Уэбстера: Пер. с англ. — М.: Мир, 1983. – 544 с.
  6. Аифичер Э. С.. Джервис Б. У. Цифровая обработка сигналов: практический подход / Пер. с англ. — М.: ИД «Вильямс», 2004. — 992 с.
  7. Рабинер Л., Гоулд Б. Теория и применение цифровой обработки сигналов. — М.: Мир, 1978. — 848 с.
  8. Хемминг Р. В. Цифровые фильтры. — М.: Сов. радио, 1980. — 224 с.
  9. Микрокомпьютерные медицинские системы: Проектирование и применение / Под ред. У.Томпкинса, Дж. Уэбстера: Пер. с англ. — М.: Мир, 1983. – 544 с.
  10. Аифичер Э. С.. Джервис Б. У. Цифровая обработка сигналов: практический подход / Пер. с англ. — М.: ИД «Вильямс», 2004. — 992 с.
  11. Хемминг Р. В. Цифровые фильтры. — М.: Сов. радио, 1980. — 224 с.
  12. Микрокомпьютерные медицинские системы: Проектирование и применение / Под ред. У.Томпкинса, Дж. Уэбстера: Пер. с англ. — М.: Мир, 1983. – 544 с.
  13. Рабинер Л., Гоулд Б. Теория и применение цифровой обработки сигналов. — М.: Мир, 1978. — 848 с.
  14. Хемминг Р. В. Цифровые фильтры. — М.: Сов. радио, 1980. — 224 с.
  15. Кардиомониторы. Аппаратура непрерывного контроля ЭКГ / Л. Л. Барановский, А. Н. Калиниченко, Л. А. Манило и др.; под ред. А. Л. Барановского и А. П. Немирно. — М.: Радио и связь, 1993. — 248 с.
  16. Микрокомпьютерные медицинские системы: Проектирование и применение / Под ред. У.Томпкинса, Дж. Уэбстера: Пер. с англ. — М.: Мир, 1983. – 544 с.
  17. Дюк В., Эмануэль В. Информационные технологии в медико-биологических исследованиях. СПб. «Питер», 2003. — 528 стр.