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

Все статьи цикла:

• Часть 1. Синтез оптимальных (по Чебышеву) КИХ-фильтров программными средствами MATLAB
• Часть 2. Синтез оптимальных цифровых БИХ-фильтров программными средствами MATLAB
• Часть 3. Описание структур КИХ- и БИХ-фильтров в MATLAB
• Часть 4. Моделирование структур цифровых фильтров с фиксированной точкой программными средствами MATLAB: анализ характеристик КИХ-фильтров
• Часть 5. Моделирование структур цифровых фильтров с фиксированной точкой программными средствами MATLAB: анализ характеристик БИХ-фильтров
• Часть 6. Моделирование структур цифровых фильтров с фиксированной точкой программными средствами MATLAB: квантование воздействия и вычисление реакции
• Часть 7. Моделирование цифровых фильтров средствами программ GUI MATLAB: GUI FDATool
• Часть 8. Моделирование цифровой фильтрации средствами программ GUI MATLAB: GUI SPTool
• Часть 9. Моделирование цифровых преобразователей Гильберта и дифференциаторов программными средствами и средствами GUI FDATool MATLAB

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

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

(1)

и при (N−1)≤(M−1) (по умолчанию) имеет порядок R = M−1.

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

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

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

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

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

• частоту дискретизации f_д (Гц);
• граничные частоты полос пропускания (ПП) и полос задерживания (ПЗ), для которых введены условные обозначения:
  • f_χ — граничная частота ПП для ФНЧ и ФВЧ;
  • f_k — граничная частота ПЗ для ФНЧ и ФВЧ;
  • f_−χ, f_χ — левая и правая граничные частоты ПП для ПФ и РФ;
  • f_−k, f_k — левая и правая граничные частоты ПЗ для ПФ и РФ;
• допустимые отклонения от Â(f) (дБ) (см. формулу (6) в [6]):
  • a_max (дБ) — максимально допустимое затухание в ПП;
  • a_min (дБ) — минимально допустимое затухание в ПЗ.

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

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

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

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

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

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

1. Задание требований к АЧХ ЦФ.
2. Выбор метода синтеза.
3. Формирование требований к АЧХ АФП.
4. Выбор типа аппроксимирующей функции. Четырем типам аппроксимирующих функций соответствуют четыре разновидности аналоговых (и цифровых) фильтров:

• Баттерворта (Butterwhorth) — с АЧХ, максимально плоской в ПП и монотонной в ПЗ;
• Чебышева I рода (Chebyshov Type I) — с АЧХ, равноволновой в ПП и монотонной в ПЗ;
• Чебышева II рода (Chebyshov Type II) — с АЧХ, максимально плоской в ПП и равноволновой в ПЗ;
• Золотарева – Кауэра (эллиптические фильтры) (Eleptic) — с АЧХ, равноволновой в ПП и ПЗ.

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

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

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

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

Здесь R — порядок АФП; Wn — вектор частот среза в шкале ω = 2πf (рад/с), содержащий один элемент — для ФНЧ и ФВЧ и два — для ПФ и РФ (частотами среза называют частоты, на которых нормированная АЧХ АФП Â(f) равна 1/√2≈0,707, а Â(f) (дБ) — 3 дБ); rp, rs — соответственно максимально и минимально допустимые затухания a_max (дБ) в ПП и a_min (дБ) в ПЗ для Â(f) (дБ); ftype — параметр, указывающий тип избирательности и принимающий значения: ‘high’ — для ФВЧ; ‘stop’ — для РФ; по умолчанию (если значение параметра не задано явно) — для ФНЧ и ПФ; ‘s’ — признак аналогового фильтра, при его отсутствии по умолчанию подразумевается ЦФ; bs, as — векторы коэффициентов числителя и знаменателя передаточной функции АФП H_a(p) в порядке возрастания степеней p; as(1) = 1.

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

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

[R,Wn]=buttord(Wp,Ws,rp,rs,'s') 
[R,Wn]=cheb1ord(Wp,Ws,rp,rs,'s') 
[R,Wn]=cheb2ord(Wp,Ws,rp,rs,'s') 
[R,Wn]=ellipord(Wp,Ws,rp,rs,'s')

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

Пример 1

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

 <  <  ft=1000; fk=1500; 
 <  <  Wp=2.*pi.*ft; Ws=2.*pi.*fk; 
 <  <  rp=0.4455; rs=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]=ellipord(Wp,Ws,rp,rs,'s'); 
 <  <  [bs1,as1]=butter(R1,Wn1,'s'); 
 <  <  [bs2,as2]=cheby1(R2,rp,Wn2,'s'); 
 <  <  [bs3,as3]=cheby2(R3,rs,Wn3,'s'); 
 <  <  [bs4,as4]=ellip(R4,rp,rs,Wn4,'s');

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

 <  <  R=[R1 R2 R3 R4] 
R = 
	15   7   7   5

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

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

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

