Быстрое рекурсивное вычисление одномерных и двумерных конечных сверток Текст научной статьи по специальности «Математика»

БЫСТРОЕ РЕКУРСИВНОЕ ВЫЧИСЛЕНИЕ ОДНОМЕРНЫХ И ДВУМЕРНЫХ КОНЕЧНЫХ СВЕРТОК

А.В. Чернов

Самарский государственный аэрокосмический университет Аннотация

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

Введение

Обработка одно- и двумерных сигналов в режиме «скользящего окна» заключается в преобразовании дискретизированного сигнала линейной системой с постоянными параметрами (ЛПП-системой) с конечной импульсной характеристикой - КИХ-фильтром [1]. Значения сигнала на выходе КИХ-фильтра являются результатом цифровой свертки входного сигнала с импульсной характеристикой (ИХ) фильтра и могут быть найдены взвешенным суммированием входных отсчетов в пределах окна обработки. Однако такое вычисление свертки («прямая» реализация КИХ-фильтра) имеет практический смысл лишь для короткой импульсной характеристики, поскольку объем вычислений здесь пропорционален числу ненулевых отсчетов последней. Для больших окон (в задачах фильтрации и восстановления сигналов, вычисления признаков, корреляционного обнаружения и т.д.) прямое вычисление свертки оказывается чрезмерно трудоемким. В этой связи представляется целесообразным применение алгоритмов, воплощающих идею рекурсивной реализации КИХ-фильтров, развитую в работах [2, 3, 4, 5]. Основная идея этих алгоритмов -приближение импульсной характеристики ЛПП-системы семейством функций специального вида, вычисление свертки с которыми допускает рекурсивную реализацию с помощью разностных схем.

В данной работе предлагается обобщенный подход, позволяющий адаптивно подобрать базис разложения среди функций, удовлетворяющих ЛРС Я-того порядка, и построить наилучшее приближение заданной конечной импульсной характеристики. Как показано в [2], вычислительная сложность нахождения выходного отсчета для такой реализации зависит только от порядка рекуррентности Я и не зависит от размеров окна обработки N. Практически для Я<^ эту задачу можно рассматривать как задачу минимизации ошибки вычисления выходного сигнала при заданной пользователем вычислительной сложности обработки. Необходимо также заметить, что переход к разностным уравнениям - это единственный способ радикального снижения вычислительной сложности по сравнению со сложностью вычисления прямой свертки или ее реализации с помощью дискретных ортогональных преобразований.

Общие сведения из теории линейно рекуррентных соотношений

Приведем необходимые в дальнейшем общие сведения из теории линейно-рекуррентных соотношений [6, 7].

Линейно-рекуррентным соотношением порядка Я называется последовательность, удовлетворяющая соотношению:

Ьп, при 0 < п < Я

Ь(п) IV иг л > Я (1)

[V агп(п-г) при п> Я

I г=1

аг е Я, г = 1,..,Я называют коэффициентами ЛРС.

Ь е Я, г = 0,..,Я -1 называют начальными условиями ЛРС.

Последовательность, удовлетворяющую на заданном отрезке ЛРС (1), будем для краткости называть рекуррентной последовательностью.

ЛРС (1) полностью определяется совокупностью ее коэффициентов и начальных условий. Многочлен

хЯ - а1хЯ-1 -... - аЯ = 0 (2)

называется характеристическим уравнением ЛРС (1).

Если все а1,...,аЯ - вещественные или комплексные корни характеристического уравнения (2)

- имеют кратность единица, то общее решение ЛРС (1) записывается в виде:

и(п) = £егагп , (3)

І=1

где Сг определяются только начальными условиями:

V Сга = Ьк, к = 1,...Я

г =1

Если же среди корней а1,...,аЯ присутствует а у степени к> 1, то в сумму (3) он входит с учетом произведений на многочлены соответствующей степени в следующем виде:

И(п) = с1а1п +... + С:1 а'’ + с :2 п а" +... +

к (4)

+ с jkJn у а ” +... + Ся а Я.

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

Частными случаями ЛРС являются следующие семейства функций.

1) Многочлены степени Я, удовлетворяющие ЛРС порядка Я+1, характеристическое уравнение для которых представимо в виде

(х -1)Я+' = 0 и имеет один корень степени

(Я+1), равный единице. Коэффициенты рекурсии

аг = Сг + являются биномиальными коэффициентами, коэффициенты многочлена определяются начальными условиями.

Постоянные значения к(п)=Ь, удовлетворяющие ЛРС первого порядка к(п)=к(п-1) с начальным условием И(0)=Ь.

