A new discrete fractional Fourier transform based on constrained eigendecomposition of DFT matrix by Largrange multiplier method

A NEW DISCRETE FRACTIONAL FOURIER TRANSFOIRM BASED O N

CONSTRAINED EIGENDECOMPOSITION

OF DFT MATRIX BY LARGRA.NGE

MULTIPLIER METHOD

Soo-Chang Pea Chien-Cheng Tseng

Department of Electrical Engineering, National Taiwan University, Taipei, Taiwan, R. 0. C. Email address: peiOcc.ee.ntu.edu.tw

ABSTRACT

This paper is concerned with the definition of the discrete fractional Fourier transform (DFRFT). First, an eigende- composition of the discrete Fourier transform (DFT) matrix is derived by sampling the Hermite Gauss functions which are eigenfunctions of the continuous Fourier transform and by performing a novel error removal procedure. Then, the result of the eigendecomposition of the DFT matrix is used to define a new DFRFT. Finally, a numerical example is

illustrated to demonstrate the proposed DFRFT is a better approximation to the continuous fractional Fourier trans- form than the conventional defined DFRFT.

1. INTRODUCTION

In recent years, many researchers have paid attention to the investigation of a new signal processing tool called frac- tional Fourier transform (FRFT). This transform has found many applications in the solution of differential equation, quantum mechanics and quantum optics, and optical sys- tems and optical signal processing, swept-frequency filter, time-variant filtering and multiplexing,. pattern recognition, and study of time-frequency distribution [l]. Besides, the FRFT has been proved to relate to other signal analysis tools, such as Wigner distribution, neural network, wavelet transform and various chirprelated operations [2]. Sever- al useful properties of FRFT are currently under study in signal processing community [3].

So far, many methods for implementing FRFT has been developed. However, most of them are to utilize the optical instruments or numerical integration. Because the FRFT is a potentially useful tool for signal processing, the direct computation of FRFT in digital computer has become an important issue. Basically, the computation of the discrete fractional Fourier transform (DFRFT) needs to obey ad- ditivity property and similarity condition. The additivity property means that application of the transform with an- gular parameter a followed by an application of the trans- form with angular parameter

p

is equivalent to the appli- cation of the transform with angular parameter Q

+

p.

The

similarity condition means that the transform results of D- FRFT must be similar to those of the continuous FRFT. In [4], a method for digital computing FRFT was proposed, but their method does not obey the additivity property and the signal can not be recovered from its transform result-

s. In [5], another DFRFT is defined, but this definition does not provide the similar transform results as those of

continuous case. The purpose of this paper is to present a new DFRFT which obey additivity property and similarity condition simultaneously.

0-8186-7919-0/97 $10.00 0 1997 IEEE

2. EIGENDECOMPOSITION OF THE

DFT

MATRIX

2.1 The eigenvalues and eigrenvectors of

DFT

matriz Now, we briefly review the properties of the t?igenvalues and eigenvectors of the DF'r matrix F whose elements de- fined by

' p ) ) O < n , k < N - l

F,,&

=

-

1 (cos(-) N 27rkn

-

j sin(-

JT

~~

(1) From the results in [6][7], the properties can be summarized as the following two facts:

Fact 2 The eigenvalues of

F

are { 1, -1, j , - j } and its mul- tiplicities are listed below:

N

4m I m + l I rn

[ M ul. of 1

I

Mul

m f 1 m + l m Fact 2 Let w =

9

and matirix

S

be

0

l

1

...

1

;

2cos(w) . .

. - .

0

...

1

...

i

1 0 * - e ~ c o s ( ( N

-

1 ) ~ : ) then it can be shown that

FS

=

SF.

Because matrix

S,

with distinct eigenvalues, commutes with

F,

the eigenvectors of

S

will be the desired set of eigenvec- tors of F. Note that S is a real and symmetlic matrix,

so its eigenvectors will be leal and orthogonal. Although Fact 2 can help us to find a real orthogonal eigenvector set of the matrix F, this solutxon is not unique because any

linear combination of the eigenvetors which correspond to the same eigenvalue is also im eigenvector. Thus, there ex- ist infinite eigendecomposition forms of the DFT matrix. If we use the eigendecompoition of the DFT matrix

F

to define the discrete fractional Fourier transform (DFRFT), then we have infinite choice. However, under the condition that transform results of DFRFT needs to be similar to those of continuous FRFT, the eigendecomposition of DFT matrix must be found trickly. In the following, we w i l l de- rive an eigendecomposition form by sampling the Hermite Gauss functons which are ithe eigenfunctions of the con- tinuous Fourier transform and by performing a novel error removal propcedure. Using the proposed decomposition to

define DFRFT, the transform results w i l l obey simialarity condition.

2.2 An eigendecomposition of

DFT

matriz

The usual continouos Fourier transform pair is defined as

.

rea

.

rea

It can be shown that the eigenfunctions of the Fourier trans- formation operator are Hermite Gauss function H , ( t ) e - g ,

where H , ( t ) are the Hermite polynomials of order m. We thus have

Now, we will use this equation to derive an approximate eigendecomposition of the DFT matrix.

Our

derivation is mainly based on the following two facts:

Fact 9: If the sequence gm(n) is obtained by sampling the Hermite Gauss function H , ( t ) e T with sampling interval T =

fi,

i.e.,

- i 2

+

gm(n) = H,(nT)e (4)

then it can be shown that

7 r-1

for sufficiently large N .

Because the degree of Hermite polynomial Hm(t) is m, the decay rate of the Hermite Gauss function H , ( t ) e T is proportional to tme-" for sufficiently large t. And, the larg- er order m is, the slower decay rate Hermite Gauss function has. Thus, when order m becomes large, the approximation in eq(5) become worse.

Fact 4: If the sequence gm(n) defined in the range [0, N- 13 is obtained by shifting Hermite Gauss samples gm(n) de- fined in the range

[+,

%

-

l] in the following way:

-1'

then it can be shown that the DFT of the &(n) can be approximated by (-j)'"#,(k), i.e.,

7 N-1

n = O

for sufficiently large N .

From the Fact 4, it is clear that &(n) are the approxi- mate eigenfunctions of the discrete Fourier transform. Be- cause the Hermite Gauss functions are orthogonal each oth- er for different orders, the sequences gml(n) and g,2(n) are approximately orthogonal for m l

#

m2, i.e.,

N-1

n=O

Let us define the vectors vm as follows:

then eq(7) means that

where V, =

&

is normalized version of the vector vm.

Thus, V, is an approximate eigenvector of the DFT matrix

F corresponding to the eigenvalue (-j),. Although the approximate expression in eq(l0) is valid for any order m, the DFT matrix F with size N x

N

only has N eigenvectors

whose eigenvalues need to satisfy the multiplicity property in Fact 1. Thus, we are required to select N orders denoted by the set 0 = { m l , m z , . - - , m N } (ml

<

m2

<

.-.

<

mN) to construct an eigendecomposition of the matrix F. Two rules of the selection in this paper is listed as follows: ( I ) The set {(-j)m1,(-j)m2,.--,(-j)"'N} formed by eigenvalues must satisfy the multiplicity property in

Fact

1.

(2) The approximation error Il(-j)"iiimi -FVmi

11

must be less than the error Il(-j)"V,

-Fo,ll

if m is not in the set

e.

Because the approximation error 1[(-j)Wm

-

FV.,II be-

comes large when order m increases, a suitable choice of set

Q which obeys two rules is described in the following table:

N I \Y

-

-

{ m l , . . . , m N } 4n I 0.1.2..

..

.4n

-

2.4n l 4 n + 1 I 0.1.2 . . . e . 4 n - 1 1 . 4 n I 4n+2

I

4n+3

I

0,1,2,... , 4 n

+

1,4n

+

2 0, 1 , 2 , .

.

-

,4n, 4n

+

2

Based on this choice, an approximate eigendecomposition of the DFT matrix

F

is given by

N

F x ~ ( - j ) m i V , , t ~ i (11)

i = l

In order to remove the error in this decomposition, an eigen- vector calibration procedure is developed as follows. As- be corrected into the eigenvector set {uml, ,',U

-

,UmN}

and the vectors from uml to u,,-~ have been obtained. Then, the eigenvector ,,U is found by minimizing the squared error (U,,

-

om,):

subjected to two prescribed constraints which are the eigenvector constraint F'u,, =

(-j)"ku,, and the orthogonal constraint ufumr = 0 for

i = m l ,

- -

, m k - ~ . After some maniputation, two con- straints can be rewritten as matrix form below:

sumed that the eigenvector set {Vml,

vma,

-

1 ,B"} will

C m k - l U m k = (12)

where the matrix C,,-, is given by Real(F

-

( -j),, I)

Img(F

-

( - j ) * k I )

The notation Real(.) and Img(.) denote the real part and imaginary part of a complex matrix, and I is identity mar- tix. Using the QR decomposition, the matrix

c,,-,

can be rewritten as

Substitute eq(14) into eq(12), the eq(l2) reduces to

Rmk-,umk = 0 (15)

If the rank of matrix C,,-, is equal to r, the size of the ma- trix R,,-, is r x

N.

Now, using the well-known Largrange multiplier method, the solution of this constrained opti- mization problem is given by

Umk = (I

-

~ k _ l ( R m k - l R k ~ _ , ) - l R m k _ l ) vmk (16) Finally, the entire eigenvector calibration procedure is sum-

marized as follows: Given DFT matrix F and the ap- proximate eigenvector set

{vml,

vma,.

. ,

v T N }

we take the following steps to compute the exact agenvector set Step 1: Let matrix C m l be [I

-

Real(F)',Img(F)'It and use eq(16) to find the vector uml. Note that we normalize uml to unit norm. Set k = 2.

Step 2.- Perform the following two computations:

{Uml U,,, *

(1) Use the eq(13)(14) to compute the matrix Rmk. (2) Use the eq(16) to calculate the vector Umk and nor- malize it to unit norm.

Step 3: Let k = k

+

1. If

k

>

N, stop the procedure. Oth- erwise go to Step 2.

After this calibration, the exact eigendecomposition of the DFT matrix F is given by

N

i = l

The unique feature of this eigendecomposition is that the shape of the eigenvector is similar to the shape of the Her- mite Gauss functions which is the eigenfunction of the con- tinuous Fourier transform. In the next section, we will use this decomposition to define a discrete fractional Fourier transform.

3. NEW D E F I N I T I O N O F DFRFT

The DFRFT of the data vector x is defined by 2p

Ra[x] = F

,,

x

Since $th power of the DFT matrix F can be calculated from its eigendecomposition by taking the F t h power for its eigenvalues, the matrix F* is given by

N

i = l

Because ( - j ) m i % = e-jmia

,

the eigenvalues of the new

transform matrix F* are consistent with those of the con- tinuous FRFT. Moreover, the eigenvectors umi are obtained by sampling Hermite Gauss functions with an error removal procedure, so the eigenvectors of new DFRFT are similar

to those of the continuous E'RFT. Due to these two agree- ments, the transform result of our DFRFT will be similar to that of continuous FRFT. In order to demonstrate the advantage of our DFRFT, we consider the FRFT of impulse function 6(t). The continuous FRFT of this special signal has the closed form formula given by

Fig.1 shows the continuous l?RFT of the impulse signal for various angles a. For compzrision, we examine the DFRFT of the pnit sample function (defined by

1 for n = 0

0 otherwise 6(n) =

Fig.2 shows the trasnform result of the DFRFT defined by Santhanam and McCleUan [.5], and Fig.3 depicts the result of our DFRFT for N = 36. It is clear our result is more similar to that of continuous case than the result of the conventional DFRFT.

4. C O N C L U S I O N S

In this paper, a new definition of the discrete fractional Fourier transform (DFRFT) based on an eigendecomposi- tion of DFT matrix has been presented. The eigencomposi- tion of the DFT matrix is derived by sampling the Hermite Gauss functions which are eigenfunctions of the continuous fractional Fourier transform imd by performing a novel error removal procedure. A numerical numerical example is illus- trated to demonstrate the proposed DFRFT is a better a p proximation to the continuous fractional Fourier transform than the conventional defined DFRFT. However, the com- plexity for implementing DI'RFT is O ( N Z ) which is same as that of DFT. Thus, it is interesting to develop a fast algorithm to compute DFRFT.

REFERENCES

H.M. Ozaktas and

B.

Barshan, "Convolution, filtering, and multiplexing in fractional Fourier domains, and their relationship to chirp and wavelet transform", J . Opt. Soc. Amer. A, vol. 11, pp. 547-559, 1994.

S.Y.

Lee and H.H. Szu, "Fractional Fourier transforms, wavelet transform, and adaptive neural network", Opti- cal Engineering, vo1.33, pp.2326-2330, 1994.

X.G. Xia, "On bandlimited signals with fractional Fouri- er transform

",

IEEE Signal Processing Letter, vo1.3, H.M. Ozaktas, O.Arikan, A. Kutay and Gozde: Bozdagi, "Digital computation of the fractional FourieI transfor- m", IEEE Trans. Signal

process.,

vo1.44, pp.2141-2150, Sept. 1996

B.

Santhanam and J.H. AkClellan, "The DRFT

-

A r w tation in time frequency space", Proceedings o j ICASSP, pp. 921-924 May 1995.

J.H. McClellan and T.W. Parks, "Eigenvalue and eigen- vector decomposition of tlhe discrete Fourier transform",

IEEE Trans. Audio and Electroacoustics., vol-20, pp. 66-74,

Mar.

1972.

B.W. Dickinson and

K.

Steiglitz, "Eigenvectors and functions of the discrete Fourier transform", IEEE Tmn-

s. Acoust., Speech, and Signal Process., vol. 31D, pp. 25- 31, Feb. 1982.

pp. 72-74, March 1996.

0 IN I 1N I n , <" 0 00 IN I P 3968