An FPT algorithm with a modularized structure for computingtwo-dimensional discrete Fourier transforms

2148 IEEE TRANSACTIONS ON SIGNAL PROCESSING. VOL. 39, NO. 9. SEPTEMBER 1991 for the vector U . Now, if this equation is premultiplied by U = U ,

or U,, then U also satisfies

UCV = Ud. (21)

With the change of variables U = U*q, substitution into (21) results in

UCU*q = Ud. (22) Consequently, the solution vector U of (20) can be obtained from the solution vector q of (22) or q can be obtained from U .

Note that if a system of real equations is Toeplitz-plus-Hankel ( T

+

H ) , where T i s symmetric Toeplitz and H i s skew-centrosym- metric Hankel, then the equations may be transformed into Her- mitian Toeplitz and solved with 1 .25n2

+ O ( n ) complex multiplies

or 3.752,

+

O ( n ) real multiplies. This is a significant improvement in complexity over the approach of [8] which requires 12n2

+

O ( n ) real multiplies, and is slightly lower in complexity than the ap- proach found in [9] which uses an entirely different development and requires 6n2

+

O ( n ) real multiplies.

111. CONCLUSION

In this correspondence, we have shown that constant unitary ma- trices exist which transform a Hermitian Toeplitz matrix into a real Toeplitz-plus-Hankel structure. As a consequence of this property, some real symmetric Toeplitz-plus-Hankel matrices may be con- verted to Hermitian Toeplitz matrices and vice versa. The unitary matrices presented are different from the one given in [2] and have two advantages: i) they transform Hermitian Toeplitz matrices into real matrices also possessing a special form, i.e., Toeplitz-plus- Hankel, (the result of [ 2 ] does not) and i i ) the unitary matrices are simple to express mathematically, thus making it easier to manip- ulate and interpret results based on them for analytical purposes. A consequence of the unitary transformation

U ,

presented in this correspondence is a relationship between the real and imaginary parts of eigenvectors of Hermitian Toeplitz matrices. This rela- tionship also differs from the eigenvector symmetry property often given in the literature, e . g . , [ I ] .

As a practical observation on eigenspace decomposition consider

the following. If all eigenvalues and eigenvectors are required, i t is more efficient to transform a Hermitian Toeplitz matrix to a real matrix (using either the result of this paper o r that of [2]) and use a standard algorithm such as the QR algorithm o r Jacobi rotations

[ 5 ] . If only a few eigenvaludeigenvector pairs are needed, how-

ever, it will be more efficient to work on with the Hermitian Toe- plitz matrix using a fast algorithm such a s the Toeplitz eigensystem solver (TESS) [IO] that exploits the special structure and can find specific eigenpairs.

REFERENCES

111 M. J. Goldstein, “Reduction of the eigenproblem for Hermitian per- symmetric matrices,” Math. Computation, vol. 28, no. 125, pp. 237- 238, Jan. 1974.

[2] A. Lee, “Centrohermitian and skew-centrohermitian matrices,” L i n - earAlg. Appl., vol. 29, pp. 211-216, 1980.

[3] N . Levinson, “The Wiener rms (root-mean-square) error criterion in filter design and prediction,” J . Marh. P h s . , vol. 25, pp. 261-278. [4] H . Krishna and S . D. Morgera, “The Levinson recurrence and fast algorithms for solving Toeplitz systems of liner equations,” IEEE Trans. Acoust., Speech, Signal Processing, vol. A S P - 3 5 , no. 6 , pp. 839-848, June 1987.

Baltimore, MD: John Hopkins University Press, 1983.

[5] E. H . Golub and C . Van Loan, Marrix Compurarions.

J . H. Wilkinson, The Algebraic Eigenvalue Problem. New York: Clarendon, 1965.

A. Cantoni and P. Butler, “Eigenvalues and eigenvectors of sym- metric centrosymmetric matrices.” Linear Alg. Its Appl., pp. 275- 288, 1976.

G . A. Merchant and T. W . Parks, “Efficient solution of a Toeplitz-

plus-Hankel coefficient matrix system of equations,” IEEE Trans. Acoust., Speech, Signal Processing, vol. ASSP-30, no. 1, pp. 40- 44, Feb. 1982.

