Метод синтеза цифровых полиномиальных фильтров с помощью базисных частотных функций

УДК 681.5.015.52 + 004.932.4

М. А. Щербаков, С. В. Сорокин

МЕТОД СИНТЕЗА ЦИФРОВЫХ ПОЛИНОМИАЛЬНЫХ ФИЛЬТРОВ С ПОМОЩЬЮ БАЗИСНЫХ ЧАСТОТНЫХ ФУНКЦИЙ

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

Введение

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

Традиционным подходом к изучению нелинейных систем является рассмотрение одного или нескольких классов таких систем и разработка методов для их анализа и проектирования. Одним из таких классов нелинейных систем являются полиномиальные фильтры [1, 2], основанные на представлении выходного сигнала в виде функционального полинома Вольтерра. Для синтеза полиномиальных фильтров оказывается полезным их представление в частотной области, что позволяет рассматривать нелинейную фильтрацию как процесс преобразования гармонических составляющих входного сигнала во множество комбинационных составляющих выходного сигнала, определяемых частотными свойствами и степенью нелинейности фильтра.

1. Цифровая полиномиальная фильтрация

Цифровой полиномиальный фильтр порядка М определяется дискретным функциональным полиномом Вольтерра вида

М М х х т

у

т=0 т=0 и =—х пт =-х 1=1

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

Нт

Отдельные составляющие фильтра (1), определяемые сверткой вида

х х т

Ут

П1=—х пт =—х 1=1

называются однородными полиномиальными фильтрами порядка т .

При т = 1 выражение (3) представляет собой обычную линейную свертку, определяющую линейный дискретный фильтр с импульсной характеристикой ^(«). При т > 1 нелинейная свертка определяет однородный фильтр т -го порядка с ядром Нт(«1,..., пт).

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

Нт

п =—^ п =—^

/«1 /«т

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

Таким образом, физически реализуемый устойчивый полиномиальный фильтр порядка М можно представить в виде

М N—1 N—1 т

У

т=0 «1=0 «т =0 ¿=1

Многомерный полиномиальный фильтр определяется следующим

функциональным полиномом:

М

У

т=0

где х(П) и У(П) обозначают г -мерные входной и выходной сигналы,

П = [«1, «2, ..., «г]; Нт[х(П)] - однородный г -мерный дискретный функ-

ционал т -го порядка, определяющий выходной сигнал ут (П) однородного нелинейного фильтра т -го порядка

т

Ут

П1 Пт ;=1

где I обозначает г -кратное суммирование по всем элементам вектора П .

П

В дальнейшем однородный фильтр вида (8) будем называть гт -фильтром.

Следует отметить, что имеет место тесная взаимосвязь между многомерной линейной и полиномиальной фильтрацией [2]. Как известно [3, 4], линейная многомерная фильтрация сигнала и(щ, ..., ns) описывается многомерной линейной сверткой вида

У

;1=0 ¿8 =0

Допустим, что 5 -мерный входной сигнал является сепарабельной функцией, т.е. представим его в виде произведения т сигналов меньшей размерности г = 5 / т . Используя векторные аргументы, запишем

т

и

} =1

Для данного воздействия выражение (9) принимает вид

т

У

*1 ¡т У=1

Выделяя из выходного сигнала лишь диагональные блоки размерности г , т.е. полагая П1 = П2 =... = Пт = П , получим

т

Ут

*1 ¡т ^=1

Данное выражение есть не что иное, как нелинейная свертка, характеризующая гт -фильтр.

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

2. Описание многомерных полиномиальных фильтров в частотной области

Выходной сигнал гт -фильтра связан с выходом 5 -мерного линейного фильтра-прототипа соотношением

Ут

а входной сигнал представляет собой произведение вида

,

п=п, =...=п„

И(Пц, ..., П1Г, ..., Пт1, ..., птг) = П Х(ПД, ..., п]г).

]=1

Зависимость между входом и выходом 5 -мерной линейной системы определяется известным соотношением

У

