Моделирование цифровой обработки сигналов ЦОС в Matlab. Часть 2. Синтез оптимальных БИХ-фильтров программными средствами Matlab Текст научной статьи по специальности «Математика»

Алла СОЛОНИНА

Моделирование цифровой обработки сигналов в MATLAB.

Часть 2. Синтез оптимальных БИХ-фильтров программными средствами MATLAB

Основные этапы проектирования цифровых фильтров (КИХ и БИХ) были рассмотрены в [6]. В этой статье автор знакомит читателей с синтезом оптимальных БИХ-фильтров — первым этапом их проектирования.

Свойства БИХ-фильтров

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

N-1

— (1) 1+^акг~к к=\

и при (М-1)<(М-1) (по умолчанию) имеет порядок Я = М-1.

Сложность БИХ-фильтра определяется порядком Я передаточной функции (1).

БИХ-фильтры характерны следующими особенностями:

• нелинейной ФЧХ;

• необходимостью проверки на устойчивость.

Задание требований к частотным характеристикам БИХ-фильтров

При синтезе частотно-избирательных БИХ-фильтров с существенно нелинейной ФЧХ последняя обычно не контролируется, и требования задаются к АЧХ. Они не отличаются от требований к АЧХ КИХ-фильт-ров, за тем исключением, что для рассматриваемого далее метода синтеза значение АЧХ в полосе пропускания не превышает единицы. Кроме того, для БИХ-фильтров требования задаются к характеристике затухания — АЧХ (дБ) (см. формулу (6) в [6]) и включают в себя:

• частоту дискретизации /д (Гц);

• граничные частоты полос пропускания (ПП) и полос задерживания (ПЗ), для которых введены условные обозначения:

- / — граничная частота ПП для ФНЧ иФВЧ;

- /к—граничная частота ПЗ для ФНЧ и ФВЧ;

- /у, / — левая и правая граничные частоты ПП для ПФ и РФ;

- /-к, /к — левая и правая граничные частоты ПЗ для ПФ и РФ;

• допустимые отклонения от А(/) (дБ)

(см. формулу (6) в [6]):

- атах (дБ) — максимально допустимое затухание в ПП;

- яш;п (дБ) — минимально допустимое затухание в ПЗ.

В статье рассматривается синтез оптимальных БИХ-фильтров методом билинейного 2-преобразования на основе аналоговых фильтров-прототипов (АФП).

Идея синтеза БИХ-фильтров на основе АФП возникла из желания воспользоваться давно известными и хорошо себя зарекомендовавшими методами синтеза аналоговых фильтров. Обоснование такой возможности вытекает из следующих положений:

• передаточные функции АФП и БИХ-филь-тров — дробно-рациональные;

• импульсные характеристики АФП и БИХ-фильтров — бесконечные.

Для того чтобы подчеркнуть контраст типа фильтра (аналоговый или цифровой), будем использовать аббревиатуры АФП и ЦФ, по умолчанию подразумевая под ЦФ БИХ-фильтр.

Процедура синтеза БИХ-фильтра

Процедура синтеза ЦФ на основе АФП включает в себя [4]:

1. Задание требований к АЧХ ЦФ.

2. Выбор метода синтеза.

3. Формирование требований к АЧХ АФП.

4. Выбор типа аппроксимирующей функции. Четырем типам аппроксимирующих функций соответствуют четыре разновидности аналоговых (и цифровых) фильтров:

- Баттерворта (Butterwhorth) — с АЧХ, максимально плоской в ПП и монотонной в ПЗ;

- Чебышева I рода (Chebyshov Type I) — с АЧХ, равноволновой в ПП и монотонной в ПЗ;

- Чебышева II рода (Chebyshov Type II) — с АЧХ, максимально плоской в ПП и равноволновой в ПЗ;

- Золотарева - Кауэра (эллиптические фильтры) (Eleptic) — с АЧХ, равноволновой в ПП и ПЗ.

5. Расчет передаточной функции АФП.

6. Преобразование передаточной функции АФП в передаточную функцию ЦФ.

Для лучшего понимания синтеза в MATLAB ЦФ на основе АФП коротко познакомимся с синтезом АФП.

Синтез аналоговых фильтров

