Быстрое рекурсивное вычисление одномерных и двумерных конечных сверток

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

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

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

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

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

h ( n ) =

b_n , при 0 <  n <  R

R

^ aih(n — i) при n > R ’ i=1

a i e R , i = 1,.., R называют коэффициентами ЛРС.

b i e R , i = 0,.., R — 1 называют начальными условиями ЛРС.

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

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

^R R - — ¹

x — a i x — ... — a R = о

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

Если все a 1 ,..., a _R - вещественные или комплексные корни характеристического уравнения (2) – имеют кратность единица, то общее решение ЛРС (1) записывается в виде:

h ( n ) = !№ , i = 1

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

^ciak = bk, k = 1,...R i=1

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

h ( n ) = c ₁ a n ' + ... + c j-1 a ” + c j ₂ n a ” + ... + k _nn

+ c jk , n j a j + ... + c_R a _R .

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

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

1) Многочлены степени R, удовлетворяющие ЛРС порядка R+1, характеристическое уравнение для которых представимо в виде (x — 1)R+1 = 0 и имеет один корень степени (R+1), равный единице. Коэффициенты рекурсии ai = CR+1 являются биномиальными коэффициентами, коэффициенты многочлена определяются начальными условиями.

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

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

Cos ( n ) = 2 Cos (1) Cos ( n - 1) - Cos ( n - 2) ^

h ( n ) = 2 Cos (1) h ( n - 1) - h ( n - 2).

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

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

Пусть h(n) - КИХ-фильтра, удовлетворяющая ЛРС порядка R с начальными условиями { b i } R = 0 ¹ и коэффициентами { a i } i ^R 1 , на отрезке [0,N-1]:

( x ( k )    '

x ( k + 1)

= A^k - ^R

^b о b i

_V x ( k + R - 1) J        ( b_{R - 1} ,

( 0

1 ... 0 A

где матрица

0    ..

...    0

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

V ^aR ... a 2 ^a 1 J провождающей матрицей ЛРС (3).

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

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

Если известны коэффициенты ЛРС и его R последовательных значений ( x ( k - R ), x ( k - R + 1), ...,

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

а ( - 1) ak

- a R - k , k = 1, R - 1, a_R

( - 1)- ¹

R =

aR

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

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

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

M + N - 1

y ( n ) = x ( n )* h ( n ) = X h ( k ) x ( n - k ) , (7) k = M

b n , при 0 <  n <  R ,

R h (n) =

0, при n< 0, n > N.

Подставляем (8) в (7).

N-1

y (n) = X h (k) x (n - k) = X h (k) x (n - k) + k=0

N-1 ( R           A

+ X I X a i h (k - i) Ix (n - k) = X h (k) x (n - k) + k=r V i=1            J

R ( N-1-i                   A R-1                   R

+X ai I X h (k) x (n - k - i) I=X h (k) x (n - k) +X aiy (n - i) -i=1 V k=R-i                     J k=0                    i=1

R (R-1-i                   A R ( N-1

-Xai I X h(k)x(n -k -i) I-Xai I X h(k)x(n-k- i) i=1 V k=0                     J i=1 V k=n -i

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

R-1                  R               R-1

У(n) = X h(k)x(n - k) +X aiy(n - i) -X x(n - k)I Xaih(k - i) k =0                   i=1               k=0

R-1              ( R

-X x(n - N -k)I X aih(N - (i -k)) I = Xaiy(n -i) - k=0                 V i= k+1

R-1         ( k           A R-1

-X x (n - k )I X aih (k - i) I-X x (n - N - k )I X aih (N + k - i)

k=0          V i=0           J k=0                V i= k+1

₍n (.A

Заметим, что значения d^l) =IXaih(k-i) |и V i=0

(-) (.A dk^ =I X aih(N + k-i) | можно рассчитать зара-

V i=k+1

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

R               R-1

У(n) = Xaiy(n - i) -X dkl)x(n - k) - X dkr)x(n - N - k) (9) i=1                k=0

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

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

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

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

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

(-1)                ¹

aR) = aR=aR, откуда aR=i, ak 1 = ak = -aR k, откуда ak = - aR _k, k = 1, R-1,

I                aR а для нечетной ИХ aV1 = - aR = Л"’ Откуда aR = -1, a a^1 = -ak = - aR-k , откуда ak = -aR-k, k = 1, R -1. aR

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

R dR-1-k =   Z   aih (N + (R - k -1) - i) = i=( R-1-k )+1

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

N-1

е²= Z 1^(m)-H(m)| ^min, m=0                              ai, bi

N -1

где H(m) =Z h(n) ш™^п , щ = ^exPf 2nj / N) - корень N-n=0

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

Обозначив a₀=-1, запишем выражение для

H(m), используя (11),(12)

H (m) = Z h (n )щ^тп = Z bn щ™ +1La_lk Z h (n -1 )щ™ n=0           i=1   ^n=R

R-1           R    NN-1-i               \

= Zbnщ™ +Zail Z h(nK(n+,)l = n=0           i=1   V n=R - i               J

kk

Z ^aR - ih⁽N^{-1 - (}k^-i)) =Z aih^{(k -}i) = dkl)

i=0                                  i=0

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

dR-1-k =   Z   aih(N + (R - k -1) - i) = i=( R-1-k )+1

k

k

R-1          R Л             N+R-1

= Z bnNmn +Z ailH(m)щ™ + Z h (n - iKn n=0           i =1 V                n=N

R                R-1 ( R           A

= ZaiЩ^miH(m) + Z | Zaih(n - i) I ш™ .

i=1                  n=0 V i=0            J

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

R-1

Z ^aR-ih⁽N^{- 1 - (k -}i)) = ^-Z aih^{(k -}i) = ^-dkl)

i=0                                      i=0

С учетом этих соотношений (9) переписывает-

ci =

Z^aj^h(i^{- j}) Vj=⁰

i=0

ся в следующем виде. Для четной ИХ:

[ R/2+1]

y⁽n) = Z a(y⁽n^-i)^-y⁽n^{- R+}i))⁺y⁽n^{- R})^-i=1

R-1

-Z dkl) (x (n - k) + x (n - N - R + k +1)), k=0

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

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

[ R/2+1]

y(n) = Z a(y(n - i)- y(n - R+i11 y(n - R)- i =1

Легко показать, что при известных ai они линейно связаны с начальными условиями bi. Для доказательства этого факта надо выразить присутствующие в (16) h(n), n > R, через a, и b,. Из представления отсчетов в матрично-векторном виде (5) видно, что h(n), n > R, линейно выражаются через начальные

R-1

-Z dk (x(n - k) - x(n - N - R + k +1)).

k=0

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

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

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

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

Нахождение оптимальных коэффициентов ЛРС, приближающего КИХ

Как было показано выше, для импульсной характеристики вида (8), удовлетворяющей ЛРС порядка R, вычислительная сложность нахождения значение выходного сигнала зависит только от порядка рекуррентности R. Поэтому возникает задача об оптимальном приближении произвольной импульсной характеристики идеального исходного фильтра {g(n)}N^¹ функцией вида (8).

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

R-1

Z щ™ сывается в форме: H(m) = i =0R------, а минимиза-

1 -Z a^-'

i=1

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

е2

N-1 R-1

= Z Z ci®' m=0 i=0

- G (m )fa,m mi i=1

Г(m) ^ min ^ai, ^ci

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

R

1-Z ai® i=1

„ N 1- е2 = Z(g(n)-h(n)) ^min.(13)

a.-, b n=0

Заметим, что значения {h(n)} нелинейно зависят от {ai}^1, что делает невозможным прямое решение задачи с помощью «приравнивания к нулю» частных производных по ai и bi.

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

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

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

Запишем явные выражения для нахождения ai и с_i в (17), приравняв нулю частные производные по ним. Нетрудно показать, что ds²= d aj

R (N -1       ,                 N -1       .Л

£ai I £ G(m)|2 Г(m)юm(i"j) + £ G(m)|2 Г(m)юm(j"i) | + i =1 V m=0                           m=0

R-1 (N-1                      N-1 _

+ £ ci I £ G(m)Г(m)rnm(i-j) + £ G(m)Г(m)®m(j-i) |, i =0 V m=0                        m=0

ds2 = d Cj

R (N-1                     N-1 _

£ ai I £ G(m)Г(m)®m(i-j) + £ G(m)Г(m)®m( j-i) | + i=1 V m=0                        m=0

R-1 (N-1                 N-1

£ ci I £ Г(m)rnm(i-j) + £ Г(m)rnm(j-i) |.(19)

i=0 V m=0                 m=0

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

N-1

Y(n) = £г(m)rn-mn, gY (n) = g(n) ®y(n), m=0

g 2y (n) = g (n) ® g (-n) ®Y( n),               (20)

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

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

' R                   R-1                                  ____

£ aig2y⁽i^-j) ⁺£ ^CigY^(j - i) = g2y^(j) ^j = ^{1, R}

< i=1                       i=0

R               R-1                           ________

£ aigY ⁽i^{- j}) ⁺ £ ^CiY^{( j -}i) = gY ^(j) ^j = 0,^{R - 1}

. i=1                     i =0

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

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

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

Г(m) ^ 1, y(n) = 8(n) -»дельта импульс».

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

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

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

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

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

• 6)    вычисляются новые значения весов Г(m) по (18) и y(n) по (20);

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

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

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

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

• 1)    Аппроксимация постоянным значением

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

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

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

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

• 3)    Аппроксимация многочленом

Многочлен степени (R-1) дискретного аргумента удовлетворяет ЛРС порядка R. Характеристическое уравнение его ЛРС (x -1)R = 0 имеет один корень степени R, равный единице. Коэффициенты рекурсии ai = (-1)iCR являются биномиальными коэффициентами и известны заранее, поэтому при минимизации (13) определению подлежат только параметры bi. Как было отмечено ранее, отсчеты h(n) для произвольного ЛРС при известных ai линейно выражаются через начальные условия h (n) = £ а i (n) bi i

bn, при 0 < n < R, б _n • b при n > R,

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

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

R

t⁽n) = y⁽n)^-2 ay⁽ⁿ - i) = i=1