где и(о>1, ..., ю5) и У(Ю1, ..., юх) представляют собой многомерные преобразования Фурье входного и выходного сигнала; Н(а^, ..., ю5)- частотная

п

п

т

характеристика линейной системы, представляющая собой преобразование Фурье импульсного откликаН(\, ..., is).

Для воздействия вида (14) соотношение (15) принимает вид

Y

т

= H

j=1

T

Используя векторные обозначения wj = [юд, ..., юjr] , можно запи-

сать

У

М

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

преобразованию (13). С этой целью выразим выходной сигнал (13) нелиней-

ного гт -фильтра через преобразование Фурье У(Ш1, ..., Шт) 5 -мерного линейного фильтра

1 п п

Ут

(2п) ' '

47 -п -п

Преобразование Фурье У(Ш1, ..., Шт) сигнала Ут (П) имеет вид

1 п п

Ут

(2п)гт ' '

n -п -п

• Т

где вектор ш = [о>1, ..., юг]Т . Используя свойство £е~]П ш = (2п)г8(ш), получаем

1 п п

Ут (“) =----— |... | У (“Ь - «т) 8(Ш - Ш1 - ... - Шт )й«1 ... й«т =

(2п) ' '

47 -п -п

1 п п

=-----1... 1 У

(2п) -п -п

Подставляя в данное выражение формулу (17), определяющую частотный отклик многомерного линейного фильтра, получаем

1 п п

Ут (ш) =----Гт^ 1... 1 Нт (ШЬ ..., Шт-1- « - «1 - ... - “т-1) Х

(2п) -п -п

XX

или в эквивалентной форме

тт тт Í ... \

где функция Нт (Ш1, ..., шт) связана с ядром Нт(П1, ...,Пт) во временной области следующим преобразованием Фурье:

Частотный отклик У(ш) многомерного полиномиального фильтра М-го порядка в частотной области при воздействии X (ш) принимает вид

М

где первое слагаемое характеризирует линейную составляющую фильтра.

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

При синтезе цифровых фильтров для обработки изображений очень часто требуется обеспечить условие изотропности оператора F фильтрации, состоящее в том, чтобы его реакция была инвариантна к ориентации входного изображения [5].

Условие изотропности оператора F можно записать в виде

где 0[P] - оператор изменения ориентации изображения P .

В задачах фильтрации изображений условие изотропности, как правило, сводится к инвариантности результата относительно вращения изображения на углы, кратные 90°, и его зеркального отражения относительно вертикальной оси. Обозначим данные преобразования соответственно через rota и ref , где а обозначает угол поворота. Заметим, что отражение относительно горизонтальной оси может быть получено последовательным применением операторов rot90 и ref . Повороты на углы 90°i, i = 2, 3 могут быть получены последовательным i -кратным применением оператора rot90 . Определим также обратные операторы rot- и ref. Нетрудно видеть, что ref= ref и

Y (ш) = £ Ym (ш) = ^(w) X (ш) +

3. Симметрия двухмерных изотропных фильтров

rota = rot360-a = rot-,

Преобразования вращения и отражения обладают следующими очевидными свойствами:

rota (A  B) = rota(A)  rota(B); rota (A o B) = rota (A) o rota (B);

ref

где через о обозначено произведение Адамара (поэлементное произведение матриц).

Как следствие (26), имеют место равенства:

с-1

A о ref (B) = ref (ref (A) o B) = ref (ref (A) o B);

A o rota (B) = rota (roí,-1 (A) o B) = rota (rot^Q-a (A) ° B).

Воспользуемся данными свойствами для определения условий, обеспечивающих инвариантность полиномиальных фильтров в задачах обработки изображений. Используя матричное представление полиномиальных фильтров [2], двухмерный 2т -фильтр можно представить в виде

У,,, (П) = £{H„ о X<")}.

эле-

где X(m) - кронекеровская степень m вектора X ; А} - сумма всех

ментов матрицы А.

Опорная область фильтра (28) на практике обычно выбирается квадратной и определяется в виде

