Тема 10. Рекурсивные частотные цифровые фильтры благословен Госп

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

Digital signals processing. Digital recursive frequency filters.

Тема 10. РЕКУРСИВНЫЕ ЧАСТОТНЫЕ ЦИФРОВЫЕ ФИЛЬТРЫ

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

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

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

Отец Дионисий, в миру В.Лебедев. Геофизик Уральской школы, XX в.

Содержание

Введение.

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

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

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

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

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

ВВЕДЕНИЕ

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

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

10.1. НИЗКОЧАСТОТНЫЙ ФИЛЬТР БАТТЕРУОРТА /12,24/.

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

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

|H(W)|2 = H(W)H*(W) = 1/(1+W2N).

где W = ω/ωc — нормированная частота, ωc — частота среза АЧХ фильтра, на которой |H(ω)|2 = 1/2 (соответственно H(ω) = 0.707, или 3 дб), N — порядок фильтра, определяющий крутизну среза АЧХ. Функция |H(W)|2 – представляет собой энергетический спектр сигнала (спектральную плотность мощности) и не имеет фазовой характеристики, т.е. является четной вещественной, образованной произведением двух комплексно сопряженных функций H(W) и H*(W), При W → 0 коэффициент передачи фильтра стремится к 1. Учитывая, что результаты вычислений будут относиться к цифровым фильтрам и при z-преобразовании с переходом в главный частотный диапазон произойдет искажение частот, до начала расчетов фактические значения задаваемых частотных характеристик (значения ωc, ωp и ωs) следует перевести в значения деформированных частот по выражению:

ωд = (2/Δt) tg(ωΔt/2), -π/Δt<ω<π/Δt. (10.1.1)

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

K = 20 log|H(ω2)/H(ω1)|, (10.1.2)

где ω1 и ω2 — частоты с интервалом в одну октаву, т.е. ω2 = 2ω1.

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

Порядок фильтра. Принимая ω1=Wc, ω2=Ws и подставляя в (10.1.2) значения H(W) с приведенными данными, получим приближенное выражение для определения порядка фильтра по заданному значению К:

N = K/6. (10.1.6′)

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

Исходные требования к передаточной функции фильтра обычно задаются в виде значений ωp, ωs и коэффициентов неравномерности (пульсаций) Ap и As (см. рис. 10.1.1). Для определения частоты среза ωc по уровню 0.707 и порядка фильтра введем параметр δ, связанный с коэффициентом Ар следующим соотношением:

(1-Ар)2 = 1/(1+δ2).

δ = [1/(1-Ар)] EMBED Equation.3 = Ap EMBED Equation.3 /(1-Ap). (10.1.3)

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

ωdp= 2 tg(ωpΔt/2)/Δt, (10.1.4)

ωds= 2 tg(ωsΔt/2)/Δt.

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

1/(1+δ2) = 1/[1+(ωdpdc)2N], (10.1.5)

As2 = 1/[1+(ωdsdc)2N].

Отсюда:

δ2 = (ωdpdc)2N, 1/As2 — 1 = (ωdsdc)2N.

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

N = ln [δ/ EMBED Equation.3 ] / ln(ωdpds), (10.1.6)

ωdc = ωdp1/N. (10.1.7)

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

Рис. 10.1.2.

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

— Шаг дискретизации данных Δt = 0.0005 сек. Частота Найквиста fN = 1/2Δt = 1000 Гц, ωN = 6.283·103 рад.

— Граничная частота пропускания: fp = 300 Гц, ωp = 1.885·103 рад.

— Граничная частота подавления: fs = 500 Гц, ωs = 3.142·103 рад.

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

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

1. Значение δ по формуле (10.1.3): δ= 0.484.

2. Деформированные частоты по формуле (10.1.4): ωdp = 2.038·103 рад. ωds = 4·103 рад.

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

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

4. Частота среза по формуле (10.1.7): ωdc(N1) = 2.443·103 рад (389 Гц), ωdc(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 EMBED Equation.3 . (10.1.9)

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

pn = j exp(jπ(2n-1)/2N), n = 1,2, … ,2N. (10.1.10)

pn = j cos[π(2n-1)/2N] — sin[π(2k-1)/2N]. (10.1.10′)

Рис. 10.1.2.

Продолжение примера.

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

Как следует из формулы (10.1.10) и наглядно видно на рис. 10.1.2, все полюса с 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) по комплексно сопряженным полюсам, при этом для каждой пары получаем вещественные квадратичные множители:

Вm(p) = Bn(p)·BN+1-n(p) =

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

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

= p2+2p sin(π(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 EMBED Equation.3 1/Bm(p) = G EMBED Equation.3 1/(p2+amp+1), (10.1.16)

am = 2 sin(π(2m-1)/2N), m = 1,2, … ,N/2. (10.1.17)

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

H(p) = (G/p+1) EMBED Equation.3 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 = γ(1-z)/(1+z). (10.1.18)

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

γ = 2/(Δt·ωdc). (10.1.19)

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

H(z) = G EMBED Equation.3 Gm (1+z)2 /(1-bm z+cm z2). (10.1.20)

Gm = 1/(γ2 + amγ + 1). (10.1.21)

bm = 2·Gm2 — 1). (10.1.22)

cm = Gm2 — amγ + 1). (10.1.23)

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