R-1                  R-1

= -2 dkl)x(n - k) — 2 dkr)x(n - N - k) k=0                k=0

является просто R-той конечной разностью t (n) = A^R (n).

A0( n) = y (n), Aj (n) = A j-1( n)-Aj-1( n -1) j = 1,... R которые рассчитываются рекурсивно. Поэтому для получения y(n) по (9) на очередном шаге можно сначала вычислить t (n) = AR (n) =

R-1                   r-1                       ,

- 2 dkl) x (n - k) -2 dkr) x (n - N - k)

k=0                k=0

а затем, используя соотношение A^j-1(n) = A^j(n) + A^j-1(n -1) , вычислить «в обратном порядке» (см. рис. 1) все конечные разности A^j(n), j = 1,...,R -1, что требует R сложений и не требует умножений. При этом y(n) = A⁰(n).

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

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

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

R1 R 2

h⁽ⁿ1,ⁿ2 ) = 22 aijh⁽ⁿ1 ^-i’ n2 ^{- j}).

i =0 j =0

Более того, не удается реализовать этот алгоритм даже в предположении о разделимости аппроксимирующей функции h (n1, n 2) = h1(n1) ■ h 2( n 2) =

R1                 R 2

= 2 ^a1ih1⁽ⁿ1^-i) ^2 a 2 jh 2⁽ⁿ 2^{- j}).

i=0                 j=0

Однако если исходная ИХ разделима, и известно ее представление, g(n₁, n₂) = g(n₁) ■ g(n₂), то можно аппроксимировать каждую из одномерных ИХ gi (ni) функциями hi(ni) вида (1), удовлетворяющими ЛРС.

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

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

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

M_t+ N1-1 M 2 + N2-1

= 2      2  ^h(k1,^k2)^x(n1^{- k}1,ⁿ2^{- k}2)

k1 = M1    k 2 = M 2

	y(n-l)		y(n-R-l)	y(n-R)
	&‘(n-l)		^(n-R-1)
A"(n)	^'(n-l)
^(n)=t(n)

Рис. 1. Вычисление значения свертки с многочленом

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

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

R

t(n1,n2) = 2 a1it(n1 - i, n2) - i=1

R1 -1                           R1-1

• ^- 2 d 1k) x⁽ⁿ1 ^{- k}, n 2) ^- 2 dYk x⁽ⁿ1 ^{- N - k}, ⁿ 2)

k=0                      k=0

R 2

y (n1,n2) = 2 a2iy (n1,n2 - i) - i=1

R 2 -1                        R 2 -1

• ^-2 d 2 kt⁽n1, n 2^-k)^- 2 d 2 k) t⁽ⁿ1, n 2^-N^-k).

. k=0                       k=0

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

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

Рассмотрим вопрос о нахождении N₁+N₂ неизвестных отсчетов разделимой импульсной характеристики h₁ (n₁, n₂) = h₁(n₁)h₂(n₂), аппроксимирующей исходную импульсную характеристику g(n₁, n₂) :

^е2 = Z[g⁽ni,n2)-hi(ni)h2(n2)] ----->minX²⁷)

n₁,n₂

Взяв производные по   h1(j), j = 0,...N1 -1, и h2( j), j = 0,...N2 -1, получим систему уравнений:

де² "hj дее д h2( j)

_         N2 -1_ hi(j) ■ E2 - Z h2 (n2 ) ■ g( j, n2) = 0, n 2 =0

N-1_              _

Z hi⁽ni) ■ g⁽n1, j) ^-h2^{( j})■ Ei = ⁰

n =0

N 2-1 _       2        N'-1 -      2

где E2 = Z [^h2⁽ⁿ2)] , E1 = Z [^h1⁽ⁿ1)] . n 2 =0                          n_x =0

Перейдем к эквивалентной системе линейных уравнений. Умножив первые (N1-1) уравнений (28) на E1, а вторые (N2-1) уравнений на E2, получаем:

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

• ^hЛ = t ■ e Л t * ⁰ lie ЛI =¹.

Для каждого Л определим такое значение 1_Л, чтобы выполнялось второе условие (34): N1 -1             ₂N 2 -1              ₂

Z[h1(n1)]Z[h₂(n₂)] = Л, которое запишет-n1=0                 n 2 =0

• ^N1-7 -    q2 ^N 2-V - _П2

ся в виде Z [ Ле 1( n₁)] Z [ Л²( n 2)] = Л, откуда n1 =0                  n 2 =0

1/4

•

N2-1         _           _

Z g⁽j, n 2) ■ ^h 2 ⁽n 2 ) ■ E1 = h1⁽j) ■ E1E1

n 2 =0

N1 -1          _             _

Z g⁽ⁿ1, j) ■ ^h1^{( n}1) ■ E 2 = ^h 2 ⁽j) ■ E1E 2

n1=0

Л =±

Л

N1 -1 _       2 N 2-1 _        2

Z [ e 1⁽ⁿ1)] Z [ e2⁽n2)] n1 =0               n 2 =0

Заменим h₂(n₂) ■ E₁ и h₁(n₁) ■ E₂ в левой части уравнений (29) их выражениями из (28):

' N₂-1           < N₁-1 _                 '

Z g⁽j, n2) ■ Z ^h1⁽n1) ■ g⁽”1’n2)

. n =0            (n 0                     _y

N1-1           < N2-1 _

Z g⁽”1’j)■ Z ^h2⁽ⁿ2) ■ g⁽”1’ⁿ2) n1=0            ( n 2 =0

^h1⁽j ) ■ ^E1^E 2

^h2⁽j ) ■ ^E1^E 2

Очевидно, что соответствующие N1+N2 собственных векторов hл = 1_Лел будут являться решением задачи (33).

Замечание 1. Согласно (27) h₁(n₁) и h₂(n₂) будут определены с точностью до взаимного множителя. Рассчитанное значение t_Л соответствует

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

N-1 h (n1) ■ Z1 g (n, n 2) g (j, n 2) = Ehh1 (j), n1 =0          n 2 =0

N2 -1_

• Z ^h 2⁽n 2) ■Z g⁽ⁿ1, n 2) g⁽ⁿ1,^j ) = Ehh2^(j ), (32)

n 2 =0

Ni -1 _      2 N2 -1 _

Z[^h1⁽ⁿ1)] Z[h2 ⁽n2)] = Eh.

^ n1 =0               n 2 =0

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

значениям h₁(n₁) и h₂(n₂) c равной энергией (суммой квадратов отсчетов).

Замечание 2. Легко показать, что все ненулевые собственные числа матриц GG^T и G^TG совпадают. Поэтому можно ограничиться нахождением максимального собственного значения матрицы размера min(N₁,N₂).

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

KGh = Eh h

N1 -1 _      ₂N 2-1 _       ₂

Z [^h1⁽ⁿ1)] Z [^h2⁽n2)] = Eh ^, n1 =0               n 2 =0

где h = (h 1(0),...h 1(N₁-1),h2(0),...h2(N₂-1)), а кле-

точная матрица KG =

(GGT

I 0

GTG)

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

состоит из попарных скалярных произведений строк (верхняя клетка) и столбцов (нижняя клетка) исходной матрицы G = {g(n₁, n₂)} отсчетов аппроксими-

руемой двумерной функции.

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

^h(n1,ⁿ2) = Z⁸i^hT⁽ⁿ1⁾, ^h2i⁽n2),                 ⁽³⁶⁾

i определяет разложение произвольной ИХ в сумму разделимых ИХ, где 8, - собственные числа матриц GTG и GGT , а h1Ti(n1) – компонента с номером n1 соответствующего собственного вектора матрицы GTG, а h2i(n2) – компонента с номером n2 соответствующего собственного вектора матрицы GGT . Исходя из этих соображений легко показать, что в описанном выше алгоритме минимум ошибки аппроксимации будет достигаться при максимальном собственном значении. Для его нахождения удобно использовать, например, метод обратных итераций.

Полный алгоритм нахождения неизвестных отсчетов h₁(n₁) и h₂(n₂), выглядит следующим образом.

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

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

• 3)    Находится h2 = t₂e2.

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

R2

-I

Z Z g ⁿ1, ⁿ2)«1 j⁽ⁿ1^)a2i⁽ⁿ2) b2i = i =1    n., n 2 =0

дет определяться отношением

M

z

i=1,2* 2 ,

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

Пусть для исходной ИХ g(n₁,n₂) найдена по

(27) аппроксимация h₁ (n₁) ■ h₂(n₂), и для каждой из

h₁(n₁), h₂(n₂) найдены, согласно (9), оптимальные

наборы коэффициентов ЛРС {a₁i }i ^ и {a₂i}i 21. Не-

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

R₁             R₂

{b₁i }i 1 и {b₂i}i 1 - коэффициенты разложения

h₁(n₁) и h₂(n₂) по рекуррентному базису. Как было отмечено ранее, отсчеты h(n) для произвольного ЛРС при известных a_i линейно выражаются через начальные условия. Можно, конечно, выбирать начальные условия b₁= {b₁i} и b₂= {b₂i} как и в од-

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

£ 2 = Z[ g ( n1, n 2) - h1( n1) h 2( n 2)] = n1,n2

z

R1-1               R 2 -1

g⁽ⁿ1, n2) ^- Z ^a1i⁽n)^b1i ■ Z ^a2i⁽n)^b2i

n1,n2 L                 i=0                i=0

^ min ^b1i^,b2i

Взяв производные по b₁j, j = 0,...,R₁ -1

и

b2 j, j = 0,...,R2 -1, получим систему уравнений:

' . 2 N1-1, N 2 -1            _   _    __

— = Z ^[g⁽nl^{, n}2 ) - hl ⁽nl ) h (ⁿ2 )] ^j (nl ) h (ⁿ2 ) = ⁰

'b j n=0 L                         J

9g 2 Nн, N 2 -1          _   __

— = Z [ g⁽”1’ ⁿ 2) - ^h1⁽ⁿ1) h 2⁽n 2)1 «2 j⁽”1) ^h1(”1) = ⁰

^db1 j      n,, n 2 =0

С учетом обозначений

N2V   12

^E2 = Z ^[h2⁽n2)] , n 2 =0

N1 -1 _

E1 = Z [h1(n1)]   система (37) переписывается в n =0

виде.

N1-1, N2-1                     

Z g⁽”bⁿ 2^)a1j^{( n}1) h 2⁽n 2) = Z h1⁽n1^)a1j⁽ⁿ1) ■^E 2

n1, n 2 =0

' N]-1, N2-1                    _

Z g⁽ⁿ1, ⁿ 2^)a2 j⁽ⁿ1) ^h1⁽ⁿ1) = Z h 2( ⁿ 2^)a2 j⁽ⁿ 2) ■ ^E1

n1, n 2 =0                                       n 2 =0

Подставляя в предыдущее равенство выраже ние для h1(n1) и h2(n2) в виде (37), получаем:

R1

= Z

i=1

N1 -1

Z ^a1 j⁽ⁿ1^)a1i⁽ⁿ1)

n1 =0

^b1 i ■ ^E 2

R₁  N₁-1, N 2-1

Z Z g ⁿ1,ⁿ2^)a2 j⁽ⁿ2^)a1i⁽ⁿ1)

i =1 Щ, n 2 =0

R2

= Z

N2-1

Z ^a2 j⁽ⁿ2^)a2i⁽ⁿ2)

n 2 =0

^b2 i ■ ^E1

^b1i =

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

K₁b₂= A₁b₁E₂

K ₂b₁= A₂b₂E₁,

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

^K1 = ^{k1⁽i, j)^}, ^k1⁽i, j) = Z g⁽ⁿ1,ⁿ2 ^)a1 j⁽ⁿ1 ^)a2i⁽ⁿ2 ), n₁,n₂

K2 = {k2(i, j)}, k2(i, j) = Z g(n1,n2 )a2j(n2 )a1i(n1), n1,n2

A1 = {a1(i,j)}, a1(i,j) = Z a1 j(n1 )a1i(n1), n1

A2 = {a2 (i,j)}, a2 (i, j) = Z a2 j (n2 )a2i (n2 ), n2

^b1 = ^{b1i^{}, b}2 = ^{b2i^}.

Умножим первое уравнение (39) на E1, а второе – на E2, и выполним эквивалентные преобразования, используя выражения для b₂E₁ и b₁E₂ из (39).

( K₁b₂E₁= A₁b₁E₂E₁

( K ₂b₁E₂= A₂b₂E₂E₁

of A-1K1A21K2b1 = b1 E2 E1, v A21K 2 A-1K1b2 = b2 E2 E1.

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

' A-¹K1A2¹K2b1 = b1 E_h

A-¹K2A-¹K1b2 = b2Eh ,                   ⁽40)

₄Eh = (bt 1 ■ b1)(bt2 ■ b2)

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

KG =    ^{1 1 2 2}                    .

(      0        A2¹K2 A1-¹K1J

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

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

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

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

( е2 = Е g(ni,n2)-ЕhV)(ni)h2)(n2) I —>min,

In,n 2 V                    i                           )

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

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

• 1)    При наличии I звеньев. Фиксирование отсче

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

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

g(n1, n2) -Е h\1)(ni)h2i)(n2) разделимой им- i * j пульсной характеристикой общего вида h1(n1)■ h2(n2), где отсчеты h1(n1) и h2(n2) находятся по (33)-(35).

• 3)    Одномерные аппроксимации каждой из h₁(n₁),h₂(n₂) функциями вида (1) – h₁^(j) (n₁), h₂^(j) (n₂) , удовлетворяющими ЛРС соответствующего порядка R.

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

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

A(n1,n 2) = g (n1,n 2) -Е h1( i)(n1) h 2( i)( n 2 ) и i ошибки аппроксимации е2(А(n1,n2)). Если ошибка аппроксимации устраивает, то – конец алгоритма.

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

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

Заключение

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

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