Синтез и анализ цифровых фильтров с использованием программного пакета MatLab

Содержание

1. Начало работы . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .1

2. Ввод коэффициентов передаточной функции и построение графиков частотных характеристик аналогового фильтра-прототипа . . . . . . . . . . . . .1

3. Нули и полюсы аналогового фильтра-прототипа и их преобразование

в коэффициенты передаточной функции . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3

4. Расчёт аналогового фильтра-прототипа по заданным требованиям к

его характеристике затухания . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4

5. Нахождение передаточной функции цифрового фильтра по

аналоговому прототипу методом билинейного z-преобразования и

методом инвариантной импульсной характеристики . . . . . . . . . . . . . . . . . . 6

6. Структуры цифровых фильтров и соответствующие им алгоритмы цифровой фильтрации . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .9

7. Просмотр характеристик синтезированного цифрового фильтра . . . . . 19

8. Синтез цифрового фильтра с использованием программы fdatool . . . . 20

9. Исследование влияния квантования сигналов и коэффициентов

фильтра на характеристики фильтра . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23

10. Моделирование работы цифрового фильтра . . . . . . . . . . . . . . . . . . . . . . . . . 29

1. Начало работы

Произведите загрузку системы MatLab. Откроется командное окно (Command Window). Вы будете иметь возможность работы в интерактивном режиме, вводя операторы после значка >>. Кроме того, в командном окне можно ввести указанное в текущем каталоге имя той или иной программы, которой затем будет передано управление. Для очистки командного окна используется команда clc.

2. Ввод коэффициентов передаточной функции фильтра и построение графиков частотных характеристик аналогового фильтра-прототипа

Передаточная функция аналогового прототипа задаётся в виде:

, m  n. (1)

- 2 -

Для построения графиков АЧХ и ФЧХ нужно ввести векторы b и a коэффициентов передаточной функции, а затем вызвать подпрограмму расчёта комплексной частотной характеристики и построения графиков. Пусть, например, m=1, n=2, b₀=2, b₁= -1, a₀=1, a₁=3, a₂=2.5. В интерактивном режиме вводятся следующие операторы:

>> b=[2 –1];

>> a=[1 3 2.5];

>> freqs (b,a)

Элементы векторов разделяются пробелами и заключаются в квадратные скобки. После ввода вектора ставится точка с запятой, а затем следует нажать клавишу . После оператора freqs(b,a) точка с запятой не ставится, ввод тоже следует завершить нажатием клавиши . Команды запоминаются, и их можно воспроизводить, пользуясь клавишами управления курсором:  и . Будут построены графики АЧХ и ФЧХ (АЧХ в логарифмическом масштабе, но без пересчёта в децибелы, ФЧХ в градусах). По умолчанию выбираются 200 частот, логарифмически равномерно распределённых в диапазоне от 0.1 до 10.

Если нужно построить АЧХ в линейном масштабе, в ином диапазоне частот, ФЧХ в радианах и с устранением скачков на 2k радиан, то следует ввести следующие операторы (пусть, для примера, границы частотного диапазона от f1=100 Гц до f2=1000 Гц, число расчётных точек N =451):

>> f=100:(1000-100)/450:1000;

>> b=[1 0 7.971e6];

>> a=[1 7.427e2 1.501e6 5.536e8];

>> k=freqs (b,a,f*2*pi);

>> subplot (2,1,1)

>> plot (f,abs(k)/max(abs(k))), grid

>> subplot (2,1,2)

>> plot (f,unwrap(angle(k))), grid

При вводе строки subplot (2,1,1) появится графическое окно

Figure No.1. Его нужно сместить, чтобы оно не загораживало командное окно, или временно свернуть. В последнем случае, по окончании ввода операторов нужно развернуть графическое окно, чтобы увидеть графики.

График АЧХ аналогового фильтра-прототипа можно использовать для нахождения частоты дискретизации Fs. Например, может быть задано, что АЧХ на частоте Fs/2 должна достигать определённого уровня. Для детального рассмотрения нужной части графика рекомендуется использовать инструменты изменения масштаба (уменьшение, увеличение), имеющиеся

на панели инструментов графического окна.

• 3 -

3. Нули и полюсы аналогового фильтра-прототипа и их преобразование в коэффициенты передаточной функции

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