H(z) = G EMBED Equation.3 EMBED Equation.3 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 коэффициента передачи фильтра при ω = 0. Для ФНЧ при использовании вышеприведенных формул значение G равно 1.

При z=exp(-jω) главный диапазон функций H(z) от -π до π. Для получения передаточной функции в шкале физических частот достаточно вместо z в выражения (10.1.20, 10.1.24) подставить значение z=exp(-jωΔt), где Δt – физический интервал дискретизации данных, и проверить соответствие расчетной передаточной функции заданным условиям.

Рис. 10.1.3.

Продолжение примера.

8. Вычисляем значения коэффициентов Gm, bm и cm:

— N=4: γ = 1.637, G1 = 0.203, G2 = 0.149, b1 = 0.681, b2 = 0.501, c1 = 0.492, c2 = 0.098.

— N=5: γ = 1.698, G1 = 0.203, G2 = 0.151, b1 = 0.763, b2 = 0.568, c1 = 0.574, c2 = 0.171.

9. Подставляем вычисленные коэффициенты в выражения (10.1.20, 10.1.24) и вычисляем значения передаточных функций при z = exp(-jωΔt). Графики полученных функций приведены на рис. 10.1.3. На рис. 10.1.4

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

yk = xk ③ {h0(i)} ③ h1(i) ③ … ③ hМ(i), i = 0,1,2.

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

yk = Gm (xk+2xk-1+xk-2) + bm yk-1 — cm yk-2. (10.1.25)

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

y0 = (xk+xk-1)/(γ+1) + yk-1·(γ-1)/(γ+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.2. ВЫСОКОЧАСТОТНЫЙ ФИЛЬТР БАТТЕРУОРТА /12/.

Синтез фильтров методом частотного преобразования. Высокочастотные и полосовые фильтры конструируются путем частотной трансформации передаточных функций фильтров низких частот. Если обозначить аргумент передаточных функций ФНЧ через p=jW, a функций ФВЧ и ПФ через s=jw, то всегда можно найти такую функцию частотного преобразования p=F(s), которая превращает один тип фильтров в другой. Для преобразования ФНЧ → ФВЧ функция частотного преобразования имеет вид:

p = 1/s, (10.2.1)

В этом нетрудно убедиться сравнением двух видов преобразования. Как известно, передаточная функция ФВЧ может быть получена из ФНЧ разностью между широкополосным фильтром (H(ω)=1) и ФНЧ. Применяя этот метод для функции Баттеруорта, получаем:

|H(w)|2 = 1-|H(W)|2 = 1- 1/(1+W2N) = W2N/(1+W2N). (10.2.2)

С другой стороны, при W = p/j: |H(p)|2 = 1/(1-p2N). Выполняя подстановку (10.2.1) в это выражение, получаем:

|H(s)|2 = s2N/(s2N-1).

Возвратимся из последнего выражения к аргументу w с учетом принятого равенства s=jw:

|H(s)|2 = (jw)2N/((jw)2N-1) =(w)2N/(1+(w)2N),

что полностью повторяет (10.2.2) при w=W.

Подставляя p=1/s непосредственно в выражение H(p) (10.1.16) для четного значения N, получаем:

H(s) = G EMBED Equation.3 s2/(s2+am s+1). (10.2.3)

Для нечетного N:

H(s) = [G·s/(s+1)] EMBED Equation.3 s2/(s2+am s+1). (10.2.4)

После билинейного z-преобразования выражения с подстановкой s=γ(1-z)/(1+z), для четного и нечетного значений N соответственно:

H(z) = G EMBED Equation.3 γ2·Gm·(1-z)2/(1-bm z+cm z2). (10.2.5)

H(z) = G EMBED Equation.3 EMBED Equation.3 γ2·Gm·(1-z)2/(1-bm z+cm z2). (10.2.6)

Gm = 1/(γ2 + amγ + 1). (10.2.7)

bm = 2·Gm2 — 1).

cm = Gm2 — amγ + 1).

Значения коэффициентов Gm, bm, cm остаются без изменения (сравнить с (10.1.21-10.1.23)). При задании частотных параметров ФВЧ в том же виде, что и для ФНЧ, формула расчетов N и ωdc получается аналогично ФНЧ, при этом в знаменателе выражения (10.1.6) отношение ωdpds заменяется на ωdsdp:

N = ln [δ/ EMBED Equation.3 ] / ln(ωdsdp), (10.2.8)

а в (10.1.7) деление членов правой части меняется на умножение:

ωdc = ωdp·δ1/N. (10.2.9)

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

yk = γ2·Gm (xk-2xk-1+xk-2) + bm yk-1 — cm yk-2. (10.2.10)

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

y0 = γ·(xk-xk-1)/(γ+1) + yk-1·(γ-1)/(γ+1). (10.2.11)

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

Техническое задание:

— Шаг дискретизации данных Δt = 0.0005 сек. Частота Найквиста fN = 1/2Δt = 1000 Гц, ωN = 6.283·103 рад.

— Граничная частота полосы пропускания: fp = 700 Гц, ωp = 4.398·103 рад.

— Граничная частота полосы подавления: fs = 500 Гц, ωs = 3.142·103 рад.

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

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

1. δ = Ap EMBED Equation.3 /(1-Ap): δ= 0.484.

2. Деформированные частоты по формуле (10.1.4): ωdp = 7.85·103 рад. ωds = 4·103 рад.

3. Порядок фильтра по формуле (10.2.8): N = 4.483. Для расчетов принимаем N=4.

Рис. 10.2.1.

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

ωdc = 6.549·103 рад (1042 Гц),

5. Строим график функции H(w), w = ω/ωdc, (рис.10.2.1).

6. Полюса pn фильтра полностью повторяют полюса ФНЧ (рис. 10.1.2), а, соответственно, повторяются и значения коэффициентов am. Остальные коэффициенты: γ = 0.611, G1 = 0.543, G2 = 0.4, b1 = — 0.681, b2 = — 0.501, c1 = 0.492, c2 = 0.098.

Рис. 10.2.2.

При сравнении коэффициентов bm, cm и коэффициентов в числителе передаточных функций ФВЧ с соответствующими коэффициентами ФНЧ предыдущего примера можно заметить, что в данном фильтре относительно ФНЧ произошла только смена знаков коэффициентов при нечетных степенях z. Это объясняется тем, что заданные в данном примере параметры ФВЧ по частоте соответствуют частотному реверсу ФНЧ: ω’ = π-ω, что приводит к частотному реверсу передаточной функции низкочастотного фильтра и превращению его в высокочастотный фильтр. Этот способ обращения ФНЧ также может использоваться для расчетов ФВЧ.

7. Импульсная реакция фильтра, вычисленная по (10.2.10) при подаче на вход фильтра импульса Кронекера приведена на рис. 10.2.2.

Курсовая работа 15-07. Разработка программы устранения сдвига фазы сигналов при использовании фильтров Баттеруорта.

10.3. ПОЛОСОВОЙ ФИЛЬТР БАТТЕРУОРТА /12/.

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

p = s+1/s. (10.3.1)

Подставив в (10.3.1) значения p = jW и s = jw, получим:

W = [w2-1]/w,

w2-Ww-1 = 0. (10.3.2)

Корни уравнения (10.3.2):

(w)1,2 = W/2 EMBED Equation.3 EMBED Equation.3 . (10.3.3)

Расщепление спектра. При W=0 имеем w = EMBED Equation.3 1, т.е. центр полосы пропускания ФНЧ (от -Wc до +Wc) расщепляется на два (как и положено, для полосовых фильтров) и смещается в точки w = EMBED Equation.3 1. Подставив в (10.3.3) граничную частоту Wс=1 нормированного ФНЧ, определяем граничные частоты нормированного полосового фильтра в виде пары сопряженных частот:

w1 = EMBED Equation.3 0.618, w2 = EMBED Equation.3 1.618

Рис. 10.3.1. Расщепление полосы.

Сущность произведенного преобразования наглядно видна на рис. 10.3.1. Ширина полосы пропускания нормированного ПФ равна 1.

Полученное преобразование можно распространить на полосовой фильтр с ненормированными частотами ωн и ωв.

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

ωо= EMBED Equation.3 . (10.3.4)

Ширина полосы пропускания ПФ связана (см. рис.10.3.1) с граничной частотой ФНЧ соотношением:

Δω = ωвн = ωс = ωн.

В долях средней геометрической частоты:

Wн = (ωвн)/ωо = Wc. (10.3.5)

Заменяя в (10.3.4-10.3.5) значение ωв на произвольную частоту ω и подставляя в (10.3.5) значение ωн = ω·ωо2 из (10.3.4), получаем произвольную частоту W:

W = (ω-ωн)/ωо = ω/ωoo/ω. (10.3.6)

Отсюда, в выражении (10.1.1) вместо нормированной частоты W = ω/ωс можно применить функцию частоты полосового фильтра w(ω):

w(ω) = (ω2о2)/[ω(ωвн)],

или, подставляя (10.3.4) вместо ωо:

w(ω) = (ω2нωв)/[ω(ωвн)]. (10.3.7)

Тем самым передаточная функция ФНЧ выражается в единицах, которые позволяют после применения преобразования (10.3.1) использовать для задания необходимые граничные частоты ωн и ωв полосового фильтра.

Пример расчета полосового фильтра Баттеруорта.

Техническое задание:

— Шаг дискретизации данных Δt = 0.0005 сек. Частота Найквиста fN = 1/2Δt = 1000 Гц, ωN = 6.283·103 рад.

— Нижняя граничная частота полосы пропускания: fн = 340 Гц, ωн = 2.136·103 рад.

— Верхняя граничная частота полосы пропускания: fв = 470 Гц, ωв = 2.953·103 рад.

— Крутизна срезов в децибелах на октаву: Кр = 45.

Расчет параметров:



Страницы: 1 | 2 | Весь текст