На главную     ↑ Выше

Общий подход к расчёту цифровых фильтров


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



Этой формуле соответствует большое количество звеньев реализации, как пассивных LC, так и активных на операционных усилителях. По ней реализуются фильтры Баттерворта, Чебышева, Бесселя и др. Каскадное соединение звеньев 2-го порядка позволяет получить практически любые заданные амплитудно-частотные характеристики. Имеются методы расчётов и таблицы, позволяющие получить коэффициенты c0, c1, c2, d0, d1, d2 по заданным характеристикам.

Обобщенная передаточная функция цифрового звена второго порядка имеет вид:


Этой функции соответствует простая схема ЦФ второго порядка.


где Т – элемент задержки на один период дискретизации.
Задержка на Т в формуле для передаточной характеристики соответствует умножению на


Цифровой фильтр реализуется программно следующим фрагментом кода

y = a0*x + a1*x1 + a2* x2 + b1*y1 + b2*y2
x2 = x1
x1 = x
y2 = y1
y1 = y
y_out = k1*y

начальные значения x1, x2, y1, y2 принимаются равными нулю.
k1 – нормирующий коэффициент, например для обеспечения единичного усиления на нулевой или резонансной частоте; может находится как на выходе, для компенсации усиления, так и на входе для исключения переполнения цифрового вычислителя.

Таким образом, реализация ЦФ сводится к арифметическим операциям с цифровыми данными. Входное значение х представляет собой оцифрованный отсчёт аналогового сигнала, полученный на выходе аналого-цифрового преобразователя. Выходное значение y_out может быть подано на цифро-аналоговый преобразователь.
В отличие от АФ характеристики ЦФ стабильны во времени и не требуют подстройки. Кроме этого они легко перестраиваются и могут быть адаптивными.
Соединяя каскадно цифровые звенья 2-го порядка, можно получить практически любые требуемые АЧХ.

Задача проектирования цифрового фильтра на основе прототипа состоит:
1) В расчете коэффициентов аналогового фильтра-прототипа c0, c1, c2, d0, d1, d2 по заданным характеристикам. Методы расчётов приводятся в обширной литературе.
2) В пересчёте этих коэффициентов в коэффициенты цифрового фильтра a0, a1, a2, b1, b2, k1

В качестве примера рассмотрим фильтр Чебышева 2-го порядка с неравномерностью в полосе пропускания 3 дб.
Этот фильтр с весьма крутым спадом АЧХ в полосе непропускания, за что приходится «платить» неравномерностью АЧХ в полосе пропускания. Для него из таблиц [Д. Джонсон и др. Справочник по активным фильтрам, 1983 г.,с. 69] возьмём два параметра: В = 0.64490, C = 0.7007948, по которым рассчитаем коэффициенты аналогового, а затем и цифрового фильтров.

Расчёт

Ниже приведён текст программы, который можно скопировать в текстовое окно калькулятора.
Расчётный блок:



; расчет цифрового фильтра Чебышева 2-го порядка
fd = 1000 ; частота дискретизации
fc_dig = 100 ;частота среза цифрового фильтр
fc_an = (fd/pi)*tg(pi*fc_dig/fd) ; частота среза аналогового фильтра прототипа
wc_an = 2*pi*fc_an
;исходные параметры для фильтра Чебышева 3 дБ
B = 0.64490
C = 0.7007948
;расчёт коэффициентов аналогового фильтра-прототипа
c0 = C*wc_an^2
c1 = 0
c2 = 0
d0 = C*wc_an^2
d1 = B*wc_an
d2 = 1
; расчёт коэффициентов цифрового фильтра
q = 2*fd
c1 = q*c1
c2 = q^2*c2
d1 = q*d1
d2 = q^2*d2
CC = c0 + c1 + c2
DD = d0 + d1 + d2
;собственно коэффициенты
a0 = 1
a1 = 2*(c0 - c2)/CC
a2 = (c0 - c1 + c2)/CC
b0 = 1
b1 = - 2*(d0 - d2)/DD
b2 = - (d0 - d1 + d2)/DD
k1 = CC/DD
val(a0,a1,a2,k1,b1,b2): a0 = 1, a1 = 2, a2 = 1, k1 = 0.05764, b1 = 1.44292, b2 = -0.67349
; основной цикл
f = 0
[
w = 2*pi*f/fd
ch = a0^2 + a1^2 + a2^2 + 2*(a1*a2 + a1*a0)*cos(w) + 2*a0*a2*cos(2*w)
zn = 1 + b1^2 + b2^2 + 2*(b1*b2 - b1)*cos(w) - 2*b2*cos(2*w)
k = k1*sqr(ch/zn)
k_dB = 20*lg(k)
val(f,k_dB)
f = f + 10
{f>400 exit}
]




