Проектирование цифровых фильтров в MATLAB с использованием возможностей пакета Signal Processing Toolbox

3.1.0. Расчет частотной характеристики

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

freqz(b,a);

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

При расчете используются нормированные значения частот, измеряемые в радианах на отсчет (при такой нормировке частота дискретизации равна , а частота Найквиста - ). По умолчанию выбирается 512 частотных точек, равномерно распределенных в диапазоне .

При отсутствии выходных параметров функция freqz строит графики АЧХ (в децибелах) и ФЧХ (в градусах) фильтра.

Построение графика частотных характеристик фильтра с использованием функции freqz производится следующим образом:

>> b=[1 2 3 4];

>> a=[1 -0.1];

>> freqz(b,a)

Полученное в результате выполнения команд графическое окно с графиками АЧХ и ФЧХ цифрового фильтра имеет вид, приведенный на рисунке.

Рисунок - графическое окно с графиками АЧХ и ФЧХ цифрового фильтра, полученного с использованием функции freqz(b,a)

 

Если при вызове функции указаны выходные параметры, построение графика не производится:

[h,w] = freqz(b,a);

В этом случае в векторе h возвращаются рассчитанные значения комплексного коэффициента передачи, а в векторе w – использованные при расчете значения нормированных частот.

Чтобы задать количество частотных точек для расчета, используется третий входной параметр n (частоты при этом по прежнему равномерно распределяются в диапазоне ):

freqz(b,a, n);

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

>> b=[1 2 3 4];

>> a=[1 -0.1];

>> freqz(b,a,5)

Полученное в результате выполнения команд графическое окно с графиками АЧХ и ФЧХ цифрового фильтра имеет вид, приведенный на рисунке.

Рисунок - графическое окно с графиками АЧХ и ФЧХ цифрового фильтра, полученного с использованием функции freqz(b,a,5)

 

Чтобы перейти от использования нормированных частот к использованию абсолютных значений частот в Гц, используется четвертый входной параметр Fs, равный частоте дискретизации (частоты при этом равномерно распределяются в диапазоне ):

freqz(b,a, n, Fs);

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

>> b=[1 2 3 4];

>> a=[1 -0.1];

>> freqz(b,a,512,2000)

Полученное в результате выполнения команд графическое окно с графиками АЧХ и ФЧХ цифрового фильтра имеет вид, приведенный на рисунке.

Рисунок - графическое окно с графиками АЧХ и ФЧХ цифрового фильтра, полученного с использованием функции freqz(b,a,512, 2000)

 

Можно заставить функцию freqz использовать ряд частот, равномерно распределенных на интервале вплоть до частоты дискретизации. Это является необходимым при анализе фильтров с комплексными коэффициентами, частотные характеристики которых не являются симметричными относительно частоты Найквиста. Такое указание дается функции freqz с помощью четвертого входного параметра – строки ‘whole’:

freqz(b,a, n, ‘whole’);

Есть возможность задать вектор нормированных круговых частот w для расчета частотной характеристики:

freqz(b,a,w);

3.1.1. Проектирование цифровых фильтров с бесконечной импульсной характеристикой

3.1.1.1. Использование аналоговых прототипов

Для проектирования цифровых фильтров с бесконечной импульсной характеристикой (БИХ, IIF) используются функции расчета фильтров, учитывающие аналоговые прототипы [3]:

[b,a]=butter(n,Wn,’ftype’) (фильтр Баттерворта)

[b,a]=cheby1(n,Rp,Wn,’ftype’) (фильтр Чебышева I)

[b,a]=cheby2(n,Rs,Wn,’ftype’) (фильтр Чебышева II)

[b,a]=ellip(n,Rp,Rs,Wn,’ftype’)  (эллиптические фильтры).

