где r_k =Res K(p) – вычет передаточной функции аналогового прототипа в

p=p_k

полюсе p_k. Общее количество полюсов p_k равно N (предполагается, что все полюсы простые).

Вычеты r_k и полюсы p_k можно найти, используя следующий оператор MatLab:

>> [r, p, k]= residue(b, a)

- 9 -

Если не ставить в конце строки точку с запятой и нажать клавишу , то на экран монитора будут выведены вектор вычетов r, вектор полюсов p и вектор коэффициентов целой части k. Полюсы и вычеты могут быть

действительными или образовывать комплексно-сопряжённые пары. Масштабирующий коэффициент  подбирается таким, чтобы значение gц(0) было порядка единицы. При этом значения величин rk по модулю тоже порядка единицы. Если далее сумму дробей (4) привести к одной дроби, то можно получить передаточную функцию в стандартном виде (3). Полюсы цифрового фильтра в методе инвариантной импульсной характеристики связаны с полюсами аналогового прототипа стандартным z-преобразованием: zk=exp(pkT).

Синтез цифрового фильтра по методу инвариантной импульсной характеристики осуществляется в MatLab путём ввода оператора:

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

где bz, az – коэффициенты числителя и знаменателя передаточной функции цифрового фильтра (см. (3)); b, a – коэффициенты числителя и знаменателя передаточной функции аналогового прототипа (см. (1)); Fs – частота дискретизации (в Гц).

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

Одной и той же передаточной функции K(z) цифрового фильтра соответствуют различные формы реализации и разные алгоритмы преобразования отсчётов входного сигнала в отсчёты выходного. Выбор той или иной структуры цифрового фильтра имеет смысл при учёте эффектов конечной разрядности представления коэффициентов фильтра и отсчётов сигнала на входе, в промежуточных точках и на выходе. Дело в том, что дисперсия шума квантования на выходе устройства, чувствительность характеристик к точности представления коэффициентов, быстродействие, необходимый объём памяти и другие характеристики качества цифрового фильтра зависят от алгоритма работы и структуры. Естественно, что предпочтительнее та структура, которая обеспечивает наилучшие характеристики качества при заданных ограничениях. При реализации цифровых фильтров в виде специализированных вычислительных устройств имеет также значение, в какой последовательности выполняются операции, есть ли возможность производить несколько операций параллельно, что позволяет повысить быстродействие; каковы аппаратурные затраты при реализации того или иного алгоритма.

- 10 -

6.1. Прямая структура рекурсивного фильтра

Передаточной функции

соответствует структура вида

x[n] b0 y[n]

1/z 1/z

b1 - a1

x[n-1] y[n-1]

1/z 1/z

b2 -a2

x[n-2] y[n-2]

1/z bM -aN 1/z

x[n-M] y[n-N]

Алгоритм (разностное уравнение) записывается следующим образом:

y[n]= b0 x[n]+b1 x[n-1]+b2 x[n-2]++bM x[n-M]- a1 y[n-1]-a2 y[n-2]--aN y[n-N]

6.2. Каноническая структура рекурсивного фильтра

Каноническая структура получается из прямой путём разделения сумматора на две части (одна для прямых связей, другая – для обратных связей) с последующей перестановкой левой и правой частей схемы и дальнейшим слиянием параллельных цепочек элементов памяти в одну.

- 11 -

x [n] v[n] b0 y[n]

- a1 1/z b1

v[n-1]

1/z

- a2 v[n-2] b2

1/z

- aN v[n-N] bN

1/z

v[n-M] bM

На данной схеме показано, что М>N, однако это не обязательно; возможны случаи M=N или M

v[n]= x[n] - a1 v[n - 1] – a2 v[n - 2] - -aN v[n - N]

y[n]= b0 v[n]+b1 v[n - 1]+b2 v[n - 2]++bM v[n - M]

Сначала производится вычисление отсчёта сигнала v[n] в промежуточной точке (на выходе первого сумматора), а затем уже с его использованием – отсчёта выходного сигнала y[n]. Каноническая форма интересна тем, что в ней, в отличие от прямой структуры, представлена одна последовательность элементов памяти, а не две. Это позволяет экономить память. Однако абсолютные значения отсчётов промежуточного сигнала v[n] могут превосходить значения отсчётов входного и выходного сигналов, так что может потребоваться увеличенная разрядность ячеек памяти по сравнению с разрядностью регистров для ввода и вывода отсчётов входного и выходного сигналов соответственно.

- 12 -

6.3. Транспонированная структура рекурсивного фильтра

Преобразование прямой структуры, связанное с изменением порядка операций задержки и суммирования, приводит к транспонированной структуре.

x[n]

b0 b1 b2 bN bM-1 bM

