Дискретная фильтрация в MATLAB

Дискретная свертка– основа алгоритма дискретной фильтрации. В MATLAB вычисляется с помощью функции conv.

z=conv(x,y)

Обращение свертки – зная результат свертки и один из сворачиваемых векторов, можно найти второй.

[q,r]=deconv(b,a)

Здесь b – результат свертки, a – один из сворачиваемых векторов, q – искомый второй вектор свертки, r – остаток (если вектор b действительно является сверткой вектора a с чем-нибудь, остаток будет нулевым).

Функция дискретной фильтрации.Основная функция, реализующая дискретную фильтрацию в MATLAB, носит имя filter. В простейшем виде имеет следующий синтаксис:

y=filter(b,a,x)

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

Если первый элемент вектора а не равен 1, значения векторов b и а нормируются – делятся на а(1).

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

Расчет импульсной характеристики–выполняется функцией impz.

h=impz(b, a, n)

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

Расчет частотной характеристики – выполняется функцией freqz.

freqz (b,a,n)

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

Расчет групповой задержки дискретной системы –выполняется функцией grpdelay.

grpdelay(b,a)

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

Отображение нулей и полюсов фильтра.Для отображения нулей и полюсов функции передачи фильтра на комплексной плоскости предназначена функция zplane.

zplane(z, p)

zplane(b,a)

Входными параметрами являются векторы-столбцы нулей и полюсов (z,p) либо векторы-строки коэффициентов полиномов числителя и знаменателя функции передачи (b,a). Различение этих двух случаев производится именно по ориентации передаваемых функции векторов.

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

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

tf2zp – преобразование коэффициентов полиномов функции передачи в нули и полюсы функции передачи;

zp2tf — преобразование нулей, полюсов и коэффициента усиления системы в коэффициенты полиномов функции передачи;

tf2ss — преобразование коэффициентов полиномов функции передачи в параметры пространства состояний;

ss2tf — преобразование параметров пространства состояний в коэффициенты полиномов функции передачи;

ss2zp — преобразование параметров пространства состояний в нули и полюсы функции передачи;

zp2ss — преобразование нулей, полюсов и коэффициента усиления системы в параметры пространства состояний;

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

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

[A, B, C, D]=tf2ss(b,a);

[b, a]=ss2tf(A, B, C, D);

[z, p, k]=ss2zp(A, B, C, D);

[A, B, C, D]=zp2ss(z, p, k);

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

[b1, a1]= eqtflength(b,a);

    Что касается разложения на простые дроби, здесь преобразование дискретной системы несколько отличается от аналогового случая, поэтому предусмотрена специальная функция reziduez. Данная функция умеет выполнять преобразование в обе стороны. Требуемый тип преобразования определяется количеством выходных параметров:

[r, p, k]=residuez(b, a)

[b, a]=residuez(r, p, k)

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

Разработка фильтра с помощью FDATool

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

» fdatool

После этого появится окно программы, показанное на рис. 3.

Рис. 3. Интерфейс программы fdatool

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

Расчет фильтра

Расчет фильтра начинается с задания требуемых параметров на вкладке Design Filter (см. рис. 3). Тип синтезируемой АЧХ выбирается с помощью переключателя Filter Туре. Возможны следующие варианты: Lowpass (ФНЧ), Highpass (ФВЧ), Bandpass (полосовой фильтр), Bandstop (режекторный фильтр). Выбор пятого, самого нижнего, положения переключателя позволяет использовать раскрывающийся список, в котором перечислены более сложные варианты: Differentiator (дифференцирующий фильтр), Hilbert Transformer (преобразователь Гильберта), Multiband (многополосный фильтр), Arbitrary Magnitude (произвольная АЧХ) и Arbitrary Group Delay (произвольная групповая задержка).

Выбрав категорию синтезируемой АЧХ, следует выбрать тип синтезируемого фильтра, установив переключатель, расположенный в разделе Design Method, в положение IIR (рекурсивный) или FIR (нерекурсивный). Каждому положению переключателя соответствует список возможных методов синтеза. Состав этого списка меняется в зависимости от выбранного типа АЧХ. Например, при синтезе фильтра с произвольной зависимостью групповой задержки от частоты (Arbitrary Group Delay) переключатель автоматически установится в положение IIR, а в списке будет доступен всего один метод – метод минимизации p-нормы ошибки (Constrained Least Pth Norm). В случае синтеза АЧХ четырех простейших типов набор возможных методов синтеза значительно шире:

□ нерекурсивные фильтры (FIR). Здесь доступны следующие методы:

· Equiripple – синтез фильтров с равномерными пульсациями АЧХ методом Ремеза;

· Least–Squares – минимизация среднеквадратического отклонения АЧХ от заданной;

· Window – синтез с использованием весовых функций (окон);

□ рекурсивные фильтры (IIR). Здесь доступны четыре варианта синтеза по различным аналоговым прототипам методом билинейного z-преобразования:

· Butterworth – синтез фильтра Баттерворта;

· Chebyshev Type I – синтез фильтра Чебышева первого рода;

· Chebyshev Type II – синтез фильтра Чебышева второго рода; Elliptic – синтез эллиптического фильтра.

Замечание      

Набор доступных методов синтеза зависит от того, установлен ли пакет Filter Design.

Далее необходимо выбрать порядок фильтра в разделе Filter Order. Во многих случаях помимо явного указания порядка в поле ввода Specify order возможен автоматический выбор порядка путем установки переключателя в положение Minimum order.

Наконец, необходимо задать числовые параметры этой АЧХ в разделах Frequency Specifications и Magnitude Specifications (при выборе типа АЧХ из дополнительного раскрывающегося списка эти два раздела объединяются под общим названием Frequency and Magnitude Specifications). Содержимое этих областей окна меняется в зависимости от выбранного типа АЧХ. На рис. 3 показаны поля ввода параметров, соответствующие синтезу ФНЧ. При заполнении полей ввода можно ссылаться на переменные, существующие в данный момент в рабочей области памяти MATLAB.

Числовые параметры, которые необходимо задать, иллюстрируются графиком, выводимым в разделе Filter Specifications. Вид этого графика также меняется в зависимости от выбранного типа АЧХ, на рис.3 он соответствует синтезу ФНЧ.

Выбрав метод синтеза и задав характеристики фильтра, щелкните на расположенной в нижней части окна кнопке Design Filter. Программа fdatool вызовет нужную функцию синтеза, передав ей указанные вами спецификации фильтра.

Замечание      

Если установлен пакет Filter Design, в программе fdatool можно установить флажок Turn quantization on, чтобы сделать доступной вкладку Set Quantization Parameters, позволяю­щую задавать параметры квантования коэффициентов фильтра.