p₁=-0.2+0.5i и p₂=-0.2-0.5i, а также нули z₁=0.1i и z₂=-0.1i. Осуществляем ввод:

>> p=[-0.2+0.5i -0.2-0.5i]’;

>> z=[0.1i -0.1i]’;

>> k=5;

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

Далее производим преобразование полюсов и нулей в коэффициенты передаточной функции:

>> [b,a]= zp2tf (z,p,k);

Здесь b – строка коэффициентов числителя, a – строка коэффициентов знаменателя передаточной функции фильтра.

Диаграмму нулей и полюсов можно вывести, если задать следующие команды:

>> plot (p,’x’)

>> hold on

>> plot (z,’o’)

>> hold off

>> axis equal

>> axis ([-1 1 -1 1])

- 4 -

Последняя команда позволяет задать граничные значения для действительной (первые два числа) и мнимой (следующие два числа) осей. Возможно и обратное преобразование коэффициентов b и a в полюсы и нули:

>> [z,p,k]= tf2zp (b,a);

Если требуется вывести значения каких-либо переменных рабочей области MatLab на монитор, то нужно просто после значка >> ввести имя переменной, завершив ввод нажатием клавиши . Например, если указать

>> z

то будет выведен вектор нулей. Вслед за оператором

>> p

будет выведен вектор полюсов. Если не ставить точку с запятой после операторов преобразования, то выходные параметры будут выведены на монитор без дополнительных команд.

4. Расчёт аналогового фильтра-прототипа по заданным требованиям к его характеристике затухания

Если заданы граничные частоты, допустимые затухания, тип фильтра(ФНЧ, ФВЧ, ППФ, ПЗФ), класс фильтра (Баттерворта, Чебышёва, инверсный Чебышёва, эллиптический, Бесселя), то рассчитать коэффициенты передаточной функции фильтра можно, используя одну из следующих функций MatLab:

1. butter (n, w0, type, ‘s’) – расчёт фильтра Баттерворта;

2. cheby1 (n, Rp, w0, type, ‘s’) – расчёт фильтра Чебышёва;

3. cheby2 (n, Rs, w0, type, ‘s’) – расчёт инверсного фильтра Чебышёва;

4. ellip (n, Rp, Rs, w0, type, ‘s’) – расчёт эллиптического фильтра;

5. besself (n, w0, type) – расчёт фильтра Бесселя.

Здесь n – порядок ФНЧ-фильтра-прототипа; Rp – неравномерность характеристики затухания в полосе пропускания (в дБ); Rs – минимальное затухание в полосе задерживания (в дБ); w0 – скаляр или двухэлементный вектор частот среза (границ полосы пропускания ) (в рад/с) (для фильтра Баттерворта частота среза задаётся по уровню 3 дБ). Если для фильтра Баттерворта граничная частота полосы пропускания wp задана для

- 5 -

значения затухания Rp (Rp<3 дБ), то её нужно пересчитать в частоту среза

type – строковая константа, задающая тип фильтра. Строка ‘s’ говорит о том, что синтезируется аналоговый фильтр.

При синтезе ФНЧ: w0 – скаляр , параметр type отсутствует.

При синтезе ФВЧ: w0 – скаляр , type=’high’.

При синтезе ППФ: w0 – вектор [w01 w02], причём w01

При синтезе ПЗФ: w0 – вектор [w01 w02], причём w01

Если порядок ФНЧ-прототипа не задан, то его можно определить (имеется в виду минимальный порядок), применяя следующие функции MatLab:

[n,wn]=buttord (wp, ws, Rp, Rs, ‘s’)

[n,wn]=cheb1ord (wp, ws, Rp, Rs, ‘s’)

[n,wn]=cheb2ord (wp, ws, Rp, Rs, ‘s’)

[n,wn]=ellipord (wp, ws, Rp, Rs, ‘s’)

Здесь n – минимальный порядок фильтра; wn – частота среза фильтра (для фильтра Баттерворта она определяется по уровню 3 дБ); wp – граничная частота полосы пропускания (в рад/с); ws – граничная частота полосы задерживания (в рад/с).

Для ППФ и ПЗФ wp и ws –двухэлементные векторы: [wp1 wp2] и [ws1 ws2].

Должны выполняться неравенства:

для ФНЧ: wp

для ФВЧ: wp>ws;

для ППФ: ws1

для ПЗФ: wp1

