A novel description and mathematical analysis of the Fractional Discrete Fourier Transform
aa r X i v : . [ m a t h . G M ] S e p A novel description and mathematical analysisof the Fractional Discrete Fourier Transform
Evan M. ZayasPhysics Department, Massachusetts Institute of TechnologyOctober 1, 2019
Abstract
I discuss the nature of a Fractional Discrete Fourier Transform (FrDFT)described algorithmically by a combination of chirp transforms and or-dinary DFTs. The transform is shown to be consistent with a continu-ous two-dimensional rotation between the time and frequency domains.I further present a new closed-form expression for the transformationmatrix and some preliminary analysis of its properties.
The Discrete Fourier Transform (DFT) has long been established as a prin-cipal tool for frequency analysis of signals. In particular, the transform de-composes a vector ~x in the time domain to its frequency spectrum ~f accordingto the transformation matrix: f j = N − X k =0 B jk x k (1)with B jk = 1 √ N e − πiN jk for 0 ≤ j, k < N (2)where N is the size of the data vector and i ≡ √−
1. Furthermore, it is wellunderstood that performing the (forward) DFT again on frequency-domain1ata returns it to the time domain but with the ordering reversed comparedthe original vector. We may verify this explicitly with the square of the trans-formation matrix above: (cid:0) B (cid:1) jk = N − X l =0 B jl B lk = 1 N N − X l =0 e − πiN l ( j + k ) (3)= j + k ) mod N = 00 else (4)Thus, the matrix element is nonzero only when k = N − j and j >
0, orwhen instead j = k = 0: B = . . . . . . . . . . . . (5)Acting with this matrix on an input vector ~x reverses the order of allelements except the first. It is common in such discussions to periodicallyextend the finite vector ~x , so that x j + N = x j for any integer j . Then, theaction of two DFTs in series is equivalent to the parity flip operation x j → x − j .Indeed, one can easily check from the above equation that B is a square-rootof the identity, consistent with the parity operator. This also immediatelyimplies that the DFT itself is a fourth root of the identity.From these properties, we can form a consistent picture of the FourierTransform in the two-dimensional space of time and frequency. The DFT isinterpreted as a rotation by π/
2, for example from the time domain (1 ,
0) to theorthogonal frequency domain (0 , B and B which we have just discussed follow sensibly as well.A monotone signal exists at a constant frequency for all time: a horizontalline f ( t ) = f in this 2D description of the time and frequency domains. Itsfrequency spectrum is of course a delta function (in the large N limit), whichindicates that the signal is maximally localized in f ; or equivalently, that arotation of the coordinate space ( t, f ) by π/ t . This2escription of the DFT provides a natural extension to the Fractional DFT(FrDFT), where a vector is instead rotated by an arbitrary angle α . This makesaccessible a continuum of mixed time/frequency domains parametrized by α .The special property of a monotone signal is then extended to any line in thecoordinate space: f ( t ) = f + qt , called a chirp signal. An appropriate choiceof α fully localizes this signal just as the ordinary DFT does to a pure tone.The FrDFT therefore has potential applications in the study and processingof not only constant-frequency tones, but also chirped tones or locally-chirpedtones including those reconstructed in experiments such as Project 8 [1] andLIGO [2].In this work, I begin by studying an algorithm to perform the FrDFTproposed by Garcia et al [3] in the context of this 2D rotation. I furtheruse the Garcia algorithm to construct a simple closed-form expression for theFrDFT transformation matrix F ( α ). Lastly, I show that this expression for F ( α ) behaves consistently in the limiting cases of α = 0 and α = π/ The actual implementation of a FrDFT has been the subject of some consid-erable study for a variety of applications; of most interest to us is an algorithmpublished by Garcia et al [3]. The algorithm consists only of two kinds of op-erations: (a) an ordinary DFT, and (b) multiplication by a quadratic phase.This is especially intriguing for two reasons: first, it is clearly related to chirpsignals as the quadratic phase multiplication takes a plane wave exp ( − iωt ) toa chirp signal exp ( − i ( ω + qt ) t ). Second, since the quadratic phase multipli-cation is an O ( N ) operation, this algorithm can be executed just as fast asthe ordinary DFT, in O ( N log N ) time as is well established by Cooley andTukey [4]. By contrast, direct computation of all matrix elements F jk ( α ) isnecessarily O ( N ) or worse.The algorithm consists of five steps:1. Multiply by the quadratic phase: x j → exp (cid:0) − iπN q j (cid:1) x j
3. Perform the DFT3. Multiply by the quadratic phase: x j → exp (cid:0) − iπN q j (cid:1) x j
4. Perform the inverse DFT5. Multiply again by the quadratic phase: x j → exp (cid:0) − iπN q j (cid:1) x j where the chirp rates q and q are related to the rotation angle α by: q = tan ( α/
2) (6) q = sin ( α ) (7)To describe this algorithm in the ( t, f ) coordinate space, we must find a2x2 matrix that describes the action of each step; in fact, we have alreadydone so for the DFT which takes t → f and f → − t via a π/ t ′ f ′ ! = −
11 0 ! tf ! (8)The quadratic phase multiplication takes a constant tone with frequency f to a chirp tone with frequency f + qt and preserves t → t : t ′ = t (9) f ′ = qt + f (10)and from the above equations we may simply read off the desired matrixelements. With the quadratic phase shift denoted A ( q ) and the ordinary DFTdenoted B , we have: A ( q ) = q ! (11) B = −
11 0 ! (12)Note that with the above expression for B we can once again easily checkthat the DFT is a fourth root of the identity. The FrDFT from the Garciaalgorithm is: F ( α ) = A ( q ) B − A ( q ) BA ( q ) (13)4nd the claim is that with q and q related to α by Equations 6 and 7, thematrix F ( α ) is a rotation by precisely the same angle α . I will next explicitlyshow this claim to be true. F ( α ) = q ! − ! q ! −
11 0 ! q ! = − q q − q q (2 − q q ) 1 − q q ! (14)With Equation 7 we can immediately replace F ( α ) = − q with − sin ( α ).The diagonal elements both simplify to:1 − q q = 1 − tan ( α/
2) sin ( α )= 1 − sin ( α/ α/
2) 2 sin ( α/
2) cos ( α/ − ( α/ ⇒ − q q = cos ( α ) (15)and all that remains is to massage the last element F : q (2 − q q ) = q (2 q /q − q )= sin ( α ) (cid:18) α/ α/
2) cos ( α/ − tan ( α/ (cid:19) = sin ( α ) (cid:0) sec ( α/ − tan ( α/ (cid:1) ⇒ q (2 − q q ) = sin ( α ) (16)Indeed, we have shown that this representation of F ( α ) as a product oflinear shifts and π/ F ( α ) = cos ( α ) − sin ( α )sin ( α ) cos ( α ) ! (17) N -dimensionaldata space We have studied the Garcia algorithm in the coordinate space of time andfrequency, and from it recovered the 2D rotation matrix. This solidifies our5nderstanding of the DFT in this space, and provides good confidence thatthe FrDFT described by this algorithm is consistent with that understanding.Next, we return to the N -dimensional space of an input data vector to onceagain find an explicit expression for F ( α ).The DFT in this space is already given by Equation 2, and we may similarlyread off the matrix elements A jk ( q ) from the quadratic phase transformation:˜ x k = exp (cid:18) − iπN qk (cid:19) x k = N − X j =0 exp (cid:18) − iπN qjk (cid:19) x j δ jk (18) ⇒ A jk ( q ) = exp (cid:18) − iπN qjk (cid:19) δ jk (19)where δ jk is the Kronecker delta. For brevity, I now introduce the 2 N th rootof unity ζ = e − iπ/N ; the simplified expressions for A jk and B jk are: A jk ( q ) = δ jk ζ qjk (20) B jk = 1 √ N ζ jk (21)and an explicit expression for F jk ( α ) follows again from the matrix product asin Equation 13: F jk = N − X l,m,n,p =0 A jl ( q ) B − lm A mn ( q ) B np A pk ( q )= 1 N N − X l,m,n,p =0 ζ q jl ζ − lm ζ q mn ζ np ζ q pk δ jl δ mn δ pk = 1 N N − X m =0 ζ q j ζ − jm ζ q m ζ mk ζ q k = 1 N ζ q ( j + k ) N − X m =0 ζ q m +2 m ( k − j ) (22)Substituting in ζ we have the full expression: F jk = 1 N exp − iπN q ( j + k ) N − X m =0 exp − iπN (cid:16) q m + 2 m ( k − j ) (cid:17) (23)This is perhaps a simpler closed-form expression for a FrDFT matrix thanpresently exists in the literature. 6 Limiting behavior of F ( α ) Next, we will study this new expression for the FrDFT matrix in the twolimiting cases of the rotation angle to prove the following results:(a) For α = 0, the transformation matrix F jk simply reduces to the identity: F jk (0) = δ jk .(b) For α = π/
2, the fractional transform corresponds to an ordinary DFT: F jk ( π/
2) = B jk up to an overall phase.It is quite straightforward to show that with α = 0 we recover the identity;both q and q are zero in this case, and thus A jk = δ jk . The algorithmthen reduces to a DFT followed by an inverse DFT, which trivially gives theidentity. We can also arrive at this result directly from Equation 23 with the q and q terms removed: F jk = 1 N N − X m =0 ζ m ( k − j ) = δ jk (24)To prove (b), we set α = π/ q = q = 1 and simplify F jk : F jk = 1 N ζ j + k N − X m =0 ζ m +2 m ( k − j ) (25)Combining both factors of ζ , we have in the exponent: j + k + m + 2 m ( k − j ) = ( m + k − j ) + 2 jk (26)The last term is precisely what we need to isolate B jk : F jk = B jk √ N N − X m =0 ζ ( m + k − j ) (27)Thus, up to a normalization of √ N and an overall phase, we must onlyshow that the sum in the above equation reduces to 1, i.e. it is independentof j and k . To accomplish this, I will first prove and then apply a trick withthe roots of unity. 7onsider a sum of the form: S k = k + N − X s = k ζ s (28)for arbitrary integer k . The difference between S k +1 and S k is: S k +1 − S k = ζ ( k + N ) − ζ k (29)since all of the remaining terms (when k < s < N ) appear in both sums. Weexpand the exponent and recall that ζ Ns = 1 for any integer s , or equivalentlythat ζ Ns = 1 for any even integer s . Then, provided only that N is even, thedifference between the sums is: S k +1 − S k = ζ k ζ Nk ζ N − ζ k = ζ k − ζ k ⇒ S k +1 = S k (30)The result is a very general and somewhat surprising property of thesesums: every S k is in fact the same. It follows by induction from Equation 30that S k = S for every k , which explicitly removes the k dependence: k + N − X s = k ζ s = N − X s =0 ζ s for any integer k (31)This completes the trick. Stated another way, we may freely displace thequantity s in Equation 28 by any integer s → s + δs and the sum S k = S remains unchanged. It is ultimately a simple consequence of the fact that ζ s is periodic in s with period N , and thus any sequence of N terms will sum tothe same value.Returning now to Equation 27, with this trick the majority of the proofbecomes trivial. The remaining sum is exactly of the desired form with s = m + k − j , which we may simply replace with s = m : F jk = B jk √ N N − X m =0 ζ m (32)8nd the j, k dependence is now completely isolated to B jk : F jk = σ B jk with σ = S √ N (33)All that remains is to show σ has unit norm. This is easily accomplishedwith the same trick, which we note is equally valid with ζ replaced by ζ ∗ . Thenorm squared of σ is: | σ | = 1 N N − X j =0 ζ j N − X k =0 ζ − k (34)We displace the exponent in the sum over k by j : | σ | = 1 N N − X j =0 ζ j N − X k =0 ζ − ( k + j ) = N − X k =0 ζ − k N N − X j =0 ζ − jk ! = N − X k =0 ζ − k δ k ⇒ | σ | = 1 (35)Thus, up to an overall phase σ = e iθ we have recovered the ordinary DFTmatrix and obtained the desired result: F jk ( π/
2) = B jk (36)This completes the proof of (b). The work presented here has provided new insight into the nature of Frac-tional DFTs, but leaves a large amount open to potential future study. Usingthe Garcia algorithm, we have established a curious relationship between chirptransforms, rotations in two dimensions, and Fourier analysis. There is per-haps much more to be understood about the interconnected nature of these9perations. Equation 23 describes a FrDFT transformation matrix in closedform, from which each matrix element can be calculated from the sum of N terms. This may also provide a new avenue to study the practical implemen-tation of fractional transforms, and further analysis of this matrix like thatin Section 4 to answer additional questions about its properties and otherbehavior. Acknowledgements
I would like to thank Sruthi Narayanan for insightful discussions pertainingto much of this work, and in particular for an initial proof of the roots of unitytrick which was instrumental in reaching the conclusions of Section 4. Inaddition, my research has been supported in part by NSF Award No. 1806251.
References [1] A. A. Esfahani et al.
Determining the neutrino mass with cyclotron ra-diation emission spectroscopy − Project 8 . Journal of Physics G: Nuclearand Particle Physics , 5. (2017)[2] E. Chassande-Motti et al. Detection of GW bursts with chirplet-like tem-plate families . Classical and Quantum Gravity , 19. (2010)[3] J. Garcia et al. Fractional-Fourier-transform calculation through the fast-Fourier-transform algorithm . Applied Optics , 35. (1996)[4] J. W. Cooley and J. W. Tukey. An algorithm for the machine calculationof complex Fourier series . Mathematics of Computation19