分数域信号与信息处理及其应用 (18).pdf

IEEE TRANSACTIONS ON SIGNAL PROCESSING,VOL.44,NO.9,SEPTEMBER 1996 2141 Digital Computation of the Fractional Fourier Transform Haldun M.Ozaktas,Orhan Ankan,Member,IEEE,M.Alper Kutay,Student Member,IEEE,and Gozde Bozdaki,Member,IEEE Abstract-An algorithm for efficient and accurate computation of the fractional Fourier transform is given.For signals with time-bandwidth product N,the presented algorithm computes the fractional transform in O(N log N)time.A definition for the discrete fractional Fourier transform that emerges from our analysis is also discussed.I.INTRODUCTION HE fractional Fourier transform 11-4 has found many T applications in the solution of differential equations 2,31,quantum mechanics and quantum optics 5-111,optical diffraction theory and optical beam propagation(including lasers),and optical systems and optical signal processing 11,131-25,swept-frequency filters 4,time-variant filtering and multiplexing I,pattern recognition 26,and study of time-frequency distributions 27.The recently studied Radon transformation of the Wigner spectrum 281-30 is also known to be the magnitude square of the fractional Fourier transform ll,31.The fractional Fourier transform has been related to wavelet transforms l,32,neural networks 32,and is also related to various chirp-related operations l,33-35.It can be optically realized much like the usual Fourier transform l,131-17,20 and,as we will show in this paper,can be simulated with a fast digital algorithm.Other applications that are currently under study or have been suggested include phase retrieval,signal detection,radar,tomography,and data compression.In this paper,we will be concerned with the digital com-putation of the fractional Fourier transform.We are not only interested in a numerical method to compute the continuous transform but also in defining the discrete fractional Fourier transform and show how it can be used to approximate the continuous transform.More precisely,we will show that the samples of the continuous time fractional Fourier transform of a function can be approximately evaluated in terms of the samples of the original function in O(N log N)time,where N is the time-bandwidth product of the signal.In many of the above-mentioned applications,it is possible to improve performance by use of the fractional Fourier transform instead of the ordinary Fourier transform.Since in this paper we show that the fractional transform can be Manuscript received February 3,1995,revised January 9,1996.The associate editor coordinating the review of this paper and approving it for publication was Dr.Bruce Suter.The authors are with Electrical Engineering,Bilkent University,06533 Bilkent,Ankara,Turkey.Publisher Item Identifier S 1053-587X(96)06680-puted in about the same time as the ordinary transform,these performance improvements come at no additional cost.To give one concrete example,in some cases,filtering in a fractional Fourier domain,rather than the ordinary Fourier domain,allows one to decrease the mean square error in estimating a distorted and noisy signal 36.In Section 11,preliminaries about the fractional Fourier transform are given,including its relation to the Wigner distribution.Section I11 reviews some straightforward yet inefficient methods of computing the fractional Fourier trans-form.Fast computational algorithms are presented in Section IV,and simulation examples are given in Section V.Some alternate methods are discussed in Section VI to better situate the suggested algorithm among other possibilities.Section VI1 deals with the issue of defining the discrete fractional Fourier transform in some detail.The remainder of the paper constitutes concluding sections.11.PRELIMINARIES A.The Fractional Fourier Transform Let 3Tf)(x)denote the Fourier transform of f(s).Integral powers 3-7 of the operator 3 e 3 may be defined as its successive applications.Then,we have 3f(s)=f(-x)and 34f(x)=f(s).The ath-order fractional Fourier transform Ff(x)of the function f(x)may be defined for O(a(xa)=lfa(xa)l2(1 1)where R,is the Radon transform operator.724 takes the integral projection of the 2-D function Wf(z,U)onto an axis making angle=with the z axis.We will refer to this axis as the z,axis or the ath fractional Fourier domain(Fig.1).The xo axis is the usual time domain z,and the z 1 axis is the usual frequency domain U.Notice that(7)and(8)are special cases of this equation.In general,the projection of the Wigner distribution on the ath fractional Fourier domain gives the magnitude squared of the ath fractional Fourier transform of the original function.There is actually nothing special about any of the continuum of domains;the privileged status we assign to the time and frequency domains can be interpreted as an arbitrary choice of the origin of the parameter a.All of the fractional transforms,including the 0th transform(the function itself),are different functional representations of an abstract signal in different domains.The unitary transformation between these different representations is the fractional Fourier transform 11.C.Compactness in the Time Domain,the Frequency Domain,and Wigner Space A function will be referred to as compact if its support is so.The support of a function is the subset of the real axis in which the function is not equal to zero.In other words,a function is compact if and only if its nonzero values are confined to a finite interval.It is well known that a function and its Fourier transform cannot be both compact(unless they are identically zero).In practice,however,it seems that we are always working with a finite time interval and a finite bandwidth.This discrepancy between our mathematical idealizations and Authorized licensed use limited to:BEIJING INSTITUTE OF TECHNOLOGY.Downloaded on April 06,2021 at 06:35:48 UTC from IEEE Xplore.Restrictions apply.OZAKTAS et al.:DIGITAL COMPUTATION OF THE FRACTIONAL FOURIER TRANSFORM 2143 the real world is usually not a problem when we work with signals of large time-bandwidth product.The time-bandwidth product can be crudely defined as the product of the temporal extent of the signal and its(double-sided)bandwidth.It is equal to the number of degrees of freedom and the number of complex numbers required to uniquely characterize the signal among others of the same time-bandwidth product.We will assume that the time-domain representation of our signal is approximately confined to the interval-Atla,At121 and that its frequency-domain representation is confined to the interval-Af/2,Af/Z.With this statement,we mean that a sufficiently large percentage of the signal energy is confined to these intervals.For a given class of functions,this can be ensured by choosing At and Af sufficiently large.We then define the time-bandwidth product N z AtA f,which is always greater than unity because of the uncertainty relation.Let us now introduce the scaling parameter s with the dimension of time and introduce scaled coordinates x=t/s and U=f s.With these new coordinates,the time and frequency domain representations will be confined to intervals of length Atls and A f s.Let us choose s=dm so that the lengths of both intervals are now equal to the dimensionless quantity d m,which we will denote by Ax.In the newly defined coordinates,our signal can be represented in both domains with N=Ax2 samples spaced Axp1=1/fi apart.From now on,we will assume that this dimensional normal-ization has been performed and that the coordinates appearing in the definition of the fractional Fourier transform,Wigner distribution,etc.,are all dimensionless quantities.If the representation of the signal in the ath domain is confined to a certain interval around the origin,the Wigner distribution will be confined to an infinite strip perpendicular to the xa axis defined by that interval.Thus,assuming that the representation of the signal in all domains is confined to an interval of length Ax around the origin,is equivalent to assuming that the Wigner distribution is confined within a circle of diameter Ax,With this,we mean that a sufficiently large percentage of the energy of the signal is contained in that circle,For any signal,this assumption can be justified by choosing Ax sufficiently large.(Of course,it is in our interest to choose it as small as possible to reduce computational complexity.)For convenience,we will require Ax to be an integer.Signals whose energy are not concentrated around the origin of Wigner space might be treated more efficiently than by simply choosing Ax large enough to include them,but this extension to the“bandpass”case is not treated in this paper.D.The Discrete Fourier Transform The discrete Fourier transform(DFT)is a mapping RN t RN.The matrix elements of this transformation are defined as(12)The DFT is related to the continuous Fourier transform as follows 39:Assume that a function f(x)and its Fourier transform f l (U)are both confined to the interval-Ax/2,Ax/2.Then,the N=Ax2 samples of the Fourier transform may be found by taking the DFT of the N samples of the original function,where the sample spacing in both domains is l/fi=l/Ax.A more precise statement of the above belongs to a class of relations known as Poisson formulas 39:111.METHODS OF COMPUTING THE CONTINUOUS FRACTIONAL FOURIER TRANSFORM The defining integral(1)can be rarely evaluated analyti-cally;therefore,numerical integration is called for.Numerical integration of quadratic exponentials,which often appear in diffraction theory,require a very large number of samples if conventional methods are to be employed,due to the rapid oscillations of the kernel.The problem is particularly pronounced when a is close to 0 or 52.If we assume both the function and its Fourier transform to be confined to a finite interval,we can walk around this difficulty as follows:If a E 0.5,1.5 or a E 2.5,3.5,we evaluate the integral directly.If a E(-0.5,0.5)or a E(1.5,2.5),we use the property Fa=F1FT”-l,noting that in this case,the(a-1)th transform can be evaluated directly.(Essentially similar issues are discussed in 28-30.)Another method of evaluating(1)would be to use the spectral decomposition of the kernel(see(5)l,14-161.This is equivalent to first expanding the function f(x)as C?o c,$,(x),multiplying the expansion coefficients c,re-spectively,with e-1nnr/2,and summing the components.Although both ways of evaluating the fractional Fourier transform may be expected to give accurate results,we do not consider them further since they take O(N2)time.Iv.FAST COMPUTATION OF THE FRACTIONAL FOURIER TRANSFORM The fractional Fourier transform is a member of a more general class of transformations that are sometimes called linear canonical transformations or quadratic-phase transforms 20.Members of this class of transformations can be broken down into a succession of simpler operations,such as chirp multiplication,chirp convolution,scaling,and ordinary Fourier transformation.Here,we will concentrate on two particular decompositions that lead to two distinct algorithms.A.Method Z First,we choose to break down the fractional transform into a chirp multiplication followed by a chirp convolution followed by another chirp multiplication 171,391.In this approach,we assume a E-1,11.Manipulating(l),we can write(14).fa(x)=exp-inz2 tan(dP)g(x),Authorized licensed use limited to:BEIJING INSTITUTE OF TECHNOLOGY.Downloaded on April 06,2021 at 06:35:48 UTC from IEEE Xplore.Restrictions apply.2144 IEEE TRANSACTIONS ON SIGNAL PROCESSING,VOL.44,NO.9,SEPTEMBER 1996 0 0 g(x)=A4 1 exp z/?(x-d)g(d)d d,(15)(16)-oo g(x)=exp-ii7x2 tan($/2)f(x)where g(x)and g/(z)represent intermediate results,and/3=In the first step(see(16),we multiply the function f(x)by a chirp function.As we discuss in the Appendix,the bandwidth and time-bandwidth product of g(x)can be as large as twice that of f(x).Thus,we require the samples of g(x)at intervals of 1/2Ax.If the samples of f(x)spaced at l/Ax are given to begin with,we can interpolate these and then multiply by the samples of the chirp function to obtain the desired samples of g(z).There are efficient ways of performing the required interpolation 40.The next step is to convolve g(x)with a chirp function,as given in(15).To perform this convolution,we note that since g(x)is bandlimited,the chirp function can also be replaced with its bandlimited version without any effect,that is csc 4.0 0 g(x)=A4 exp i;.P(x-z)g(x)dx J-00 J-00 where A X H(v)exp(i27rvz)du(18)and where 1 H(v)-e2w/4 exp(-i;.v/P)(19)Jp is the Fourier transform of exp(irr/3x2).It is possible to express h(x)explicitly in terms of the Fresnel integral defined as(20)Now,(15)can be sampled,giving This convolution can be evaluated using a fast Fourier trans-form.Then,after performing the last step(see(14),we obtain the samples of fa(z)spaced at 1/2Ax.Since we have assumed that all transforms of f(x)are bandlimited to the interval-Ax/2,Ax/2,we finally decimate these samples by a factor of 2 to obtain samples of fa(z)spaced at l/Ax.Overall,the procedure starts with N samples spaced at l/Ax,which uniquely characterize the function f(z),and returns the same for f a (x).If we let f and fa denote column vectors with N elements containing the samples of f(x)and fa(z),the overall procedure can be represented as f a =G f,(22)F4=DAHl,AJ.(23)Here,D and J are matrices representing the decimation and interpolation operations 40.A is a diagonal matrix that corresponds to chirp multiplication,and Hz,corresponds to the convolution operation.We notice that 3:allows us to obtain the samples of the ath transform in terms of the samples of the original function,which is the basic requirement for a definition of the discrete fractional Fourier transform matrix.If we are merely interested in computing and plotting the fractional Fourier transform of a given continuous f(z),then the decimation and interpolation steps can be eliminated.Note that the described algorithm works for-1 5 a 5 1.If a lies outside this interval,we can use the properties F4f(x)=f(x)and F2f(z)=f(-z)to easily obtain the desired transform.B.Method I1 We now turn our attention to an alternative method that does not require Fresnel integrals.The defining equation for the fractional Fourier transform can be put in the form 00 e-z2Tpzx f(.)dx J L Fat(.)=A+ezTcuz2(24)where a=cot q!J and/3=csc 4.We are again assuming that the Wigner distribution of f(.)is zero outside a circle of diameter Ax centered around the origin.(This was_ discussed in detail in Section 11-C.)Under this assumption,and by limiting the order a to the interval 0.5 5 la1 5 1.5,the amount of vertical shear in Wigner space resulting from the chirp modulation is bounded by Ax/2.Then,the modulated function eznaxt2 f(x)is bandlimited to A x in the frequency domain.Thus,ezTcux2 f(d)can be represented by Shannons interpolation formula sinc 2Ax x-)(25)(2Ax where N=(Ax)2.The summation goes from-N to N since f(d)is assumed to be zero outside-Ax/2,Ax/2.By using(25)and(24)and changing the order of inte