DFT and FFT: An Algebraic View 1.9 2008/05/27 09:58:50 GMT-5 2008/09/11 16:47:12.933 GMT-5 Markus Pueschel pueschel@ece.cmu.edu Markus Pueschel pueschel@ece.cmu.edu C. Sidney Burrus csb@rice.edu Daniel Collins Williamson dwilliamson1285@gmail.com by Markus Pueschel, Carnegie Mellon UniversityIn infinite, or non-periodic, discrete-time signal processing, there is a strong connection between the z-transform, Laurent series, convolution, and the discrete-time Fourier transform (DTFT) . As one may expect, a similar connection exists for the DFT but bears surprises. Namely, it turns out that the proper framework for the DFT requires modulo operations of polynomials, which means working with so-called polynomial algebras . Associated with polynomial algebras is the Chinese remainder theorem, which describes the DFT algebraically and can be used as a tool to concisely derive various FFTs as well as convolution algorithms , , , (see also Chapter Winograd’s Short DFT Algorithms). The polynomial algebra framework was fully developed for signal processing as part of the algebraic signal processing theory (ASP). ASP identifies the structure underlying many transforms used in signal processing, provides deep insight into their properties, and enables the derivation of their fast algorithms , , , . Here we focus on the algebraic description of the DFT and on the algebraic derivation of the general-radix Cooley-Tukey FFT from Factoring the Signal Processing Operators. The derivation will make use of and extend the Polynomial Description of Signals. We start with motivating the appearance of modulo operations. The z-transform associates with infinite discrete signals X=(⋯,x(-1),x(0),x(1),⋯) a Laurent series: X ↦ X ( s ) = ∑ n ∈ Z x ( n ) s n . Here we used s=z-1 to simplify the notation in the following. The DTFT of X is the evaluation of X(s) on the unit circle X ( e - j ω ) , - π < ω ≤ π . Finally, filtering or (linear) convolution is simply the multiplication of Laurent series, H * X ↔ H ( s ) X ( s ) . For finite signals X=(x(0),⋯,x(N-1)) one expects that the equivalent of becomes a mapping to polynomials of degree N-1, X ↦ X ( s ) = ∑ n = 0 N - 1 x ( n ) s n , and that the DFT is an evaluation of these polynomials. Indeed, the definition of the DFT in Winograd’s Short DFT Algorithms shows that C ( k ) = X ( W N k ) = X ( e - j 2 π k N ) , 0 ≤ k < N , i.e., the DFT computes the evaluations of the polynomial X(s) at the nth roots of unity. The problem arises with the equivalent of , since the multiplication H(s)X(s) of two polynomials of degree N-1 yields one of degree 2N-2. Also, it does not coincide with the circular convolution known to be associated with the DFT. The solution to both problems is to reduce the product modulo sn-1: H * circ X ↔ H ( s ) X ( s ) mod ( s n - 1 ) .

Concept Infinite Time Finite Time Signal X ( s ) = ∑ n ∈ Z x ( n ) s n ∑ n = 0 N - 1 x ( n ) s n Filter H ( s ) = ∑ n ∈ Z h ( n ) s n ∑ n = 0 N - 1 h ( n ) s n Convolution H ( s ) X ( s ) H ( s ) X ( s ) mod ( s n - 1 ) Fourier transform DTFT: X ( e - j ω ) , - π < ω ≤ π DFT: X ( e - j 2 π k n ) , 0 ≤ k < n

Infinite and finite discrete time signal processing.

The resulting polynomial then has again degree N-1 and this form of convolution becomes equivalent to circular convolution of the polynomial coefficients. We also observe that the evaluation points in are precisely the roots of sn-1. This connection will become clear in this chapter. The discussion is summarized in . The proper framework to describe the multiplication of polynomials modulo a fixed polynomial are polynomial algebras. Together with the Chinese remainder theorem, they provide the theoretical underpinning for the DFT and the Cooley-Tukey FFT. In this chapter, the DFT will naturally arise as a linear mapping with respect to chosen bases, i.e., as a matrix. Indeed, the definition shows that if all input and outputs are collected into vectors X=(X(0),⋯,X(N-1)) and C=(C(0),⋯C(N-1)), then is equivalent to C = DFT N X , where DFT N = [ W N k n ] 0 ≤ k , n < N . The matrix point of view is adopted in the FFT books , .

Polynomial Algebras and the DFT In this section we introduce polynomial algebras and explain how they are associated to transforms. Then we identify this connection for the DFT. Later we use polynomial algebras to derive the Cooley-Tukey FFT. For further background on the mathematics in this section and polynomial algebras in particular, we refer to .

Polynomial Algebra An algebra A is a vector space that also provides a multiplication of its elements such that the distributivity law holds (see for a complete definition). Examples include the sets of complex or real numbers C or R, and the sets of complex or real polynomials in the variable s: C[s] or R[s]. The key player in this chapter is the polynomial algebra. Given a fixed polynomial P(s) of degree deg(P)=N, we define a polynomial algebra as the set C [ s ] / P ( s ) = { X ( s ) ∣ deg ( X ) < deg ( P ) } of polynomials of degree smaller than N with addition and multiplication modulo P. Viewed as a vector space, C[s]/P(s) hence has dimension N. Every polynomial X(s)∈C[s] is reduced to a unique polynomial R(s) modulo P(s) of degree smaller than N. R(s) is computed using division with rest, namely X ( s ) = Q ( s ) P ( s ) + R ( s ) , deg ( R ) < deg ( P ) . Regarding this equation modulo P, P(s) becomes zero, and we get X ( s ) ≡ R ( s ) mod P ( s ) . We read this equation as “X(s) is congruent (or equal) R(s) modulo P(s).” We will also write X(s)modP(s) to denote that X(s) is reduced modulo P(s). Obviously, P ( s ) ≡ 0 mod P ( s ) . As a simple example we consider A=C[s]/(s2-1), which has dimension 2. A possible basis is b=(1,s). In A, for example, s·(s+1)=s2+s≡s+1mod(s2-1), obtained through division with rest s 2 + s = 1 · ( s 2 - 1 ) + ( s + 1 ) or simply by replacing s2 with 1 (since s2-1=0 implies s2=1).

Chinese Remainder Theorem (CRT) Assume P(s)=Q(s)R(s) factors into two coprime (no common factors) polynomials Q and R. Then the Chinese remainder theorem (CRT) for polynomials is the linear mappingMore precisely, isomorphism of algebras or isomorphism of A-modules. Δ : C [ s ] / P ( s ) → C [ s ] / Q ( s ) ⊕ C [ s ] / R ( s ) , X ( s ) ↦ ( X ( s ) mod Q ( s ) , X ( s ) mod R ( s ) ) . Here, ⊕ is the Cartesian product of vector spaces with elementwise operation (also called outer direct sum). In words, the CRT asserts that computing (addition, multiplication, scalar multiplication) in C[s]/P(s) is equivalent to computing in parallel in C[s]/Q(s) and C[s]/R(s). If we choose bases b,c,d in the three polynomial algebras, then Δ can be expressed as a matrix. As usual with linear mappings, this matrix is obtained by mapping every element of b with Δ, expressing it in the concatenation c∪d of the bases c and d, and writing the results into the columns of the matrix. As an example, we consider again the polynomial P(s)=s2-1=(s-1)(s+1) and the CRT decomposition Δ : C [ s ] / ( s 2 - 1 ) → C [ s ] / ( x - 1 ) ⊕ C [ s ] / ( x + 1 ) . As bases, we choose b=(1,x),c=(1),d=(1). Δ(1)=(1,1) with the same coordinate vector in c∪d=(1,1). Further, because of x≡1mod(x-1) and x≡-1mod(x+1), Δ(x)=(x,x)≡(1,-1) with the same coordinate vector. Thus, Δ in matrix form is the so-called butterfly matrix, which is a DFT of size 2: DFT2=[ 1 1 1 -1 ].

Polynomial Transforms Assume P(s)∈C[s] has pairwise distinct zeros α=(α0,⋯,αN-1). Then the CRT can be used to completely decompose C[s]/P(s) into its spectrum: Δ : C [ s ] / P ( s ) → C [ s ] / ( s - α 0 ) ⊕ ⋯ ⊕ C [ s ] / ( s - α N - 1 ) , X ( s ) ↦ ( X ( s ) mod ( s - α 0 ) , ⋯ , X ( s ) mod ( s - α N - 1 ) ) = ( s ( α 0 ) , ⋯ , s ( α N - 1 ) ) . If we choose a basis b=(P0(s),⋯,PN-1(s)) in C[s]/P(s) and bases bi=(1) in each C[s]/(s-αi), then Δ, as a linear mapping, is represented by a matrix. The matrix is obtained by mapping every basis element Pn, 0≤n<N, and collecting the results in the columns of the matrix. The result is P b , α = [ P n ( α k ) ] 0 ≤ k , n < N and is called the polynomial transform for A=C[s]/P(s) with basis b. If, in general, we choose bi=(βi) as spectral basis, then the matrix corresponding to the decomposition is the scaled polynomial transform diag 0 ≤ k < N ( 1 / β n ) P b , α , where diag0≤n<N(γn) denotes a diagonal matrix with diagonal entries γn. We jointly refer to polynomial transforms, scaled or not, as Fourier transforms.

