Низкочастотный фильтр Баттеруорта /12,24/

ЦИФРОВАЯ ОБРАБОТКА СИГНАЛОВ

Тема 6:   РЕКУРСИВНЫЕ ЧАСТОТНЫЕ ФИЛЬТРЫ

Благословен Господь, кто содеял все нужное нетрудным, а все трудное ненужным.

Григорий Сковорода. Украинский философ, ХIII в.

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

Отец Дионисий, в миру В.Лебедев.

 Геофизик Уральской школы, XX в.

Содержание:

Введение.

6.1. Низкочастотный фильтр Баттеруорта. Передаточная функция. Крутизна среза. Порядок фильтра. Преобразование Лапласа. Билинейное преобразование.

6.2. Высокочастотный фильтр Баттеруорта. Синтез фильтров методом частотного преобразования.

6.3. Полосовой фильтр Баттеруорта. Расщепление спектра. Полосовой фильтр на s-плоскости. Передаточная функция.

6.4. Фильтры Чебышева. Фильтры первого рода. Фильтры второго рода.

6.5. Дополнительные сведения.

Литература.

Введение.

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

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

Низкочастотный фильтр Баттеруорта   /12,24/.

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

Передаточная функция. Гладкий вид амплитудно-частотной характеристики фильтра Баттеруорта  (рис. 6.1.1) задают квадратом передаточной функции вида:

|H(W)|² = H(W)H*(W) = 1/(1+W^2N).

где W = w/w_c - нормированная частота, w_c - частота среза АЧХ фильтра, на которой |H(w)|² = 1/2 (соответственно H(w) = 0.707), N - порядок фильтра, определяющий крутизну среза АЧХ. При W → 0 коэффициент передачи фильтра стремится к 1. Учитывая, что результаты вычислений будут относиться к цифровым фильтрам и при z-преобразовании с переходом в главный частотный диапазон произойдет искажение частот, до начала расчетов фактические значения задаваемых частотных характеристик (значения w_c, w_p и w_s) следует перевести в значения деформированных частот по выражению:

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

Крутизна среза. Наклон частотной характеристики фильтра при переходе от области пропускания к области подавления можно характеризовать коэффициентом крутизны среза фильтра K в децибелах на октаву:

                                                      K = 20 log|H(w₂)/H(w₁)|,                              (6.1.2)

где w₁ и w₂ - частоты с интервалом в одну октаву, т.е. w₂ = 2w₁.

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

Порядок фильтра. Принимая w₁=W_с, w₂=W_s и подставляя в (6.1.2) значения H(W) с приведенными данными, получим приближенное выражение для определения порядка фильтра по заданному значению К:

