Тема 5: РЕКУРСИВНЫЕ ФИЛЬТРЫ

Рекурсия – свойство живой природы. И амеба, и человек принимают решения на основании текущей ситуации и прошлого опыта. Всевышний может подтвердить, что это было самое удачное решение при сотворении Мира.

Писецкий. Введение в специальность.

(Уральский геофизик)

Творца и отца Вселенной и найти то трудно, а найдя нельзя показать его толпе.

Платон. Тимей.

(Греческий философ)

Содержание:5.1. Принципы рекурсивной фильтрации. Конструкция РЦФ. Каскадная форма. Параллельная форма. Устранение сдвига фазы. 5.2. Режекторные и селекторные фильтры. Комплексная z-плоскость. Режекторный фильтр постоянной составляющей сигнала. Режекторный фильтр произвольной частоты. Селекторный фильтр. 5.3. Билинейное z-преобразование. Принцип преобразования. Деформация частотной шкалы. 5.4. Типы рекурсивных частотных фильтров. Аппроксимационная задача. Передаточная функция. Виды фильтров. Литература.

5.1. Принципы рекурсивной фильтрации.

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

Конструкция РЦФ отображается в z-образе передаточной функции фильтра в виде отношения двух многочленов:

H(z) = H₀+H₁z+H₂z²+...= B(z)/[1+A(z)], (5.1.1)

где: B(z) = B₀+B₁z+B₂z²+ ... +B_Nz^N, A(z) = A₁z+A₂z²+ ... +A_Mz^M.

Естественно, что переход на РЦФ имеет смысл только в том случае, если степень многочленов A(z) и B(z) во много раз меньше степени многочлена H(z) прямого z-преобразования импульсной реакции фильтра. При z-образе входных данных Х(z), на выходе РЦФ имеем:

Y(z) = H(z)Х(z) = X(z)B(z)/[1+A(z)],

Y(z)[1+A(z)] = Y(z)+Y(z)A(z) = X(z)B(z),

Y(z) = X(z)B(z)-Y(z)A(z). (5.1.2)

При обратном z-преобразовании выражения (5.1.2) получаем уравнение рекурсивной цифровой фильтрации:

y_k = b_n x_k-n – a_m y_k-m. (5.1.3)

Рис. 5.1.1. Схема РЦФ.

Рекурсивная фильтрация требует задания начальных условий как по x_k, так и по y_k при k<0. Схема рекурсивной фильтрации приведена на рис. 5.1.1.

Как следует из выражения (5.1.3), при вычислении значения у_k текущей точки используются предыдущие вычисленные значения у_k-m, (m>0), что и определяет принцип рекурсии - фильтрации с обратной связью. Такая фильтрация обеспечивает определенную "память" по значениям предыдущих отсчетов, которая, в пределе, может быть бесконечной. С учетом этого фактора рекурсивные фильтры получили название фильтров с бесконечной импульсной характеристикой (БИХ-фильтров), в отличие от нерекурсивных фильтров, всегда имеющих конечную импульсную характеристику (КИХ-фильтры).

Другой особенностью РЦФ является их односторонность и физическая реализуемость в реальном масштабе времени. При машинной обработке данных многочлен B(z) передаточной функции фильтра может реализоваться и в двухстороннем варианте.

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

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

Практическая реализация РЦФ осуществляется в двух вариантах.

Рис. 5.1.2. Каскадная форма.

Каскадная форма. Находятся корни многочленов А(z),B(z) и производится разложение H(z):

H(z) = , (5.1.4)

где G - масштабный множитель. Это позволяет применять каскадное построение фильтров, показанное на рис. 5.1.2, в котором:

H(z) = G H₁(z) H₂(z) ..... H_N(z),

H_n(z) = B_n(z)/A_n(z).

Функции А_n(z) и B_n(z) обычно представляются в виде биквадратных блоков (фильтров второго порядка):

B_n(z) = b_n.0 + b_n.1 z + b_n.2 z²,

