Рекурсивные цифровые фильтры с импульсной характеристикой, описываемой полиномиальной функцией Текст научной статьи по специальности «Компьютерные и информационные науки»

Известия Тульского государственного университета Естественные науки. 2011. Вып. 2. С. 211-221

= ИНФОРМАТИКА =

УДК 681.5.015

Рекурсивные цифровые фильтры с импульсной характеристикой, описываемой полиномиальной функцией *

Д. А. Никитин

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

Ключевые слова: рекурсивный цифровой фильтр, импульсная характеристика, синтез ЦФ.

Введение

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

Рассмотрим, какие общеизвестные методы синтеза цифровых рекурсивных фильтров используют в качестве исходной информацию о значениях отсчётов импульсной характеристики рассчитываемого фильтра. Наиболее полная методика расчета цифровых рекурсивных фильтров (фильтров с бесконечной импульсной характеристикой, БИХ-фильтров) представлена в классической монографии [1]. Методы расчёта БИХ-фильтров можно поделить на следующие 3 вида:

• на основе уже имеющегося (рассчитанного) аналогового фильтра;

* Работа выполнена в рамках реализации ФЦП Научные кадры и научнопедагогические кадры инновационной России на 2009-2013 годы (контракт П1032 от 27.05.2010).

• без аналогового прототипа (прямые методы синтеза);

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

Во-первых, к интересуемым методам относится метод инвариантного преобразования импульсной характеристики, принадлежащий классу методов синтеза с использованием аналогового прототипа. Из-за названия может сложиться впечатление, что в качестве исходных данных используется ИХ аналогового фильтра. Однако на самом деле в этом методе для расчёта рекурсивного ЦФ используется передаточная функция аналогового прототипа. А такое название метод получил из-за того, что импульсная характеристика получаемого ЦФ представляет собой дискретизированную импульсную характеристику аналогового прототипа, а также ИХ используется при доказательстве корректности метода (только её общее комплексное выражение, а не значения отсчетов ИХ). Среди остальных методов с применением аналогового прототипа также в качестве исходных данных используется передаточная функция в s-плоскости.

Среди других методов синтеза — прямых и с использованием процедуры оптимизации — имеется 2 метода, более близких к интересующей задаче, то есть они используют для расчёта значения отсчетов ИХ. Так как импульсная характеристика БИХ-фильтра бесконечна, они на самом деле используют только некоторое количество её начальных отсчетов.

Среди прямых методов синтеза имеется метод, основанный на аппроксимации Паде. Данный вид аппроксимации предназначен для аппроксимации степенных рядов рациональной функцией. Предполагается, что значения отсчетов импульсной характеристики являются значениями степенного ряда, а передаточная функция рекурсивного ЦФ, как известно, является рациональной функцией. Таким образом, результатом данного метода являются коэффициенты передаточной функции. Как показано в [1, с. 295-297], результат данного метода может содержать значительные отклонения аппроксиманта, поэтому там же предлагается использовать коэффициенты фильтра, полученные таким методом, в качестве начальных значений при расчете БИХ-фильтров более сложными методами оптимизации.

Среди прямых методов синтеза с использованием процедуры оптимизации известен метод Прони [2]. Современный вариант метода Прони обобщен на модели, состоящие из затухающих синусоид, и модели, которые аппроксимируются экспоненциально-гармоническими суммами [3]. Метод является процедурой отыскания полюсов некоторого авторегрессионного процесса. Метод Прони позволяет находить фильтр, импульсная характеристика которого аппроксимирует заданную импульсную характеристику. Отдельные этапы метода не имеют аналитических методов

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

1. Основные результаты

В данной статье доказывается существование и единственность БИХ-фильтра, значения отсчетов импульсной характеристики которого являются значениями полинома (то есть постановка аналогична апрроксимации Паде: ИХ описывается степенным рядом, последние члены которого отброшены), а также попутно показывается, что коэффициенты передаточной функции можно найти точно (а не аппроксимируя) и без использования пошаговых приближающих методов.