Синтез частотно-избирательных АФП Баттерворта, Чебышева I и II рода и Золотарева - Кауэра выполняется соответственно с помощью функций:

[bs,as]=butter(R,Wn,ftype,y)

[bs,as]=cheby1(R,rp,Wn,ftype,V)

[bs,as]=cheby2(R,rs,Wn,ftype,'s')

[bs,as]=ellip(R,rp,rs,Wn,ftype,'s')

Здесь R — порядок АФП; Wn — вектор частот среза в шкале ю = 2nf (рад/с), содержащий один элемент — для ФНЧ и ФВЧ и два — для ПФ и РФ (частотами среза называют частоты, на которых нормированная АЧХ АФП A(f) равна 1/V2«0,707, а A(f) (дБ) — 3 дБ); rp, rs — соответственно максимально и минимально допустимые затухания amax (дБ) в ПП и «min (дБ) в ПЗ для A(f) (дБ); ftype — параметр, указывающий тип избирательности и принимающий значения: 'high' — для ФВЧ; 'stop' — для РФ; по умолчанию (если

значение параметра не задано явно) — для ФНЧ и ПФ; У — признак аналогового фильтра, при его отсутствии по умолчанию подразумевается ЦФ; Ь8, а8 — векторы коэффициентов числителя и знаменателя передаточной функции АФП На(р) в порядке возрастания степеней р; ав(1) = 1.

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

Как правило, при синтезе АФП порядок фильтра (Я) и частоты среза (”Шп) заранее неизвестны. Их можно определить по требованиям к АЧХ с помощью следующих функций, соответственно для АФП Баттерворта, Чебышева I и II рода и Золотарева - Кауэра:

[К,'Шп]=ЬиИ:оЫ('Шр,М8,гр,Г8,'8')

[К,'Шп]=сЬеЬ1оЫ(Шр,М8,Гр,Г8,,8')

[К,'Шп]=сЬеЬ2оЫ(Шр,М8,Гр,Г8,,8')

[К,'Шп]=еШроЫ('Шр,М8,гр,т8,У)

Здесь И — минимальный порядок при заданных требованиях, соответствующий оптимальному АФП; Wp, Ws — соответственно векторы граничных частот ПП и ПЗ в порядке их следования слева направо в шкале частот ю = 2п/ (рад/с). Остальные параметры были определены ранее.

Пример 1

Синтезировать оптимальные АФП Баттер-ворта, Чебышева I и II рода и Золотарева -Кауэра по заданным требованиям к АЧХ ФНЧ (см. табл. 3 в [6]). Значения атах = 0,4455 дБ и ятт = 40 дБ (гр и ге) были вычислены в [6] в примере 1:

>> Й=1000; £к=1500;

>> Шр=2.*р1.*Й; Ш8=2.*р1.*£к;

>> гр=0.4455; к=40;

>> [R1,Wn1]=buttord(Wp,Ws,rp,rs,'s');

>> [R2,Wn2]=cheb1ord(Wp,Ws,rp,rs,'s');

>> [R3,Wn3]=cheb2ord(Wp,Ws,rp,rs,'s');

>> [R4,Wn4]=eШpoгd(Wp,Ws,гp,гs,У);

>>

>> [bs2,as2]=cheby1(R2,гp,Wn2,'s');

>> [bs3,as3]=cheby2(R3,гs,Wn3,'s');

>> [bs4,as4]=eШp(R4,гp,гs,Wn4,'s');

Выведем значения порядков R1, R2, R3, и R4 соответственно оптимальных ФНЧ Бат-терворта, Чебышева I и II рода и Золотарева - Кауэра:

>> R=[R1 R2 R3 R4]

R =

15 7 7 5

ЧГц)

f(rit)

ЧГц)

f(rit)

Рис. 1. АЧХ аналоговых ФНЧ: а) Баттерворта; б) Чебышева I рода; в) Чебышева II рода; г) Золотарева — Кауэра

Наименьший порядок имеет ФНЧ Золотарева - Кауэра.

Построим графики АЧХ аналоговых ФНЧ Баттерворта, Чебышева I и II рода и Золотарева - Кауэра на густой сетке частот (выберем 1000 точек) и выведем их в основной полосе частот [0; /д/2] ЦФ при частоте дискретизации 8000 Гц (для сравнения с ним впоследствии).