A_n(z) = 1 + a_n.1 z + a_n.2 z².

Рис. 5.1.3. Параллельная форма.

Параллельная форма. Функция H(z) разлагается на элементарные дроби:

H(z) = H_o(z) B_n(z) / [1+A_n(z)],

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

Устранение сдвига фазы. Рекурсивные фильтры являются фазосдвигающими фильтрами. Если требуется обеспечить нулевой фазовый сдвиг, то операция фильтрации производится дважды, в прямом и обратном направлении числовой последовательности массива данных, при этом амплитудно-частотная характеристика (АЧХ) фильтрации будет равна |H(w)|² фильтра, что необходимо учитывать при конструировании фильтра.

5.2. Режекторные и селекторные фильтры.

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

Рис. 5.2.1. Синтез фильтров.

Комплексная z-плоскость. Простейший фильтр типа НЦФ имеет один нуль на единичной окружности в z-плоскости в точке с частотой, которую необходимо подавить. Так, например, если из входного сигнала требуется исключить постоянную составляющую (нулевая частота), то импульсная реакция фильтра НЦФ имеет вид:

H(z) = 1-z. (5.2.1)

Нуль функции (5.2.1) равен z_n1=1. Как можно видеть на рис. 5.2.1, коэффициент передачи сигнала H(w) на любой частоте w_i от 0 до w_N=p/Dt - частоты Найквиста, определяемый выражением (5.2.1), будет равен длине вектора V_n1, проведенного из нуля функции H(z) - точка n1 на действительной оси, до соответствующей частоты w_i - точки z(w_i) на единичной окружности. На частоте w_i = 0 длина этого вектора равна нулю. Амплитудно-частотная характеристика фильтра, приведенная на рисунке 5.2.2 для передаточной функции (5.2.1) пунктиром, далека от идеальной для фильтр-пробки.

Рис. 5.2.2. АЧХ фильтров.

Режекторный фильтр постоянной составляющей сигнала. Сконструируем простейший РЦФ, добавив к оператору (5.2.1) один полюс вне единичной окружности на малом расстоянии от нуля:

H_п(z) = G(1-z)/(1-az), z_p= 1/a. (5.2.2)

Допустим, что полюс помещен в точке z_p1= 1.01, при этом а=0,99. Масштабный коэффициент G получим нормировкой H(z) к 1 на частоте Найквиста. Для приведенных условий G=0.995. Отсюда, при Dt=1:

H_п(z) = 0,995(1-z)/(1-0.99z),

y_k = 0.995(x_k-x_k-1)+ 0.99y_k-1.

Отображение нуля n₁ и полюса р₁ на z-плоскости и АЧХ фильтра для исключения постоянной составляющей приведены на рис.5.2.2. Коэффициент передачи сигнала на произвольной частоте w_i равен отношению длин векторов V_n1(z) и V_p1(z) соответственно из нуля и полюса до точки z(w_i) на единичной окружности и близок к единице для всех частот, за исключением нулевой:

|H_п(z)| = G V_n1(z)/V_p1(z).

Рис. 5.2.0.

Фазочастотная характеристика фильтра приведена на рис. 5.2.0 и определяется разностью фазовых углов векторов V_n1(z) и V_p1(z):

j_п(w) = j_n1-j_p1.

Режекторный фильтр произвольной частоты. При проектировании на подавление любой другой частоты w_v нули и полюсы располагаются на соответствующем радиусе z-плоскости. Радиальный угол направления на нуль и полюс определяются выражением:

j_v = p·w_v/w_N. (5.22.3)

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

H_v(z) = G(z-z_n)(z-z_n*)/[(z-z_p)(z-z_p*)]. (5.2.4)

Нули фильтра располагаются на единичной окружности:

z_n = cos j_v + j sin j_v = Re z_n + j Im z_n. (5.2.5)

Полюсы - на полярном радиусе R:

z_p = R·cos j_v + j R·sin j_v = Re z_p + j Im z_p. (5.2.6)