Поскольку полином является расходящейся функцией, полученный ЦФ будет неустойчивым. Однако есть области, где такие фильтры могут быть использованы в целях, отличных от фильтрации [4-6]. Итак, ниже приводится доказательство существования и единственности БИХ-фильтра с импульсной характеристикой, которая описывается произвольной полиномиальной функцией, то есть отсчеты ИХ являются значениями такой функции, взятыми на равномерной сетке. Теорема даёт алгоритм синтеза таких фильтров. Также показана зависимость порядка фильтра от степени описывающего полинома и связь между коэффициентами фильтра и коэффициентами полинома.

Теорема 1. Для последовательности У длины Ь ^ 2(п + 1), члены которой являются равноотстоящими отсчетами полиномиальной функции одного переменного уг = рп(г) = а0 ■ гп + а\ ■ гп-1 + ... + ап, ао = 0, и заданы любые п + 1 подряд идущих члена ут, ут+\, ... ут+п, существует единственный рекурсивный цифровой фильтр порядка п + 1 (М = п + 1, N = п), импульсная характеристика которого совпадает с У.

Доказательство. Часть 1. Пусть т = 0, а отсчеты функции берутся с шагом равным единице. Формула цифровой фильтрации имеет вид

N М

Уп = ^2 ЬкХ(п - к) -^2 аку(п - к), (1)

к=0 к=1

где ак и Ьк — действительные коэффициенты. Большее из чисел М и N называется порядком фильтра. Примем N = М — 1 (таким образом, количество коэффициентов ак и Ьк будет равным) и построим начальную систему уравнений для М + N + 1 первых элементов последовательности У, приняв её за импульсную характеристику искомого ЦФ. Матрица системы имеет вид:

Р

1 0 0 . .0 0 0 .. .0 0

0 1 0 . .0 —Рп (°) 0 .. .0 0

0 0 1 . .0 —Рп(1) —Рп(0) .. .0 0

0 0 0 . .1 —Рп(п — 1) — Р ( п. 1 2) . —Рп(0) 0

0 0 0 . .0 —Рп(п) —Рп(п — 1) .. . —Рп(1) —Рп (0)

0 0 0 . .0 —Рп(п + 1) —Рп(п) .. . —Рп(2) —Рп(1)

0 0 0 . .0 —Рп(2п — 1) —Рп (2п — 2) .. . —Рп(п) —Рп(п — 1)

0 0 0 . .0 —Рп(2п) 1) —— п (2 (п Р — . —Рп(п + 1) —Рп(п)

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

Видно, что матрица Z является теплицевой. При приведении матрицы к треугольному виду будем производить над ней такие эквивалентные преобразования, чтобы из элементов, находящихся на одной диагонали вычитались одинаковые выражения, и таким образом матрица останется теплицевой.

Такими преобразованиями может быть, например, следующая последовательность действий (алгоритм приведения теплицевой матрицы специального вида к треугольному виду):

Будем нумеровать строки и столбцы с нуля.

1а) отнять п-й столбец из всех столбцов с нулевого по (п — 1)-й; отнять нулевую строку из всех строк с первой по п-ю; отнять (к + 1) раз (п — 1)-й столбец из всех столбцов с нулевого по (п — 2)-й, где к — разность индексов столбцов;

отнять (к + 1) раз первую строку из всех строк со второй по п-ю, где к — разность индексов строк; к + 2

отнять (к + 1)----- раз (п — 2)-й столбец из всех столбцов с нулевого

по (п — 3)-й, где к — разность индексов столбцов;

3б) отнять (к + 1) — раз вторую строку из всех строк с третьей по п-ю,

где к — разность индексов строк;

1б)

2а)

2б)

3а)

к + 2 к + 3

4а) отнять (к + 1)—2-3— раз (п — 3)-й столбец из всех столбцов с

23

нулевого по (п — 4)-й, где к — разность индексов столбцов;

4б) отнять (к + 1)

к + 2 к + 3

раз третью строку из всех строк с четвертой

23

по п-ю, где к — разность индексов строк; и так далее, п раз (до 1-го столбца и (п — 1)-й строки) Таким образом, получим

/ —Рп(п) —Рп(п — 1) ■ . . — Рп(1) —Рп(0) \

—Рп(п + 1) —Рп(п) •.. —Рп(2) —Рп(1)

Z =

—Рп(2п — 1) —Рп(2п — 2) . . . —Рп(п) —Рп(п — 1)

\ —Рп(2п) —Рп(2п — 1) . . . —Рп (п + 1) —Рп(п) /