Для построения графиков АЧX АФП используем функцию:

Ha=freqs(bs,as,W)

где Ь8, а8 — коэффициенты числителя и знаменателя передаточной функции АФП; W — вектор, задающий сетку частот в шкале ю = 2 л/ (рад/с).

Выведем значения АЧХ всех АФП в одинаковом диапазоне [0;1] по оси ординат с помощью функции у1іт([0 1]) (рис. 1):

>> %f — густая сетка частот в Гц

>> %W — густая сетка круговых частот в рад/с

>> %oHa1,Ha2,Ha3,Ha4 — передаточные функции АФП Баттерворта,

... Чебышева I и II рода и Золотарева - Кауэра

>> Fs=8000;

>> f=0:((Fs/2)/1000):Fs/2;

>> W=2.*pi.*f;

>> Ha1=freqs(bs1,as1,W); MAG1=abs(Ha1);

>> Ha2=freqs(bs2,as2,W); MAG2=abs(Ha2);

>> Ha3=freqs(bs3,as3,W); MAG3=abs(Ha3);

>> Ha4=freqs(bs4,as4,W); MAG4=abs(Ha4);

>> subplot(2,2,1),plot(f,MAG1),xlabel('f(Hz)'),grid,... ylabel('MAGNITUDE'),title('Analog Filter Butterworth'),ylim([0 1]) >> subplot(2,2,2),plot(f,MAG2),xlabel('f(Hz)'),grid,... ylabel('MAGNITUDE'),title('Analog Filter Chebyshov I'),ylim([0 1]) >> subplot(2,2,3),plot(f,MAG3),xlabel('f(Hz)'),grid,... ylabel('MAGNITUDE'),title('Analog Filter Chebyshov II'),ylim([01]) >> subplot(2,2,4),plot(f,MAG4);xlabel('f(Hz)'),grid,... ylabel('MAGNITUDE'),title('Analog Filter Eleptic'),ylim([0 1])

Синтез БИХ-фильтров методом билинейного Z-преобразования

Отображение р-плоскости в г-плоскость выполняется в соответствии с формулой билинейного 2-преобразования:

lnz = 2

Z— 1 (Z—1

—+

Z+1

Z+1

и сохранения первого члена.

Формула (2) позволяет представить передаточную функцию ЦФ Н(г) на основе передаточной функции АФП На(р):

т=щр) (3)

р=у—

1+г~1

Используя (2), выражаем г через р:

г = (у+р)/(у-р).

И подставляя г = те’ю и р = ]П, где — обо-

значение оси частот АФП (во избежание путаницы), получаем:

Y+/Q y'2arctg -■ = Є

У-jn

(4)

Р = У

1 —z

1+Z"1 ’

(2)

где у = 1/Т, получаемой из стандартного 2-преобразования г = ерТ ^ р = (1/Т)1п2, путем разложения Ыг в ряд Тейлора:

Откуда имеем нелинейные зависимости между частотами АФП и ЦФ:

г. t тТ ♦ ®

^ = Ytg—= ?tg-; (5)

Q Т

m=yarctg— . (б)

Согласно (4), ось частот ] П р-плоскости, как и при стандартном 2-преобразовании, отображается на г-плоскости в единичную окружность (радиус равен единице), однако каждому ее обороту (изменению нормированной частоты на Дю = 2п), а именно: ..., -3п<(в<-п, -п<(в<п, л<Ю<3п, ..., соответствует не отрезок оси, как при стандартном Z-преобразовании, а вся ось ]П (так как зави-

симость между частотами определяется функцией агй£: ..., -ю<^<ю, -ю<^<ю, -TO<Q<<»...

Связь между частотными характеристиками АФП и ЦФ, соответственно И(ёюТ) и Иа(]ю), имеет вид [1, 4]:

' 2л)

a>+m—

V т,

. (7)

При этом частотная характеристика АФП Ha(jQ), бесконечная в шкале частот Q, в шкале частот ( , согласно (б), сжимается в отрезок Д( = ( д, то есть становится финитной. Соответственно, частотная характеристика ЦФ H(ej (T), согласно (7), оказывается равной (с точностью до множителя 1/T) бесконечной сумме копий финитных частотных характеристик АФП длины Дю = юд, сдвинутых друг относительно друга на частоту ( д. Вследствие этого элайсинг (Aliasing) принципиально отсутствует, и АЧX ЦФ и АФП в основной полосе частот [0; юд/2] совпадают (с учетом нелинейной зависимости между частотами).

«Платой» за отсутствие элайсинга, помимо нелинейной зависимости между частотами АФП и ЦФ, будет полное несовпадение их импульсных характеристик и ФЧX.

Метод билинейного 2-преобразования реализуется по-разному, в зависимости от поставленной задачи, а именно:

• Если ЦФ синтезируется на основе имеющегося АФП (копирует его АЧX с учетом нелинейного соотношения между частотами), то в этом случае удобно воспользоваться функцией bilinear следующих основных форматов:

[b,a]=bilinear(bs,as,Fs[,Fp])

[q,p,K]=bilinear(qs,ps,Ks,Fs[,Fp])

Здесь bs, as — векторы коэффициентов числителя и знаменателя передаточной функции АФП На(р) в порядке возрастания степеней р; ав(1)=1; Fs — частота дискретизации /д (Гц); Рр — необязательный параметр — частота /(Гц), на которой значения АЧХ АФП и ЦФ должны совпадать; Ь, а — векторы коэффициентов числителя и знаменателя передаточной функции ЦФ Н(г) (1) в порядке возрастания отрицательных степеней г, а(1) = 1; ц, р, К — соответственно векторы нулей и полюсов и коэффициент усиления передаточной функции, представленной в виде произведения простейших множителей:

ад-Чїо^'

тренной ранее для АФП, более того, используются те же функции, но без параметра '8':

[Ь,а]=Ьи«ег(К^Оп,Йуре)

[b,a]=cheby1(R,rp,WDn,ftype)

[Ь,а]=сЬеЬу2(К,г8^0п,Йуре)

[b,a]=eШp(R,rp,rs,WDn,ftype)

Здесь R — порядок ЦФ; WDn — вектор нормированных частот среза в шкале / (частотами среза называют частоты, на которых нормированная АЧХ ЦФ А(/) равна 1/^2*0,707, а А(/) (дБ) — 3 дБ), содержащий один элемент для ФНЧ и ФВЧ, равный:

™п(1,=/0 = ^,

где f0 — абсолютная частота среза, и два — для ПФ и РФ, равные:

■ Лі

WDn(l)=/01 = WDn(l) =f02 ■■

Л/2’

. /о2

Л/2;

(8)

где г0к, г*к — соответственно к-е нуль и полюс передаточной функции (1), а Ь0 — коэффициент усиления; qs, ps, Ks — аналогичные параметры для передаточной функции АФП.

• Если ЦФ синтезируется непосредственно по заданным требованиям к АЧХ, то в этом случае процедура синтеза подобна рассмо-

где f01, f02 — абсолютные частоты среза; rp, rs — соответственно максимально и минимально допустимые затухания amax (дБ) в ПП и amin (дБ) в ПЗ для A(f) (дБ) (они совпадают с допустимыми отклонениями для АФП); ftype — параметр, указывающий тип избирательности и принимающий значения: 'high' — для ФВЧ; 'stop' — для РФ; по умолчанию (если значение параметра не задано явно) — для ФНЧ и ПФ; b, a — векторы коэффициентов числителя и знаменателя передаточной функции ЦФ (1) в порядке возрастания отрицательных степеней z где a(1)=1.

Примечание. Здесь и далее в обозначениях частот ЦФ добавлена буква “D” от слова “Digital”. Согласно (5-6) зависимость между частотами WDn и Wn нелинейная.

При синтезе ЦФ Баттерворта, Чебышева I и II рода и Золотарева - Кауэра по методу билинейного Z-преобразования свойство оптимальности ЦФ сохраняется: синтезируемый ЦФ, подобно АФП, при заданных требованиях к АЧХ (характеристике затухания) имеет минимальный порядок.

Расчет минимального порядка R и вектора частот среза WDn для ЦФ Баттерворта, Чебышева I и II рода и Золотарева - Кауэра выполняется с помощью тех же функций, что и для АФП, но без параметра 's':

[R,WDn]=buttord(WDp,WDs,rp,rs)

[R,WDn]=cheb1ord(WDp,WDs,rp,rs)

[R,WDn]=cheb2ord(WDp,WDs,rp,rs)

[R,WDn]=ellipord(WDp,WDs,rp,rs)

Здесь WDp, WDs — соответственно векторы граничных нормированных частот ПП и ПЗ в порядке их следования слева направо в шкале /. Остальные параметры были определены ранее.

При синтезе ЦФ с помощью данных функций в MATLAB реализуется алгоритм билинейного Z-преобразования, а именно:

• по требованиям к АЧХ ЦФ автоматически формируются требования к АЧХ АФП путем пересчета граничных частот по формуле (5);

• синтезируется АФП;

• в соответствии с (3) передаточная функция АФП Ha(p) преобразуется в передаточную функцию ЦФ H(z) (1).

Подобно функции bilinear, выходными параметрами функций butter, cheby1, cheby2 и ellip могут быть [q,p,K].

Приведем примеры синтеза оптимальных БИХ-фильтров ФНЧ и ПФ непосредственно по заданным требованиям к АЧХ (дБ) ЦФ.

Пример 2

Заданы требования к АЧХ ФНЧ (табл. 3 и пример 2 в [6]). Значения amax = 0,4455 дБ и «min = 40 дБ (rp и rs) были вычислены в [6] в примере 1. Синтезировать оптимальные БИХ-фильтры Баттерворта, Чебышева I и II рода и Золотарева - Кауэра методом билинейного 2-преобразования:

>> Fs=8000;

>> ft=1000; fk=1500;

>> ft=1000; fk=1500;;

>> [R2,WDn2]=cheb1ord(WDp,WDs,rp,rs);

>> [R3,WDn3]=cheb2ord(WDp,WDs,rp,rs);

>> [R4,WDn4]=ellipord(WDp,WDs,rp,rs);

>> [b1,a1]=butter(R1,WDn1);

>> [b2,a2]=cheby1(R2,rp,WDn2);

>> [b3,a3]=cheby2(R3,rs,WDn3);

>> [b4,a4]=ellip(R4,rp,rs,WDn4);

Выведем рассчитанные значения порядков R1, R2, R3, и R4 соответственно оптимальных ФНЧ Баттерворта, Чебышева I и II рода и Золотарева - Кауэра:

>> R=[R1 R2 R3 R4]

R =

12 7 7 5

Поскольку свойство оптимальности синтезируемых ЦФ сохраняется, их порядки совпадают с порядками соответствующих АФП (переменная R в примере 1).

Построим графики АЧХ БИХ-фильтров ФНЧ Баттерворта, Чебышева I и II рода и Золотарева - Кауэра на густой сетке частот (выберем 1000 точек) в основной полосе [0; /д/2] и одинаковом диапазоне [0;1] по оси ординат, установленном с помощью функции уНш([0 1]). АЧХ рассчитывается с помощью функции (рис. 2):

>> %f — густая сетка частот

>> %Ha1,Ha2,Ha3,Ha4 — передаточные функции АФП Баттерворта, ... Чебышева I и II рода и Золотарева - Кауэра >> Fs=8000;

>> f=0:((Fs/2)/1000):Fs/2;

>> Ha1=freqz(b1,a1,f,Fs); MAG1=abs(Ha1);

>> Ha2=freqz(b2,a2,f,Fs); MAG2=abs(Ha2);

>> Ha3=freqz(b3,a3,f,Fs); MAG3=abs(Ha3);

>> Ha4=freqz(b4,a4,f,Fs); MAG4=abs(Ha4);

>> subplot(2,2,1),plot(f,MAG1),xlabel('f(Hz)'),grid,... ylabel('MAGNITUDE'),title('Digital Filter Butterworth'),ylim([0 1]) >> subplot(2,2,2),plot(f,MAG2),xlabel('f(Hz)'),grid,... ylabel('MAGNITUDE'),title('Digital Filter Chebyshov I'),ylim([0 1])

и

0,5

1

0 1000 2000 3000 4000 0 1000 2000 3000 4000

f(rit)

f(ru)

0,5

V

0 1000 2000 3000 4000 0 1000 2000 3000 4000

ПГЦ) f(ru)

Рис. 2. АЧХ БИХ-фильтров ФНЧ:

а) Баттерворта; б) Чебышева I рода; в) Чебышева II рода; г) Золотарева — Кауэра

0,5

AJVVi

/ V

0 1000 2000 3000 4000 0 1000 2000 3000 4000

f (Г ц)

«Гц)

0,5

'(Гц)

'(Гц)

Рис. 3. АЧХ БИХ-фильтров ПФ:

а) Баттерворта; б) Чебышева I рода; в) Чебышева II рода; г) Золотарева — Кауэра