Тригонометрические функции (косинусы и синусы дискретного аргумента). Удовлетворяют ЛРС второго порядка. Например,

Cos(n) = 2Со^(1)Со5(п -1) - Св$(п - 2) ^

к(п) = 2Cos(1)k(n -1) - к(п - 2).

Иногда для удобства рассуждений используют представление отсчетов ЛРС в матрично-векторном виде

2)

3)

Г х(к) А Г Ь А Ь0

х( к +1) = Ак -Я Ь1 , к > Я , (5)

V х(к + Я - 1)у V ЬЯ-1 у

где матрица А =

Г 0 0 0

аЯ

... 0 А

.. 0

0 1

а2 а

называется со-

уля ... и.2 1 у

провождающей матрицей ЛРС (3).

Из общего вида решения линейно-

рекуррентных соотношений (3), (4) видно, что сумма двух ЛРС порядка Я1 и Я2 является ЛРС порядка не больше, чем (Я1+Я2).

Нетрудно показать [3], что общий вид 2-преобразования семейства функций, удовлетворяющих ЛРС, является дробно-рациональной функцией от г, и соответствующее семейство совпадает с множеством КИХ-фильтров, допускающих реализацию с помощью разностных схем с Я-звеньями.

Если известны коэффициенты ЛРС и его Я последовательных значений (х(к - Я), х(к - Я +1), ..., х(к -1)), то по этой информации можно восстановить не только «будущее» последовательности, но и прошлое. При этом коэффициенты ЛРС «обратного к данному» определяются следующим выражением:

-,(-!) = .

а

Я-к к = 1, Я -1, 1

а

Я

(5)

а(-1) =-аЯ

а

Я

Вычисление одномерной свертки с конечной импульсной характеристикой, удовлетворяющей ЛРС

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

Как известно [1], для ЛПП-систем с конечной импульсной характеристикой преобразование входного сигнала х(п) в выходную последовательность у(п) описывается соотношением «конечной свертки»:

М+N-1

у (п) = х(п)* к(п) = V к (к)х (п - к) , (7)

к =М

где к(п) - импульсная характеристика фильтра, равная нулю вне интервала [мМ + N -1], параметр М задает положение окна обработки относительно формируемого выходного отсчета, N - размер окна. Без ограничения общности дальнейших рассуждений можно считать М=0.

Пусть к(п) - КИХ-фильтра, удовлетворяющая

ЛРС порядка Я с начальными условиями {Ьг }Я_()1 и коэффициентами {аг , на отрезке [0^-1]:

Ьп, при 0 < п < Я, к(п) = <^ агк(п - г), при Я < п < N -1,

г=1

0, при п < 0, п > N. Подставляем (8) в (7).

N-1 Я-1

У (п) = к(к)х(п - к) = V к(к)х(п - к) +

к=0 к=0

N-1 Г я А я-1

+х| V а1к(к - г) |х(п - к) = V к(к)х(п - к) +

(8)

к=Я V 1=1

к=0

Заметим, что значения ) = |Х агк(к - г)

V 1=0

^к > =| V агк(N + к - г)| можно рассчитать зара

Я Г N-1-1 А Я-1 Я

+2 аг | V к(к)х(п - к - г) I = V к(к)х(п - к) +V ау(п - г) -г=1 V к= Я-1 У к=0 г=1

Я Г Я-1-1 А Я Г N-1 А

- V аг | V к(к)х(п - к - г) I- V аг | V к(к)х(п - к - г)

г=1 V к=0 у г=1 V к=N -1 у

Меняя порядок суммирования в третьей и четвертой суммах и введя обозначение а0= -1, получаем

У(п) = V к(к)х(п - к) + V агУ(п -г) - V х(п - к) | V агк(к -г) I к=0 г=1 к=0 V 1=1

-V х(п - N - к) | V агк(N - (г - к)) | = £ агУ(п -г) -к =0 V 1=к+1 У г=1

Я-1 Г к А Я-1 Г К

-Vх(п-к)IVак(к-г) |-^х(п-И-к)I V а-к^ + к-г) к=0 V 1=0 У к=0 V 1=к+1

V 1=к+1 у

нее, а для вычисления

у(п) = Vагу(п - г) -V ё(к‘)х(п - к) - V ^(кг)х(п - N - к) (9)

г=1 к=0 к=0

необходимо выполнить 3Я операций умножения и 3(Я-1)+2=3Я-1 операции сложения.

Соотношение (9) является основной расчетной формулой для вычисления линейной свертки входного сигнала с рекуррентной последовательностью.

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

В ряде задач фильтрации сигналов, импульсные характеристики фильтров обладают свойством центральной симметрии к(п) = к(N -1 - п) или асимметрии к(п) = -к(N -1 - п). Для функций данного класса удается дополнительно сократить количество умножений при вычислении (9).

Используя (8) и свойства четности/нечетности сигнала, получим, что для четной ИХ:

а,

4 1) = а4 = -1, откуда а4 = 1,

а для нечетной ИХ

„(-1) _

откудаак = -а4-к, к = 1,4 -1,

1

= -а4 = -—, откуда а4 = -1,

,(-1) _

= -ак = -

лК-к

откуда ак = -а4-к, к = 1,4 -1.

Получим связь между коэффициентами ёк ) и й?кг) в выражении (9). Исходя из тех же соображений для четной ИХ

4-1-к = V ак(N + (Я - к -1) - г) =

г=( Я-1-к )+1 к к

V ая-гк(N -1 - (к - г)) = V агк(к - г) = ё(к1) г=0 г=0

и, соответственно, для нечетной ИХ

/г) = у

иЯ-1-к _ ^

(10)

агк(N + (Я - к -1) - г) =

г=( Я-1-к )+1 к к

V аЯ-гк( N -1 - (к - г)) = - V а,к(к - г) = - $)

г=0 г=0

С учетом этих соотношений (9) переписывается в следующем виде. Для четной ИХ:

[ Я/2+1]

у(п) = V аг (у(п - г) - у(п - Я + г))+ у(п - Я) -

г =1

Я-1

- V ) ((п - к) + х(п - N - Я + к +1)+,

к=0

где [Я/2+1] обозначает целую часть числа.

Для нечетной ИХ:

[ Я/2+1]

у(п) = V аг (у(п - г) - у(п - Я + г))- у(п - Я) -

(11)

4-1

-у

к=0

(12)

п.—і

У dlk ((п - к) - х(п - N - 4 + к +1)).

Вычисление по (11) и (12) требует (34/2) умножений и (34) сложений на отсчет.

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

фильтра { (п)}^1 функцией вида (8).

Будем минимизировать квадрат отклонения

N-1 2

е2 = У((п) - к(п)) ^ тіп. (13)

п=о аЬ

Заметим, что значения {к(п)} нелинейно зависят от {аг- 4, что делает невозможным прямое

решение задачи с помощью «приравнивания к нулю» частных производных по а и Ьі.

Согласно равенству Парсеваля, задача минимизации (13) эквивалента среднеквадратичной минимизации отклонения дискретных спектров:

N-1

= У Р(ш) - Н(ш)|

ш=0 N-1

^ тіп,

а,ь

(14)

гдеН(т) =^ к(п) штп , щ = ехр(2л/ /^ - корень N-

п=0

ой степени из единицы, _Р= -1.

Обозначив а0=-1, запишем выражение для Н (т), используя (11),(12)

Я-1 Я (N-1 А

Н(т) = Vк(п)щтп = V ЬпШ”" + V«г I У к(п - гК" I =

п=0 г=1 V п=Я у

4-1

4

N-1-і

= У ЬпЩтп +Уаі І У М(п)Щт(п+і) 1 =

п=0 і=1 \ п=4-і

4-1 4 ( N+4-1

= У Ьпщтп + У аі І Н(ш)щті + У к(п - і)щтп | =

п=0 і=1 V n=N

= У аг-щтіН(т) +У і'^аік(п - і) 1 щтп . і=1 п=0 V і=0 /

Введем новые переменные

4-1

(15)

(

У а]к(і - ])

7=0

(16)

і=0

Легко показать, что при известных аі они линейно связаны с начальными условиями Ьі. Для доказательства этого факта надо выразить присутствующие в (16) к(п), п > 4, через аі и Ьі. Из представления отсчетов в матрично-векторном виде (5) видно, что к(п), п > 4, линейно выражаются через начальные

условия к(п) = бп • Ь, где бп - последняя строка (п-4+1)-ой степени сопровождающей матрицы.

Поэтому система уравнений (16) эквивалентна линейной относительно Ьі системе уравнений, которая легко решается.

В новых обозначениях равенство (15) перепи-

У с щті

сывается в форме: Н (ш) = —1^-------------

а минимиза-

1 -У ащ1

і=1

ция функционала (14) эквивалентна минимизации функционала

N-1

= У

4-1 4

У сют1 - с(т)У аюи

і=0 і=1

с «весами» Г(ш) =

1

Г(ш) ^ тіп (17)

а с

(18)

1-У аг ю”

г=1

Задача отыскания минимума функционала (17) по-прежнему остается нелинейной, ибо веса Г(т)

зависят от неизвестных коэффициентов а1.

Для получения неизвестных коэффициентов в

(17) можно использовать различные квазиоптималь-ные методы типа «аппроксимации Паде» [8], используемые для построения фильтров с бесконечной

2

4

е

і=1

2

2

е

импульсной характеристикой. В нашем случае КИХ-фильтра, задав начальные значения Г(т), разумно организовать итерационный вычислительный процесс попеременного нахождения а1 и с1 из соотношения (17) и вычисления Г(т) по (18).

Запишем явные выражения для нахождения а1 и с1 в (17), приравняв нулю частные производные по ним. Нетрудно показать, что

5е2 =

5 ау

я Г N-1 2 N-1 2

Vа11 V |О(т)| Г(т)ют(г-1) + У |О(т)| Г(т)ют(1-1) г=1 V т=0 т=0

Я-1 (N-1

N-1

+ У сг I V О(т)Г(т)ют(г 1 + V 0(т)Г(т)ю

г=0 V т=0 т=0

2

т( ]-г)

5е

5 СУ

Я / N-1

г=1 V т=0

Я-1 / N-1

Vаг | V О(т)Г(т)ют(г 1) +У О(т)Г(т)ют( 1 1) | +

(19)

V С I V Г(т)ют(г 1) + V Г(т)ю

1=0 V т=0 т=0

Введем обозначения

т( у-1 )

У(п) = V Г(т)ю тп, g1 (п) = я(п)®у(п)

т=0

§2У (п) = Я(п) ® Я (-п) ® у(п)

(20)

где ® обозначает циклическую свертку.

Вычисляя в явном виде суммы (19) через обратное ДПФ, окончательно получаем СЛУ порядка 2Я для нахождения а1 и с1.

Я Я-1 _____

V а1<?2у (1 - у) + V С1§у (1 - 1) = Я2у (1) 1 = 1, Я

1=1 1=0

Я

Я-1

(21)

V агЯу (г - Л + V СгУ(1' - г) = £у (Л 1' = 0, Я - 1 1=1 1=0

Соотношения (16), (20) и (21) являются основными расчетными формулами для нахождения неизвестных коэффициентов и начальных условий ЛРС.

Общий итерационный процесс нахождения параметров ЛРС выглядит следующим образом.

1) Задаются начальные значения

Г(т) = 1, у(п) = 5(п) -»дельта импульс».

На каждом шаге итерации:

2) вычисляются яу и я2у по формуле (20);

3) вычисляются а1 и с1 по формуле (21);

4) вычисляются Ь1 по формуле (16);

5) вычисляются значения отсчетов к( п) и ошибка аппроксимации е2 (17);

6) вычисляются новые значения весов Г(т) по

(18) и у(п) по (20);

7) если е2 удовлетворительна, или векторы параметров а1 и Ь1 «мало меняются» по сравнению с предыдущей итерацией, то выход, иначе переход к п.2 на следующий шаг итерации.

Замечание. Для аппроксимации четной или нечетной импульсной характеристики необходимо использовать данный алгоритм для аппроксимации «на половине интервала».

Вычисление одномерных сверток с рекуррентными последовательностями специального вида

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

1) Аппроксимация постоянным значением Постоянное значение к(п)=Ь удовлетворяет ЛРС первого порядка к(п)=к(п-1) с начальным условием к(0)=Ь. Оптимальное значение Ь вычисляется как «среднее значение» отсчетов к(п) и для реализации свертки с такой импульсной характеристикой необходима одна операция умножения и 3 сложения: у(п)=Ьу(п-1)+х(п)-х(п-М).

Аппроксимация произвольной импульсной характеристики семейством функций из прямоугольного базиса подробно рассмотрена в [2].

2) Аппроксимация тригонометрическими функциями

Тригонометрические функции (косинусы и синусы дискретного аргумента) удовлетворяют ЛРС второго порядка. Для отыскания возможного вида совокупности базисных функций для аппроксимации я(п) можно воспользоваться разложением я(п) по тригонометрическому базису (то есть нахождения О(т) с помощью ДПФ) и отбором симметричных трансформант (частот) с максимальной энергией. Количество частот выбирается в соответствии с желаемой сложностью реализации одномерной свертки. Оценки вычислительной сложности для косинусного базиса и Фурье-базиса подробно рассмотрены в [9], [2].

3) Аппроксимация многочленом Многочлен степени (Я-1) дискретного аргумента удовлетворяет ЛРС порядка Я. Характеристическое уравнение его ЛРС (х -1)я = 0 имеет один корень степени Я, равный единице. Коэффициенты рекурсии аг = (-1)1 Ся являются биномиальными коэффициентами и известны заранее, поэтому при минимизации (13) определению подлежат только параметры Ьг. Как было отмечено ранее, отсчеты И(п) для произвольного ЛРС при известных аг линейно выражаются через начальные условия

к(п) = Уаг (п)Ьг =

Ьп, при 0 < п < Я,

б п • Ь при п > Я,

(22)

где б п - последняя строка (п-Я+1)-ой степени сопровождающей матрицы. Для нахождения коэффициентов Ьг из условия минимума функционала (13) решается простая система линейных уравнений порядка Я.

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

что биномиальные коэффициенты аг = (-1)1 Ся

удовлетворяют свойству «треугольника Паскаля», и

Я

*(п) = у(п) -Vагу(п -1) =

г =1

Я-1 Я-1

= - V ёк)х(п - к) - V ёкг) х(п - N - к) к=0 к=0

является просто Я-той конечной разностью *(п) = ДЯ (п) .

Д°(п) = у(п), Д1 (п) = Д1-1 (п)-Д1-1(п-1) / = 1,...Я

которые рассчитываются рекурсивно. Поэтому для получения у(п) по (9) на очередном шаге можно сначала вычислить

Ґ(п) = Д4 (п) =

4-1 4-1 ,

-У dk )х(п - к) - У dkr')х(п - N - к)

к=0 к=0

(23)

а затем, используя соотношение

Д1-1 (п) = Д1 (п) + Д1-1 (п -1), вычислить «в обратном порядке» (см. рис. 1) все конечные разности

Д1 (п), / = 1,...,Я -1, что требует Я сложений и не

требует умножений. При этом у(п) = Д0 (п).

у(п) у(п-1) у(п-Я-1) у(п-Я)

А ‘(п) ге- А.'(п-1) А‘(п-К-1)

Ак,(п) А"'(п-1)

Ая(п)=і(п)

Рис. 1. Вычисление значения свертки с многочленом Окончательно для вычисления одномерной свертки с произвольным многочленом степени (Я-1) требуется и*(Я) = 2Я +1 операции умножения и и+ (Я) = 3Я + 3 операции сложения, что практически совпадает с оценками для параллельнорекурсивного вычисления свертки с многочленом, приведенные в [5]. Заметим, что рекурсивное вычисление свертки с многочленами общего вида является самостоятельной задачей и широко используется, например, для вычисления моментных инвариантов.

Для вычисления свертки с многочленами четных и нечетных степеней по (10) и (11) также требуется сокращенное число умножений порядка Я на отсчет.

Расчет рекурсивных двумерных КИХ-фильтров

К сожалению, не удается обобщить одномерный алгоритм на двумерный случай аппроксимации Я(п1,п2) двумерной разностной схемой общего вида

Я1 Я2

к(пъп2) = УУ а1ук(пх -1,п2 - /).

1 =0 1 =0

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

к(щ, п2) = кх(щ) • к2(п2) =

Я1 Я2 (24)

= V а11к1 (п1 - г) •V а21к2 (п2 - Л .

1=0 1=0

Однако если исходная ИХ разделима, и известно ее представление, я(п1, п2) = я (п1) • я (п2), то можно аппроксимировать каждую из одномерных ИХ (пг) функциями кг (пг) вида (1), удовлетво-

ряющими ЛРС.

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

Вычисление двумерной свертки с разделимой рекуррентной последовательностью

В предположении разделимости исходной ИХ, используя соотношение (24), удается «каскадно» реализовать вычисление конечной двумерной свертки

у (п1, п2) = х(п1, п2) **к(п1, п2) =

М1+N1 -1 М2 +N2-1 (25)

= V V к(кх,к2)х(п -кх,п2-к2)

к! =М1 к 2 =М 2

в следующей рекуррентной форме

Ґ (пі, щ) = У аиґ (щ - і, п2) -

і=1

-У dl(Jk) х(п1 - к, п2) - У d1(k) х(п1 - N - к, п2)

к=0

к=0

(26)

Ж , п2) = У а2і>’(п1, п2 - і) -

і =1

-У d2k)ґ (пь п2 - к) - £ d2k)ґ (пь п2 - N - к). к=0 к=0

Вычисление у(п],п2) по формуле (26) требует и*=3(41+42+1) умножений и U+=3(Я1+Я2+3) сложений, где 41 и 42 - порядки рекуррентности соответствующих одномерных ЛРС.

Аппроксимация неразделимой ИХ разделимыми ИХ общего вида