+ 1/z + 1/z + 1/z + 1/z + 1/z

v1[n] v2[n] vN[n] vM-1[n] vM[n-1]

-a1 -a2 -aN

y[n]

Алгоритм для транспонированной структуры:

y[n]= b0 x[n]+v1[n-1]

v1[n]= b1 x[n] – a1 y[n]+v2[n-1]

v2[n]= b2 x[n] – a2 y[n]+v3[n-1]





vN[n]= bN x[n] – aN y[n]+vN+1[n-1]





vM-1[n]= bM-1 x[n]+vM[n-1]

vM[n]= bM x[n]

Разумеется, возможны случаи М=N или M

- 13 -

6.4. Каскадная (последовательная) структура

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

K(z) = K1(z) K2(z) KL(z).

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

Каскадная структура изображается следующим образом:

Алгоритм:

v[n]= b₁₀ x[n]+ b₁₁ x[n-1] – a₁₁ v[n-1]

y[n]= b₂₀ v[n]+ b₂₁ v[n-1]+ b₂₂ v[n-2] – a₂₁ y[n-1] – a₂₂ y[n-2]

В данной схеме каскады выполнены в виде прямых структур. Возможно их реализовать и в виде канонических форм:

- 14 -

Алгоритм в данном случае содержит три уравнения:

v[n] = x[n] – a₁₁ v[n-1]

w[n] = b₁₀ v[n]+ b₁₁v[n-1] – a₂₁w[n-1] – a₂₂ w[n-2]

y[n] = b₂₀ w[n]+ b₂₁ w[n-1]+ b₂₂ w[n-2]

Можно реализовать каскады и в виде транспонированных структур:

- 15 -

Теперь алгоритм включает в себя пять уравнений:

w[n]= b₁₀ x[n] +v₁[n-1]

v₁[n]= b₁₁ x[n] – a₁₁ w[n]

y[n]= b₂₀ w[n]+u₁[n-1]

u₁[n]= b₂₁ w[n] – a₂₁ y[n]+ u₂[n-1]

u₂[n]= b₂₂ w[n] – a₂₂ y[n]

Анализ показывает, что каскадная форма рекурсивного фильтра обладает меньшей чувствительностью частотной характеристики цифрового фильтра к точности представления коэффициентов фильтра, чем прямая, каноническая и транспонированная формы. Это означает, что при одинаковой разрядности округлённых коэффициентов b и a частотная характеристика фильтра в каскадной форме будет в меньшей степени отличаться от расчётной, чем характеристики некаскадных форм.

Для получения коэффициентов передаточных функций для каскадного представления фильтра можно использовать оператор MatLab преобразования фильтра в соединение секций второго порядка (sos – second-order sections):

>> [sos,g] = tf2sos (b, a);

где b,a – векторы-строки коэффициентов передаточной функции цифрового фильтра, sos – шестистолбцовая матрица, каждая строка которой соответствует одной секции и имеет структуру

[ b₀ b₁ b₂ 1 a₁ a₂ ],

что соответствует передаточной функции секции вида

В частном случае какая-либо из секций может иметь первый порядок; тогда соответствующие элементы строки матрицы sos будут нулевыми:

[ b₀ b₁ 0 1 a₁ 0]. Передаточная функция, соответствующая такой строке:

- 16 -

g – масштабный коэффициент, который реализуется в схеме в виде отдельного умножителя либо учитывается путём умножения на него коэффициентов b одного из каскадов.

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

Для вывода рассчитанной матрицы sos на монитор достаточно не ставить точку с запятой в конце оператора или набрать имя матрицы после значка >>.

Обратное преобразование матрицы sos в коэффициенты b и a фильтра осуществляется оператором:

>> [b, a]= sos2tf (sos, g);

6.5. Параллельная структура

Представим передаточную функцию K(z) в виде суммы дробей:

, (5)

где r₁ , r₂ , …r_N - вычеты, p₁, p₂ , … p_N – полюсы, k₀, k_1, …k_MN – константы.

Последние появляются в разложении, если MN, то есть прямых связей в структуре фильтра не меньше, чем обратных. Подобное разложение осуществляется в MatLab с использованием оператора residuez. После ввода или расчёта коэффициентов системной функции, представленных векторами b и a, нужно задать:

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

Если не ставить точку с запятой в конце строки и нажать клавишу , то будут выведены значения вектора вычетов r, вектора полюсов p и вектора коэффициентов k. При действительных b и a значения вычетов и полюсов могут образовывать комплексно-сопряжённые пары или быть действительными. Дроби с комплексно-сопряжёнными значениями вычетов и полюсов нужно объединить в одну дробь второго порядка. Ей будет соответствовать прямая, каноническая или транспонированная структура

- 17 -

