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

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

H(z) = H₀+H₁z+H₂z²+...= B(z)/[1+A(z)], (9.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). (9.1.2)

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

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

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

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

Как следует из выражения (9.1.3), при вычислении значения у_k текущей точки используются предыдущие вычисленные значения y_k-m, (m>0), что и определяет принцип рекурсии - фильтрации с обратной связью. Другой особенностью РЦФ является их односторонность и физическая реализуемость в реальном масштабе времени. При машинной обработке данных многочлен B(z) передаточной функции фильтра может реализоваться и в двухстороннем варианте.

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

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

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

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

H(z) = , (9.1.4)

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

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

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

Функции А_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².

В принципе, порядок расположения блоков в каскадной форме, равно как и порядок множителей B(z) и A(z) в числителе и знаменателе функции (9.1.4), значения не имеет. Однако следует учитывать, что полюса знаменателя, близкие к единичной окружности на z-плоскости (близкие по модулю к 1) формируют большие коэффициенты усиления на соответствующих частотах в блоках, в которых они находятся, и при обработке сигналов могут вызывать переполнение разрядов числовых ячеек этих блоков, если их разрядность ограничена. С учетом этого при формировании каскадов желательно объединять в пары B_i(z)/A_i(z) нули и полюса, близкие по модулю к 1, и располагать их в концевые блоки каскадной формы. Такое комбинирование полезно также с позиций наилучшего отношения сигнал/шум.

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

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

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

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

Стандартные блоки рекурсивных фильтров обычно реализуются биквадратными звеньями в канонической форме, которая имеет минимальное количество элементов задержки. Уравнения звена:

v(k) = x(k) –a(n) v(k-n), y(k) = b(n) v(k-n). (9.1.6)

Рис. 9.1.4. Каноническая форма

Функциональная схема реализации звена приведена на рис. 9.1.4.

Вторая форма реализации – по уравнению (9.1.5) в прямой форме, приведенная на рис. 9.1.5:

y(k) = b(n) x(k-n) –a(n) y(k-n). (9.1.7)

Рис. 9.1.5. Прямая форма

При определенных условиях прямая форма лучше канонической с точки зрения шумовых характеристик.

При нулевых значениях коэффициентов a₂ и b₂ звенья второго порядка превращаются в звенья первого порядка.

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

9.2. Разработка Рекурсивных цифровых фильтров [43].

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

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

Этапы разработки рекурсивных фильтров включают:

1. Задание частотной характеристики или передаточной функции фильтра.

2. Аппроксимация и расчет коэффициентов b(n) и a(m) передаточной функции фильтра (9.1.3). Этот этап может выполняться четырьмя методами:

· Метод размещения нулей и полюсов на комплексной z-плоскости.

· Метод инвариантного преобразования импульсной характеристики.

· Согласованное z-преобразование.

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

3. Выбор структуры реализации фильтра – параллельная или каскадная, блоками второго и/или первого порядка.

4. Программное или аппаратное обеспечение реализации фильтра.

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

|H(w)| = U_i /V_j, (9.2.1)

Каждой точке z_s = exp(-jw_sDt) может быть поставлен в соответствие вектор (z_s – n_i) на n_i -нуль, модуль которого U_i = |(z_s – n_i)| отображает расстояние от z_s до i-нуля, а равно и вектор (z_s – p_j) на p_j-полюс с соответствующим расстоянием V_j = (z_s – p_j). Наибольшее влияние на изменение АЧХ по частоте оказывают нули и полюсы, расположенные ближе к единичной окружности. При расположении нуля непосредственно на окружности гармоника w_s в этой точке полностью обнуляется (коэффициент передачи фильтра равен нулю). И, наоборот, при перемещении w_s к полюсу, близкому к единичной окружности, происходит резкое нарастание коэффициента усиления системы.

В качестве иллюстрации метода выполним расчет фильтра со следующими параметрами:

1. Полная режекция сигнала на частотах 0 и 250 Гц.

2. Полоса пропускания с центром на f_p = 125 Гц с шириной полосы по уровню 3 дб Dp =10 Гц.

3. Частота дискретизации данных f_D = 500 Гц.

При частоте дискретизации 500 Гц интервал временной дискретизации Dt = 1/f_D, а частота Найквиста f_N = 1/2Dt = 250 Гц. Соответственно, нули передаточной функции располагаются в точках n₁ = exp(-j2p 0 Dt) = 1 и n₂ = exp(-j2p 250 Dt) = -1. Угол из начала координат z-плоскости на полюс p с учетом его сопряженности для получения действительных коэффициентов ±180^{o .} f_p/f_N = ±90^o. Значение радиуса r до полюса определяет ширину полосы пропускания и в первом приближении (при 0 < r < 1.1) оценивается по выражению:

r » 1 + (Dp/f_D)p. r » 1.063.

Передаточная функция:

H(z) = (z-1)(z+1) / [(z-r exp(jp/2) (z-r exp(-jp/2)] = (z² -1)/(z² +r²) = Y(z)/X(z).

z²Y(z)+r²Y(z) = z²X(z)-X(z). Y(z) = [z²X(z)-X(z)-z²Y(z)] / r² .

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

y(k) = [x(k-2) – x(k) – y(k-2)] / r².

При использовании символики z^-1 полюс располагается внутри единичной окружности на том же радиусе со значением (при r < 1):

r » 1 - (Dp/f_D)p. r » 0.937.

Передаточная функция и алгоритм фильтра:

H(z) = (z-1)(z+1) / [(z-r exp(jp/2) (z-r exp(-jp/2)] = (z² -1)/(z² +r²) = (1-z^-2) / (1+r² z^-2).

y(k) = x(k) – x(k-2) – r² y(k-2).

Характеристики фильтров приведены на рис. 9.2.1. Индексы h1, H1, y1 относятся к первому фильтру с полюсом за пределами единичной окружности, индексами h2, H2, y2 – внутри окружности (символика z^-1). Импульсные отклики фильтров получены подачей на их входы импульса Кронекера, частотные характеристики вычислены по импульсным откликам. Значение r первого фильтра подобрано по АЧХ под равный коэффициент усиления гармоники f_p со вторым фильтром, после чего коэффициенты фильтров нормированы по коэффициенту усиления к 1 на частоте f_p.

Рис. 9.2.1.

Как следует из рисунков, изменение многочлена по степеням z на 1/z хотя и изменило коэффициенты разностного уравнения фильтра, но практически не повлияло на его амплитудно-частотную характеристику. Однако при этом произошло изменение области сходимости фильтра с соответствующим изменением фазовых углов направления на полюсы из всех точек единичной окружности, что отразилось на фазово-частотной характеристике и отсчетах импульсного отклика фильтра изменением фазы на p.

Синтез систем непосредственно в z-области применяется, в основном, только для режекторных и селекторных фильтров и более детально рассматривается ниже.

Метод инвариантного преобразования импульсной характеристики применяется для получения из подходящей аналоговой передаточной функции H(s) с помощью преобразования Лапласа импульсной характеристики h(t), которая затем дискретизируется и подвергается z-преобразованию,

Допустим, имеется простая аналоговая система с передаточной функцией:

H(s) = C / (s-p).

Выполняем обратное преобразование Лапласа функции H(s) и дискретизируем результат преобразования с определенной постоянной времени Dt:

h(t) = TL^-1[H(s)] = C exp(pt) → C exp(pnDt).

Выполняем z-преобразование и формируем передаточную функцию H(z):

H(z) =h(nDt) zⁿ =C exp(pnDt) zⁿ = C / (1-z exp(pDt)).

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

Согласованное z-преобразование применяется для преобразования аналоговых фильтров в эквивалентные цифровые непосредственным переводом всех полюсов и нулей с s-плоскости в z-плоскость:

(s-a) → z exp(aDt).

Большинство полюсов и нулей являются комплексно сопряженными и реализуются фильтрами второго порядка:

(s-a)(s-a*) → 1 – 2z exp(Re(a)Dt) cos(Im(a)Dt) + z² exp(Re(a)Dt).

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

Билинейное z-преобразование является основным методом получения коэффициентов рекурсивных БИХ-фильтров и использует следующую замену:

s = g (1-z) / (1+z), g = 1 или 2/Dt,

при этом ось jw s-плоскости отображается в единичную окружность z-плоскости, правая половина s-плоскости – внутрь единичной окружности, а левая половина с полюсами устойчивых аналоговых фильтров – снаружи единичной окружности. Аналогичная замена при отрицательной символике z^-1 с соответствующей сменой отображения:

s = g (z-1) / (z+1).

Для фильтров верхних и нижних частот порядок фильтра H(z) равен порядку фильтра H(s). Для полосовых и заградительных фильтров порядок H(z) вдвое больше порядка H(s). Для сохранения частотных характеристик фильтра при нелинейном сжатии частотной шкалы аналогового фильтра (переход от ∞ к w_N) предварительно выполняется деформация частотной шкалы аналогового фильтра. Более подробно эти вопросы рассмотрены ниже.