Программу легко обобщить для фильтров любого типа. Отличие будет в расчёте коэффициентов прототипа c0, c1, c2, d0, d1, d2. Их расчёт специфичен для каждого типа звеньев. Необходимо учитывать, что частота дискретизации fd должна быть в разы (как минимум в два с лишним раза) больше верхней частоты fc полосы пропускания. фильтра. Чем она выше, тем точнее совпадение АЧХ с прототипом.

Копируем эту программу в текстовое окно калькулятора и нажимаем кнопку Вычислить. Получаем текст с расчитанными коэффицентами.
В том же текстовом окне получаем таблицу АЧХ ЦФ:


f	k_dB
0	0
10	0.08153
20	0.32634
30	0.73291
40	1.28936
50	1.9514
60	2.59549
70	2.95683
80	2.67065
90	1.56833
100	-0.1247
110	-2.05501
120	-3.99524
130	-5.8505
140	-7.5951
150	-9.2313
160	-10.77066
170	-12.2268
180	-13.61272
190	-14.94004
200	-16.21897
210	-17.45841
220	-18.6662
230	-19.84932
240	-21.0141
250	-22.16637
260	-23.31164
270	-24.45522
280	-25.60232
290	-26.75821
300	-27.92834
310	-29.11844
320	-30.33469
330	-31.58394
340	-32.87384
350	-34.21321
360	-35.61233
370	-37.0835
380	-38.64166
390	-40.3054
400	-42.09839




В графическом окне калькулятора будет показан график данной АЧХ.


С помощью калькулятора можно провести моделирования полученного фильтра во временной области, вставив следующий текст программы:


; моделирование цифрового фильтра 2-го порядка
;
;-------------- коэффициенты фильтра ----------------
a0 = 1
a1 = 2
a2 = 1
b0 = 1
b1 = 1.44292
b2 = - 0.67349
k1 = 0.05764
;
fd = 1000 ; частота дискретизации
f = 70 ; частота входного сигнала
w = 2*pi*f/fd
x1 = 0
x2 = 0
y1 = 0
y2 = 0
;-------------- основной цикл -----------------
t = 0 ; дискретное время (номер отсчёта)
[
x = sin(w*t) =; входной сигнал
y = a0*x + a1*x1 + a2*x2 + b1*y1 + b2*y2
x2 = x1
x1 = x
y2 = y1
y1 = y
y_out = k1*y ; выходной сигнал после нормировки
val(t,x,y_out)
t = t + 1
{t>100 exit}
]




В том же текстовом окне получим таблицу отсчётов входного и выходного сигналов:
В графическом окне будет выведена диаграмма входного x и выходного y_out сигналов:



Видно, что на частоте 70 Гц выходной сигнал имеет усиление 1.41 и заметно отстает по фазе от входного сигнала.
Виден также переходный процесс продолжительностью примерно в один период. При выбранной частоте дискретизации 1000 Гц одному отсчёту соответсвует 1/1000 = 1 мс. Там же в графическом окне можно посмотреть спекты входного x и выходного y_out сигналов.

Видны две линии на частоте 7 гармоники. Т.к. период наблюдения составляет 100 мс, то частота первой гармоники для разложения в спектр равна 1/100 мс = 10 Гц, 7-ой гармоники - соответственно 70 Гц. В идеале синусоидальный сигнал должен давать одну линию, но из-за переходного процесса в спектре возникают и другие гармоники.

Выполнение расчётов под заказ.