Пример положения нулей (n2 и n2*) и полюсов (р2 и р2*) приведен на рис.5.2.1. Подставляя (5.2.5-5.2.6) в (5.2.4), получаем:

H_v(z) = , (5.2.7)

G = [1+(1+2Re z_p)/R²] / (2+2Re z_n). (5.2.8)

При приведении уравнения (5.2.7) в типовую форму:

H_v(z) = , (5.2.7')

b₀ = 1, b₁ = -2·Re z_n, b₂ = 1. (5.2.9)

a₁ = - (2·Re z_p)/R², a₂ = 1/R².

Соответственно, алгоритм вычислений:

y_k = G·(x_k+b₁·x_k-1+x_k-2) – a₁·y_k-1 – a₂·y_k-2. (5.2.10)

В качестве примера проведем расчет режекторного фильтра на сетевую частоту питания приборов f_s = 50 Гц, которая очень часто попадает в измеренные данные. При шаге дискретизации данных Dt = 0.001 сек радиальный угол на нули и полюса фильтра в z-плоскости:

f_N = 1/2Dt = 500 Гц, j = p·f_s/f_N = 0.1π.

Радиус полюса фильтра примем равным R = 1.01. Значения нуля и полюса:

z_n = cos j + j sin j = 0.951 + 0.309 j, z_p = R·cos j_v + j R·sin j_v = 0.961 + 0.312 j.

Рис. 5.2.3.

Значение масштабного множителя G по (5.2.8):

G = 0.99.

Значения коэффициентов передаточной функции:

b₁ = -2·Re z_n = -1.902,

a₁ = - (2·Re z_p)/R² = -1.883, a₂ = 1/R² = 0.98.

При подстановке коэффициентов в уравнение (5.2.7') и замене z = exp(-jω) может быть получена непосредственно частотная передаточная функция фильтра, которая приведена на рис. 5.2.3:

0.99[1-1.902·exp(-jω)+exp(-2jω)]

H(w) = ---------------------------------------------------------

1-1.883·exp(-jω)+0.98·exp(-2jω)

Алгоритм фильтра:

y_k = 0.99·(x_k - 1.902·x_k-1 + x_k-2) + 1.883·y_k-1 – 0.98·y_k-2.

На рис. 5.2.4 приведен модельный входной сигнал фильтра, состоящий из суммы двух равных по амплитуде гармоник с частотой 50 и 53 Гц, и сигнал на выходе фильтра (смещен вверх). Справа на рисунке приведены спектры входного и выходного сигналов. Спектр выходного сигнала зарегистрирован после интервала установления реакции фильтра, который хорошо заметен на начальной части графика выходного сигнала. После установления сигнал на выходе фильтра практически полностью освобожден от гармоники 50 Гц.

Рис. 5.2.4.

Рис. 5.2.5.

При R → 1 ширина полосы подавления фильтра становится все более узкой, но при этом увеличивается длительность импульсной реакции фильтра и, соответственно, время установления фильтра при изменении спектра входного сигнала. В первом приближении значимая часть импульсной реакции режекторных фильтров равна (4÷5)/(R-1). Пример импульсной реакции для фильтра, вычисленного выше, приведен на рис. 5.2.5. Отклик фильтра получен при подаче на вход РЦФ импульса Кронекера. Для наглядности реакции на графике не показан начальный пик отклика (отсчет на нулевой точке), амплитуда которого равна значению G.

Селекторный фильтр. Если в уравнении (5.2.4) опустить нули, то получим селекторный фильтр, выделяющий сигналы одной частоты ω_s – частоты селекции, с передаточной функцией:

H_s(z) = G/[(z-z_p)(z-z_p*)], (5.2.11)

H_s(z) = , (5.2.11')

Рис. 5.2.6.

Характер передаточной функции (5.2.11) можно представить непосредственно по z-плоскости (рис. 5.2.1). При расположении полюсов фильтра за пределами единичного круга (например, в точках р2 и р2*) значение коэффициента передачи фильтра на произвольной частоте ω на единичной окружности будет обратно пропорционально величине векторов из этих точек окружности на полюса фильтра. При изменении ω от нуля до ±π (движение по единичной окружности на z-плоскости по или против часовой стрелки) один из векторов (на полюс противоположной полуплоскости) изменяется в достаточно небольших пределах (не превышая значения 2), в то время как второй из векторов (на полюс в своей полуплоскости) будут сначала уменьшаться, достигает минимума при расположении ω на полярном радиусе полюса (на частоте селекции ω_s), а затем снова начинает увеличиваться. Соответственно, значение H_s(ω) максимально на частоте селекции ±ω_s и при R → 1 может быть очень высоким. Пример передаточной функции (при G1=1) приведен на рис. 5.2.6. При необходимости фильтр может быть пронормирован к 1 на частоте селекции определением значения G1 по условию H_s(ω) = 1 при ω = ω_s, т.е.:

G1 = 1+a₁ z(w_s)+a₂ z(w_s)².

Рис. 5.2.7.

Фильтр (5.2.11) в принципе не может иметь нулевого коэффициента передачи на других частотах главного диапазона. Если последнее является обязательным, то фильтр выполняется методом обращения режекторного фильтра H_v(z):

H_s(z) = 1-H_v(z).

H_s(z) = . (5.2.12)

с₀ = 1-G, c₁ = a₁-Gb₁, c₂ = a₂-G.

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

Рис. 5.2.8. Фильтрация сигнала селекторным РЦФ.

(!!!KP30- Исследование возможности дополнения интегрирующих фильтров Симпсона и прямоугольников режекторными фильтрами на частоту Найквиста).

(!!!KP12- Разработка программы расчета ФЧХ режекторных и селекторных РЦФ и сдвига фазы выходных сигналов).

(!!!KP13- Исследование характеристик и разработка программы проектирования полосовых РЦФ на базе последовательных селекторных РЦФ).

(!!!KP14- Исследование характеристик и разработка программы проектирования гладких полосовых РЦФ на базе многократного селекторного РЦФ).

5.3. Билинейное z-преобразование.

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

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

Принцип преобразования. При стандартном z-преобразовании передаточной функции используется замена переменной вида:

z = exp(-pDt), (5.3.1)

где Dt - шаг дискретизации данных, p – комплексная переменная, р = s+jw.

Уравнение (5.3.1) можно записать в виде ln z = -pDt и разложить ln z в ряд:

ln z = -2[(1-z)/(1+z)+(1-z)³/(3(1-z)³)+ ....], z > 0.

Первый член этого разложения и представляет собой билинейное z- преобразование:

p = (2/Dt)(1-z)/(1+z). (5.3.2)

По сути, оно представляет собой отображение точек комплексной p-плоскости в точки комплексной z-плоскости, и наоборот. В общем виде:

p = g(1-z)/(1+z), (5.3.3)

z = (g-p)/(g+p). (5.3.4)

Значение множителя g не меняет формы преобразования, в связи с чем обычно принимают g = 1. Подставим p = jw в (5.3.4) и выразим z в показательной форме:

z = r exp(jj(w)), r = |z| = 1.

j(w) = 2 arctg(w/g),_­

Рис. 5.3.1.

При изменении w от - ¥ до ¥ фазовый угол j(w) монотонно изменяется от -p до p (см. рис. 5.3.1), т.е. мнимая ось p-плоскости (p = jw, -¥ < w < ¥) отображается в единичную окружность z-плоскости. В частности:

w = 0, z = exp(j0) = 1,

w = ¥, z = exp( jp) = -1

Деформация частотной шкалы. Реальное отображение передаточных функций фильтров является непрерывным (в силу своей физической сущности) и для упрощения дальнейших расчетов обычно задается в аналитической форме в комплексной р-плоскости по частотному аргументу ω от -¥ до +¥. При билинейном z-преобразовании происходит нелинейное искажение шкалы частот: полный частотный диапазон от - ¥ до ¥ непрерывных функций в р-плоскости сжимается до главного частотного диапазона от -p/Dt до p/Dt дискретных функций в z-плоскости. При задании уравнений непрерывных передаточных функций в частотной области это должно сопровождаться соответствующей обратной деформацией частотной шкалы, которая будет скомпенсирована при билинейном z-преобразовании. Подставляя в (5.3.2) z = exp(-jwDt) и умножая числитель и знаменатель правой части полученного уравнения на exp(jwDt/2), получим:

p = (2/Dt)[exp(jwDt/2)-exp(-jwDt/2)] / [exp(jwDt/2)+exp(-jwDt/2)],

p = (2/Dt) th(jwDt/2). (5.3.5)

Обозначим шкалу частот в р-области через индекс w_д (деформированная) и, полагая p = jw_д , с учетом тождества th(x) = - jtg(jx), получаем:

w_д = (2/Dt) tg(wDt/2) = g tg(wDt/2), -p/Dt<w<p/Dt. (5.3.6)

Рис. 5.3.2. Деформация частоты.

Выражение (5.3.6) позволяет осуществлять переход от фактических частот w главного частотного диапазона, которым должен соответствовать оператор РЦФ, к деформированным частотам w_д комплексной p-плоскости, на которой можно задавать требуемую форму передаточной функции проектируемого фильтра, при этом аппроксимация передаточных функций, учитывая область существования w от -¥ до ¥, может производиться многочленами и рациональными функциями. Связь частот приведена на рис. 5.3.2 (в начальной части p пространства деформированных частот).

5.4. Типы рекурсивных частотных фильтров.

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

Синтез рекурсивных фильтров, как и НЦФ, выполняется на базе фильтров низких частот (ФНЧ). Другие типы фильтров (ФВЧ - высоких частот, ПФ - полосовые, РФ - режекторные) образуются на основе ФНЧ путем частотного преобразования.

Рис. 5.4.1. Частотная характеристика ФНЧ.

Аппроксимационная задачанизкочастотного фильтра. В качестве основных исходных данных для решения аппроксимационных задач принимаются граничные частоты w_p - полосы пропускания, и w_s – начала полосы подавления сигнала. Как правило, задаются также допуски А_p - на максимальное значение неравномерности в полосе пропускания, и А_s – на максимальное отклонение АЧХ от нуля в полосе подавления (уровень шума фильтра). Разность между граничными частотами w_p и w_s будет определять ширину переходной зоны. Типичный пример задания формы АЧХ приведен на рис. 5.4.1. В допустимой зоне передаточной функции условно показана возможная форма АЧХ, удовлетворяющая заданным условиям.

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

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

|H(W)|² = H(W)·H*(W) = 1/(1+A_n(W)), (5.4.1)

где А_n(W) - многочлен n-го порядка, W - нормированная частота (например, W = w/w_p). Вид многочлена А_n(W) выбирается таким образом, чтобы выполнялось условие: А_n(W) << 1 при 0<W<1, что обеспечивает |H(W)|² ® 1, и A_n(W) >> 1 при W>1, соответственно |H(W)|² ® 0. Крутизна переходных зон фильтра устанавливается величиной порядка фильтра (чем больше значение n, тем больше крутизна переходных зон).

По знаменателю правой части выражения (5.4.1) достаточно просто могут быть определены комплексные полюса передаточной функции в p-области преобразования Лапласа и соответствующим комбинированием и объединением комплексно-сопряженных полюсов получены передаточные функции в виде биквадратных блоков при четном порядке, и с одним линейным блоком при нечетном порядке:

H(p) = G В_n(p), n-четное, (5.4.2)

H(p) = В_n(р), n-нечетное, (5.4.3)

где В_n(р) выражается в форме:

В_n(p) = 1/[(p-p_n)(p-p_n*)] = 1/(p²-2 a_np+b_n). (5.4.4)

Рис. 5.4.2. АЧХ фильтра Баттеруорта.

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