^2 = {

Матрица Hm фильтра обладает свойством восьмиугольной симметрии

вида

Hm = rot90

Вследствие такого рода симметрии матрица Hm полностью определяется одним из своих восьми сегментов.

Для определения групп равных элементов в матрице Hm рассмотрим характер преобразования индексов опорной области (29) rm -фильтра при операциях вращения и отражения на примере опорной области фильтра размерности 3 х 3

^0 1 2'

= 3 4 5 .

6 7 8

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

.

"0 1 2 3 4 5 6 7 00 " 0 1 2 3 4 5 6 7 8Л

II 1 О) V2 5 8 1 4 7 0 3 6, ; Qref V 2 1 0 5 4 3 8 7 6,

Каждая подстановка характеризуется порядком, представляющим собой наименьшую степень, возведение в которую приводит ее к тождественной подстановке, играющей роль единицы в данных группах. Порядок подстановки |2га?| равен 4, порядок \0г^\ равен 2. Введенные подстановки допускают разложения на следующие независимые циклы:

Под циклами подразумеваются подстановки, в которых каждый элемент замещается следующим по порядку, а последний переходит в первый [6].

Пусть Бт =^2 х^2 Х".х^2 - декартова степень порядка т опорной области ^2 фильтра (маска фильтра). Назовем две точки (¿1, ¿2, ..., ¿т)е Ет и (л, ]2, ..., ]т)е °т эквивалентными: (¿1, ¿2, ..., ¿т)~(Л, к, ..., ]т), если

они могут быть преобразованы одна в другую путем перестановки индексов и выполнения перестановок 0го1 и 0ге^ . Например, точка (1,2) эквивалентна

точке (2,3), т.к. имеют место следующие преобразования:

0-го1 0Гв/ 0-го1 Рег

(1,2) — (3,6) — (1,4) — (3,2) — (2,3),

где Рег обозначает операцию перестановки индексов.

Очевидно, что введенное бинарное отношение ~ обладает всеми свойствам отношения эквивалентности: рефлексивностью, симметричностью и транзитивностью [7]. Поэтому множество Бт допускает разбиение на непе-ресекающиеся классы эквивалентности

Dm = ^тД и Вт,2 и ... и £>т,д .

Множество классов эквивалентности Бт1, ¿ = 1, ..., д образует так называемое фактор-множество относительно отношения ~ и обозначается Бт / ~ . Каждый класс Бт г- в формуле (34) полностью определяется одним

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

0-ГО1 0-го1 0-го1

—— а1 —— а2 —— а3;

й0 (І1, •••, Іт) '

С) о о о

^¿ге/ ^¿гоґ ^¿гоґ ^¿гоґ

—— а4 —— Й5 —— —— ау •

Выбирая, например, т = 2, а в качестве представителя элемент (0,1), в результате выполнения подстановок (35) получаем следующие восемь элементов:

Сгої Сгої Сгої

— (2,5) — (8,7) — (6,3);

Сге/ Сгої Сгої Сгої

— (2,1) — (5,8) — (7,6) — (0,3).

Если среди индексов элемента-представителя (¿1, i2, im) имеется s групп, содержащих по mk, к = 1, ..., s, равных между собой индексов, то с

учетом всех возможных перестановок cm...lm максимальное число элементов в классе составит

i Qrot ii Qref i cm-im = .

Таким образом, количество классов эквивалентности определяется числом возможных расположений m точек в опорной области ^2 таких, что ни одно из них не может быть преобразовано в другое путем применения операций rot и ref к ^. Такие комбинации точек являются представителями классов. Количество элементов в каждом классе будет равно s !K, где K - число возможных расположений s точек на ^2.

Элементы классов эквивалентности для фильтра второго порядка (m = 2) сведены в таблицу 1, где подчеркнутые элементы - представители классов.

Таблица 1

Разбиение коэффициентов квадратичного фильтра на классы эквивалентности

Номер класса Элементы класса эквивалентности D2 i

1

2

3

4

5

6

7

8

9

10

11

Таким образом, из таблицы видно, что для описания изотропного фильтра второго порядка потребуется лишь 11 коэффициентов, вместо 45.

4. Синтез двухмерных изотропных фильтров в частотной области