Фильтры Баттерворта отличаются максимально гладкой и плоской АЧХ в полосе пропускания, однако спады и подъемы частотной характеристики этих фильтров менее крутые, чем у фильтров Чебышева и эллиптических фильтров, особенно при малых порядках [4]. Фильтры Чебышева I рода характеризуются равновеликими пульсациями АЧХ в полосе пропускания и монотонным изменением ослабления в полосе задержания. Фильтры Чебышева II рода характеризуются равновеликими пульсациями АЧХ в полосе задержания и монотонным изменением ослабления в полосе пропускания. Они не обеспечивают такого же крутого спада АЧХ, как фильтры I рода, однако в полосе пропускания их АЧХ меняются монотонно, т.е. отсутствуют пульсации. Эллиптические фильтры характеризуются более крутыми спадами и подъемами АЧХ, чем фильтры Баттерворта и Чебышева. Их АЧХ имеют равновеликие пульсации в полосе пропускания и в полосе задержания.

Параметр Wn задает частоту среза ФНЧ, если это скаляр. Следует помнить, что при обращении к указанным функциям при проектировании цифровых фильтров частота среза Wn задается в единицах нормированной частоты в диапазоне от 0 до 1, где 1 соответствует половине частоты дискретизации или частоте Найквиста. Частота среза Wn соответсвует значению амплитуды, равному . Если его задать как двухкомпонентный вектор Wn=[w1,w2], то функция будет возвращать характеристики полосового фильтра. При отсутствии параметра ’ftype’ получаются характеристики ФНЧ или полосового фильтра. Параметр ’ftype’ имеет два значения, что позволяет проектировать еще два важных типа фильтров: high (ФВЧ), stop (режекторный с заданием Wn=[w1,w2]).

Параметр Rp задает неравномерность АЧХ в децибелах в полосе пропускания. Параметр Rs задает минимальное затухание АЧХ в децибелах в полосе задержания.

3.1.1.2. Использование функций maxflat и yulewalk  

В группу функций проектирования фильтров с БИХ и КИХ включены еще две функции [3]:

[b,a]=maxflat(nb,na,Wn)  

[b,a]=yulewalk(n,f,m)  

Функция maxflat задает проектирование обобщенного фильтра Баттерворта, у которого передаточная функция представлена полиномами числителя и знаменателя с разной степенью: nb и na соответственно. Частота среза этого фильтра Wn задается на уровне – 3 дБ. Она находится между 0 и 1, где 1 соответствует частоте Найквиста.

Расчёт цифрового фильтра с использованием функции maxflat осуществляется следующим образом:

 >> nb=7;

>> na=2;

>> Wn=0.5;

>> [b,a]=maxflat(nb,na,Wn,'plots')

b = 0.0642 0.3009 0.5436 0.4514 0.1407  -0.0156 -0.0097 0.0021

a = 1.0000 0.0000 0.4776

Полученное в результате выполнения команды графическое окно с графиками АЧХ синтезированного цифрового фильтра имеет вид, приведенный на рисунке.

Рисунок – графическое окно с графиками АЧХ цифрового фильтра ФНЧ, синтезированного с использованием функции [b,a]=maxflat(nb,na,Wn,'plots')

 

Функция yulewalk служит для проектирования рекурсивного цифрового фильтра Юла-Уокера порядка n с заданной АЧХ, которая представляется векторами частот f и значений m. Частоты должны быть заданы в интервале , где 1 соответствует частоте Найквиста (частота, равная половине частоты дискретизации). Вектор частот должен иметь первый элемент 0 и последний 1, а в промежутке частоты указываются в возрастающем порядке. Если надо задать скачок АЧХ, то частота скачка указывается подряд дважды, а в векторе m указываются значения АЧХ до скачка и после скачка.

Рисунок – графическое окно с графиками АЧХ цифрового фильтра, синтезированного с использованием функции [b,a]=yulewalk(10,f,m)

Использование функции yulewalk иллюстрируется следующим примером:

>> f=[0 0.2 0.2 0.6 0.6 1];