( Хп(А) гп-і(А) ... ¿1(А) х0(А) \

()А о ... х2(а) х1(а)

, где А — ОД, &1-, ... ап

0 0 ... гп(А) ^п-і(А))

V 0 0 . . . 0 Хп(А) /

То же самое справедливо и для системы с большим количеством уравнений:

( Хп(А) Хп-і(А) ... х1(А) хо(А) \

0 Хп(А) ... х2 (а) х1(а)

Z — 0 0 ... Хп(А) Хп-1(А))

0 0 . . . 0 Хп(А)

0 0 . . . 0 0

0 0 . . . 0 0

То есть матрица является верхнетреугольной, так как она теплицева и х^о

— 0, і > 0.

Рассмотрим теперь расширенную матрицу такой же системы уравнен

Аналогично будем рассматривать только четвертый квадрант матрицы:

( —Рп(п) —Рп(п — 1) ... —Рп(1) —Рп (0) Рп(п +1) \

—Рп (п + 1) —Рп( п) . .. —Рп (2) —Рп(1) Рп(п + 2)

Z —

1) 1 п (п Р — п (2 (п Р — — 2) . . . —Рп(п) —Рп(п — 1) Рп(2п)

\ —Рп(2п) п (2 (п Р — — 1) . . . —Рп(п + 1) —Рп(п) Рп(2п + 1) у

Переместим сначала последний столбец на первое место и умножим его на минус один, чтобы матрица стала теплицевой, а после этого произведём

над ней точно такие же преобразования, что и над матрицей Z. Получим следующий результат:

Ъ

1 Рп (и + 1) Рп (и) —Рп (и — 1) . . . —Рп(1) —Рп(0) \

—Рп(и + 2) ~Рп(и + 1) —Рп(1 и) . . . —Рп (2) —Рп(1)

—Рп(2и) 1) 1 и 2 (п Р —Рп(2и — 2) . . . —Рп (и) Рп (и — 1)

\ —Рп(2и + 1) —Рп(2и) —Рп(2и — 1) . . . —Рп(и + 1) —Рп(и) /

( о Хп ( А) гп -1(А) . • • Х1(А) го(А) \

0 0 г п(А) . •• х2(а) х1(а)

0 0 0. •• гп(А) Хп-1(А))

0 0. .. 0 *п(А)

Действительно, так как результирующая матрица является теплицевой, значит все элементы Х^о =0, г ^ 0, как и все поддиагональные элементы матрицы Z.

То же самое справедливо и для системы с большим количеством уравнений:

0 Хп (А) Хп-1(А) . • Х1(А) хо(А) \

0 0 Хп(А) . • Х2 ( а) х1(а)

0 0 0 . • Хп(А) Хп-1(А))

0 0 0 . • 0 Хп(А)

0 0 0 . •0 0

0 0 0 . •0 0

То есть ранги матриц системы равны. Из вышеприведенного алгоритма нетрудно посчитать, что Хп(А) = —и!а0. Таким образом, система обладает полным рангом при ао = 0, что соответствует условию теоремы. Следовательно, она имеет единственное решение: коэффициенты цифрового фильтра порядка и + 1.

Пока что доказана справедливость теоремы только для случая, когда У является последовательностью отсчетов полинома степени и, взятыми с шагом равным единице. Покажем, что она справедлива и в случае любого постоянного шага:

Уг = ^ ау ■ (г ■ вгвр)п- = ^2а ■ ^ерпЧ) ■ (г)п- = ^

а ■гп-у ,

г = о,ь — 1,] = о, и, в1ер = о.

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

Часть 2. Пусть т > 0. Для любых и + 1 фиксированных точек существует единственный многочлен степени и, проходящий через них. Следовательно, зная у г, г = т,т + и, можно построить систему уравнений

^ ау ■ тп-3 = ут з

^ а3 ■ (т + 1)п-3 = Ут+1 < 3

^2 ау ■ (т + и)п-3 = ут+п^

, 3

Данная система имеет единственное решение относительно ау для фиксированного т > 0. Значит, если фиксированы уг, г = т,т + и, то фиксированы коэффициенты полинома ау. А если фиксированы коэффициенты полинома, то фиксированы уг, г = 0,и, так как они однозначно выражаются через коэффициенты полинома. Следовательно, доказательство в части 1 данной теоремы справедливо и при т > 0.