I. Gohberg and I. Koltracht, “Efficient algorithm for Toeplitz-plus- Hankel matrices,” Integral Equarions Oper. Theory, vol. 12, no. I , Y. H . Hu and S.-Y. Kung, “Toeplitz eigensystem solver,” IEEE Trans. Acoust., Speech, Signal Processing, vol. ASSP-33, no. 4, pp. pp. 136-142, 1989.

1264-1271, Oct. 1985.

An FPT Algorithm with a Modularized Structure for

Computing Two-Dimensional Discrete Fourier

Transforms

Ja-Ling Wu and Yuh-Ming Huang

Abstract-In this correspondence, the fast polynomial transform

(FPT) for computing two-dimensional (2-D) discrete Fourier trans-

forms (DFT’s) is modularized into identical modules. In this new

method, only FPT’s and FFT’s of the same length are required. As a

consequence, the architecture is more regular and naturally suitable for multiprocessor and VLSI implementations.

I. INTRODUCTION

Polynomial transforms, defined in rings of polynomials, have been shown to give efficient algorithms for the computation of 2-D convolutions [ l ] , [2] and to map multidimensional DFT’s into 1-D DFT’s [3]. In this approach, the Z N - 1 polynomial ring is first factorized into several irreducible cyclotomic polynomials and then operated on by the polynomial transforms. Thus, one has to use FFT’s and FPT’s with different lengths through the F P T stages. This fact complicates the control of multiprocessor and/or VLSI FPT.

Recently, Truong et al. [4] proposed the modularized F P T struc- ture for computing 2-D cyclic convolutions. In this approach, the I - D cyclic polynomial convolution is decomposed into cyclic con- volutions of polynomials, all of the same length. The regularity of the new algorithm makes it of great practical value. Based on the ideas in [ 4 ] , in this correspondence, we modularized the F P T al- gorithm for computing 2-D DFT’s. In this new method, only FPT’s of length N , and FFT’s of length N2/2‘ are required.

11. T H E MODULARIZED F P T ALGORITHM FOR 2-D DFT’S

Without loss of generality, let us consider the 2-D D F T of the

N : N 2 complex sequence x,,,,~ defined by N I - l N , - I

where N , = 2”, and W , = f o r i = 1, 2 , and t , 5 t z . Manuscript received February 6, 1989; revised February 12, 1991. The authors are with the Department of Computer Science and Infor- mation Engineering, National Taiwan University, Taipei, Taiwan, 10764. Republic of China.

IEEE Log Number 9101527.

1053-587X/91/0900-2148$01 .OO 0 1991 IEEE

1

IEEE TRANSACTIONS ON SIGNAL PROCESSING. VOL. 39. NO. 9. SEPTEMBER 1991 2149

-

In order to compute X k l , k 2 by the polynomial transform (PT) [3], we represent once again (1) by a set of three polynomial equations

N I - I ,,I = o - X k l ( Z ) =

c