>> subplot(2,2,3),plot(f,MAG3);xlabel('f(Hz)'),grid,... ylabel('MAGNITUDE'),title('Digital Filter Chebyshov II'),ylim([0 1]) >> subplot(2,2,4),plot(f,MAG4),xlabel('f(Hz)'),grid,... ylabel('MAGNITUDE'),title('Digital Filter Eleptic'),ylim([0 1])

Пример 3

Заданы требования к АЧХ ПФ (табл. 4 и пример 3 в [6]). Значения атах = 0,4455 дБ и аШш = 40 дБ (гр и ге) были вычислены в [6] в примере 1. Синтезировать оптимальные БИХ-фильтры Баттерворта, Чебышева I и II рода и Золотарева - Кауэра методом билинейного 2-преобразования.

Параметры WDp и WDs представляют собой векторы из двух элементов:

>> Рэ=8000;

>> £к1=1000; Ш=1400; 02=2000; Ас2=2400;

>> Й=[Й1 £12]; £к=[£к1 £к2];

>> WDp=ft./(Fs/2); WDs=£k./(Fs/2);

>> гр=0.4455; к=40;

>> [R1,WDn1]=buttoгd(WDp,WDs,гp,гs);

>> [R2,WDn2]=cheb1oгd(WDp,WDs,гp,гs);

