DFT and FFT: An Algebraic Viewm16331DFT and FFT: An Algebraic View1.142008/05/27 09:58:50 GMT-52009/09/18 16:08:41.940 GMT-5MarkusPueschelMarkus Pueschelpueschel@ece.cmu.eduMarkusPueschelMarkus Pueschelpueschel@ece.cmu.eduC.SidneyBurrusC. Sidney Burruscsb@rice.eduDanielCollinsWilliamsonDaniel Williamsondcwill@cnx.orgMarkusPueschelMarkus Pueschelpueschel@ece.cmu.eduMathematics and Statisticsenby 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 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∈Zx(n)sn.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=0N-1x(n)sn,and that the DFT is an evaluation of these polynomials. Indeed, the
definition of the DFT in Winograd’s Short DFT Algorithms shows thatC(k)=X(WNk)=X(e-j2πkN),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*circX↔H(s)X(s)mod(sn-1).ConceptInfinite TimeFinite TimeSignalX(s)=∑n∈Zx(n)sn∑n=0N-1x(n)snFilterH(s)=∑n∈Zh(n)sn∑n=0N-1h(n)snConvolutionH(s)X(s)H(s)X(s)mod(sn-1)Fourier transformDTFT:X(e-jω),-π<ω≤πDFT:X(e-j2πkn),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
Winograd’s Short DFT Algorithms is equivalent toC=DFTNX,whereDFTN=[WNkn]0≤k,n<N.The matrix point of view is adopted in the FFT books
, .Polynomial Algebras and the DFTIn 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 AlgebraAn 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 setC[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, namelyX(s)=Q(s)P(s)+R(s),deg(R)<deg(P).Regarding this equation modulo P, P(s) becomes zero, and we getX(s)≡R(s)modP(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)≡0modP(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 rests2+s=1·(s2-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)modQ(s),X(s)modR(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]/(s2-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=[111-1].Polynomial TransformsAssume 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 isPb,α=[Pn(αk)]0≤k,n<Nand 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 transformdiag0≤k<N(1/βn)Pb,α,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 TransformWe show that the DFTN is a polynomial transform for A=C[s]/(sN-1) with basis b=(1,s,⋯,sN-1). Namely,sN-1=∏0≤k<N(x-WNk),which means that Δ takes the formΔ:C[s]/(sN-1)→C[s]/(s-WN0)⊕⋯⊕C[s]/(s-WNN-1),X(s)↦(X(s)mod(s-WN0),⋯,X(s)mod(s-WNN-1))=(X(WN0),⋯,X(WNN-1)).The associated polynomial transform hence becomesPb,α=[WNkn]0≤k,n<N=DFTN.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:Pb,α=[WN(k+1/2)n]0≤k,n<N=DFT-3N,The DFTs of type 2 and 4 are scaled polynomial
transforms .Algebraic Derivation of the Cooley-Tukey FFTKnowing 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 NotationWe denote the N×N identity matrix with IN, and diagonal
matrices withdiag0≤k<N(γk)=γ0⋱γN-1.The N×Nstride
permutation matrix is defined for N=KM by the permutationLMN:iK+j↦jM+ifor 0≤i<K,0≤j<M. This definition shows that LMN
transposes a K×M matrix stored in row-major order.
Alternatively, we can writeLMN:i↦iMmodN-1,for0≤i<N-1,N-1↦N-1.For example (· means 0),L26=1·······1·······1··1·······1·······1.LN/2N is sometimes called the perfect shuffle.Further, we use matrix operators; namely the direct sumA⊕B=ABand the Kronecker or tensor productA⊗B=[ak,ℓB]k,ℓ,forA=[ak,ℓ].In particular,In⊗A=A⊕⋯⊕A=A⋱Ais block-diagonal.We may also construct a larger matrix as a matrix of
matrices, e.g.,ABBA.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(AB)T=BTAT,(AB)-1=B-1A-1,(A⊕B)T=AT⊕BT,(A⊕B)-1=A-1⊕B-1,(A⊗B)T=AT⊗BT,(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 FFTThe DFT decomposes A=C[s]/(sN-1) with basis b=(1,s,⋯,sN-1) as shown in . We assume N=2M.
Thens2M-1=(sM-1)(sM+1)factors and we can apply the CRT in
the following steps:C[s]/(sN-1)→C[s]/(sM-1)⊕C[s]/(sM+1)→⊕0≤i<MC[s]/(x-WN2i)⊕⊕0≤i<MC[s]/(x-WM2i+1)→⊕0≤i<NC[s]/(x-WNi).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 areB=IM*IM*,where the entries * are determined next.
For the base elements sM+n, 0≤n<M, we havesM+n≡snmod(sM-1),sM+n≡-snmod(sM+1),which yields the final resultB=IM-IMIM-IM=DFT2⊗IM.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 isDFT2M=LMN(DFTM⊕DFT-3M)(DFT2⊗IM).To obtain a better known form, we use DFT-3M=DFTMDM,
with DM=diag0≤i<M(WNi), which is evident from
. It yieldsDFT2M=LMN(DFTM⊕DFTMDM)(DFT2⊗IM)=LMN(I2⊗DFTM)(IM⊕DM)(DFT2⊗IM).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:DFT2M=(DFT2⊗IM)(IM⊕DM)(I2⊗DFTM)L2N.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 FFTTo algebraically derive the general-radix FFT, we use the decomposition property of sN-1. Namely, if N=KM thensN-1=(sM)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]/(sN-1)=C[x]/((sM)K-1)→⊕0≤i<KC[s]/(sM-WKi)→⊕0≤i<K⊕0≤j<MC[s]/(x-WNjK+i)→⊕0≤i<NC[s]/(x-WNi).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 computesn=sℓM+m=(sM)ℓsm≡Wkℓmsmmod(sM-WKi).This shows that the matrix for
is given by DFTK⊗IM.In step , each C[s]/(sM-WKi) is completely
decomposed by its polynomial transformDFTM(i,K)=DFTM·diag0≤i<M(WNij).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:LMN⊕0≤i<KDFTM·diagj=0M-1(WNij)(DFTk⊗IM)=LMN(IK⊗DFTM)TMN(DFTk⊗IM).The matrix TMN is diagonal and usually called the
twiddle matrix. Transposition using
yields the corresponding decimation-in-time version:(DFTk⊗IM)TMN(IK⊗DFTM)LKN.Discussion and Further ReadingThis 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 AlgorithmsAs 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 TheoryThe 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 modelInfinite timeFinite timeA∑n∈ZH(n)sn∣(⋯,H(-1),H(0),H(1),⋯)∈ℓ1(Z)C[x]/(sn-1)M∑n∈ZX(n)sn∣(⋯,X(-1),X(0),X(1),⋯)∈ℓ2(Z)C[s]/(sn-1)ΦΦ:ℓ2(Z)→MΦ:Cn→Mdefined in defined in

Infinite and finite time models as defined in ASP.

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.Auslander, L. and Feig, E. and Winograd, S.Abelian Semi-simple Algebras and Algorithms for the Discrete Fourier TransformAdvances in Applied Mathematics1984531–55Bürgisser, P. and Clausen, M. and Shokrollahi, M. A.Algebraic Complexity TheorySpringer1997Beth, Th.Verfahren der Schnellen Fouriertransformation [Fast Fourier Transform Methods]Teubner1984Britanak, V. and Rao, K. R.The fast generalized discrete Fourier transforms: A unified approach to the discrete sinusoidal transforms computationSignal Processing199979135–150Egner, S. and Püschel, M.Automatic Generation of Fast Discrete Signal TransformsIEEE Transactions on Signal Processing20014991992–2002Fuhrman, Paul A.A Polynomial Approach to Linear AlgebraSpringer Verlag1996New YorkMaslen, D. and Rockmore, D.Generalized FFTs – A survey of some recent resultsProceedings of IMACS Workshop in Groups and Computation199528182–238Nicholson, P. J.Algebraic Theory of Finite Fourier TransformsJournal of Computer and System Sciences19715524–547Nussbaumer, H. J.Fast Fourier Transformation and Convolution AlgorithmsSpringer19822ndOppenheim, A. V. and Schafer, R. W. and Buck, J. R.Discrete-Time Signal ProcessingPrentice Hall19992ndPüschel, M. and Moura, J. M. F.Algebraic Signal Processing Theoryavailable at http://arxiv.org/abs/cs.IT/0612077Püschel, M. and Moura, J. M. F.Algebraic Signal Processing Theory: Cooley-Tukey Type Algorithms for DCTs and DSTsIEEE Transactions on Signal Processing20085641502-1521a longer version is available at http://arxiv.org/abs/cs.IT/0702025Püschel, M. and Moura, J. M. F.Algebraic Signal Processing Theory: Foundation and 1-D TimeIEEE Transactions on Signal Processing20085683572-3585Püschel, M. and Moura, J. M. F.Algebraic Signal Processing Theory: 1-D SpaceIEEE Transactions on Signal Processing20085683586-3599Püschel, M. and Moura, J. M. F.The Algebraic Approach to the Discrete Cosine and Sine Transforms and their Fast AlgorithmsSIAM Journal of Computing20033251280–1316Rockmore, D.Some applications of generalized FFT'sProceedings of DIMACS Workshop in Groups and Computation199528329–370Tolimieri, R. and An, M. and Lu, C.Algorithms for Discrete Fourier Transforms and ConvolutionSpringer19972ndVan Loan, C.Computational Framework of the Fast Fourier TransformSiam1992Voronenko, Y. and Püschel, M.Algebraic derivation of general radix Cooley-Tukey algorithms for the real discrete Fourier transformProc. International Conference on Acoustics, Speech, and Signal Processing (ICASSP)20063876-879Voronenko, Y. and Püschel, M.Algebraic Signal Processing Theory: Cooley-Tukey Type Algorithms for Real DFTsIEEE Transactions on Signal Processingto appearWinograd, S.On Computing the Discrete Fourier TransformMathematics of Computation197832175–199Winograd, S.On the Multiplicative Complexity of the Discrete Fourier TransformAdvances in Mathematics19793283–117Winograd, S.Arithmetic Complexity of ComputationSiam1980