второго порядка. Например, применение оператора residuez привело к следующему результату:

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

r=

-0.5000 - 0.1000i

-0.5000 + 0.1000i

2.1000

p=

0.4000 - 0.7000i

0.4000 + 0.7000i

0.6000

k=

1.5000

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

>> [b1, a1]= residuez (r(1:2), p(1:2), [ ] )

Здесь функция residuez в обратном направлении: объединяет две суммируемые дроби в одну и вычисляет коэффициенты полиномов числителя и знаменателя этой дроби. [ ] – символ пустой матрицы (не задаём коэффициент k). Будет выведено:

b1=

-1.0000 0.2600

a1=

1.0000 -0.8000 0.6500

что соответствует дроби

Покажем теперь, как изображается параллельная структура и как записать для неё алгоритм. Пусть, например,

- 18 -

Такому разложению соответствует схема

Алгоритм:

v[n]= b₁₀ x[n] – a₁₁ v[n-1]

w[n]= b₂₀ x[n]+b₂₁ x[n-1] – a₂₁ w[n-1] – a₂₂ w[n-2]

y[n]= k₀ x[n]+v[n]+w[n]

Разумеется, отдельные части схемы можно реализовывать в виде канонических или транспонированных структур (см. подраздел 6.4 ). Так же как и каскадная форма, параллельная обеспечивает меньшую чувствительность частотных и временных характеристик фильтра к точности представления коэффициентов по сравнению с прямой, канонической и транспонированной структурами. Это позволяет при реализации цифрового фильтра в виде специализированного вычислительного устройства обеспечивать заданные допуски на отклонение характеристик от расчётных при меньшем количестве двоичных разрядов, используемых для представления коэффициентов фильтра.

- 19 -

6.6. Нерекурсивный фильтр

Структура нерекурсивного не содержит обратных связей. Значит, все коэффициенты a_k равны нулю, кроме a₀=1. Передаточная функция такого фильтра

K(z)=b₀ + b₁z ^–1 + b₂z ^–2 +…+b_Mz ^–M

Схема:

Алгоритм цифровой фильтрации:

y[n]= b₀x[n]+b₁x[n-1]+b₂x[n-2]+…+b_Mx[n-M].

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

фильтра.

Для просмотра частотных и временных характеристик синтезированного цифрового фильтра используют оператор

>> fvtool (b, a)

После задания функции fvtool (filter visualization tool) нужно нажать клавишу . Аргументы функции – коэффициенты числителя и знаменателя передаточной функции, расположенные в порядке возрастания отрицательных степеней z (см. выражение (3)). Для просмотра характеристик нескольких фильтров одновременно нужно указать несколько пар векторов в списке входных параметров функции fvtool:

>> fvtool (b1, a1, b2, a2, b3, a3)

- 20 -

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

• амплитудно-частотную характеристику (АЧХ, magnitude response),

• фазочастотную характеристику (ФЧХ, phase response),

• АЧХ и ФЧХ совместно,

• xарактеристику группового времени запаздывания (group delay),

• импульсную характеристику (impulse response),

• переходную характеристику (step response),

• диаграмму полюсов и нулей (pole/zero plot),

• коэффициенты фильтра (filter coefficients).

Для детального просмотра участков графика используется кнопка zoom in панели инструментов. Щёлкнув по этой кнопке, нужно затем подвести курсор мыши к той точке графика, которая должна оказаться в центре увеличенного изображения и нажать левую кнопку мыши. Можно щёлкать левой кнопкой мыши многократно, что будет приводить ко всё большему увеличению. Возврат к исходному состоянию осуществляется щелчком по правой кнопке мыши, что соответствует нажатию кнопки zoom out панели инструментов.

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

Расчёт коэффициентов передаточной функции цифрового фильтра при заданных требованиях к частотной характеристике может быть произведён с использованием тех же самых функций, которые применяются для расчёта аналоговых фильтров (см. раздел 4). Отличие заключается в том, что строковый параметр ‘s’ вводить не нужно. Кроме того, все частоты (w0, wn, wp, ws) указываются в долях от частоты Найквиста (Fs/2), то есть лежат в интервале от 0 до 1. Есть и специальный пакет программ, где собраны многие из функций, рассмотренных нами по отдельности. Он носит название fdatool (filter design & analysis tool). Этот пакет содержит удобный пользовательский интерфейс и позволяет производить расчёт передаточной функции рекурсивных и нерекурсивных фильтров разнообразными методами синтеза, просматривать характеристики фильтра, анализировать изменение характеристик при квантовании коэффициентов фильтра, отсчётов входного сигнала и результатов промежуточных вычислений. Возможна работа с различными структурами фильтра. Вызов пакета осуществляется путём ввода его имени в командном окне MatLab: