A family of computationally efficient algorithms for multichannel signal processing—A tutorial review

Signal Processing 5 (1983) 5-19 North-Holland Publishing Company

5

A FAMILY OF COMPUTATIONALLY EFFICIENT ALGORITHMS FOR MULTICHANNEL SIGNAL PROCESSING- A TUTORIAL REVIEW N. K A L O U P T S I D I S

Division of Electronics, University of Athens, 104 Solonos St., Athens, Greece D. MANOLAKIS

Division of Electronics, University of Athens and Ionospheric Institute, National Observatory of Athens, Theseum, Athens (306), Greece G. CARAYANNIS

Division of Electronics, National Technical University of Athens, 42 Patission Ave., Athens, Greece Received 15 July 1981 Revised 8 February 1982 and 7 June 1982

Abstract. This paper is concerned with the solution of block linear systems arising in multichannel digital signal processing. First the general problem of block matrix inversion and linear system solution is considered and a corresponding algorithm, recursive in nature, is developed, together with a block form of triangularization theorem. Subsequently these general schemes are specialized to block Toeplitz structures yielding most existing efficient algorithms. To this respect, matrices related to multichannel Wiener fltering, AR and ARMA modelling as well as s + 1 steps ahead prediction are examined and related algorithms are described. Finally block banded Toeplitz systems are considered and two new efficient algorithms are presented, one recursive in nature and the other FFT based. They constitute natural extensions of methods already available for the single channel case and their derivation is simple due to the unified approach introduced in this paper. Zusammenfassung. Dieser Beitrag besch~iftigt sich mit Bl6cken linearer Gleichungssysteme, wie sie in der Signalverarbeitung bei verwendung yon Mehrkanalsystemen auftreten. Zuniichst wird das allgemeine Problem des Investierens yon Blockmatrizen im Zusammenhang mit der L6sung linearer Gleichungssysteme behandelt. Hierzu wird ein rekursiver Algorithmus angegeben, auJ~erdem wird die L/R-Zerleg~mg der Matrixrechnung auf den Fall der Blockmatrix erweitert. Im folgenden werden diese allgemeinen Methoden dann auf Blockmatrizen mit Toeplitzstruktur angeweudet, wobei sich in den meisten F~illen schnelle Algorithmen einsetzen lassen. Diskutiert werden insbesondere Mehrkanal-Wienerfilter, autoregressive und ARMA-Modelle sowie ein lineares Pr~idiktionsverfahren, bei welchem s + 1 Werte gleichzeitig vorhergesagt werden; die entsprechenden Algorithmen weren beschrieben. Fiir Gleichungssysteme mit Toeplitzstruktur weden zwei neue effiziente Algorithmen angegeben; der eine ist rekursiv, der andere basiert auf der schnellen Fouriertransformation. Diese Algorithmen stellen natiirliche Erweiterungen von Methoden dar, die fiir einkanalige Systeme bereits zur Verfiigung stehen; durch die vereinheitlichte Darstellung, wie sie im vorliegenden Beitrag vorgeschlagen wird, lassen sie sich verh~iltnismii~ig einfach herlieten.

R~sum~. Dans'cet article les auteurs s'int6ressent au probl6me de la r6solution de syst~mes-blocslin6aires, que l'on rencontre en traitement du signal num6rique multicanal. Dans un premier temps, le probl~me g6n~ral de l'inversion d'une matrice-bloc et de la r6solution de syst~me-block lin~aire est considerS; un algorithme de nature r6cursive est d~velopp~, et en m~me temps une forme-bloc du th6or~me de triangularisation est donn~e. Ensuite ces r6sultats g~n~raux sont specialists ~ des structures blocs de Toeplitz, permettant de retrouver des algorithmes efficaces existants. De ce point de rue, les cas de matrices li~es au filtrage de Wiener multicanal, h la mod61isation Ar et ARMA ainsi qu'au probl~me de la pr6diction (s + 1) pas en avant, sont consid6r6s et des algorithmes sont d~crits pour r6soudre ces probl~mes. Enfin des syst~mes-blocs de Toeplitz avec bandes sont consid6r6s et deux nouveaux algorithmes sont present,s, l'un 6tant de nature r~cursive et l'autre bas~ sur la FFT. Ces algorithmes constituent une extension naturelle de m~thodes d~jh disponibles pour le cas mono-canal, et leur d6rivation est obtenue de mani~re simple ~t partir de l'approche unitize introduite darts cet article.

Keywords. Multichan'nel signal processing (MSP), fast algorithms, Wiener filtering, parametric signal modelling, Toeplitz and banded Toeplitz matrices. 0165-1684/83/0000-0000/$03.00

© 1983 North-Holland

6

N. Kalouptsidis et al. / Fast algorithms for MSP

1. Introduction Multichannel and multidimensional signals arise in many fields as geophysics, image processing, control etc. Seismograph array systems used in earth resource exploration, earthquake studies and nuclear surveillance programs, antenna array systems in radar and radio communications, meteorological and air-pollution stations, as well as hydrophone array systems, provide some wellknown examples of multichannel signal sources. On the other hand multichannel EEGs and ECGs are well known representatives of this type of signals. Although processing on a per channel basis is a possible alternative, multichannel signal processing techniques are often preferred since they take into account the interrelations among the individual channels [1, 2]. Multidimensional signals, e.g. pictures, can be processed either as one-dimensional ones or in one-shot as a whole to get advantage of the correlation among the various dimensions. This last approach is computationally very tedious, and the multichannel approach is often taken as a compromising point of view [3-5]. For example in digital image processing this approach can be effected by dividing a picture in strips and then considering each strip as a multichannel signal. Multichannel linear modelling is a powerful tool capable to handle most of these types of signals. The model is determined by using an optimization procedure leading to the solution of block systems of equations. As is well known the structure of the block matrix, as well as the structure of each individual matrix element, have important consequences on the computational complexity of the various alternative algorithms leading to the solution. Besides linear system solution, matrix inversion and triangularization [10] are two other closely related topics of particular significance in signal processing. These subjects are also discussed here, though the main emphasis is focused on the solution of equations. It is also pointed out Signal Processing

that the subsequent algorithmic development concerns only the possible 'rich structure' of the block matrix and in no case that of the individual matrix elements. The resulting operations count refers to the number of matrix by matrix multiplications and matrix inversions. Fast algorithms for these operations can also be developed, however this task is not further pursued here. Section 2 gives a general recursive scheme. It provides the theoretical background for the subsequent development. Section 3 is devoted to block Toeplitz matrices, which is the most frequent case in signal processing. Previous results given by Akaike [6] and Whittle [7] are included and discussed in this general context. Section 4 is dealing with special cases of Toeplitz matrices appearing in various applications as multichannel Wiener filtering, multichannel AR and ARMA modelling and (s + 1)-steps ahead prediction. Special emphasis is given to the solution of systems with block banded Toeplitz matrices, encountered in image processing and tomography applications. Two efficient schemes extending the single channel problem are developed.

2. Recursive algorithms for block matrices This section is concerned with the development of a general recursive scheme for the solution of block linear equations and the inversion of block matrices. Despite its simplicity, this scheme will play a profound role in deriving fast algorithms when additional structure is present, as is the case in many multichannel signal processing applications. The approach followed here extends in a natural and nontrivial way related work dealing with the scalar case [8]. Making use of the general algorithm, a block form of the triangularization theorem is proved and subsequently utilized to improve the computational efficiency achieved thus far. Special consideration is given to cases where the matrix is featured by various symmetries.

IV. Kalouptsidis et al. / Fast algorithms for MSP

Finally the general scheme studied in this section can be proved useful in cases where submatrices are observed to have special structure. An efficient strategy will be to isolate the 'structured' blocks by an adequate block partitioning• The partitioned matrix inversion a n d / o r system solution may be performed with considerable gain in the computational effort when compared to the non-partitioned approach. 2.1. A general recursive s c h e m e

Let (1)

Ra =-d

be a general block system of equations, where R is a p × p block matrix consisting of elements which are l × l matrices, i.e.

R

=

[

Let us consider now the following partitioning of R and d ,

R,n+l = tZm T

dm+x =

r12

• • •

rlp l

r21

r22

" " "

?'2

l <~m <~p - 1 ,

Ym

(4)

[dV~d i n + l ]

with Rp = R and d, = d. Thus R.,+I has block order m + 1 and is built from the "old" block matrix R,. by adding on a column block vector xm, a row block vector z ~ and a matrix ym. Applying the matrix inversion by partitioning lemma [9], we can compute the inverse as follows: Xm

R~,L, = I-Z,n

Ym

[R~I =

rll

rp ,

7

-1 T

-I-WmClm V m

-1 v

Ot rn V m

OLm

T

(2)

w.~l] -1

J

T

wm = - R .1 x.,,

rppj

T

T

(5) (6)

o~,. = ym + z .~wm = y., + v reX.,, -I

rv;._ . . .

,

(7)

--1

v,. = - z , . R . , .

(8)

Thus R is a pl x pl scalar matrix, whereas

a=

ae

and

d=

p

d2

(3)

p

with ai, di, i = 1, 2 . . . . . p being I x k matrices.' The symbol T will denote interchange of rows with columns without affecting the matrix elements, whereas t will denote the usual transposition. Thus for example a T = [ a , a2 " ' ' ap], a t=[atla~

The above relations constitute an inversion procedure for R. Let us next apply the above relations to obtain a recursive scheme for solving (4). To this end, let

Using (5) and (9), (10) becomes a.,+l =

a., - w . , ~ . ,

-1

T

( v . , d , . +d,.+~) l

T -c~m (v,.dm +din+l) -1

J

(11)

which can be written more concisely as follows:

[

Wm

where km+l = - a ~1¢~,~,

' Scalar matrices are denoted by small letters, block matrices by capital letters and block vectors by bold face small letters.

(10)

Rm+,a.,+l = -d.,+l•

.., at].

The well known properties of t-transposition are not in general valid for the T-operation.

(9)

Rma., = -d..

~m

T = v m d m +din+,,

(13) (14) Vol. 5, No. 1, January 1983

8

N. Kalouptsidis et al. / Fast algorithms for M S P

1 is the l × l identity matrix and 0 is the l × l zero matrix. Finally using (8), (14) becomes

where

[wm ] W,n

T

/3"` = Z m a m

+dm+l

Win+l=

(15)

0T

1

'

(22) which provides an alternative way for the calculation of/3.,. From (5) and (6) we also get det a., = (det R.,+l)/(det Rm).

(16)

and A"`+l=diag{ao O~1 " , '

The partitioning we selected in the previous description, originates from the upper left corner and proceeds as indicated in (4). Other "block order" updating algorithms can similarly be derived if different partitioning strategies are adopted. When a linear system is to be solved the proper choice of partitioning is an important step, since algorithms resulting from various types of partitioning may differ substantially in computational effort. This result becomes important when matrices of special structure are considered. If the unknown matrices a~ of the block a multiply the elements of matrix R from the left, then the following system of equations results

aTR = - d T .

(17)

Proceeding exactly as for system (1) we obtain the following algorithm a ma-+ l

= [ a T i0]+k.,+a[v T i 1],

(18)

k.,+l = -/3.,a ~,1,

(19)

/3,, = a ,~x., + d., +1.

(20)

The system arrangement given by (17) appears more naturally than (1) in optimum multichannel signal filtering and prediction.

This paragraph demonstrates how the previous algorithm can be used to give the "block" triangularization theorem which states that any block matrix R can be decomposed as: Rml+l = Win+lAin1+1V"`+I Signal Processing

(21)

(23)

Indeed it is easy to see that it holds for m = 1. Assuming it is true for m > 1 and using (5) we obtain

[w"`a~lv"`+w,.a~lv~

-1

R"`+I :

a m-l V mT

-1

= Wm+lAm+l

Vm+l.

WmO/"`1] Olin1 _1

(24)

Using the above theorem, eqs. (7) and (8) become

w ,. = - W.,A ~nI V"`xm,

(25)

v.,T = - z . ,TW " ` A " ` 1 V"`.

(26)

The computation of w,, and v"` via (25) and (26) requires fewer l × l matrix multiplication operations than (7) and (8) in conjunction with (5).

2.3. "Symmetry" and its implications In the block case and contrary to the scalar one, a distinction has to be made between three types of symmetry. Thus a matrix R will be called 'block symmetric' if R = R T,

i.e.,

rii = r~i.

(27)

R will be named symmetric if it is symmetric in the usual sense R = Rt,

2.2. Recursive algorithms and the triangularization theorem

Ogre},

i.e.,

rij = r~i.

(28)

Finally R will be called "strongly symmetric" if it is symmetric and block symmetric. In this case we have t

rq = rii,

rii = rii

(29)

which implies that the elements of R are symmetric matrices.

9

N. Kalouptsidis et al. / Fast algorithms for MSP

Next we restrict ourselves to the symmetric case. Block symmetry as well as strong symmetry are of interest when matrices with additional structure are considered. If R is symmetric, (7) and (8) easily imply that t

T

w ,~ = v ,1

(30)

yielding then a considerable computational reduction. M o r e o v e r (30) simplifies the triangularization t h e o r e m if we note that W t = Vm. Thus (21) becomes RT.I+I

=

, -1 V,n+l V,n+lAm+l

=

Wm+,A~I+I

t Wm+

(31)

1.

In a similar way we can easily develop Cholesky type algorithms for the block case, following the methods described in [10, 8]. It is worthwhile to mention at this point that such types of symmetric systems are encountered in multichannel signal prediction using a 'covariance' type formulation as well as in multi-input multi-output system identification in a nonstationary environment. Time recursive schemes similar to the well known recursive least-squares algorithm can also be developed [11, 12].

tion of w and v which considerably reduce the overall complexity of the inversion procedure. M o r e o v e r it will be seen that the recursions yielding w and v are of the same type as the one obtained for a in Section 2. Let J denote the block matrix whose elements at the antidiagonal are identity l x ! matrices and the remaining are zero matrices. J satisfies the following properties: j 2 = I,

j = jT,

j = jt.

(34)

Then any block Toeplitz matrix R has the following persymmetry property: JRmJ = R ~ ,

(35)

l <~m<<-p

which easily implies the following:

R., =JR,~], ( R ~ ) -I =