>> [R3,WDn3]=cheb2oгd(WDp,WDs,гp,гs);

>> [R4,WDn4]=eШpoгd(WDp,WDs,гp,гs);

>> [b1,a1]=butteг(R1,WDn1);

>> [b2,a2]=cheby1(R2,гp,WDn2);

>> [b3,a3]=cheby2(R3,гs,WDn3);

>> [b4,a4]=eШp(R4,гp,гs,WDn4);

Выведем рассчитанные значения порядков R1, R2, R3, и R4 соответственно оптимальных ПФ Баттерворта, Чебышева I и II рода и Золотарева - Кауэра:

>> R=[R1 R2 R3 R4]

R =

5 5 tv 4

Наименьший порядок имеет ФНЧ Золотарева - Кауэра.

Графики АЧХ БИХ-фильтров (рис. 3) ПФ Баттерворта, Чебышева I и II рода и Золотарева - Кауэра строятся так же, как для ФНЧ (пример 2).

Анализ БИХ-фильтра

В состав MATLAB входит программа GUI FVTool (Filter Visualization Tool — средства визуализации фильтра), предназначенная для анализа характеристик синтезированных ЦФ в окне Figure : Filter Visualization Tool, обращение к которой производится с помощью функции fvtool:

fvtool(b,a)

Здесь b, a — векторы коэффициентов передаточной функции БИХ-фильтра. ■

Литература

1. Ingle V., Proakis J. Digital Signal Processing Using MATLAB. Second Edition — Thomson.

2. Оппенгейм А., Шафер Р. Цифровая обработка сигналов. М.: Техносфера, 2006.

3. Сергиенко А. Б. Цифровая обработка сигналов, 2-е изд. СПб.: ПИТЕР, 2006.

4. Солонина А. И., Улахович Д. А., Арбузов С. М., Соловьева Е. Б. Основы цифровой обработки сигналов. 2-е изд. СПб.: БХВ-Петербург, 2005.

5. Солонина А. И., Арбузов С. М. Цифровая обработка сигналов. Моделирование в MATLAB. СПб.: БХВ-Петербург, 2008.

6. Солонина А. Моделирование цифровой обработки сигналов в MATLAB. Часть 1. Синтез оптимальных (по Чебышеву) КИХ-фильтров программными средствами MATLAB // Компоненты и технологии. 2008. № 11.