Рассмотренные свойства симметрии изотропных фильтров во временной области имеют соответствующие аналоги в частотной области, где опорную область ^2 удобно определить в виде

Г N -1 N -1 1

^2 = т, «2):------2“ “ П] ^ ] = 1,2 Г

Рассмотрим свойства нелинейной частотной характеристики двухмерного фильтра, определяемой выражением

Нт

= ... ^т (n11, n12,..., nm1, ит2) Х

(«11, «12)е^2 (ит1, «т2)е^2

хе~ 1

Учитывая свойства симметрии многомерной импульсной характеристики Нт(иц, «12,..., «т1, пт2) и используя векторные обозначения для арТ Т

гументов П1 = [«11, «12] , ш 1 = [ю 11, ю 12] , выражение (38) может быть записано в виде

Нт

Пт Х'^т /~

где (П1, ..., Пт)г- обозначает элемент-представитель г -го класса Бт г- эквивалентности, а Нт (П1, ..., Пт)г- - соответствующее ему значение импульсной характеристики.

Суммирование в (39) выполняется по всем представителям (П1, ..., Пт)г- классов эквивалентности, принадлежащим фактор-множеству Бт /~. Частотные функции ¥т г (Ш1, ..., шт) определяются для каждого класса Бт г- эквивалентности в виде

гт, г г

Vт,г (Ы1 ,•••, Шт ) = —2— {сте(иПю11 + ••• + ит1ют1)сте(и12ю11 + ••• + пт2ют2 ) +

+ с°5

где («11, «12, ..., «т1, пт2) = (П1, ..., Пт)г- и обозначает представителя г -го

класса Бт г-; {•} - операция симметризации относительно перестановки

пространственных частот Ш1, ..., шт.

Для квадратичного фильтра частотные функции (40) будут выглядеть следующим образом:

г2 ¿1

V 2,1 (Ш1, ш2) = [^(ЯцЮц + ^^^(я^Ю^ + ... + Я22Ю22) +

+ C0S(Иl2Юll + «22Ю21)с0$(Я11Ю12 + ... + Я21Ю22) +

+ с0^и21ю11 + и11ю21)с0^и22ю12 + ••• + и12ю22) +

+ с0^и22ю11 + и12ю21)с0^и21ю12 + ••• + и11ю22)].

Таким образом, каждому классу эквивалентности Бт г- может быть поставлена в соответствие базисная частотная функция ¥т г (ш1, •••, шт)• Эти

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

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

В качестве примера рассмотрим квадратичный фильтр с маской 3 х 3. Лексикографическое упорядочивание элементов данной маски показано на рисунке 1.

Рис. 1 Нумерация элементов опорной области Ж 2 размером 3 х 3

Найдем, например, базисную частотную функцию ¥2 7(ш1, шт) для

седьмого класса эквивалентности. Согласно данным таблицы 1 и нумерации на рисунке 1, представителем данного класса является элемент ^2 (1, 3) = ^2 (0, 1, -1, 0). Подставляя в выражение (41) значения пц = 0,

«12 = 1, П21 =-1, П22 = 0, |^2 71 = 8, получаем

¥2 7

Приближенное представление о частотных свойствах нелинейного фильтра можно получить на основе исследования его интегральной частотной характеристики для синусоидального воздействия. Для двухмерного случая синусоидальный сигнал с пространственной частотой Л = [А^, А2] определяется выражением

х

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

x

/2 2

где A = ^Aj +А2 , a = arctg(A1/А2).

Реакция двухмерного изотропного фильтра на синусоиду вида (44) остается неизменной для углов а± k90°, k = 1, 2, ... Поэтому при анализе частотных свойств фильтра достаточно исследовать его реакцию в диапазонах 0 <а<п/2 и 0 <А<л/2л.

П

Рис. 2 Пространственная синусоида

5. Синтез нелинейного оператора для повышения четкости изображения

Рассмотрим синтез изотропного нелинейного фильтра, обладающего свойством круговой симметрии, в классе квадратичных фильтров с маской (3 х 3), имеющих вид