RT.' =J(R T) i j, (36)

mlf.

JR

To simplify the exposition the block order of matrix J is not explicitly indicated. Note that in general (RT) -1 does not coincide with ( R - l ) T. Next we develop recursive relations for Win, Vm. Indeed, notice first that using (36), Wm+l and T Vm+l may be written -1

Wm+l = -Rm+lXm+x = -RT.l+lJJxm+l

3. Block Toeplitz matrices

T

=-J(Rm+l)

The results of the previous section are applied here to derive a fast scheme for the inversion and linear system solution of block Toeplitz matrices. A p x p block matrix R with elements l x l matrices is called block Toeplitz if (32)

rii = r i + k , ] + g .

T T 10m+ 1 =--Zm+lgm+

T Xrn+l

=

[Xm+l

T x.,],

T Zm+l

= [zm+l

Z mT ] ,

(33) It will be demonstrated below that relations (33) lead to simple recursive formulas for the computa-

(37)

Jxm+l,

-1 1

T

T

= - z m+lJ(R,,+l)

1

(38)

J.

To compute the inverse of R ~ + I we write T

Rm+l

Taking into consideration (4), one observes that the columns and rows of block Toeplitz matrices obey the following reproduction pattern:

--1

[

= . Rv~ LX T

(39)

Yo "

We put Yo instead of ym, since for a Toeplitz matrix ym = yo for rn = 1, 2 . . . . . p and apply the matrix inversion by partitioning lemma:

T

(R,.+I

)-,

T

[(R~)

= [

-1

--1

T

+l~6-,g~

6-1 T rn g m

fm 8:.' l aT.' J (40)

Vol. 5, No. 1, January 1983

N. Kalouptsidis et al. / Fast algorithms for MSP

10

where

and finally

f,.

= - ( R "T)

T

g,. = -x,.

-1

z,.,

T/nT

x-1

t~x ,.)

(41) ,

Jf,.+l:[J~"]+[

(51)

l Jk"+l

(42)

T

•,. = yo+x~d',. = Y0 + g , .Tz , . .

(43)

where

Using (33) and (40), (37) becomes Jw,.+l = r ( R T ) - I J x , . +[,. 8mlgTmJx,. +[,.

k~+l = - a ,.-1 13,., r

(52)

/3~ = v ~ z , . + z,.+l.

(53)

In a similar way we obtain the following recursion for g "

~mlx,.+l]

~ ~lg ~rX,,~ + ~ ~nlx"+l

--[

[ W ,.'I f

J"

Hence we conclude

T

g

g,.+tJ=[gr..JiO]+km+l[v,,

T

i l]

(54)

where

where

/3~ = g ~ x , . + x,.+l,

(45)

k.~+l = -

(46)

,./J,..

g g -1 k,.+l = -/3ma,. ,

(55)

/3g = xT"JW" + Xm+x .

(56)

a " , 6 " can be computed from relations (6), (43) respectively. An alternative, more efficient way is described below. Combining (6) and (44) we obtain w

In exactly the same way we obtain the following recursion for v,.: T

T

v

T

v , . + l J = [OmJ i O ] - { - k " + l [ g " i

O~m+l =

yo + Z m T W" +[zTj

.[ f "kkW"++l l J] iZm+lj[

T

1]

w

= a " + (z,.J[,. + z " + A k , . + .

(47)

Thus

where v

v

o

-1

k " + l = -/3,. t$,.,

(48)

/30,. = z r~Yf,. + z,.+l.

(49)

The above expressions indicate that recursive procedures for computing w,. and v,. can be obtained, provided that recursions of the same type can be found for f,. and g,.. This is accomplished next for f,,. Indeed from (41) we have f"+l

=

- ( R , T. + I )

w

a " + l = am +/3"k,.+1.

Similarly, using (43), (51), (56) we obtain the following recursion for 3,. g f 8,.+1 = 6,. +/3,.k,.+1.

/3.~=/3~

and

/3v=/3~

Indeed, from (42), (45), (56), it follows

Jf"+x = -RL½1Jz"+I.

/37.= g,.Jx,. T + x,.+l (50)

Application of the matrix inversion by partitioning lemma (5) gives

[R m'Jz" + w " a m1 (vTmJ.~" "Jr-Z , . + I ) ] I/"+1 = - [ ,~;I(vT, j z " + z"+O J Signal Processing

(58)

A further simplification results if we observe that

-- l z " + 1

or, using ( 3 6 )

(57)

T ( R T ) - I J x , . + X,.+I

= x~JR/.ix,. + x,.+l T = x , . J w , . +X,.+l=Bgm.

In the same manner we get/3~ =/3~.

(59)

11

N. Kalouptsidis et al. / Fast algorithms for MSP

Finally notice that det a.~ = det &..

(60)

(60) is a consequence of (16), since

Table 1

det am = (det Rm+l)/(det Rm)

Algorithm for block Toeplitz matrix inversion

T d e t J det Rm-1 d e t J - detJdetRTdetJ =detSm"

~fm = ~ v = 12Tjzm + Zm+l = zT Jfm + Zm+l

For block Toeplitz system solution O(3p 2) matrix multiplications and 2p matrix inversions are needed. Some supplementary operations are necessary for block Toeplitz matrix inversion. Indeed the data specified so far yield the first and the last block columns and rows of Rm+l. In fact by making use of (36), (40) we obtain:

n G1 = J ( R ~ + I ) - 9 r [ ( R T~.) 1 +f,,,&,,- 1 g.~T =

-1

~[

fm~m.l --1

T

~mgm

=dfrnt~ll

t

87~

[

(~mI

= [jf~8~

Table I summarizes the portion of the algorithm related to matrix inversion only.

t~ mI

JJ

t~nlgTj

k ~ + l = --CCmlfl~

k~+~= - - ~ n

1

+[ l j km+l t~T+l J : [IJTJ ', 0]q- k~n+l [gT i 1] gL+,J=~LJ lo]+g~+,tv T 11] am÷l = a,~ +/~k,~+l det a,. 8m+1= &~+/3~k~+1 = a,. + k,~,+d3~ = det t~,. =6,.÷k~+ff3~ (R.I)11 = ~pl_ 1

(Rp-1)lj+l

(Rpl)i+l,1

= 3p_,gp lJ /'=1,2 . . . . . p - 1

(gT)-lJ+fmt~rnlgTq

8;.'g~J

k,~+l = -8~'fl~, k~+, = -/3~ 62,',

J ]

R~, ~ +JfmS~gT.,jJ"

Thus the first block row and column of R ~ ~ are determined by gp_~, fp-1, 8p-1. To compute the remaining elements we observe from the above relation: (Rp-1 )i+1,/+1 -'~ ( R p 1 , )ij q- (Jfp-1 ~ p ! l g T - 1 J ) i j , i,j=1,2 ..... p-1

i=1,2 . . . . . p - 1 (R~l)i+l,i+l = (R-pq)q+(jfp_.l t~o 1lgpT 1 J - w p l°ttJ1lVpT 1)q i,/=1,2 . . . . . p - 1 (R~ l )ij = (Rp 1)p+b j.p+l-i (persymmetry)

Table 2 summarizes computations required for block Toeplitz system solutions. Finally we note that the initialization of the various algorithms described so far can be easily done by considering the linear systems defining a,., w,., vm, g , . , f m for rn = 1.

which combined with the partitioning (5) gives (Rp 1)i+1,i+1 = (Rp 1)q + ( J f p - I 8 p-1- l g p T- l J - W p - l O t p - l V p1

T 1)q.

Note that for the system solution R a = d ( a T R = d T) the recursions for win, fro,/3 w,/3~, am, &, (v T, T v gin, tim, fl~, am, &,) are sufficient. The algorithm for block Toeplitz matrix inversion is derived in a different way in [6], whereas the linear system solution algorithm was first given in [13]. The triangularization aspect alone is examined in [10] for block Hankel and block Toeplitz matrices.

4. Special cases appearing in multichannel signal processing Block matrices with special structure appear quite often in multichannel signal processing. Thus the topics of block matrix inversion and the corresponding linear system solution via efficient recursive algorithms, are of special importance in such diverse fields as geophysical and seismic signal processing, digital image processing, tomography, economic forecasting etc. Vol. 5, No. 1, January 1983

12

N. Kalouptsidis et al.

Table 2

Fast algorithms for M S P

equations [1, 15]

Algorithms for block Toeplitz linear systems solution of Ra = - d , aTR = - d r

aTRxx =--r~x

(62)

where Ra = -d,

aXR = - d T

13,, = z ~ a , , +d,.+a,

T •m = a mxm +

kin+ 1 = --ot ma [3m,

k m + l = --~mOtm 1

t32 = xLYw.,+x,,+,, k~,+a = -8,,at3,~,

t3~ =z~Yf,. + zm+a,

a=[ao

(63)

rax=[ra~(O) rd~(1) • .. tax(p)] T,

~ =g~x,,,+x.,+~

[ r~x(0)

r~(1)

•..

Rx~= rx~(-1)

rxx(0)

...

k~,+~ = -t3~ 87.,~ ~,=v~]z.,

al . . . ap] T,

d,~ +l

+ z,.+l

k~+l = -~'.,c~,1 r~(-p)

+[ 1 ]k-+,,

[v~dio]

r~x (i - ] ) = E { x (n - ])x'(n - i)},

(66)

+k=+,[g~ i 1]

rax(i) = E ( d ( n )x' (n - i)}.

(67)

gL+,] [g~ i o] =

+kGl[VL i 1] 6,~+1= 6,. +/3~,k~+, = 6m+ kg.,+,/3~. /3~ =/3~, detain = det &,, /3f =/3~,

= rtxx ( - i ) .

(68)

Thus according to the definitions of Section 2.3, Rxx is symmetric but not block symmetrix. Relation (68) with the use of (30) and (45) gives then (B~)' =/3~.

(69)

This last result is due to Burg [3] and leads to an additional saving in computations.

4. I. Multichannel Wiener filters The input-output relationship for an FIR multichannel digital filter with /-input channels and k-output channels is P

(61)

i=o

where the filter coefficients ao, al . . . . . ap are k x l scalar matrices. Let d ( n ) be a desired response signal consisting of k channels. If the signals x(n), d ( n ) are jointly stationary the optimum Wiener filter, in the sense of minimizing the MSE E { e T ( n ) e ( n ) } , e ( n ) = d ( n ) - y ( n ) , can be easily found by using the orthogonality principle [14]. It is given by the solution of the following block system of normal Signal Processing

The algorithms found in Section 3 apply thus here since Rxx is block Toeplitz. Furthermore, Rx~ being an autocorrelation matrix implies

r~x (i)

The purpose of this section is to discuss block matrix structures arising in such practical applications and the resulting simplifications on the algorithms described so far.

y(n) = - Y aix(n - j )

1/1,

rx~i0) 3 (65)

+ k~+l[vT., i 1] v~+,J =

rx~(p) 1

rxx(-p + 1) . . .

a~+, = [,,~ io]

J-m+,=t o j

(64)

4.2. Some other applications In certain practical situations R appears in even simpler form for example it is strongly symmetric. This for instance occurs in the digital processing of arrays of underwater sound sources [16]. If R is strongly symmetric, combining the results of Sections 2.3 and 3 we deduce the following simplifications: wm = vm = f,, = gin, k.~ = k ~ = k ~ =k°~,

(70)

O~rtl ~ ~rrt*

The reader can easily write down the algorithms for the inversion and linear system solution by

N. Kalouptsidis et al. / Fast algorithms for MSP

13

simple inspection of Table 1 and Table 2. The resulted inversion algorithm was obtained also in [17] under the additional assumption that the ! × l matrix elements of R are persymmetric.

due to the special form of rq+l.q+p we have

4.3. Multichannel A R a n d A R M A

This relation leads to a further saving since the recursion for gm can be omitted. Let us next define a * as follows:

modeling

This paragraph is concerned with the treatment of equations involved in the determination of the AutoRegressive (AR) parameters of an AutoRegressive-Moving Average (ARMA) model. The special case of pure A R modeling is also considered. The further simplifications that will be established in this case result solely from the form of the right-hand side vector and apparently they concern only the system solution mechanism and in no way affect the inversion process. Let s ( n ) be a k-channel stochastic process governed by the following relation p

s(n)=-Y

q

ais(n-j)+

(71)

Y. b ~ ( n - j )

i=1

i=0

where

T

T

for n # j .

(72)

The above scheme is known as a (p, q) order multichannel A R M A model. As is well known [7], the A R - p a r a m e t e r s can be obtained from the following modified set of block Y u l e - W a l k e r equations T

a qr~Rqo = - r q

(73)

+ l , q +P"

where

RqT=

rq

rq-1

rq+ 1

rq

[ r,+p-i r qT+ l . q + p

T

a~T

1

~ [Fq+l

• • •

rq-p+l

"'"

rq-P+2

••• rq+2

'''

,

(74)

rq

rq+p].

(75)

The matrix Rap is non-symmetric block Toeplitz and the algorithm developed in Section 3 can be used to solve the system (73) recursively. H o w e v e r

T

=gin.

T = v,.J.

(76)

(77)

Then the algorithm for the solution of (73) becomes a mT+ l

=[a~iO]+k,.+~[a*..wJi 1],

(74)

a m.T +l

= [ a . T i O ] + k . , +, l [ a , . jTi

(75)

1],

km+l = - / 3 , , a 1 ,

(76)

13,, = a~/rq+l.q+m +rq+m+l,

(77)

*

*

*

km+l = --~.n(am)

-1

,

(78)

. Jrq t.q-,. + r , - m - 1 , 13*~ = a , ,T

(79)

am+l

(80)

= a m q" k * ~ + l f l m , --OLrn

+k,,+l/3m,

det o~., = det a *,,.

E{e(n)eV(j)} = 0

T

-1

= -rq+l,q+mJ(Rq,,)

~rn+l

E { e (n)} = 0,

T

a ,.J = - r q + l.~ + , . R q,.J

(81) (82)

Besides the solution of (73) the above recursions yield the solution of the following system as well a,Tn

T

T

qp lxqp = - r q

1,q-p.

(83)

The case q = 0 corresponds to a pure A R model and in contrast to the scalar case (single channel) [8] no additional simplifications can be achieved. M o r e o v e r for q = 0 the system (83) gives the parameters of a pure A R model running backwards. Finally, a,,, a , , are the covariance matrices of the forward and backward residuals respectively. Since det a,, = det a * , the generalized variances of these residuals are equal. These interpretations are no longer valid, if q # 0 [7]. 4.4. Prediction s + 1 steps a h e a d

This is a special case of the problem considered in 4.1 where the desired signal d ( n ) is equal to x (n + s + 1). The prediction coefficients then are Vol. 5. No, 1, January 1983

N. Kalouptsidis et al. / Fast algorithmsfor MSP

14 given by (62) where now

Now notice that nT

r d x = [ r ~ ( s + l ) r~(s + 2) • '" rxx(s+p+l)] T. (84) For s = 0 the solution can be obtained by the algorithm of Section 4.3 corresponding to the A R modelling, whereas for any other s the more complex scheme of Table 2 may be used. A more efficient technique can be developed however for relatively small values of s and this is described next, following closely [8] where the scalar case has been treated. The systems we are dealing with have the form

~TR ov =

a

T

T

rv+~].

(86)

The proposed algorithm consists of the following steps: (i) Form the system 0T

T

ap+sRo,p+s = -rl.p+s.

T

T

T

en T

-1 1

= - r p + s , l + n ( R O,p+s-tl )

.-1

= --Fp+s,p+s-m+ 1 [ K Orn )

= h,,

for p + s - n -= m.

(93)

Hence, substituting (93), (92) into (90) gives (88). Clearly the coefficients kn need not be computed since (88) runs backwards.

r

n+lT

Block banded Toeplitz matrices appear in various signal processing applications including image processing and tomography. In all these applications fast algorithms are needed to effect the corresponding matrix inversion and linear system solution. A block matrix R is called banded with band (M, N ) if:

(87)

and obtain its solution by recursions (74)-(81). (ii) Using ap+ff ox as initial condition, available from step (i) apply the following recursion backwards to obtain a psTJ, the desired solution, in s steps v

rli=O

xp-1

(94)

[0]

=

JXM

,

Zp-1 =

[01 JZN

(95)

where

0T

(88)

To establish (88), we apply the algorithm of Section 2 on the system T

forf>M+i,i>N+f.

If in addition R is Toeplitz from (33) it follows:

ap+~ ~J =[ap+s-,-1J i O]+k,[ap+s_n-lJ i 1].

T

h Ro,p+s =-rp+~,l

general

XM

=Ix1

X2

"'"

XM] T,

ZN

=[Zl

Z2

"'"

ZN] T.

0 will denote the zero vector of proper dimension. More generally,

(89)

to get

x,,+l=[0 Tix~r] T form>M and

T h.,+, = [ h T i 0 ] + k ~ + , [ g ~ i 1]

(90)

where h~

T

= - r l+.,p + J ( R o,v+s-. )

(85)

rs+l,v+,=[r,+l r.+2 " "

T

-1

4.5. Block banded Toeplitz systems - r s +T l , p + s .

Rop is given by (74) for q = 0 and

nT

T

ap+s J = - r l + . , p + s - . + . R o , v + s - J

1-

= --rp+s,p+s

m+l

,nT

,1

(IKOm)

(91)

and by (76)

gT =aOmTj. Signal Processing

(92)

grn+l

= [0 T i

z T J ] T for m > N .

(96)

In this section extending previous work on the single channel case reported in [18] and [19], it is shown that for banded block Toeplitz matrices the presence of a highly sparse part can be exploited to speed up the algorithms of Section 3. Then an

15

N. Kalouptsidis et al. / Fast algorithms for MSP

FFT based scheme is outlined, suitable for parallel processing.

4.5.1. A recursive algorithm Systems with banded matrices share the following property: Suppose that the last N components of the unknown vector a have been computed. Then the remaining components of a can be found by a procedure of back substitution and the required matrix multiplications are approximately 2pM - 2 M 2, (M ~ N ) . Because of the above property it is sufficient to find efficient means for computing the last N components of a. This in turn can be achieved as follows. Consider the block matrices:

F=[IMIO],

L=[Oils],

G = [LF]. (97)

F has block dimensions M x p and when operates on a vector x, it gives its first M components. L is N x p and acting on x gives its last N components. We multiply the vector resursion of the algorithm of Table 2 by the matrices F, L, G from the left accordingly and we obtain

Relations (98f, g, h) follow easily from (14), (45), (53). The set of formulas (98) involves only the last N components of a, w, f and the first M components of w, f. Thus at each step only O(p) operations are required• The algorithm starts from aN which is obtained by the scheme of Table 2 and gives Lap (the last components of a), the initial condition of the back substitution filter• By similar means we can establish an efficient inversion procedure for banded Toeplitz matrices.

4.5.2. A n F F T based algorithm The algorithms developed so far are recursive in nature. Next we briefly discuss an extension of Jain's method [18] for the block case. Its main feature is that for the most part it is nonrecursive and affords parallel processing architecture making it suitable for real time applications• The banded structure of R allows us to write R as follows: R = R c - R~

(99)

where

vo xl x2 "''XM 0 0 ' ' ' O ZN ''" Z2 Z: \ 7.

\\

~

z2

\

Wm

ZN

• •

XM

f

•

\

f

-1

km+lW

(

km+l = -c~m Bin, f am+l=am+~mkm+l, w

\ ~.

x

\

\

\ \\

\

x

\

\

\

Xl

\

(100)

x.\\ \ "-, \ ,~, X 1 x

\ XM

.

x

\

x x\

" " "

~,

\ XM

\

\ \ \ X2

xZN

x \

x\

\

= __,~"1 ~m,'~'

\

\

0

0 ,,,\ \

x.

\

X2

km+l = --a~nlflm,

" x \

\

(98b) Rc =

.

\ \ • x

•

Gw,,,+~=G[ 0 ] + G [ 1 ]k~+l, w" tJ[,.

\\

" " "

ZN

ZN

" " '

" "" Z2

Z1

YO

(98d) "'"

w f

6,.+1 = 6,. +Bmk,.+l,

(98e)

0

Z2

Z1-

~ZN

Rs

[3,. = z ~ J L a ~ +d,.+l,

(98f)

[3~ =x~JFwm + Xm+l,

(98g)

[3~ = z T L j f , , + z,.+I,

(98h)

det a,, = det 6,..

(98i)

(101)

XM

Xl

"'"

XM

• " "

0

Rc is the block circular extension of R and Rs is a block sparse matrix. Vol. 5, No. 1, January 1983

16

N. KalouptsMis et al. / Fast algorithms for MSP

Then relation (1) is equivalent to: a = b +SRsa,

(102a)

Rob = -d,

(102b)

S = R e 1.

(102c)

The solution of (102b) as well as the inversion of (102c) can be carried out by block FFT. Indeed Rc is block circular, hence it can be block diagonalized by the block FFT (see Appendix). Next let H1, H2 denote the following block triangular matrices of Toeplitz form.

H1

I

Ix -. 0 1 ° -z j

... IMJ

with the system (107) is block Toeplitz. Therefore Ul, u2 can be obtained using the algorithm of Section 3. Knowing Ul, u2, a is obtained from (105). Employing the above method the major part of the computations is carried out by FFT allowing thus parallel processing, while only a small part is executed recursively, that is, the solution of (107). A similar method with slightly more loaded recursive portion results if ul, u2 instead of being computed from (107), are determined from (104), where L a , F a are computed by the algorithm of Section 4.5.1.

5. Conclusion

(103) and ul =H1La,

u2=H2Fa

(104)

where the F and L operators are defined in (97). Then eq. (102a) becomes

a=b+S

0

(105)

.

LU2J

Multiplying the above equation by G (see relation 97) and partitioning S in an obvious way, we obtain Fa = Fb + S~ul

q-S2u2,

(106)

L a = L b + S a u i + S 4 u 2.

Taking into account (104) it follows that

([.021 H01a]--[S: SlJ]lfJ1j (107) The above system yields u 1, u2. Notice that since H1,//2 are block triangular Toeplitz matrices their inverses have also the same structure. Moreover the inverse of a block circular is block circular as well (see Appendix). Taking into account these remarks it easily follows that the matrix associated

SignalProcessing

This paper was concerned with the solution of block type equations that typically arise in multichannel signal processing. Focused on the computational complexity, it described a family of algorithms, starting with the general case of block matrix inversion and system solution and proceeding by progressively imposing more structure on the parameters of the problem. The algorithm for the general case is recursive in nature and its performance can be improved by invoking a suitable block triangularization theorem stated and proved here. Subsequently this general scheme was applied to block Toeplitz matrices and most existing algorithms were derived as special cases. Multichannel Wiener filtering and other digital processing applications were examined in the above context. The multichannel Autoregressive and ARMA modelling were considered and a new algorithm was given for the s + 1 steps ahead prediction problem. The important structure for signal processing application of banded block Toplitz matrices was examined and two algorithms extending single channel methods were presented in the general framework introduced here. Most of the theory presented in the previous sections can be extended to c~-Toeplitz matrices, c~ being the displacement rank of the matrix [20, 21] whereas formulations and interpretatio_ns

17

N. Kalouptsidis et al. / Fast algorithms for MSP

based on the ' p o l y n o m i a l ' a p p r o a c h can be f o u n d

for instance in [22]. Apparently a continuing effort towards reduction of computatonal complexity is necessary since among others, it may lead to important practical realizations. Thus further research is currently c o n d u c t e d on efficient nonrecursive algorithms based on DFT.

References [1] S. Treitel, "Principles of digital multichannel filtering", Geophysics, Vol. 35, No. 5, October 1970, pp. 785-811. [2] J.F. Claerbout, Fundamentals of Geophysical Data Processing, McGraw-Hill, New York, 1976. [3] E.A. Robinson and S. Treitel, "'Digital signal processing in geophysics", in: Alan V. Oppenheim, ed., Applications of Digital Signal Processing, Prentice-Hall, Englewood Cliffs, NJ, 1978. [4] B.R. Hunt, "Digital image processing", proc. IEEE, Vol. 63, 1975, pp. 693-708. [5] B.R. Suresh and B.A. Shenoi, "New results in twodimensional Kalman filtering with applications to image restoration", IEEE Trans. Circuits Systems, Vol. 28, 1981, pp. 307-319. [6] H. Akaike, "Block Toeplitz matrix inversion", SIAMJ. Appl. Mach., Vol. 24, 1973, pp. 234-241. [7] P. Whittle, "On the fitting of multivariate autoregressions and the approximate canonical facorization of the spectral density matrix", Biometrica, Vol. 50, 1963, pp. 129-134. [8] G. Carayannis, N. Kalouptsidis and D. Manolakis, "Fast recursive algorithms for a class of linear equations", IEEE Trans. Acoust. Speech Signal Process., April 1982. [9] F.R. Gantmacher, Theorie des Matrices, Tome I, Dunod, Paris, 1966. [10] J. Rissanen, "Algorithms for triangular decompositions of Block Hankel and Toeplitz matrices with applications to factoring positive matrix polynomials", Math. Computation, Vol. 27, 1973, pp. 147-154. [11] T.C. Hsia, System Identification, Lexington Books, 1978. [12] P. Eykhoff, System Identification, Wiley, New York, 1974. [13] R.A. Wiggins and E.A. Robinson, "Recursive solutions to the multichannel filtering problem", J. Geophys. Res., Vol. 70, 1965, pp. 1885-1891. [14] A. Papoulis, Signal Analysis, McGraw-Hill, New York, 1977. [15] E.A. Robinson, Multichannel Time Series Analysis with Digital Computer Programs, Holden Day, San Francisco, 1967. [16] R.L. Pritchard, "Mutual acoustic impedance between radiators in an infinite rigid plane", J. Acoust. Soc. Amer., Vol. 32, 1960, pp. 730-737.

[17] G.A. Watson, "An algorithm for the inversion of block matrices of Toeplitz form", JACM, Vol. 20, 1973, pp. 409-415. [18] B. Dickinson, "Efficient solution of linear equations with banded Toeplitz matrices", IEEE Trans. ASSP, Vol. 27, No. 4, 1979. [19] A.K. Jain, "Fast inversion of banded Toeplitz matrices by circular decompositions", IEEE Trans. ASSP, Vol. 26, No. 2, 1978. [20] T. Kailath, S.Y. Kung and M. Morf, "Displacement ranks of matrices and linear equations", J. Math. Anal. Appl., Vol. 68, 1979, pp. 395-409. [21] B. Friedlander, M. Morf, T. Kailath and L. Ljung, "New inversion formulas for matrices classified in terms of their distance from Toeplitz matrices", Linear Algebra Appl., Vol. 27, 1979, pp. 31-60. [22] P. Desarte, Y. Genin and Y. Kemp, "'Orthogonal polynomial matrices on the unit circle", IEEE Trans. CAS, Vol. CAS-25, 1978, pp. 149-160. [23] D. Manolakis, N. Kalouptsidis and G. Carayannis, "Fast algorithms for discrete-time Wiener filters with optimum lag", IEEE Trans. Acoust. Speech Signal Process., to appear. [24] N. Kalouptsidis, G. Carayannis and D. Manolakis, "'On block matrices with elements of special structure", proc. IEEE Int. Conf. ASSP, Paris, 3-5 May 1982, pp. 17441747.

Appendix This appendix gives a brief outline of the problems of diagonalization, inversion and linear system solution for block circular matrices by using the Discrete Fourier Transform (DFT). The use of existing fast algorithms (FFT, WFTA) for the DFT allows the efficient solution of the above problems in a finite number of steps. The Block D F T (BDFT) of an N - p o i n t s e q u e n c e of 1 × l matrices

{Xn}I~

1 is defined as

N-1 =

x.wN,

WN = e

j2rr/N

k = O , 1. . . . . N - - 1

(1)

where and

j = 4-~.

Relation (1) can be written in matrix notation as (2)

~? = W x

where X = [ X o X1 " ' "

XN-1] T,

(3) J~ = [ X o .~'1 " " " 3~N

1] T

Vol. 5. No. 1. January1983

N. Kalouptsidis et al. / Fast algorithms for MSP

18

and W is the N x N B D F T matrix W = [ w ~ " 1],

and

k,n=0,1 ..... N-I,

(4)

W A = [ w o go 1¢1c1

"'"

WN--ICN-1].

(12)

1 being the l x 1 identity matrix. Note that W is different from the l N × l N D F T matrix. If I is the l N × l N identity matrix, then it follows easily that

Taking into account that wuu±k = w~ and after a simple rearrangement of terms (see [14, p 86]) it can be shown that

(5)

which implies C W = W A and (8) easily follows• Thus a block circular matrix can beblock diagonalized by using the BDFT. Based on this fact the inversion and linear system solution for block circular matrices can be done very efficiently. Indeed from (8) and (5) it immediately follows that

1 ww, N

= i,

w-1 = L w* N

where * denotes complex conjugate. Hence the inverse B D F T (IBDFT) is given by x = 1I~¢*~.

(6)

In the following notation £ = B D F T { x } and x = IBDFT{i} will be used for (2) and (6) respectively. Now let C be a block circular matrix of the form

c=

IcNO C1 C2 ''" 1 Co

•

,

LC1

CN-I1

cl.-•.

Cz

C,

~

(7)

,

"''"'CoJ

W*

(8)

(9)

To prove (8) we represent W by its block columns, i.e. w1

"'"

WN-1]

(10)

where w k = [ 1 w k l w2kl . . .

Comparison of (8), (14) shows that C -1 is also block circular with first row the I B D F T of the vector St = [ S O 1 S 1 1

"•"

C~N- 1- 1 IT•

(15)

(16)

Ca = - d

can be obtained by a=-C

(17)

ld=-W[A-I(1w*d)].

Inversion Step 1. Step 2. Step 3.

of a block Compute Compute Compute

circular matrix C ~ = BDFT{c }• ~' = [~o 1 ¢1 --1 ' ' " c' = IBDFT{~'}.

C W = [Cwo C w l • • • CW~-l]

(11)

~,/I__1 I T•

e is the first row of C and c' the first row of which is also block circular.

W (N 1)kl]"

By using (9) and (10) we obtain

Signal Processing

(14)

C -1 = I w A - 1 W * .

~ = BDFT{c },

~- [1~0 C1 " " " I ~ N - 1 ] T•

W=[!il,'o

(13)

The above algorithms can be summarized as follows:

where A is a block diagonal matrix of the form A --- diag{~},

k=0,1 ..... N-1

The solution of the linear system

where Cj are I × l arbitrary matrices, with scalar elements, not necessarily circular. It is next shown that W diagonalizes C, i.e. C = ~1w a

C W k = W k Sk,

Solution of the linear system Ca = - d Step 1. Compute d = IBDFT{d}. Step 2. Compute ~ = BDFT{c}.

C -1

N. Kalouptsidis et al. / Fast algorithms for MSP

a

Step 3. C o m p u t e ~' = [ g o 1 C~1 1 • • • ~ N l l ] T . Step 4. Set A - l = d i a g { g ' } and compute d = la~.

Step 5. C o m p u t e a = IBDFT{ti}. The computational complexity of the above algorithms can be easily determined in terms of the n u m b e r of B D F T ' s and the l × ! matrix inversions required. Since the B D F T is a point-wise D F T its computational effort depends on the specific implementa-

19

tion of the fast algorithm used for the D F T and the structure of the block elements. Thus for arbitrary ci l 2 D F T s are required whereas for T o e plitz or circular block elements l DFTs are enough. The complexity of the required inversions depends on the structure of the block elements as well. Thus in case the block elements ci are circular matrices, due to the linearity of DFT, gi are also circular and their inversion can be done by using the above described techniques (for N = l and l = 1) in O(l log l) operations.

Vol. 5, No. I, January 1983