Рассмотрим вопрос о нахождении ^+^ неизвестных отсчетов разделимой импульсной характеристики к1 (п1, п2) = к1 (п1 )к2 (п2), аппроксимирующей исходную импульсную характеристику g(nl, п2) :

е2 = У [я(п1,п2)-к1(п1)к2(п2)] --------->тіп.(27)

п1,п2

Взяв производные по к1(І), І = 0,...N1 -1, и к2(}),} = 0,...^ -1, получим систему уравнений:

_ N2 -1_

МІ) • ^2 - У к^) • я(у, п2) = 0,

де

дМІ)

де2 дк2( І)

п2 =0

N1 -1_ _

У к1(п1) • я (п1, І) - к2 (І) • Е = 0

п1 =0 N2-1

(28)

N1-1.

(29)

где Е2 = V [к2 (п2) ] , Е1 = V [к1 (п1 )] .

п2=0 п1 =0

Перейдем к эквивалентной системе линейных уравнений. Умножив первые (Ы1-1) уравнений (28) на Е1, а вторые (N^1) уравнений на Е2, получаем:

У я(І, п2 ) • к2(п2) • Е1 = к1(У ) • Е1Е2

У Я (п1, І') • к1 (п1 ) • Е2 = к2 (І') • Е1Е2

п =0

(30)

Заменим к2 (п2) • Е1 и к1(п1) • Е2 в левой части уравнений (29) их выражениями из (28):

Чг-1 ^ N1 -1 _

У £(Л п2) • У к1 (п1) • ЯІЩ, п2)

п7 =0

N1 -1

V п1=0 ( N.-1

= М І) • Е1Е2

У Я(п1 =У) • У к2(п2) • Я(п1>п2)

пі=0 V п, =0

= к2 (І) • Е1Е2

Изменив порядок суммирования в левой части, и обозначив Ек= Е1Е2, окончательно получим систему уравнений для нахождения коэффициентов

У к1(п1) •У Я (пЬ п2 ) Я (У, п2 ) = Екк1( У X

п1 =0 п=0

N2 -1_ N1 -1 _

У к2(п2) •У Я (п1, п2 ) Я (п1, у) = Екк2 (І), (32)

п2 = 0

N1 -1

п1 =0 ^N2 -1

У[ к1(п1) ] У [ к2 (п2 ) ] = Ек.

п1 =0 п2 = 0

Перепишем (32) в матричном виде:

N1-1

КСЬ = Ек Ь

П2 %-1

У[ к1(п1)] У [к2(п2)] = Ек

(33)

где И = (кl(0),...кl(N1 -1),к2(0),...к2^2 -1)), а кле-

точная матрица Ке =

( ООт

0

0 ото

размера (N+N2)

состоит из попарных скалярных произведений строк (верхняя клетка) и столбцов (нижняя клетка) исходной матрицы О = {я (п1, п2)} отсчетов аппроксимируемой двумерной функции.

Система уравнений (33) по-прежнему не является линейной, так как зависит от И. Рассмотрим

вспомогательную задачу на нахождение собственных значений КеИ = ЛЬ. Существует лишь конечное множество собственных значений Л,...,Л , каждому из которых соответствует множество собственных векторов

ЬЛ =Ґ • eл,

ґ Ф 0 еЛ = 1.

Для каждого Л определим такое значение ҐЛ,

условие (34):

чтобы выполнялось

второе

2

У [к1 (п1 )] У[ к2 (п2 ) ] =Л, которое запишет-

п1 =0 п2 =0

N1 -1 _ 2 N2-1 _ 2

ся в виде У [ґдЄ1(п1)] У [ґдЄ2(п2)] = Л , откуда

п1 =0 (

,1/4

ҐЛ =±

N -1 _ 2 N-1 _ 2

У [ Ы^)] У [ Є2(п2)]

п2 =0

(35)

(31)

Очевидно, что соответствующие N^+N2 собственных векторов Ил = *лел будут являться решением задачи (33).

Замечание 1. Согласно (27) к1(п1) и к2(п2) будут определены с точностью до взаимного множителя. Рассчитанное значение *л соответствует значениям к1(п1) и к2(п2) с равной энергией (суммой квадратов отсчетов).

Замечание 2. Легко показать, что все ненулевые собственные числа матриц ООт и ОтО совпадают. Поэтому можно ограничиться нахождением максимального собственного значения матрицы размера тт^ьЩ.

Замечание 3. Нетрудно показать, что сингулярное разложение матрицы к(п1, п2)

к(пЪ п2 ) = У 8г (п1Х Ь2і (п2 ) ,

(36)

определяет разложение произвольной ИХ в сумму разделимых ИХ, где 5г- - собственные числа матриц

ОТО и ООт , а (п1) - компонента с номером п1

соответствующего собственного вектора матрицы ОТО , а Ъ21 (п2) - компонента с номером п2 соответствующего собственного вектора матрицы ООт . Исходя из этих соображений легко показать, что в описанном выше алгоритме минимум ошибки аппроксимации будет достигаться при максимальном собственном значении. Для его нахождения удобно использовать, например, метод обратных итераций.

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

1) Находится максимальное по модулю собственное число Л симметричной матрицы Ко и соответствующий ему нормированный собственный вектор ел.

п =0

N -1

2) Находится нормирующий коэффициент *л из условия (35).