N = K/6.                                           (6.1.6')

Так, для гарантированного ослабления сигнала в полосе подавления в 100 раз (40 децибел) порядок фильтра N = 7. В среднем, при изменении N на единицу коэффициент подавления сигнала изменяется на 6 децибел.

Исходные требования к передаточной функции фильтра обычно задаются в виде значений w_p, w_s и коэффициентов неравномерности (пульсаций) A_p и A_s (см. рис. 6.1.1). Для определения частоты среза w_c по уровню 0.707 и порядка фильтра введем параметр d, связанный с коэффициентом А_р следующим соотношением:

(1-А_р)² = 1/(1+d²).

                                                      d = [1/(1-А_р)]· .                          (6.1.3)

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

 w_dp= 2·tg(w_p·Dt/2)/Dt,                         (6.1.4)

w_ds= 2·tg(w_s·Dt/2)/Dt.

При нормированной частоте W = w/w_dc, где w_dc соответственно также деформированная частота, на границах переходной зоны выполняются равенства:

                                                      1/(1+d²) = 1/[1+(w_dp/w_dc)^2N],                        (6.1.5)

 A_s² = 1/[1+(w_ds/w_dc)^2N].

       Отсюда:

d² = (w_dp/w_dc)^2N,   1/A_s² - 1 = (w_ds/w_dc)^2N.

Решая эти два уравнения совместно, находим:

                                 N = ln [d/ ] / ln(w_dp/w_ds),                  (6.1.6)

                                                                  w_dc = w_dp/d^1/N.                                   (6.1.7)

Пример расчета фильтра низких частот Баттеруорта.

Рис. 6.1.2.

Начиная с этого параграфа, будем сопровождать рассмотрение теории последовательным расчетом конкретного фильтра низких частот с применением приводимых формул. Для расчета примем следующие исходные параметры фильтра:

- Шаг дискретизации данных Dt = 0.0005 сек.

Частота Найквиста f_N = 1/2Dt = 1000 Гц,  ω_N = 6.283·10³ рад.

- Граничная частота полосы пропускания:  f_p = 300 Гц,   w_p = 1.885·10³ рад.

- Граничная частота полосы подавления: f_s = 500 Гц, w_s = 3.142·10³ рад.

- Коэффициенты неравномерности:  А_р = А_s = 0.1.

Расчет дополнительных параметров:

1. Значение d по формуле (6.1.3) или по ее эквиваленту

d = A_p /(1-A_p): d= 0.484.

2. Деформированные частоты по формуле (6.1.4):

w_dp = 2.038·10³ рад.

w_ds = 4·10³ рад.

3. Порядок фильтра по формуле (6.1.6): N = 4.483.  

Для пояснения дальнейшего порядка расчетов при четном и нечетном порядке фильтра, принимаем N1=4, N2=5.

4. Частота среза фильтра по формуле (6.1.7):

w_dc(N1) = 2.443·10³ рад (389 Гц),     

w_dc(N2) = 2.356·10³ рад (375 Гц).  

5. По формуле H(w) = , w = ω/ω_dc, для контроля строим графики передаточных функций (рис.6.1.2).

Функция |H(W)|² – представляет собой энергетический спектр сигнала (спектральную плотность мощности) и не имеет фазовой характеристики, т.е. является четной вещественной, образованной из двух комплексно сопряженных функций H(W) и H*(W), при этом порядок фильтра N определяет число полюсов функции H(W) и комплексно сопряженных с ними полюсов функции H*(W).

Преобразование Лапласа. Переводим функцию |H(W)|² на координатную ось пространства преобразования Лапласа при p = jW, для чего достаточно подставить W = p/j:

                                                      |H(р)|² = 1/[1+(p/j)^2N].                                    (6.1.8)

Полюсы функции находятся в точках нулевых значений знаменателя:

                                                      1+(p/j)^2N = 0, p = j .                        (6.1.9)

Отсюда следует, что полюсы располагаются на единичной окружности в p-плоскости, а их местоположение определяется корнями уравнения (6.1.9). В полярных координатах:

p_n = j exp(jp(2n-1)/2N),    n = 1,2,...,2N.           (6.1.10)

p_n = j cos[p(2n-1)/2N] - sin[p(2k-1)/2N].           (6.1.10')

Рис. 6.1.2.

Продолжение примера расчета фильтра.

6. Вычисляем значения полюсов фильтра по формуле (6.1.10). Значения полюсов и их расположение на р-плоскости приведены на рис. 6.1.2. Положение первого полюса отмечено. Нумерация полюсов идет против часовой стрелки.

Как следует из формулы (6.1.10) и наглядно видно на рис. 6.1.2, все полюса с n ³ N являются комплексно сопряженными с полюсами n<N.  Устойчивую минимально-фазовую передаточную функцию фильтра образуют полюса левой половины р-плоскости:

H(p) = G/B(p),     (6.1.11)

где G - масштабный множитель, B(p) - полином Баттеруорта:

B(p) = B₁(p) B₂(p)... B_N(p),  (6.1.12)

B_n(p) = p-p_n.               (6.1.13)

Практическая реализация фильтра Баттеруорта при четном значении N производится в виде последовательной каскадной схемы биквадратными блоками, т.е. составными фильтрами второго порядка. Для этого множители B(p) в (6.1.12) объединяются попарно с обоих концов ряда по n (от 1 до N) по комплексно сопряженным полюсам, при этом для каждой пары получаем вещественные квадратичные множители:

В_m(p) = B_n(p)·B_N+1-n(p) =

= [p+j exp(jp(2n-1)/2N)][p+j exp(jp(2(N+1)-2n-1)/2N)] =

= [p+j exp(jp(2n-1)/2N)][p-j exp(jp(2n-1)/2N)] =

                              = p²+2p sin(p(2m-1)/2N)+1, n = 1,2,..., N/2;  m = n. (6.1.14)

Общее количество секций фильтра M=N/2. При нечетном N к членам (6.1.14) добавляется один линейный множитель с вещественным полюсом p_(N+1)/2 = -1, пример положения которого на р-плоскости можно видеть на рисунке 6.1.2 для N=5:

                                                      В_(N+1)/2(p)= p+1.                                         (6.1.15)

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

Таким образом, передаточная функция ФНЧ Баттеруорта в p-области при четном N:

                              H(p) = G 1/B_m(p) = G 1/(p²+a_mp+1),                  (6.1.16)

                       a_m = 2 sin(p(2m-1)/2N), m = 1,2,...,N/2.                  (6.1.17)

       При нечетном N:

                                H(p) = (G/p+1) 1/(p²+a_mp+1),                           (6.1.16')

Продолжение примера расчета фильтра.

7. Вычисляем значения коэффициентов a_m  по формуле (6.1.17):

- N=4: a₁ = 0.765, a₂ = 1.848.

- N=5: a₁ = 0.618, a₂ = 1.618.  

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

p = g·(1-z)/(1+z).                                         (6.1.18)

       С учетом автоматического возврата к нормальной (недеформированной) шкале частот в главном частотном диапазоне z-преобразования значение коэффициента g:

     g = 2/(Dt·ω_dc).                                    (6.1.19)

После перехода в z-область и приведения уравнения передаточной функции в типовую форму, для четного N получаем передаточную функцию из М=N/2 биквадратных блоков:

                                          H(z) = G G_m (1+z)² /(1-b_m z+c_m z²).                    (6.1.20)

                                          G_m = 1/(g² + a_mg + 1).                                               (6.1.21)

                                          b_m = 2·G_m (g² - 1).                                           (6.1.22)

                                                   c_m = G_m (g² - a_mg + 1).                       (6.1.23)

       При любом нечетном N добавляется один постоянный линейный блок первого порядка, который можно считать нулевым блоком фильтра (m=0):

                              H(z) = G G_m (1+z)² /(1-b_m z+c_m z²),   (6.1.24)

при этом, естественно, в выражении (6.1.24) используются значения коэффициентов G_m, b_m и c_m, вычисленные по (6.1.21-6.1.23) для данного нечетного значения N.

       При z=exp(-jw) главный диапазон функций H(z) от -p до p. Для получения передаточной функции в шкале физических частот достаточно вместо z в выражения (6.1.20, 6.1.24) подставить значение z=exp(-jwDt), где Dt – физический интервал дискретизации данных, и проверить соответствие расчетной передаточной функции заданным условиям.

Рис. 6.1.3. Рис. 6.1.4.

Продолжение примера расчета фильтра.

8. Вычисляем значения коэффициентов G_m, b_m и c_m:

- N=4: g = 1.637, G₁ = 0.203, G₂ = 0.149, b₁ = 0.681, b₂ = 0.501, c₁ = 0.492, c₂ = 0.098.

- N=5: g = 1.698, G₁ = 0.203, G₂ = 0.151, b₁ = 0.763, b₂ = 0.568, c₁ = 0.574, c₂ = 0.171.

9. Подставляем вычисленные коэффициенты в выражения (6.1.20, 6.1.24) и вычисляем значения передаточных функций при z = exp(-jwDt). Графики полученных функций приведены на рис. 6.1.3. На рис. 6.1.4 приведена фазочастотная характеристика фильтра (сплошная кривая) и групповое время задержки (пунктир) при N=4. Нелинейность ГВЗ в полосе пропускания, в принципе, не так велика, но начинает увеличиваться при увеличении порядка фильтра.

 

Во временной области фильтрация выполняется последовательной сверткой входного сигнала с операторами ячеек фильтра:

y_k = x_{k *} {h₀(i)} _* h₁(i) _* … _* h_М(i), i = 0,1,2.

Уравнение рекурсивной фильтрации для m-го оператора фильтра:

                                 y_k = G_m (x_k+2x_k-1+x_k-2) + b_m y_k-1 - c_m y_k-2.            (6.1.25)

Уравнение рекурсивной фильтрации для дополнительного h₀(i) линейного оператора фильтра при нечетном N:

                                 y₀ = (x_k+x_k-1)/(g+1) + y_k-1·(g-1)/(g+1)                     (6.1.26)

Продолжение примера расчета фильтра.

Рис. 6.1.5.

10. Каждый оператор фильтра имеет определенную передаточную функцию, что можно видеть на рис. 6.1.5. Порядок последовательной свертки сигнала с операторами фильтра значения не имеет.

11. Для оценки длительности импульсной реакции фильтра подаем на вход фильтра импульс Кронекера на отсчете k = 3, и начинаем фильтрацию со второго отсчета (что обеспечивает начальные условия фильтрации на точках k=0 и k=1). Сигналы на выходе первой и второй секции фильтра приведены на рис. 6.1.6.

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

12. Коэффициент усиления дисперсии шумов (сумма квадратов значений импульсного отклика) равен 0.341 при N=5, и 0.278 при N=4.

Рис. 6.1.6.

Значение множителя G в общем случае находится нормировкой к 1 коэффициента передачи фильтра при w = 0. Для ФНЧ и ФВЧ при использовании вышеприведенных формул значение G равно 1.