Найденный порядок n далее используется в функциях расчёта коэффициентов передаточной функции фильтра, описанных выше.

Рассмотрим пример. Пусть требуется найти коэффициенты передаточной функции эллиптического полосно-пропускающего фильтра, для которого Rp=1 дБ, Rs=40 дБ, wp1=1000 рад/с, wp2=1100 рад/с, ws1=900 рад/с, ws2=1250 рад/с. Вводим операторы:

>> [n, wn]=ellipord([1e3 1.1e3], [0.9e3 1.25e3], 1, 40, ‘s’);

>> [b,a]=ellip(n, 1, 40, [1e3 1.1e3], ‘s’);

- 6 -

Если требуется определить полюсы и нули, то вместо последнего оператора вводят:

>> [z,p,k]= ellip (n, 1, 40, [1e3 1.1e3], ‘s’);

Альтернатива: использование функции преобразования:

>> [z,p,k]= tf2zp (b,a);

5. Нахождение передаточной функции цифрового фильтра по аналоговому прототипу методом билинейного

z-преобразования и методом инвариантной импульсной характеристики

5.1. Билинейное z-преобразование

Это преобразование задаётся путём замены переменных:

, (2)

где p – комплексная частота аналогового фильтра-прототипа, z – комплексная переменная, от которой зависит передаточная функция цифрового фильтра, Т – интервал дискретизации (Т=1/Fs, где Fs – частота дискретизации). Достаточно в передаточной функции аналогового прототипа сделать подстановку (2), и будет найдена передаточная функция цифрового фильтра К(z):

(3)

Если а0 1, то нужно разделить все коэффициенты числителя и знаменателя на а0, чтобы знаменатель записывался в виде:

- 7 -

Альтернативным способом нахождения передаточной функции K(z) методом билинейного z-преобразования является пересчёт полюсов и нулей аналогового прототипа в полюсы и нули цифрового фильтра по формуле:

Затем осуществляется преобразование нулей и полюсов в коэффициенты фильтра (см. разд.3; функции MatLab, описанные здесь для аналоговых фильтров, годятся и для цифровых фильтров). Если количество полюсов аналогового ФНЧ-прототипа превышает количество его нулей, то возникают также дополнительные нули, так что общее количество нулей и полюсов цифрового фильтра, синтезированного по методу билинейного z-преобразования, оказывается равным (об этом см. в специальной литературе).

Указанные преобразования применимы, если аналоговый фильтр-прототип и цифровой фильтр – это фильтры одинакового типа (оба ФНЧ или оба ФВЧ и т.д.). Значит, если задан ФНЧ-прототип, а его нужно преобразовать в цифровой фильтр другого типа (ФВЧ, ППФ, ПЗФ), то нужно прежде найти передаточную функцию аналогового фильтра этого типа, а затем применить билинейное z-преобразование. Соответствующее преобразование типов фильтра осуществляется в MatLab операторами:

>> [b1,a1]=lp2hp (b, a, w0);

- преобразование ФНЧ в ФВЧ; w0 – граничная частота ФВЧ (рад/с).

>> [b1,a1]=lp2bp (b, a, w0, Bw);

- преобразование ФНЧ в ППФ; w0 – средняя геометрическая частота

полосы пропускания (рад/с) (w0=w1*w2); Bw – полоса пропускания (рад/с)

(Bw=w1-w2).

>> [b1,a1]=lp2bs (b, a, w0, Bw);

- преобразование ФНЧ в ПЗФ; w0 – средняя геометрическая частота полосы

з адерживания (рад/с) (w0=w1*w2 ); Bw – ширина полосы задерживания

(рад/с) (Bw=w1-w2).

Само билинейное преобразование в MatLab осуществляется следующими операторами:

>> [bz, az]= bilinear (b, a, Fs);

- 8 -

или

>> [zz, pz, kz]= bilinear (z, p, k, Fs);

Здесь b, a, z, p, k – коэффициенты передаточной функции, нули, полюсы и масштабный коэффициент передаточной функции аналогового фильтра-прототипа. bz, az, zz, pz, kz – соответствующие параметры цифрового фильтра. Векторы b и a должны задаваться как векторы-строки, z и p – как векторы-столбцы. Преобразование строки в столбец и наоборот осуществляется путём постановки символа ‘ (апостроф) после имени вектора. Например, операция

>> d=d’;