3) Находится Ъл = *лел.

Качество аппроксимации, в общем случае, бул

дет определяться отношением ----М аХ------.

V л

1=1,Л’^Лшах

Аппроксимация двумерной неразделимой ИХ произведением разделимых ИХ, удовлетворяющих ЛРС

Пусть для исходной ИХ я(п1, п2) найдена по

(27) аппроксимация к1 (п1) • к2 (п2), и для каждой из

к1(п1), к2(п2) найдены, согласно (9), оптимальные

наборы коэффициентов ЛРС {а1г }= и {а21 } = . Нетрудно заметить (из общего решения ЛРС (2)), что эти коэффициенты определяют базисные функции, а

{Ь1г 11 и {Ь21 21 - коэффициенты разложения

к1(п1) и к2(п2) по рекуррентному базису. Как было отмечено ранее, отсчеты к(п) для произвольного ЛРС при известных аг линейно выражаются через начальные условия. Можно, конечно, выбирать начальные условия Ь1 = {Ь1г} и Ь2 = {Ь21} как и в одномерном случае, однако существует возможность подобрать их значения из условия минимизации отклонения двумерной функции

е2 = У [8(п^ п2) - к1(п1)к2(п2)]

-ч * 42—1

8(nl,п2) - У а (п)Ь1і • У а (п)Ь2

^ тіп

Ь1і ,Ь2і

Взяв производные по Ь11, і = 0,...,41 -1 и

Ь2і, і = 0,...,42 -1, получим систему уравнений:

д

= У [ё(п1> п2) - к1(п1)к2(п2)]а1і(п1)к2(п2) = 0

,2 N1 -1,^7-1

дЬ

= У [ 8 (п1> п2) - к1(п1)к2(п2)]а2і(п1)к1(п1) = 0

1І =0

.2 N1 -1, N2-1

(37)

дЬ.

1 і п ,п2=0

N2 -1

С учетом обозначений Е2 = У [к2(п2)]

п2 =0

N14 _ 2

Е1 = V [к1(п1)] система (37) переписывается в

п1 =0

виде.

N1 -1,N2-1 ____ N-1 -1___

V Я(п1, п2)а1 / (п1)к2(п2) = V к1(п1)а1 / (п1) • Е2

п1,п2 = 0 п1 = 0

N1 -1,N2-1 __ N2-1___

V Я(п1,п2)а2/(п1)к1(п1) = V к2(п2)а2/(п2) • Е1

п^п = 0 п2 =0

Подставляя в предыдущее равенство выражение для к1(п1) и к2(п2) в виде (37), получаем:

~ 2

У

і =1

=У

і=1

41

У

і=1

42

= У

У 8(n1, п2)а11 (п1 )а2і (п2)

т=0

N1 -1

У а1 і (п1 )а1і (п1)

Ь2і =

Ь1і • Е2

N1 -1, N2-1

У 8 (п1, п2)«2 І (п2)«1і (п1)

п1,п2 =0

Ь1і =

2

У а2І (п2 )а2і (п2 )

Ь2і • Е1

(39)

Перепишем (38) в матрично-векторном виде

Г К1Ь2 = А1Ь1Е2

V К 2Ь1 = А 2Ь2 Е1

используя обозначения

К = {к (і, І)}, к (і, і) = У 8(п, п )«1 у (п )аъ п),

п ,п2

К = {к2(і,І)}, к2(і,у) = У 8(пьп2)«2у (ща(п1),

п1,п2

А1 = {^1 (і, І)}, Й1 (і, І) = У «11 (п1 )«1і (п),

п1

А2 = {^2 (і, І)}, ^2 (і, І) = У«2 І (п2)«2г (п2),

п2

Ь1 = {Ь1і }, Ь2 = {Ь2і }.

Умножим первое уравнение (39) на Е1, а второе - на Е2і и выполним эквивалентные преобразования,

используя выражения для Ь2 Е1 и Ь1Е2 из (39).

Г К1Ь2Е1 = А1Ь1Е2Е1 Г К 2Ь1Е2 = А 2Ь2 Е2 Е1

( А-1К1А-1К 2Ь1 = Ь1Е2 Е1, А-1К 2 А-1К1Ь2 = Ь 2 Е2 Е1.

Окончательно получим систему нелинейных уравнений для нахождения неизвестных коэффициентов Ь1, Ь2

Г *-1

А-1К1А-1К 2 Ь1 = ЬЕ А-1 К 2 А-1К1Ь2 = Ь2 Ек

(40)

1 Ек = (Ь 1 • Ь1)(Ь2 • Ь2)

где Ек=Е1Е2. Для решения (40) можно использовать алгоритм, описанный в предыдущем параграфе (см. (33)-(35) ), использующий нахождение собственных чисел и векторов клеточной матрицы

Ко =

( А-К А^К2 0

0

л

А-1К2 А-К

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

удовлетворяющим ЛРС заданного порядка (многочленам, тригонометрическим функциям и пр.).

Аппроксимация неразделимой ИХ семейством разделимых ИХ, удовлетворяющих ЛРС

В общем случае при аппроксимации исходной ИХ семейством разделимых ИХ

е1 = VI" я(п1, п2) -V к )(п1)к21 )(п2) 1 -> ш1п,

п1,п2 V 1 )

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

В общем виде алгоритм выглядит следующим образом.

1) При наличии I звеньев. Фиксирование отсчетов всех звеньев, кроме одного (звена _0. Для выбранного звена выполняются шаги 2-5.

2) Аппроксимация функции

Я (п1, п2) - V к[') (п1 )к2) (п2) разделимой им-

1 * ]

пульсной характеристикой общего вида к1 (п1) • к2 (п2), где отсчеты к1 (п1) и к2 (п2) находятся по (33)-(35).

3) Одномерные аппроксимации каждой из к1(п1), к2(п2) функциями вида (1) -

к\1 )(п1), к21 )(п2) , удовлетворяющими ЛРС соответствующего порядка Я.