X , , ( Z ) WS:’ mod (ZN’ - 1) ( 2 4 - - X A ~ . L ? X k ,

( z )

mod

(z

-

w$).

(2C) Based on [4], if the polynomial Z” - 1 is decomposed as

2‘- I

1 =

11

(ZW?’ - a h ) (3)

z N ? - A = O

where 1 5 r 5 t? and a = e’211’2‘, i.e., the 2‘th root of unity. Then the index k2 can be partitioned into 2‘ sections by

k2 = 2 ‘ . d

+

t (4)

w h e r e 0 5 t 5 2‘ - 1 a n d 0 5 d 5 N2/2‘ - 1. Furthermore, it can easily shown that

( 5 )

z

- W 2 ‘ . d + r ( z N ? / 2 ‘ - a 2 ‘ - l

N?

where “alb” means a is divided by b. Let

XLI(Z) = n 2 = 0

2

x:,,,,,,Z“’

=

X,?,(Z) mod (ZN?’” - a * ’ - ‘ ) .

N 2 / 2 ’ - I

(6) Derivations and discussions of the F P T algorithm for 2-D DFT’s are given in [3] and are not repeated here. From [3] and (3)-(6) it follows that, for k2 odd (k, = 2‘

.

d

+

t , t : odd), (2) can be refor- mulated as

N I - I -

X;+,(Z)

= X k , ( Z ) (ZN”N1)’llkl mod (ZN2/*‘ - a’’-‘) (7a) (7b) Notice that, since (k2, N I ) = 1 k2kl mod NI is just only a permu- tation of k , . Now with k2 = 2‘

.

d

+

t and let

nl = O - X k 2 k , . k z = %;?kI(Z) mod

(z

-

w$?).

Nz/2’- I - X ; 2 , . d + l ) k l ( Z ) =

c

Y;,,,Z” ( 8 4 n = O then N z l 2 ‘ - I - -

2

y ; (W$?‘.d+r),l X ( 2 ‘ ~ d + r ) k ~ . ( 2 ‘ . d + r ) - , , = 0 i.11 N?/2‘- I =

C

( Y : ~ , , ? *

ws,)

W $ / P . (8b)

Thus, for fixed k , , the evaluation of (7b) becomes 2‘- N2/2‘-point 1-D DFT’s as described in (8). Further, since k2 is odd here, (8) can be computed very efficiently in a FFT-like manner by using the generalized F F T algorithm developed in [ 5 ] and [6]. And the

n = O

calculations of remainders mod (ZN2/*‘ - a*’-‘), required in (6) and (7a), can be accomplished by a butterfly-like structure as de- scribed in [4]. Now, let us change our attention to the case o f k2 is

even. For k2 even (k2 = 2‘

.

d

+

f, f : even). Since

N2

z

-

wk:+I

1

z N ? / z r - a2‘-‘-1

the same as (7), for k2 :even, (2) can also be reformulated as

N I - I

nl = O

.

mod ( z N ? / 2 ’ - a ? ‘ - f - l

-

Xi,? + ,)kl ( Z ) =

c

[XtIl ( Z )

w;;’]

(ZN21N1

1

‘IUk1

1

(94

- -

X ( k ~ + l ) k l . k 2 = x ; k . + l ) A ~ ( z ) mod

(z

-

w%+’),

k2:2‘

.

d

+ t , t

even. (9b) In this case, if follows that

k 2 + 1 = 2 ‘ . d + t + 1 = 2 ‘ . d + t ’

where 1 5 ? ’ : o d d 5 2‘ - 1 , and

z

- ~ ; ? . d + l ’

I

Z N ? / ? ‘ - a 2 ‘ - l ’

Hence by comparing (7) with (9), it is easy to conclude that, for k, even, the F P T algorithm for 2-D DFT’s can also be modularized following the same procedures for the case of k2 odd.

111. CONCLUSIONS A N D DISCUSSIONS

In highly parallel computational environments, multiprocessor, or VLSI systems, the system’s throughput rate depends not only on the number of total arithmetic operations but also on the regu- larity of the algorithm. The modularized F P T algorithms developed in [4] and extended in this correspondence can simplify the prob- lems of control, memory management, load balancing, etc., al- though more arithmetic operations are needed than with the origi- nal ones.

REFERENCES

H . J . Nussbaumer, “Fast polynomial transform algorithms for digital

convolution,” IEEE Trans. Acousr., Speech, Signal Processing, vol. ASSP-28, pp. 205-215, Apr. 1980.

T. K . Truong et a l . , “On the application of a fast polynomial’trans- form and the Chinese remainder theorem to compute a two-dimen- sional convolution,” IEEE Truns. Acousr., Speech, Signal Process- ing., vol. ASSP-29. pp. 91-99, Feb. 1981.

H . J. Nussbaumer and Quandalle, “Fast computation of discrete Fou- rier transforms using polynomial transforms,” IEEE Trans. Acoust., Speech, Signal Processing, vol. ASSP-27, pp. 169-181, Apr. 1979. T. K . Truong er al., “An FPT algorithm with a modularized structure for computing 2-D cyclic convolutions,” IEEE Trans. Acousr., Speech, Signal Processing, vol. ASSP-36, pp. 1540-1542. Sept. 1988. P. Corsimi and G . Frosini, “Properties of the multidimensional gen- eralized discrete Fourier transform,” IEEE Truns. Cornput., vol. C-28, pp. 819-830, 1979.

1. S . Reed et a l . , “An improved FPT algorithm for computing two- dimensional cyclic convolution,” IEEE Trans. Acousr., Speech, Sig- nal Processsing, vol. ASSP-31, pp. 1048-1050, Aug. 1983.

-T-