DFT as a Polynomial Transform We show that the DFTN is a polynomial transform for A=C[s]/(sN-1) with basis b=(1,s,⋯,sN-1). Namely, s N - 1 = ∏ 0 ≤ k < N ( x - W N k ) , which means that Δ takes the form Δ : C [ s ] / ( s N - 1 ) → C [ s ] / ( s - W N 0 ) ⊕ ⋯ ⊕ C [ s ] / ( s - W N N - 1 ) , X ( s ) ↦ ( X ( s ) mod ( s - W N 0 ) , ⋯ , X ( s ) mod ( s - W N N - 1 ) ) = ( X ( W N 0 ) , ⋯ , X ( W N N - 1 ) ) . The associated polynomial transform hence becomes P b , α = [ W N k n ] 0 ≤ k , n < N = DFT N . This interpretation of the DFT has been known at least since , and clarifies the connection between the evaluation points in and the circular convolution in . In , DFTs of types 1–4 are defined, with type 1 being the standard DFT. In the algebraic framework, type 3 is obtained by choosing A=C[s]/(sN+1) as algebra with the same basis as before: P b , α = [ W N ( k + 1 / 2 ) n ] 0 ≤ k , n < N = DFT -3 N , The DFTs of type 2 and 4 are scaled polynomial transforms .

Algebraic Derivation of the Cooley-Tukey FFT Knowing the polynomial algebra underlying the DFT enables us to derive the Cooley-Tukey FFT algebraically. This means that instead of manipulating the DFT definition, we manipulate the polynomial algebra C[s]/(sN-1). The basic idea is intuitive. We showed that the DFT is the matrix representation of the complete decomposition . The Cooley-Tukey FFT is now derived be performing this decomposition in steps as shown in . Each step yields a sparse matrix; hence, the DFTN is factorized into a product of sparse matrices, which will be the matrix representation of the Cooley-Tukey FFT.

Basic idea behind the algebraic derivation of Cooley-Tukey type algorithms

This stepwise decomposition can be formulated generically for polynomial transforms , . Here, we consider only the DFT. We first introduce the matrix notation we will use and in particular the Kronecker product formalism that became mainstream for FFTs in in , . Then we first derive the radix-2 FFT using a factorization of sN-1. Subsequently, we obtain the general-radix FFT using a decomposition of sN-1.

Matrix Notation We denote the N×N identity matrix with IN, and diagonal matrices with diag 0 ≤ k < N ( γ k ) = γ 0 ⋱ γ N - 1 . The N×N stride permutation matrix is defined for N=KM by the permutation L M N : i K + j ↦ j M + i for 0≤i<K,0≤j<M. This definition shows that LMN transposes a K×M matrix stored in row-major order. Alternatively, we can write L M N : i ↦ iM mod N - 1 , for 0 ≤ i < N - 1 , N - 1 ↦ N - 1 . For example (· means 0), L 2 6 = 1 · · · · · · · 1 · · · · · · · 1 · · 1 · · · · · · · 1 · · · · · · · 1 . LN/2N is sometimes called the perfect shuffle. Further, we use matrix operators; namely the direct sum A ⊕ B = A B and the Kronecker or tensor product A ⊗ B = [ a k , ℓ B ] k , ℓ , for A = [ a k , ℓ ] . In particular, I n ⊗ A = A ⊕ ⋯ ⊕ A = A ⋱ A is block-diagonal. We may also construct a larger matrix as a matrix of matrices, e.g., A B B A . If an algorithm for a transform is given as a product of sparse matrices built from the constructs above, then an algorithm for the transpose or inverse of the transform can be readily derived using mathematical properties including ( A B ) T = B T A T , ( A B ) - 1 = B - 1 A - 1 , ( A ⊕ B ) T = A T ⊕ B T , ( A ⊕ B ) - 1 = A - 1 ⊕ B - 1 , ( A ⊗ B ) T = A T ⊗ B T , ( A ⊗ B ) - 1 = A - 1 ⊗ B - 1 . Permutation matrices are orthogonal, i.e., PT=P-1. The transposition or inversion of diagonal matrices is obvious.