Ha=freqs(bs,as,W)

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

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

 <  <  %f — густая сетка частот в Гц 
 <  <  %W — густая сетка круговых частот в рад/с 
 <  <  %Ha1,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([0 1]) 
 <  <  subplot(2,2,4),plot(f,MAG4);xlabel('f(Hz)'),grid,... ylabel('MAGNITUDE'),title('Analog Filter Eleptic'),ylim([0 1])

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

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

(2)

где γ = 1/T, получаемой из стандартного Z–преобразования z = e^pT → p = (1/T)lnz, путем разложения lnz в ряд Тейлора:

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

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

(3)

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

И подставляя z = re^jˆω  и p = jΩ, где jΩ — обозначение оси частот АФП (во избежание путаницы), получаем:

(4)

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

		(5)
		(6)

Согласно (4), ось частот jΩ p–плоскости, как и при стандартном Z–преобразовании, отображается на z–плоскости в единичную окружность (радиус равен единице), однако каждому ее обороту (изменению нормированной частоты на Δˆω  = 2π), а именно: …, −3π < ˆω  < −π, −π < ˆω  < π, π < ˆω  < 3π, …, соответствует не отрезок оси, как при стандартном Z–преобразовании, а вся ось jΩ (так как зависимость между частотами определяется функцией arctg: …, −∞ < Ω < ∞, −∞ < Ω < ∞, −∞ < Ω < ∞…

Связь между частотными характеристиками АФП и ЦФ, соответственно H(e^jΩT) и H_a(jΩ), имеет вид [1, 4]:

(7)

При этом частотная характеристика АФП H_a(jΩ), бесконечная в шкале частот Ω, в шкале частот ω, согласно (6), сжимается в отрезок Δω = ω_д, то есть становится финитной. Соответственно, частотная характеристика ЦФ H(e^jΩT), согласно (7), оказывается равной (с точностью до множителя 1/T) бесконечной сумме копий финитных частотных характеристик АФП длины Δω = ω_д, сдвинутых друг относительно друга на частоту ω_д. Вследствие этого элайсинг (Aliasing) принципиально отсутствует, и АЧХ ЦФ и АФП в основной полосе частот [0; ω_д/2] совпадают (с учетом нелинейной зависимости между частотами).

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

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

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

[b,a]=bilinear(bs,as,Fs[,Fp]) 
[q,p,K]=bilinear(qs,ps,Ks,Fs[,Fp])

Здесь bs, as — векторы коэффициентов числителя и знаменателя передаточной функции АФП H_a(p) в порядке возрастания степеней p; as(1)=1; Fs — частота дискретизации f_д (Гц); Fp — необязательный параметр — частота f (Гц), на которой значения АЧХ АФП и ЦФ должны совпадать; b, a — векторы коэффициентов числителя и знаменателя передаточной функции ЦФ H(z) (1) в порядке возрастания отрицательных степеней z; a(1)=1; q, p, K — соответственно векторы нулей и полюсов и коэффициент усиления передаточной функции, представленной в виде произведения простейших множителей:

(8)

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

• Если ЦФ синтезируется непосредственно по заданным требованиям к АЧХ, то в этом случае процедура синтеза подобна рассмотренной ранее для АФП, более того, используются те же функции, но без параметра ‘s’:

[b,a]=butter(R,WDn,ftype) 
[b,a]=cheby1(R,rp,WDn,ftype) 
[b,a]=cheby2(R,rs,WDn,ftype) 
[b,a]=ellip(R,rp,rs,WDn,ftype)

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

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

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

При синтезе ЦФ Баттерворта, Чебышева 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 — соответственно векторы граничных нормированных частот ПП и ПЗ в порядке их следования слева направо в шкале ^ˆf . Остальные параметры были определены ранее.

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

• по требованиям к АЧХ ЦФ автоматически формируются требования к АЧХ АФП путем пересчета граничных частот по формуле (5);
• синтезируется АФП;
• в соответствии с (3) передаточная функция АФП H_a(p) преобразуется в передаточную функцию ЦФ H(z) (1).

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

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

Пример 2

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

 <  <  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; f_д/2] и одинаковом диапазоне [0;1] по оси ординат, установленном с помощью функции ylim([0 1]). АЧХ рассчитывается с помощью функции freqz (рис. 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])
 <  <  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]). Значения a_max = 0,4455 дБ и a_min = 40 дБ (rp и rs) были вычислены в [6] в примере 1. Синтезировать оптимальные БИХ–фильтры Баттерворта, Чебышева I и II рода и Золотарева – Кауэра методом билинейного Z–преобразования.

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

 <  <  Fs=8000; 
 <  <  fk1=1000; ft1=1400; ft2=2000; fk2=2400; 
 <  <  ft=[ft1 ft2]; fk=[fk1 fk2]; 
 <  <  WDp=ft./(Fs/2); WDs=fk./(Fs/2); 
 <  <  rp=0.4455; rs=40; 
 <  <  [R1,WDn1]=buttord(WDp,WDs,rp,rs); 
 <  <  [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 = 
	7   5   5   4

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

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

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

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

fvtool(b,a)

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

Часть 3. Описание структур КИХ- и БИХ-фильтров в MATLAB

Литература

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.