4) Поиск начальных условий Ь1г и Ь21 (коэффициенты аппроксимации базисными функциями) по (40).

5) Вычисление разностного сигнала

Д(п1, п2) = Я (п, п2) - V к1(г) (п1 )к2(г) п) и 1

ошибки аппроксимации е2(Д(п1,п2)). Если ошибка аппроксимации устраивает, то - конец алгоритма.

6) Если ошибка уменьшается слабо относительно предыдущих шагов, то переход к шагу 1 с количеством звеньев 1:=1+1 (при этом шаг 2 начинается с последнего звена с номером 1+1, а найденные на предыдущем шаге I звеньев

используются в качестве начальных приближений). Иначе - переход к шагу 2 с фиксированием всех звеньев, кроме (j+1) mod I для получения следующих элементов семейства

h°+>i), h2j+1)(«2).

Заключение

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

Благодарность

Работа выполнена при поддержке Министерства образования РФ, Администрации Самарской области и Американского фонда гражданских исследований и развития (CRDF Project SA-014-02) в рамках российско-американской программы «Фундаментальные исследования и высшее образование» (BRHE), а также гранта Президента РФ № НШ-1007.2003.01.

Литература

1. Оппенгейм А.В., Шафер Р.В. Цифровая обработка сигналов // М.: Мир, 1979. - 416с.

2. Сергеев В.В. Параллельно-рекурсивные КИХ-фильтры в задачах обработки изображений // Радиотехника. 1990. N8. С. 38-41.

3. Сергеев В.В. Параллельно-рекурсивные КИХ-фильтры для обработки изображений // Компьютерная оптика. М.: МЦНТИ, 1992. Вып.10-11. С. 186-201.

4. Ярославский Л.П. О возможности параллельной и рекурсивной организации цифровых фильтров // Радиотехника. 1984. N 3. С. 87-91.

5. Glumov N.I., Myasnikov V.V., Sergeyev V.V.: Polynomial Bases for Image Processing in a Sliding Window// Pattern Recognition and Image Analysis, 1994. No.4. Р. 408-413.

6. Холл Г. Комбинаторика // М., Мир, 1970.

7. Гельфонд А.О. Исчисление конечных разностей // 2-е изд.,доп. М.: Гос. изд-во физ.-мат. лит-ры,1959. 398 с.

8. Каппелини В., Константинидис А.Дж., Эми-лиани П. Цифровые фильтры и их применение // М.: Энергоатомиздат. 1983.

9. Методы компьютерной обработки изображений // Под ред В.А. Сойфера. М.: Физматлит, 2001. 784 с.