Radix-2 FFT The DFT decomposes A=C[s]/(sN-1) with basis b=(1,s,⋯,sN-1) as shown in . We assume N=2M. Then s 2 M - 1 = ( s M - 1 ) ( s M + 1 ) factors and we can apply the CRT in the following steps: C [ s ] / ( s N - 1 ) → C [ s ] / ( s M - 1 ) ⊕ C [ s ] / ( s M + 1 ) → ⊕ 0 ≤ i < M C [ s ] / ( x - W N 2 i ) ⊕ ⊕ 0 ≤ i < M C [ s ] / ( x - W M 2 i + 1 ) → ⊕ 0 ≤ i < N C [ s ] / ( x - W N i ) . As bases in the smaller algebras C[s]/(sM-1) and C[s]/(sM+1), we choose c=d=(1,s,⋯,sM-1). The derivation of an algorithm for DFTN based on - is now completely mechanical by reading off the matrix for each of the three decomposition steps. The product of these matrices is equal to the DFTN. First, we derive the base change matrix B corresponding to . To do so, we have to express the base elements sn∈b in the basis c∪d; the coordinate vectors are the columns of B. For 0≤n<M, sn is actually contained in c and d, so the first M columns of B are B = I M * I M * , where the entries * are determined next. For the base elements sM+n, 0≤n<M, we have s M + n ≡ s n mod ( s M - 1 ) , s M + n ≡ - s n mod ( s M + 1 ) , which yields the final result B = I M - I M I M - I M = DFT 2 ⊗ I M . Next, we consider step . C[s]/(sM-1) is decomposed by DFTM and C[s]/(sM+1) by DFT-3M in . Finally, the permutation in step is the perfect shuffle LMN, which interleaves the even and odd spectral components (even and odd exponents of WN). The final algorithm obtained is DFT 2 M = L M N ( DFT M ⊕ DFT -3 M ) ( DFT 2 ⊗ I M ) . To obtain a better known form, we use DFT-3M=DFTMDM, with DM=diag0≤i<M(WNi), which is evident from . It yields DFT 2 M = L M N ( DFT M ⊕ DFT M D M ) ( DFT 2 ⊗ I M ) = L M N ( I 2 ⊗ DFT M ) ( I M ⊕ D M ) ( DFT 2 ⊗ I M ) . The last expression is the radix-2 decimation-in-frequency Cooley-Tukey FFT. The corresponding decimation-in-time version is obtained by transposition using and the symmetry of the DFT: DFT 2 M = ( DFT 2 ⊗ I M ) ( I M ⊕ D M ) ( I 2 ⊗ DFT M ) L 2 N . The entries of the diagonal matrix IM⊕DM are commonly called twiddle factors. The above method for deriving DFT algorithms is used extensively in .

General-radix FFT To algebraically derive the general-radix FFT, we use the decomposition property of sN-1. Namely, if N=KM then s N - 1 = ( s M ) K - 1 . Decomposition means that the polynomial is written as the composition of two polynomials: here, sM is inserted into sK-1. Note that this is a special property: most polynomials do not decompose. Based on this polynomial decomposition, we obtain the following stepwise decomposition of C[s]/(sN-1), which is more general than the previous one in –. The basic idea is to first decompose with respect to the outer polynomial tK-1, t=sM, and then completely : C [ s ] / ( s N - 1 ) = C [ x ] / ( ( s M ) K - 1 ) → ⊕ 0 ≤ i < K C [ s ] / ( s M - W K i ) → ⊕ 0 ≤ i < K ⊕ 0 ≤ j < M C [ s ] / ( x - W N j K + i ) → ⊕ 0 ≤ i < N C [ s ] / ( x - W N i ) . As bases in the smaller algebras C[s]/(sM-WKi) we choose ci=(1,s,⋯,sM-1). As before, the derivation is completely mechanical from here: only the three matrices corresponding to – have to be read off. The first decomposition step requires us to compute snmod(sM-WKi), 0≤n<N. To do so, we decompose the index n as n=ℓM+m and compute s n = s ℓ M + m = ( s M ) ℓ s m ≡ W k ℓ m s m mod ( s M - W K i ) . This shows that the matrix for is given by DFTK⊗IM. In step , each C[s]/(sM-WKi) is completely decomposed by its polynomial transform DFT M ( i , K ) = DFT