Процесс проектирования рекурсивного частотного фильтра обычно заключается в задании необходимой передаточной характеристики фильтра в частотной области и ее аппроксимации с определенной точностью какой-либо непрерывной передаточной функцией, с последующим z-преобразованием для перехода в z-область. Первые две операции хорошо отработаны в теории аналоговой фильтрации сигналов, что позволяет использовать для проектирования цифровых фильтров большой справочный материал по аналоговым фильтрам. Последняя операция является специфичной для цифровых фильтров.
Для алгебраического преобразования непрерывной передаточной функции в многочлен по z используется билинейное преобразование, известное в теории комплексных переменных под названием дробно-линейного преобразования.
Низкочастотный фильтр Баттеруорта /12,24/.
Рис. 10.1.1. АЧХ фильтра Баттеруорта.
Передаточная функция. Гладкий вид амплитудно-частотной характеристики фильтра Баттеруорта (рис. 10.1.1) задают квадратом передаточной функции вида:
|H(W)|2 = H(W)H*(W) = 1/(1+W2N).
где W = w/wc - нормированная частота, wc - частота среза АЧХ фильтра, на которой |H(w)|2 = 1/2 (соответственно H(w) = 0.707, или 3 дб), N - порядок фильтра, определяющий крутизну среза АЧХ. Функция |H(W)|2 – представляет собой энергетический спектр сигнала (спектральную плотность мощности) и не имеет фазовой характеристики, т.е. является четной вещественной, образованной произведением двух комплексно сопряженных функций H(W) и H*(W), При W → 0 коэффициент передачи фильтра стремится к 1. Учитывая, что результаты вычислений будут относиться к цифровым фильтрам и при z-преобразовании с переходом в главный частотный диапазон произойдет искажение частот, до начала расчетов фактические значения задаваемых частотных характеристик (значения wc, wp и ws) следует перевести в значения деформированных частот по выражению:
wд = (2/Dt) tg(wDt/2), -p/Dt<w<p/Dt. (10.1.1)
Крутизна среза. Наклон частотной характеристики фильтра при переходе от области пропускания к области подавления можно характеризовать коэффициентом крутизны среза фильтра K в децибелах на октаву:
K = 20 log|H(w2)/H(w1)|, (10.1.2)
где w1 и w2 - частоты с интервалом в одну октаву, т.е. w2 = 2w1.
Длительность импульсной реакции фильтра в пределах ее значимой части также зависит от крутизны среза: чем больше крутизна, тем больше длительность импульсного отклика фильтра.
Порядок фильтра. Принимая w1=Wc, w2=Ws и подставляя в (10.1.2) значения H(W) с приведенными данными, получим приближенное выражение для определения порядка фильтра по заданному значению К:
N = K/6. (10.1.6')
Так, для гарантированного ослабления сигнала в полосе подавления в 100 раз (40 децибел) порядок фильтра N = 7. В среднем, при изменении N на единицу коэффициент подавления сигнала изменяется на 6 децибел.
Исходные требования к передаточной функции фильтра обычно задаются в виде значений wp, ws и коэффициентов неравномерности (пульсаций) Ap и As (см. рис. 10.1.1). Для определения частоты среза wc по уровню 0.707 и порядка фильтра введем параметр d, связанный с коэффициентом Ар следующим соотношением:
(1-Ар)2 = 1/(1+d2).
d = [1/(1-Ар)] = Ap /(1-Ap). (10.1.3)
Для учета деформации частотной шкалы в процессе билинейного преобразования при переходе в дальнейшем к полиномам по Z, выполняем расчет деформированных частот wdp и wds по формулам:
wdp= 2 tg(wpDt/2)/Dt, (10.1.4)
wds= 2 tg(wsDt/2)/Dt.
При нормированной частоте W = w/wdc, где wdc соответственно также деформированная частота, на границах переходной зоны выполняются равенства:
1/(1+d2) = 1/[1+(wdp/wdc)2N], (10.1.5)
As2 = 1/[1+(wds/wdc)2N].
Отсюда:
d2 = (wdp/wdc)2N, 1/As2 - 1= (wds/wdc)2N.
Решая эти два уравнения совместно, находим:
N = ln [d/ ] / ln(wdp/wds), (10.1.6)
wdc = wdp/d1/N. (10.1.7)
Пример расчета фильтра низких частот Баттеруорта.
Рис. 10.1.2.
Начиная с этого параграфа, будем сопровождать рассмотрение теории последовательным расчетом фильтра низких частот с применением приводимых формул. Для расчета примем следующие исходные параметры фильтра:
- Шаг дискретизации данных Dt = 0.0005 сек. Частота Найквиста fN = 1/2Dt = 1000 Гц, ωN = 6.283·103 рад.
3. Порядок фильтра по формуле (10.1.6): N = 4.483.
Для пояснения порядка расчетов при четном и нечетном порядке фильтра, принимаем N1=4, N2=5.
4. Частота среза по формуле (10.1.7): wdc(N1) = 2.443·103 рад (389 Гц), wdc(N2) = 2.356·103 рад (375 Гц).
5. По формуле H(w)= [1/(1+w2N)]1/2, w = ω/ωdc, строим графики передаточных функций (рис.10.1.2).
Преобразование Лапласа.Переводим функцию |H(W)|2 на координатную ось пространства преобразования Лапласа при p = jW, для чего достаточно подставить W = p/j:
|H(р)|2 = 1/[1+(p/j)2N]. (10.1.8)
Полюсы функции находятся в точках нулевых значений знаменателя:
1+(p/j)2N = 0, p = j . (10.1.9)
Отсюда следует, что полюсы располагаются на единичной окружности в p-плоскости, а их местоположение определяется корнями уравнения (10.1.9). В полярных координатах:
pn = j exp(jp(2n-1)/2N), n = 1,2, ... ,2N. (10.1.10)
6. Вычисляем значения полюсов фильтра по формуле (10.1.10). Значения полюсов и их расположение на р-плоскости приведены на рис. 10.1.2. Положение первого полюса отмечено. Нумерация полюсов идет против часовой стрелки.
Как следует из формулы (10.1.10) и наглядно видно на рис. 10.1.2, все полюса с n ³ N являются комплексно сопряженными с полюсами n<N. Устойчивую минимально-фазовую передаточную функцию фильтра образуют полюса левой половины р-плоскости:
H(p) = G/B(p), (10.1.11)
где G - масштабный множитель, B(p) - полином Баттеруорта:
B(p) = B1(p) B2(p) ... BN(p), (10.1.12)
Bn(p) = p-pn. (10.1.13)
Практическая реализация фильтра Баттеруорта при четном значении N производится в виде последовательной каскадной схемы биквадратными блоками, т.е. составными фильтрами второго порядка. Для этого множители B(p) в (10.1.12) объединяются попарно с обоих концов ряда по n (от 1 до N) по комплексно сопряженным полюсам, при этом для каждой пары получаем вещественные квадратичные множители:
= p2+2p sin(p(2m-1)/2N)+1, n = 1,2, ..., N/2; m = n. (10.1.14)
Общее количество секций фильтра M=N/2. При нечетном N к членам (10.1.14) добавляется один линейный множитель с вещественным полюсом p(N+1)/2 = -1, пример положения которого на р-плоскости можно видеть на рисунке 10.1.2 для N=5:
В(N+1)/2(p)= p+1. (10.1.15)
Машинное время фильтрации на один оператор фильтра первого или второго порядка практически не отличаются, поэтому использование операторов первого порядка можно не рекомендовать и при установлении порядка фильтра по формуле (10.1.6) округлять расчетное значение N в сторону большего четного числа, что создает определенный запас по крутизне среза частотной характеристики.
Таким образом, передаточная функция ФНЧ Баттеруорта в p-области при четном N:
H(p) = G 1/Bm(p) = G 1/(p2+amp+1), (10.1.16)
am = 2 sin(p(2m-1)/2N), m = 1,2, ... ,N/2. (10.1.17)
При нечетном N:
H(p) = (G/p+1) 1/(p2+amp+1), (10.1.16')
Продолжение примера.
7. Вычисляем значения коэффициентов am по формуле (10.1.17):
- N=4: a1 = 0.765, a2 = 1.848.
- N=5: a1 = 0.618, a2 = 1.618.
Билинейное преобразование. Для перевода передаточной функции фильтра в z-область производится билинейное преобразование, для чего в выражение (10.1.16) подставляется параметр р:
p = g(1-z)/(1+z). (10.1.18)
С учетом автоматического возврата к нормальной шкале частот в главном частотном диапазоне z-преобразования значение коэффициента g:
g = 2/(Dt·ωdc). (10.1.19)
После перехода в z-область и приведения уравнения передаточной функции в типовую форму, для четного N получаем передаточную функцию из М=N/2 биквадратных блоков:
H(z) = G Gm (1+z)2 /(1-bm z+cm z2). (10.1.20)
Gm = 1/(g2 + amg + 1). (10.1.21)
bm = 2·Gm (g2 - 1). (10.1.22)
cm = Gm (g2 - amg + 1). (10.1.23)
При нечетном N добавляется один линейный блок первого порядка, который можно считать нулевым блоком фильтра (m=0):
H(z) = G Gm (1+z)2 /(1-bm z+cm z2), (10.1.24)
при этом, естественно, в выражении (10.1.24) используются значения коэффициентов Gm, bm и cm, вычисленные по (10.1.21-10.1.23) для нечетного значения N.
Значение множителя G в общем случае находится нормировкой к 1 коэффициента передачи фильтра при w = 0. Для ФНЧ при использовании вышеприведенных формул значение G равно 1.
При z=exp(-jw) главный диапазон функций H(z) от -p до p. Для получения передаточной функции в шкале физических частот достаточно вместо z в выражения (10.1.20, 10.1.24) подставить значение z=exp(-jwDt), где Dt – физический интервал дискретизации данных, и проверить соответствие расчетной передаточной функции заданным условиям.
9. Подставляем вычисленные коэффициенты в выражения (10.1.20, 10.1.24) и вычисляем значения передаточных функций при z = exp(-jwDt). Графики полученных функций приведены на рис. 10.1.3. На рис. 10.1.4
Во временной области фильтрация выполняется последовательной сверткой входного сигнала с операторами ячеек фильтра:
yk = xk ③ {h0(i)} ③ h1(i) ③ … ③ hМ(i), i = 0,1,2.
Уравнение рекурсивной фильтрации для m-го оператора фильтра:
Уравнение рекурсивной фильтрации для дополнительного h0(i) линейного оператора фильтра при нечетном N:
y0 = (xk+xk-1)/(g+1) + yk-1·(g-1)/(g+1) (10.1.26)
Продолжение примера.
Рис. 10.1.4.
10. Каждый оператор фильтра имеет определенную передаточную функцию, что можно видеть на рис. 10.1.4. Порядок последовательной свертки сигнала с операторами фильтра значения не имеет, но с учетом разрядности ячеек памяти звено H1(f) целесообразно реализовать за H2(f).
11. Для оценки длительности импульсной реакции фильтра подаем на вход фильтра импульс Кронекера на отсчете k = 3, и начинаем фильтрацию со второго отсчета (что обеспечивает начальные условия фильтрации на точках k=0 и k=1). Коэффициент усиления дисперсии шумов (сумма квадратов значений импульсного отклика) равен 0.341 при N=5, и 0.278 при N=4.
Рис. 10.1.1. АЧХ фильтра Баттеруорта.
= Ap 
/(1-Ap). (10.1.3)
] / ln(wdp/wds), (10.1.6)
Рис. 10.1.2.
. (10.1.9)
Рис. 10.1.2.
1/Bm(p) = G 
1/(p2+amp+1), (10.1.16')
Рис. 10.1.3.
Рис. 10.1.4.