DFT and FFT: An Algebraic View

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

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 by 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.

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 , . 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 M · diag 0 ≤ i < M ( W N i j ) . At this point, C[s]/(sN-1) is completely decomposed, but the spectrum is ordered according to jK+i, 0≤i<M, 0≤j<K (j runs faster). The desired order is iM+j. Thus, in step , we need to apply the permutation jK+i↦iM+j, which is exactly the stride permutation LMN in . In summary, we obtain the Cooley-Tukey decimation-in-frequency FFT with arbitrary radix: L M N ⊕ 0 ≤ i < K DFT M · diag j = 0 M - 1 ( W N i j ) ( DFT k ⊗ I M ) = L M N ( I K ⊗ DFT M ) T M N ( DFT k ⊗ I M ) . The matrix TMN is diagonal and usually called the twiddle matrix. Transposition using yields the corresponding decimation-in-time version: ( DFT k ⊗ I M ) T M N ( I K ⊗ DFT M ) L K N .

Discussion and Further Reading This chapter only scratches the surface of the connection between algebra and the DFT or signal processing in general. We provide a few references for further reading.

Algebraic Derivation of Transform Algorithms As mentioned before, the use of polynomial algebras and the CRT underlies much of the early work on FFTs and convolution algorithms , , . For example, Winograd's work on FFTs minimizes the number of non-rational multiplications. This and his work on complexity theory in general makes heavy use of polynomial algebras , , (see Chapter Winograd’s Short DFT Algorithms for more information and references). See for a broad treatment of algebraic complexity theory. Since C[x]/(sN-1)=C[CN] can be viewed a group algebra for the cyclic group, the methods shown in this chapter can be translated into the context of group representation theory. For example, derives the general-radix FFT using group theory and also uses already the Kronecker product formalism. So does Beth and started the area of FFTs for more general groups , . However, Fourier transforms for groups have found only sporadic applications . Along a related line of work, shows that using group theory it is possible that to discover and generate certain algorithms for trigonometric transforms, such as discrete cosine transforms (DCTs), automatically using a computer program. More recently, the polynomial algebra framework was extended to include most trigonometric transforms used in signal processing , , namely, besides the DFT, the discrete cosine and sine transforms and various real DFTs including the discrete Hartley transform. It turns out that the same techniques shown in this chapter can then be applied to derive, explain, and classify most of the known algorithms for these transforms and even obtain a large class of new algorithms including general-radix algorithms for the discrete cosine and sine transforms (DCTs/DSTs) , , , . This latter line of work is part of the algebraic signal processing theory briefly discussed next.

Algebraic Signal Processing Theory The algebraic properties of transforms used in the above work on algorithm derivation hints at a connection between algebra and (linear) signal processing itself. This is indeed the case and was fully developed in a recent body of work called algebraic signal processing theory (ASP). The foundation of ASP is developed in , , . ASP first identifies the algebraic structure of (linear) signal processing: the common assumptions on available operations for filters and signals make the set of filters an algebraA and the set of signals an associated A-module M. ASP then builds a signal processing theory formally from the axiomatic definition of a signal model: a triple (A,M,Φ), where Φ generalizes the idea of the z-transform to mappings from vector spaces of signal values to M. If a signal model is given, other concepts, such as spectrum, Fourier transform, frequency response are automatically defined but take different forms for different models. For example, infinite and finite time as discussed in are two examples of signal models. Their complete definition is provided in and identifies the proper notion of a finite z-transform as a mapping Cn→C[s]/(sn-1). Signal model Infinite time Finite time A ∑ n ∈ Z H ( n ) s n ∣ ( ⋯ , H ( - 1 ) , H ( 0 ) , H ( 1 ) , ⋯ ) ∈ ℓ 1 ( Z ) C [ x ] / ( s n - 1 ) M ∑ n ∈ Z X ( n ) s n ∣ ( ⋯ , X ( - 1 ) , X ( 0 ) , X ( 1 ) , ⋯ ) ∈ ℓ 2 ( Z ) C [ s ] / ( s n - 1 ) Φ Φ : ℓ 2 ( Z ) → M Φ : C n → M defined in defined in

ASP shows that many signal models are in principle possible, each with its own notion of filtering and Fourier transform. Those that support shift-invariance have commutative algebras. Since finite-dimensional commutative algebras are precisely polynomial algebras, their appearance in signal processing is explained. For example, ASP identifies the polynomial algebras underlying the DCTs and DSTs, which hence become Fourier transforms in the ASP sense. The signal models are called finite space models since they support signal processing based on an undirected shift operator, different from the directed time shift. Many more insights are provided by ASP including the need for and choices in choosing boundary conditions, properties of transforms, techniques for deriving new signal models, and the concise derivation of algorithms mentioned before.