Следствие 1. Для синтеза рекурсивного ЦФ, импульсная характеристика которого описывается полиномом степени и, достаточно 2(и + 1) подряд идущих элементов ИХ. Количество коэффициентов и прямой и обратной связи получаемого фильтра при этом равно и + 1.

При попытке доказательства были рассмотрены также случаи, когда количества коэффициентов аи и Ъи не равна, то есть N = М — 1. Однако, если М < и + 1 или N < и, то решение не существует. А если М ^ и + 1 и N ^ и, то решение существует, но все последние коэффициенты аи или Ъk с порядковым номером свыше и (если нумеровать с нуля) оказываются равными нулю. То есть такое решение совпадает с доказанным выше.

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

Предположение 1. Коэффициенты обратной связи цифрового рекурсивного фильтра с импульсной характеристикой, описываемой полиномом и-ой степени, не зависят от коэффициентов полинома и совпадают с элементами (и + 1)-ой строки треугольника Паскаля с

Таблица 1

Зависимость коэффициентов Ьь от параметров описывающего полинома

Степень полинома Формула г-го члена ИХ Коэффициенты Ьк

1 Уі = ао • г + аі Ьо = ао + аі Ь1 = —а1

2 •2 Уі = ао • г2 + аі • г + а2 Ьо = ао + аі + а2 Ьі = ао — аі — 2а2 Ь2 = а2

3 Уі = Е ак • гк к=0 Ьо = ао + аі + а2 + аз Ьі = 4ао — 2а2 — 3а3 Ь2 = ао — аі + а2 + 3аз Ьз = —аз

4 Уі = Е ак • гк к=о Ьо = ао + аі + а2 + аз + а4 Ьі = 11ао + 3аі — а2 — 3аз — 4а4 Ь2 = 11ао — 3аі — а2 + 3аз + 6а4 Ьз = ао — аі + а2 — аз — 4а4 Ь4 = а4

отброшенным первым элементом, умноженными на (—1)к , где к — порядковый номер коэффициента

ак = (—1)к&п, к = 1, и + 1 (2)

Например, если элементы импульсной характеристики описываются полиномом 2-го порядка, то коэффициенты ОС соответствующего ЦФ равны -3, 3, -1. Для полинома 3-ей степени: -4, 6, -4, 1. Для полинома 4-й степени: -5, 10, -10, 5, -1. И так далее. Это свойство было использовано в других работах автора при разработке алгоритмов сжатия данных и анализа последовательностей в целях поиска участков данных, которые можно описать полиномом, и сокращения размера данных после обработки

[4-6].

Предположение 2. Коэффициенты обратной связи цифрового рекурсивного фильтра с импульсной характеристикой, описываемой полиномом п-ой степени, взаимно однозначно связаны с параметрами описывающего полинома. Можно из а^, і = 0, п, получить йк или Ьк, не производя громоздких вычислений, приведённых выше при доказательстве, и обратно, можно по коэффициентам синтезированного фильтра восстановить вид полинома, описывающего элементы его импульсной характеристики.

Это свойство также было использовано в [4-6].

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

Таблица 2

Зависимость параметров полинома от коэффициентов Ьь соответствующего

фильтра

Степень полинома Формула г-го члена ИХ Выражения для вычисления параметров

1 Уг = ао • г + а1 ао = Ьо + Ь1 а1 = —Ь1

2 Уг = ао • г2 + а1 • г + а2 Ьо + Ь1 + Ь2 ао = 2 Ьо — Ь1 — 3Ь2 а1 = 2 а2 = Ь2

3 Уг = Е ак • гк к=0 Ьо + Ь1 + Ь2 + Ьз ао = с Ьо - Ь2 - 2Ьз а1 = 2 2Ьо — &1 + 2Ь2 + 11 Ьз а2 = с 6 аз = -Ьз

4 Уг = Е ак • гк к=о Ьо + Ь1 + Ь2 + Ьз + Ь4 ао = 24 ЗЬо + Ь1 — 02 — 3Ьз — 5Ь4 11Ьо — Ь1 —^> + 11 Ьз + 3564 а2 = 24 3Ьо — Ь1 + Ь2 — 3Ьз — 25Ь4 аз = 12 а4 = Ь4

