by Markus Pueschel, Carnegie Mellon University
In 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) 10. 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 6. 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 9, 21, 22, 1 (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 13, 14, 11, 12. 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:
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
Finally, filtering or (linear) convolution is simply the multiplication of Laurent series,
For finite signals X=(x(0),⋯,x(N–1)) one expects that the equivalent of Equation 8.1 becomes a mapping to polynomials of degree N–1,
and that the DFT is an evaluation of these polynomials. Indeed, the definition of the DFT in Winograd’s Short DFT Algorithms shows that
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 Equation 8.3, 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:
Concept | Infinite Time | Finite Time |
Signal | X ( s ) = ∑n ∈ Z x ( n ) sn | |
Filter | H ( s ) = ∑n ∈ Z h ( n ) sn | |
Convolution | H ( s ) X ( s ) | |
Fourier transform |
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 Equation 8.5 are precisely the roots of sn–1. This connection will become clear in this chapter.
The discussion is summarized in Table 8.1.
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 to
where
The matrix point of view is adopted in the FFT books 18, 17.
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 6.
An algebra A is a vector space that also provides a multiplication of its elements such that the distributivity law holds (see 6 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
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
Regarding this equation modulo P, P(s) becomes zero, and we get
We read this equation as “X(s) is congruent (or equal) R(s) modulo P(s).” We will also write to denote that X(s) is reduced modulo P(s). Obviously,
As a simple example we consider , which has dimension 2. A possible basis is b=(1,s). In A, for example, , obtained through division with rest
or simply by replacing s2 with 1 (since s2–1=0 implies s2=1).
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 mapping[1]
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
As bases, we choose . Δ(1)=(1,1) with the same coordinate vector in c∪d=(1,1). Further, because of and , Δ(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: .
Assume P(s)∈C[s] has pairwise distinct zeros . Then the CRT can be used to completely decompose C[s]/P(s) into its spectrum:
If we choose a basis in C[s]/P(s) and bases bi=(1) in each , 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
and is called the polynomial transform for A=C[s]/P(s) with basis b.
If, in general, we choose as spectral basis, then the matrix corresponding to the decomposition Equation 8.16 is the scaled polynomial transform
where denotes a diagonal matrix with diagonal entries γn.
We jointly refer to polynomial transforms, scaled or not, as Fourier transforms.