>> m=[0 0 1 1 0 0];

>> [b,a]=yulewalk(10,f,m);

>> [h,w]=freqz(b,a,128);

>> plot(f,m,w/pi,abs(h))

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

3.1.1.3. Выбор минимального порядка фильтров с БИХ

В пакет Signal Processing Toolbox включены 4 функции [3], которые обеспечивают вычисление минимального порядка и граничной частоты Wn фильтров с БИХ по следующим заданным параметрам: Wp – частота пропускания (среза), Ws – частота задержания,  Rp – максимально допустимые пульсации в области пропускания в децибелах и  Rs – минимальное затухание в полосе задержания в децибелах. К этим функциям относятся:

[n,Wn]=buttord(Wp,Ws,Rp,Rs)

[n,Wn]=cheb1ord(Wp,Ws,Rp,Rs)

[n,Wn]=cheb2ord(Wp,Ws,Rp,Rs)

[n,Wn]=ellipord(Wp,Ws,Rp,Rs)

Для ВЧ фильтров должно соблюдаться условие: Wp > Ws. Для полосовых и режекторных фильтров Wp и Ws должны задаваться двухэлементными векторами, элементы которых задают граничные частоты полос. При этом первой указывается более низкая частота. В этом случае параметр Wn также будет двухэлементным вектором.

Величина Wp нормирована и должна находиться в диапазоне , где 1 соответсвует частоте Найквиста или половине частоты дискретизации [2]. Величина Ws находится в том же диапазоне, что и Wp. Кроме минимального порядка также вычисляется частота среза Wn, которая далее может быть использована в качестве исходной величины для расчета параметров фильтров.

Расчет ФНЧ Баттерворта, обеспечивающего затухание не более 3 дБ в диапазоне от 0 до 100 Гц и не менее 15 дБ в диапазоне от 150 Гц до частоты Найквиста (частота дискретизации 1000 Гц) выполняется следующим образом >> Wp=100/500;

>> Ws=150/500;

>> [n,Wn]=buttord(Wp,Ws,3,15)

n =4

Wn = 0.2042

>> [b,a]=butter(n,Wn);

>> freqz(b,a,512,1000):

Полученное в результате выполнения команды графическое окно с графиками АЧХ и ФЧХ синтезированного цифрового фильтра имеет вид, приведенный на рисунке.

Рисунок – графическое окно с графиками АЧХ и ФЧХ цифрового фильтра, синтезированного с использованием функции [b,a]=butter(n,Wn)

 

3.1.1.4. Методы преобразования аналоговых фильтров в цифровые

Для преобразования аналоговых фильтров в цифровые используются функции[3]:

[numd,dend]=bilinear(num,den,fs) (метод билинейного преобразования)

[numd,dend]=bilinear(num,den,fs,Fp) (метод билинейного преобразования)

[bz,az]=impinvar(b,a,fs,tol]) (метод инвариантности имп. характеристики)

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

,

где  fs – частота дискретизации.

Параметр Fp является дополнительным и задает предискажения. Он определяет частоту в Гц, для которой значения АЧХ до и после преобразования совпадают. В этом случае билинейное преобразование видоизменяется следующим образом:

.

При методе инвариантного импульсного преобразования осуществляется преобразование аналогового фильтра в цифровой фильтр, имеющий ту же импульсную характеристику, что и у аналогового фильтра с точностью до множителя 1/f_s.

 

3.1.2. Проектирование цифровых фильтров с конечной импульсной характеристикой

3.1.2.1. Метод Ремеза

Для проектирования цифровых фильтров с КИХ, имеющих равные пульсации как в полосе пропускания, так и в полосе задержания, по методу ремеза служит функция [3]:

b=cremez(n,f,’fresp’)

Эта функция вычисляет вектор b, содержащий (n+1) коэффициентов фильтра с КИХ, АЧХ которого наилучшим образом представляет некоторую функцию fresp. Вектор f должен содержать пары частот, граничных для областей пропускания и задержания. Частоты должны быть в диапазоне от -1 до +1, причем +1 соответствует частоте Найквиста (половине частоты дискретизации). Частоты должны идти в строго нарастающем порядке. Параметр fresp может принимать различные значения: lowpass (ФНЧ), highpass (ФВЧ), bandpass (полосовой), bandstop (задерживающий), differentiator (дифференцирующий), hilfilt (Гильберта).

Расчёт фильтра ФНЧ 20-го порядка по методу Ремеза осуществляется следующим образом:

>> b=cremez(20,[-1 -0.4 -0.3 0.6 0.7 1],'lowpass');

>> freqz(b,1,512,'whole');

Полученное в результате выполнения второй команды графическое окно с графиками АЧХ и ФЧХ синтезированного цифрового фильтра ФНЧ имеет вид, приведенный на рисунке. Для вывода значений коэффициентов фильтра необходимо убрать знак; после первой команды, либо дополнительно ввести команду b.

Рисунок – графическое окно с графиками АЧХ и ФЧХ цифрового фильтра ФНЧ, синтезированного с использованием функции b=cremez(20,[-1 -0.4 -0.3 0.6 0.7 1],'lowpass')

 

3.1.2.2. Метод взвешивания

 

Для проектирования фильтров с КИХ методом взвешивания служит функция [3]:

b=fir1(n,Wn,’ftype’,window)

В простейшей форме b=fir1(n,Wn) функция создает вектор из n+1 коэффициентов ФНЧ с КИХ и линейной АЧХ для заданной частоты среза Wn (от 0 до 1). Групповое время задержки таких фильтров равно n/2. При проектировании используется весовое окно Хемминга. Для проектирования полосового фильтра надо задать Wn как вектор из двух частот [w1,w2], причем w1 < w2. Параметр ‘ftype’ задает тип фильтра: high (ФВЧ), stop (заграждающий). Вектор window размером n+1 позволяет задавать отсчеты окна.

Рисунок – графическое окно с графиками АЧХ и ФЧХ цифрового полосового фильтра, синтезированного с использованием функции b=fir1(32,[0.3 0.7])

 

Использование функции fir1 для проектирования цифрового полосового фильтра 32-ого порядка с граничными частотами 0.3 и 0.7 осуществляется следующим образом:

>> b=fir1(32,[0.3 0.7]);

>> freqz(b,1,512)

Полученное в результате выполнения второй команды графическое окно с графиками АЧХ и ФЧХ синтезированного цифрового фильтра имеет вид, приведенный на рисунке. Из рисунка видно, что ФЧХ фильтра в полосе пропускания действительно носит строго линейный характер.

 

3.1.2.3. Метод взвешивания для фильтра с произвольной АЧХ

 

Для проектирования фильтров с КИХ и произвольной частотной характеристикой методом взвешивания с использованием окон служит функция [3]:

b=fir2(n,f,m,npt,lap,window)

Эту функцию не рекомендуется применять для построения обычных фильтров: полосовых, режекторных, ФНЧ и ФВЧ.

Функция b=fir2(n,f,m) возвращает вектор с (n+1) коэффициентами ФНЧ порядка n. АЧХ фильтра задается векторами f (значения частот от 0 до 1) и m (отсчеты АЧХ для заданных в векторе f частот). Для ступенчатых АЧХ задаются два подряд отсчета частоты для ступеньки и два значения АЧХ. Параметр npt (по умолчанию 512) задает число отсчетов сетки для интерполяции, параметр lap (по умолчанию 25) задает число точек в области скачков АЧХ, а вектор window задает тип окна (по умолчанию окно Хемминга).

Использование функции fir2 для проектирования цифрового заграждающего фильтра 40-ого порядка с прямоугольной АЧХ и скачками на частотами 0.3 и 0.6 осуществляется следующим образом:

>> f=[0 0.3 0.3 0.6 0.6 1];

>> m=[1 1 0 0 1 1];

>> b=fir2(40,f,m);

>> [h,w]=freqz(b,1,128);

>> plot(f,m,W/pi,abs(h));

>> plot(f,m,w/pi,abs(h));

На рисунке показаны АЧХ идеального фильтра и фильтра, построенного с помощью функции fir2.

Рисунок – графическое окно с графиками АЧХ цифрового полосового фильтра, синтезированного с использованием функции b=fir2(40,f,m)

 

3.1.2.4. Метод наименьших квадратов

 

Метод наименьших квадратов в пакете Signal Processing Toolbox реализуют три функции:

- fircls;

- fircls1;

- firls.

Для проектирования фильтров с КИХ и произвольной АЧХ служит функция:

b=fircls(n,f,amp,up,lo)

Функция задает проектирование фильтра порядка n, заданного (n+1) коэффициентами вектора b с АЧХ, которая задается своими отсчетами частоты f и значениями amp. Векторы up и lo имеют то же число элементов, что и вектор amp, и задают верхние и нижние границы ЧХ для каждой из полос фильтра.

Для проектирования фильтров ФНЧ и ФВЧ с КИХ методами наименьших квадратов функция:

b=fircls1(n,W0,dp,ds,wt, ‘high’)

В простейшей форме b=fircls1(n,W0,dp,ds) функция создает вектор коэффициентов ФНЧ порядка n с частотой среза W0 между 0 и 1, где 1 соответствует половине частоты дискретизации, максимально допустимой функцией ошибки АЧХ в полосе пропускания dp и максимальной функцией ошибки в полосе задержания  ds. Параметр ‘high’ задает проектирование ФВЧ. Задается дополнительный параметр wt, выше которого или ниже которого выполняются требования к АЧХ проектируемого фильтра. Если wt в полосе пропускания, то , где  - функция ошибки. Если wt вне полосы, то пропускания . Векторы up и lo имеют то же число элементов, что и вектор amp, и задают верхние и нижние границы ЧХ для каждой из полос фильтра.

Использование функции fircls1 для проектирования цифрового ФВЧ 60-ого порядка имеет следующий вид:

>> n=60;

>> w0=0.4;

>> dp=0.03;

>> ds=0.01;

>> b=fircls1(n,w0,dp,ds,'high','plots')

Полученная при вычислении АЧХ приведена на рисунке. Из рисунка наглядно видно, что АЧХ укладывается в заданные нормы отклонений dp и ds.

Рисунок – графическое окно с графиками АЧХ цифрового ФВЧ, синтезированного с использованием  функции b=fircls1(n,w0,dp,ds,'high','plots')

 

Функция b=fircls(n,f,a,w, ‘ftype’) также использует метод наименьших квадратов для построения фильтра с КИХ и линейной ФЧХ для заданной АЧХ самой разнообразной формы. АЧХ задается векторами частот f  и значений нужной АЧХ a. В векторе w можно задать вес для аппроксимации АЧХ, причем число элементов вектора w равно половине числа элементов вектора f. Параметр ‘ftype’ позволяет задать тип фильтра: hilbert (фильтр Гильберта типов III и IV с линейной ФЧХ) и  differentiator (фильтр с АЧХ, линейно нарастающей с частотой).

Использование функции fircls для проектирования цифрового фильтра Гильберта 25-ого порядка с трапецеидальной АЧХ, отсчеты которой заданы в векторах частот f и значений a имеет следующий вид:

>> A=[0 0 1 1 0 0];

>> b=firls(25,F,A,'hilbert');

>> [H,f]=freqz(b,1,512,2);

>> plot(f,abs(H));

Полученная при вычислении АЧХ приведена на рисунке.

Рисунок – графическое окно с графиками АЧХ фильтра Гильберта, синтезированного с использованием  функции b=firls(25,F,A,'hilbert')