by Markus Pueschel, Carnegie Mellon University
In infinite, or non-periodic, discrete-time signal processing, there
is a strong connection between the zz-transform, Laurent series,
convolution, and the discrete-time Fourier transform (DTFT)
Entry 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
Entry 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 Entry 9, Entry 21, Entry 22, Entry 1
(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
Entry 13, Entry 14, Entry 11, Entry 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 zz-transform associates with infinite discrete signals X=(⋯,x(-1),x(0),x(1),⋯)X=(⋯,x(-1),x(0),x(1),⋯) a Laurent series:
X
↦
X
(
s
)
=
∑
n
∈
Z
x
(
n
)
s
n
.
X
↦
X
(
s
)
=
∑
n
∈
Z
x
(
n
)
s
n
.
(1)
Here we used s=z-1s=z-1 to simplify the notation in the following.
The DTFT of XX is the evaluation of X(s)X(s) on the unit circle
X
(
e
-
j
ω
)
,
-
π
<
ω
≤
π
.
X
(
e
-
j
ω
)
,
-
π
<
ω
≤
π
.
(2)
Finally, filtering or (linear) convolution is simply the multiplication
of Laurent series,
H
*
X
↔
H
(
s
)
X
(
s
)
.
H
*
X
↔
H
(
s
)
X
(
s
)
.
(3)
For finite signals X=(x(0),⋯,x(N-1))X=(x(0),⋯,x(N-1)) one expects that the
equivalent of Equation 1 becomes a mapping to polynomials of
degree N-1N-1,
X
↦
X
(
s
)
=
∑
n
=
0
N
-
1
x
(
n
)
s
n
,
X
↦
X
(
s
)
=
∑
n
=
0
N
-
1
x
(
n
)
s
n
,
(4)
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
,
C
(
k
)
=
X
(
W
N
k
)
=
X
(
e
-
j
2
π
k
N
)
,
0
≤
k
<
N
,
(5)
i.e., the DFT computes the evaluations of the polynomial X(s)X(s)
at the nnth roots of unity.
The problem arises with the equivalent of Equation 3,
since the multiplication H(s)X(s)H(s)X(s) of two polynomials of degree N-1N-1
yields one of degree 2N-22N-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-1sn-1:
H
*
circ
X
↔
H
(
s
)
X
(
s
)
mod
(
s
n
-
1
)
.
H
*
circ
X
↔
H
(
s
)
X
(
s
)
mod
(
s
n
-
1
)
.
(6)
The resulting polynomial then has again degree N-1N-1 and this form of
convolution becomes equivalent to circular convolution of the
polynomial coefficients. We also observe that the evaluation points in
Equation 5 are precisely the roots of sn-1sn-1. This connection will
become clear in this chapter.
The discussion is summarized in Figure 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))X=(X(0),⋯,X(N-1)) and C=(C(0),⋯C(N-1))C=(C(0),⋯C(N-1)), then
(Reference) is equivalent to
C
=
DFT
N
X
,
C
=
DFT
N
X
,
(7)
where
DFT
N
=
[
W
N
k
n
]
0
≤
k
,
n
<
N
.
DFT
N
=
[
W
N
k
n
]
0
≤
k
,
n
<
N
.
(8)
The matrix point of view is adopted in the FFT books
Entry 18, Entry 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 Entry 6.
An algebra AA is a vector space that also provides a
multiplication of its elements such that the distributivity law holds
(see Entry 6 for a complete definition). Examples include
the sets of complex or real numbers CC or RR, and the sets of
complex or real polynomials in the variable ss: C[s]C[s] or R[s]R[s].
The key player in this chapter is the polynomial algebra. Given
a fixed polynomial P(s)P(s) of degree deg(P)=Ndeg(P)=N, we define a
polynomial algebra as the set
C
[
s
]
/
P
(
s
)
=
{
X
(
s
)
∣
deg
(
X
)
<
deg
(
P
)
}
C
[
s
]
/
P
(
s
)
=
{
X
(
s
)
∣
deg
(
X
)
<
deg
(
P
)
}
(9)
of polynomials of degree smaller than NN with addition and
multiplication modulo PP. Viewed as a vector space, C[s]/P(s)C[s]/P(s)
hence has dimension NN.
Every polynomial X(s)∈C[s]X(s)∈C[s] is reduced to a unique polynomial
R(s)R(s) modulo P(s)P(s) of degree smaller than NN. R(s)R(s) is computed
using division with rest, namely
X
(
s
)
=
Q
(
s
)
P
(
s
)
+
R
(
s
)
,
deg
(
R
)
<
deg
(
P
)
.
X
(
s
)
=
Q
(
s
)
P
(
s
)
+
R
(
s
)
,
deg
(
R
)
<
deg
(
P
)
.
(10)
Regarding this equation modulo PP, P(s)P(s) becomes zero, and we get
X
(
s
)
≡
R
(
s
)
mod
P
(
s
)
.
X
(
s
)
≡
R
(
s
)
mod
P
(
s
)
.
(11)
We read this equation as “X(s)X(s) is congruent (or equal) R(s)R(s)
modulo P(s)P(s).” We will also write X(s)modP(s)X(s)modP(s) to denote
that X(s)X(s) is reduced modulo P(s)P(s). Obviously,
P
(
s
)
≡
0
mod
P
(
s
)
.
P
(
s
)
≡
0
mod
P
(
s
)
.
(12)
As a simple example we consider A=C[s]/(s2-1)A=C[s]/(s2-1), which has
dimension 2. A possible basis is b=(1,s)b=(1,s). In AA, for example,
s·(s+1)=s2+s≡s+1mod(s2-1)s·(s+1)=s2+s≡s+1mod(s2-1), obtained
through division with rest
s
2
+
s
=
1
·
(
s
2
-
1
)
+
(
s
+
1
)
s
2
+
s
=
1
·
(
s
2
-
1
)
+
(
s
+
1
)
(13)
or simply by replacing s2s2 with 1 (since s2-1=0s2-1=0 implies s2=1s2=1).
Assume P(s)=Q(s)R(s)P(s)=Q(s)R(s) factors into two coprime (no common factors)
polynomials QQ and RR. Then the Chinese remainder theorem (CRT) for
polynomials is the linear mapping
Δ
:
C
[
s
]
/
P
(
s
)
→
C
[
s
]
/
Q
(
s
)
⊕
C
[
s
]
/
R
(
s
)
,
X
(
s
)
↦
(
X
(
s
)
mod
Q
(
s
)
,
X
(
s
)
mod
R
(
s
)
)
.
Δ
:
C
[
s
]
/
P
(
s
)
→
C
[
s
]
/
Q
(
s
)
⊕
C
[
s
]
/
R
(
s
)
,
X
(
s
)
↦
(
X
(
s
)
mod
Q
(
s
)
,
X
(
s
)
mod
R
(
s
)
)
.
(14)
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)C[s]/P(s) is equivalent to computing in parallel
in C[s]/Q(s)C[s]/Q(s) and C[s]/R(s)C[s]/R(s).
If we choose bases b,c,db,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 bb with ΔΔ,
expressing it in the concatenation c∪dc∪d of the bases cc and dd,
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)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
)
.
Δ
:
C
[
s
]
/
(
s
2
-
1
)
→
C
[
s
]
/
(
x
-
1
)
⊕
C
[
s
]
/
(
x
+
1
)
.
(15)
As bases, we choose b=(1,x),c=(1),d=(1)b=(1,x),c=(1),d=(1). Δ(1)=(1,1)Δ(1)=(1,1) with the same coordinate vector in c∪d=(1,1)c∪d=(1,1). Further,
because of x≡1mod(x-1)x≡1mod(x-1) and x≡-1mod(x+1)x≡-1mod(x+1), Δ(x)=(x,x)≡(1,-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
]DFT2=[
1 1
1 -1
].
Assume P(s)∈C[s]P(s)∈C[s] has pairwise distinct
zeros α=(α0,⋯,αN-1)α=(α0,⋯,αN-1). Then the CRT
can be used to completely
decompose C[s]/P(s)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
)
)
.
Δ
:
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
)
)
.
(16)
If we choose a basis b=(P0(s),⋯,PN-1(s))b=(P0(s),⋯,PN-1(s)) in C[s]/P(s)C[s]/P(s) and
bases bi=(1)bi=(1) in each C[s]/(s-αi)C[s]/(s-αi), then ΔΔ, as a linear
mapping, is represented by a matrix. The matrix is obtained by mapping
every basis element PnPn, 0≤n<N0≤n<N, and collecting the results in
the columns of the matrix. The result is
P
b
,
α
=
[
P
n
(
α
k
)
]
0
≤
k
,
n
<
N
P
b
,
α
=
[
P
n
(
α
k
)
]
0
≤
k
,
n
<
N
(17)
and is called the polynomial transform
for A=C[s]/P(s)A=C[s]/P(s) with basis bb.
If, in general, we choose bi=(βi)bi=(βi) as spectral basis, then the
matrix corresponding to the decomposition Equation 16 is the scaled
polynomial transform
diag
0
≤
k
<
N
(
1
/
β
n
)
P
b
,
α
,
diag
0
≤
k
<
N
(
1
/
β
n
)
P
b
,
α
,
(18)
where diag0≤n<N(γn)diag0≤n<N(γn) denotes a diagonal matrix
with diagonal entries γnγn.
We jointly refer to polynomial transforms, scaled or not, as Fourier
transforms.
We show that the DFTNDFTN is a polynomial transform for A=C[s]/(sN-1)A=C[s]/(sN-1) with basis b=(1,s,⋯,sN-1)b=(1,s,⋯,sN-1). Namely,
s
N
-
1
=
∏
0
≤
k
<
N
(
x
-
W
N
k
)
,
s
N
-
1
=
∏
0
≤
k
<
N
(
x
-
W
N
k
)
,
(19)
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