систему линейных уравнений, поэтому не составляет труда вычислить обратные зависимости — параметров полинома от коэффициентов фильтра (табл. 2).

Приведённая теорема расширяет доказательную базу для синтеза БИХ-фильтров только по импульсной характеристике, а кроме того даёт алгоритм для расчёта таких фильтров по импульсной характеристике, если известно, что её отсчеты точно подчиняются произвольной полиномиальной зависимости известной степени. Данный алгоритм в отличие от классических аналогов не является численным итерационным процессом поиска аппроксиманта и всегда даёт точное единственное решение. Для ясности повторим основные шаги алгоритма:

(1) Пусть отсчеты ИХ искомого БИХ-фильтра являются значениями

алгебраического полинома степени п, взятыми на равномерной сетке. Возьмём любые 2(п + 1) подряд идущих отсчета.

(2) Составим из них, используя уравнение (1), систему линейных алгебраических уравнений с 2(n + 1) неизвестными (M = n + 1, N = n).

(3) Далее решаем СЛАУ. Решать можно любым классическим способом, либо методом Гаусса с приведением к треугольному виду (прямой ход), используя описанный выше алгоритм приведения тёплицевой матрицы специального вида к треугольному виду.

(4) Решение системы содержит 2(n + 1) значений. Из них первые n + 1 — коэффициенты прямой связи фильтра (bk), а оставшиеся n + 1 — коэффициенты обратной связи (ak).

Полученный ЦФ будет неустойчивым, так как полином — расходящаяся функция. Тем не менее, есть возможные приложения, не связанные с использованием полученного фильтра для фильтрации. Это анализ данных, сжатие, и другие, описанные в [6]. Кроме полиномиальных функций подобные теоремы доказаны автором для синусоидальных и показательных функций, а также для случая, когда отсчеты импульсной характеристики являются значениями любой периодической функции, взятыми с постоянным шагом при условии, что период функции кратен расстоянию между взятыми отсчетами. Формулировки теорем, касающихся этих видов описывающих функций, можно найти в [7].

Список литературы

1. Рабинер Л., Гоулд Б. Теория и применение цифровой обработки сигналов. М.:Мир, 1978. 848 с.

2. Марпл мл. С.Л. Цифровой спектральный анализ и его приложения. М.: Мир, 1990. 547 с.

3. Завьялов М.Н., Елизаров Д.В. Комбинированный вариант гармонического разложения Прони // Журн. СФУ. Сер. Матем. и физ. 2008. №1(4). C.443-452.

4. Ханов В.Х., Никитин Д.А. Алгоритм анализа числовых последовательностей // Вестник СибГАУ. 2006. №6(13). С.11-15.

5. Никитин Д.А. Сжатие с использованием поиска функциональных зависимостей // Технологии Microsoft в теории и практике программирования: конф.-конкурс. НГУ. Новосибирск, 2007. С.133-135.

6. Никитин Д.А. Приложения алгоритма синтеза рекурсивных цифровых фильтров по импульсной характеристике // Цифровая обработка сигналов. 2009. №4. C.8-15.

7. Никитин Д.А. Теоремы о существовании и порядках цифровых рекурсивных фильтров с импульсными характеристиками определённой формы // Информационные технологии и математическое моделирование (ИТММ-2009): матер. VIII Всероссийской научно-практической конф. с междун. участием.Томск: ТГУ, 2009. Ч.2. С.144-146.

Никитин Дмитрий Александрович ([email protected]), старший преподаватель, кафедра безопасности информационных технологий, Сибирский государственный аэрокосмический университет им. акад. М.Ф. Решетнёва, Красноярск.

The recursive digital filters with impulse response defined by

polynomial function

D. A. Nikitin

Abstract. In this paper proved the existence and uniqueness of the digital recursive filter with impulse response, which is described by an arbitrary polynomial function. The theorem also provides a synthesis algorithm of such filters not using the approximating step-by-step methods. The nature of relationship between the coefficients describing a polynomial and coefficients of corresponding DF is described.

Keywords : recursive digital filter, impulse response, DF synthesis.

Nikitin Dmitry ([email protected]), senior lecturer, department of information technologies security, Siberian State Aerospace University named after academician M.F. Reshetnev, Krasnoyarsk.

Поступила 07.06.2011