у(п1, «2) = ¿2(4,4)х2(и1, «2) +

+¿2 (0,8) [х(«1 - 1, «2 - 1)х(«1 + 1, «2 + 1) + х(«1 -1, «2 + 1)х(«1 + 1, «2 - 1)] +

+¿2

Используя формулы (39)-(41) и нумерацию на рисунке 1, частотное ядро можно определить в виде трех базисных функций:

Н 2 (ю1, ю2) = 2 ¿2 (0,8) V 2,9 (ю1, ю2) +

+ "2 ¿2

Пусть фильтр должен обеспечивать высокочастотные свойства по отношению к постоянной составляющей и основной гармонике, полностью подавляя вторую гармонику. Так как коэффициент передачи по второй гармонике определяется сечением Н2(А,А), необходимо, чтобы Н2(А,А) = 0. На основании (46) это условие во временной области будет выглядеть следующим образом:

2Н2

Заметим, что из (47) также следует Н2 (0,0) = 0 . Поэтому частотный отклик фильтра У (ю) будет включать в себя постоянную составляющую

У0 = 2А2Н2(-А,А)8(ю) = 2А [4Л2(0,8)ф1(2А, а) + 2й2(1,7)ф2(2А, а)]8(ю), (48) и первую гармонику

У1(ю) = 2 А0 А1Н 2(0-, А)8(ю -А) =

= 2А0А1 [4й2

где функции фг- (X, a), i = 1,2 определяются выражениями

Фх

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

среднеквадратичная ошибка по постоянной составляющей }q определяется

следующим выражением:

йл/2

E = І I [2(-X,X)-R(X)]2dXda,

о о

где О - верхняя граничная частота рабочего диапазона фильтра.

Подставляя в (51) выражение для Н 2(-А, А) из (48) и приравниваю нулю частные производные по коэффициентам фильтра, получаем

йл/2

^(0,8) = І І Фі(2А, a)R(X)dX da

0 0 I

Qn/2

Й2(1,7) = І І Ф2(2X, a)R(X)dX da

0 0

Qn/2

2 1 1 Фі^, a)dX d a;

0 0 йп/2

0 0

Для обеспечения дифференцирующих свойств фильтра в диапазоне частот [0, п/2] положим R(Х) = |Х| иО = п/2. В результате численного интегрирования по формулам (52) были получены следующие коэффициенты фильтра: ^(4,4) = 2,3; ^(0,8) = -1,9; ^(1,7) = -0,4. На рисунке 3 приведен пример улучшения изображения за счет повышения четкости.

Рис. 3 Улучшение качества изображения с помощью нелинейного фильтра: а - исходное изображение; б - результат фильтрации

На рисунке 3,а представлено исходное изображение. Результат фильтрации с помощью нелинейного дифференцирующего фильтра (45) представлен на рисунке 3,б. Итоговое изображение на рисунке 3,б является более контрастным (различается, например, номер дома) и субъективно воспринимается лучше, чем первоначальное.

Выводы

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

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

1. Mathews, V. J. Polynomial signal processing / V. J. Mathews, G. L. Sicuranza. - A Wiley-Interscience Publication, 2000. - 445 p.

2. Щерб аков, М. А. Теория цифровой полиномиальной фильтрации и ее приложения : дисс. ... д-ра тех. наук / М. А. Щербаков. - Пенза, 1998. - 408 с.

3. Даджион, Д. Цифровая обработка многомерных сигналов / Д. Даджион, Р. Мерсеро. - М. : Мир, 1988. - 488 с.

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

5. Рампони, Дж. Расчет изотропных характеристик квадратичных фильтров методом двухимпульсной характеристики / Дж. Рампони // ТИИЭР. - 1990. - T. 78. -№ 4. - С. 96-108.

6. Фрид, Э. Элементарное введение в абстрактную алгебру / Э. Фрид. - М. : Мир, 1979. - 260 с.

7. Кострикин, А. И. Введение в алгебру. Основы алгебры : учебник для вузов / А. И. Кострикин. - М. : Физматлит, 1994. - 320 с.