Polynomials are important in digital signal processing
because calculating the DFT can be viewed as a polynomial evaluation
problem and convolution can be viewed as polynomial multiplication
[1], [5]. Indeed, this is the basis for the important
results of Winograd discussed in Winograd’s Short DFT Algorithms. A length-N signal
x(n)x(n) will be represented by an N-1N-1 degree polynomial X(s)X(s)
defined by
X
(
s
)
=
∑
n
=
0
N
-
1
x
(
n
)
s
n
X
(
s
)
=
∑
n
=
0
N
-
1
x
(
n
)
s
n
(1)
This polynomial X(s)X(s) is a single entity with the coefficients
being the values of x(n)x(n). It is somewhat similar to the use of
matrix or vector notation to efficiently represent signals which
allows use of new mathematical tools.
The convolution of two finite length sequences, x(n)x(n) and
h(n)h(n), gives an output sequence defined by
y
(
n
)
=
∑
k
=
0
N
-
1
x
(
k
)
h
(
n
-
k
)
y
(
n
)
=
∑
k
=
0
N
-
1
x
(
k
)
h
(
n
-
k
)
(2)
n=0,1,2,⋯,2N-1n=0,1,2,⋯,2N-1 where h(k)=0h(k)=0 for k<0k<0. This is
exactly the same operation as calculating the coefficients when
multiplying two polynomials. Equation Equation 2 is the same as
Y
(
s
)
=
X
(
s
)
H
(
s
)
Y
(
s
)
=
X
(
s
)
H
(
s
)
(3)
In fact, convolution of number sequences, multiplication of
polynomials, and the multiplication of integers (except for the
carry operation) are all the same operations. To obtain cyclic
convolution, where the indices in Equation 2 are all evaluated
modulo NN, the polynomial multiplication in Equation 3 is done
modulo the polynomial P(s)=sN-1P(s)=sN-1. This is seen by noting that
N=0N=0 mod NN, therefore, sN=1sN=1 and the polynomial modulus is
sN-1sN-1.
Residue reduction of one polynomial modulo another is
defined similarly to residue reduction for integers. A polynomial
F(s)F(s) has a residue polynomial R(s)R(s) modulo P(sP(s) if, for a given
F(s)F(s) and P(s)P(s), a Q(S)Q(S) and R(s)R(s) exist such that
F
(
s
)
=
Q
(
s
)
P
(
s
)
+
R
(
s
)
F
(
s
)
=
Q
(
s
)
P
(
s
)
+
R
(
s
)
(4)
with degree{R(s)}<degree{P(s)}degree{R(s)}<degree{P(s)}. The notation that will
be used is
R
(
s
)
=
(
(
F
(
s
)
)
)
P
(
s
)
R
(
s
)
=
(
(
F
(
s
)
)
)
P
(
s
)
(5)
For example,
(
s
+
1
)
=
(
(
s
4
+
s
3
-
s
-
1
)
)
(
s
2
-
1
)
(
s
+
1
)
=
(
(
s
4
+
s
3
-
s
-
1
)
)
(
s
2
-
1
)
(6)
The concepts of factoring a polynomial and of primeness are an
extension of these ideas for integers. For a given
allowed set of coefficients (values of x(n)x(n)), any polynomial has a unique factored representation
F
(
s
)
=
∏
i
=
1
M
F
i
(
s
)
k
i
F
(
s
)
=
∏
i
=
1
M
F
i
(
s
)
k
i
(7)
where the Fi(s)Fi(s) are relatively prime. This is analogous to the
fundamental theorem of arithmetic.
There is a very useful operation that is an extension of
the integer Chinese Remainder Theorem (CRT) which says that if the
modulus polynomial can be factored into relatively prime factors
P
(
s
)
=
P
1
(
s
)
P
2
(
s
)
P
(
s
)
=
P
1
(
s
)
P
2
(
s
)
(8)
then there exist two polynomials, K1(s)K1(s) and K2(s)K2(s), such that any
polynomial F(s)F(s) can be recovered from its residues by
F
(
s
)
=
K
1
(
s
)
F
1
(
s
)
+
K
2
(
s
)
F
2
(
s
)
mod
P
(
s
)
F
(
s
)
=
K
1
(
s
)
F
1
(
s
)
+
K
2
(
s
)
F
2
(
s
)
mod
P
(
s
)
(9)
where F1F1 and F2F2 are the residues given by
F
1
(
s
)
=
(
(
F
(
s
)
)
)
P
1
(
s
)
F
1
(
s
)
=
(
(
F
(
s
)
)
)
P
1
(
s
)
(10)
and
F
2
(
s
)
=
(
(
F
(
s
)
)
)
P
2
(
s
)
F
2
(
s
)
=
(
(
F
(
s
)
)
)
P
2
(
s
)
(11)
if the order of F(s)F(s) is less than P(s)P(s). This generalizes to any
number of relatively prime factors of P(s)P(s) and can be viewed as a
means of representing F(s)F(s) by several lower degree polynomials,
Fi(s)Fi(s).
This decomposition of F(s)F(s) into lower degree polynomials is
the process used to break a DFT or convolution into several simple
problems which are solved and then recombined using the CRT of
Equation 9. This is another form of the “divide and conquer" or
“organize and share"
approach similar to the index mappings in Multidimensional Index Mapping.
One useful property of the CRT is for convolution. If cyclic
convolution of x(n)x(n) and h(n)h(n) is expressed in terms of
polynomials by
Y
(
s
)
=
H
(
s
)
X
(
s
)
mod
P
(
s
)
Y
(
s
)
=
H
(
s
)
X
(
s
)
mod
P
(
s
)
(12)
where P(s)=sN-1P(s)=sN-1, and if P(s)P(s) is factored into two
relatively prime factors P=P1P2P=P1P2, using residue reduction of
H(s)H(s) and X(s)X(s) modulo P1P1 and P2P2, the lower degree residue
polynomials can be multiplied and the results recombined with the
CRT. This is done by
Y
(
s
)
=
(
(
K
1
H
1
X
1
+
K
2
H
2
X
2
)
)
P
Y
(
s
)
=
(
(
K
1
H
1
X
1
+
K
2
H
2
X
2
)
)
P
(13)
where
H
1
=
(
(
H
)
)
P
1
,
X
1
=
(
(
X
)
)
P
1
,
H
2
=
(
(
H
)
)
P
2
,
X
2
=
(
(
X
)
)
P
2
H
1
=
(
(
H
)
)
P
1
,
X
1
=
(
(
X
)
)
P
1
,
H
2
=
(
(
H
)
)
P
2
,
X
2
=
(
(
X
)
)
P
2
(14)
and K1K1 and K2K2 are the CRT coefficient polynomials from
Equation 9. This allows two shorter convolutions to replace one
longer one.
Another property of residue reduction that is useful in DFT
calculation is polynomial evaluation. To evaluate F(s)F(s) at s=xs=x,
F(s)F(s) is reduced modulo s-xs-x.
F
(
x
)
=
(
(
F
(
s
)
)
)
s
-
x
F
(
x
)
=
(
(
F
(
s
)
)
)
s
-
x
(15)
This is easily seen from the definition in Equation 4
F
(
s
)
=
Q
(
s
)
(
s
-
x
)
+
R
(
s
)
F
(
s
)
=
Q
(
s
)
(
s
-
x
)
+
R
(
s
)
(16)
Evaluating s=xs=x gives R(s)=F(x)R(s)=F(x) which is a constant. For
the DFT this becomes
C
(
k
)
=
(
(
X
(
s
)
)
)
s
-
W
k
C
(
k
)
=
(
(
X
(
s
)
)
)
s
-
W
k
(17)
Details of the polynomial algebra useful in digital signal
processing can be found in [1], [4], [5].
The Z-transform of a number sequence x(n)x(n) is defined as
X
(
z
)
=
∑
n
=
0
∞
x
(
n
)
z
-
n
X
(
z
)
=
∑
n
=
0
∞
x
(
n
)
z
-
n
(18)
which is the same as the polynomial description in Equation 1 but
with a negative exponent. For a finite length-N sequence Equation 18
becomes
X
(
z
)
=
∑
n
=
0
N
-
1
x
(
n
)
z
-
n
X
(
z
)
=
∑
n
=
0
N
-
1
x
(
n
)
z
-
n
(19)
X
(
z
)
=
x
(
0
)
+
x
(
1
)
z
-
1
+
x
(
2
)
z
-
2
+
·
+
x
(
N
-
1
)
z
-
N
+
1
X
(
z
)
=
x
(
0
)
+
x
(
1
)
z
-
1
+
x
(
2
)
z
-
2
+
·
+
x
(
N
-
1
)
z
-
N
+
1
(20)
This N-1N-1 order polynomial takes on the values of the DFT of
x(n)x(n) when evaluated at
z
=
e
j
2
π
k
/
N
z
=
e
j
2
π
k
/
N
(21)
which gives
C
(
k
)
=
X
(
z
)
|
z
=
e
j
2
π
k
/
N
=
∑
n
=
0
N
-
1
x
(
n
)
e
-
j
2
π
n
k
/
N
C
(
k
)
=
X
(
z
)
|
z
=
e
j
2
π
k
/
N
=
∑
n
=
0
N
-
1
x
(
n
)
e
-
j
2
π
n
k
/
N
(22)
In terms of the positive exponent polynomial from Equation 1, the DFT is
C
(
k
)
=
X
(
s
)
|
s
=
W
k
C
(
k
)
=
X
(
s
)
|
s
=
W
k
(23)
where
W
=
e
-
j
2
π
/
N
W
=
e
-
j
2
π
/
N
(24)
is an NthNth root of unity (raising WW to the NthNth power gives one).
The NN values of the DFT are found from
X(s)X(s) evaluated at the N NthNth roots of unity which are equally
spaced around the unit circle in the complex ss plane.
One method of evaluating X(z)X(z) is the so-called Horner's
rule or nested evaluation. When expressed as a recursive
calculation, Horner's rule becomes the Goertzel algorithm which has
some computational advantages especially when only a few values of
the DFT are needed. The details and programs can be found in
[7], [2] and The DFT as Convolution or Filtering: Goertzel's Algorithm (or A Better DFT Algorithm)
Another method for evaluating X(s)X(s) is the residue reduction
modulo (s-Wk)(s-Wk) as shown in Equation 17. Each evaluation requires
NN multiplications and therefore, N2N2 multiplications for the NN
values of C(k)C(k).
C
(
k
)
=
(
(
X
(
s
)
)
)
(
s
-
W
k
)
C
(
k
)
=
(
(
X
(
s
)
)
)
(
s
-
W
k
)
(25)
A considerable reduction in required arithmetic can be achieved if
some operations can be shared between the reductions for different
values of kk. This is done by carrying out the residue reduction in
stages that can be shared rather than done in one step for each kk
in Equation 25.
The NN values of the DFT are values of X(s)X(s) evaluated at s
equal to the NN roots of the polynomial P(s)=sN-1P(s)=sN-1 which are
WkWk. First, assuming NN is even, factor P(s)P(s) as
P
(
s
)
=
(
s
N
-
1
)
=
P
1
(
s
)
P
2
(
s
)
=
(
s
N
/
2
-
1
)
(
s
N
/
2
+
1
)
P
(
s
)
=
(
s
N
-
1
)
=
P
1
(
s
)
P
2
(
s
)
=
(
s
N
/
2
-
1
)
(
s
N
/
2
+
1
)
(26)
X(s)X(s) is reduced modulo these two factors to give two residue
polynomials, X1(s)X1(s) and X2(s)X2(s). This process is repeated by
factoring P1P1 and further reducing X1X1 then factoring P2P2 and
reducing X2X2. This is continued until the factors are of first
degree which gives the desired DFT values as in Equation 25. This is
illustrated for a length-8 DFT. The polynomial whose roots are
WkWk, factors as
P
(
s
)
=
s
8
-
1
P
(
s
)
=
s
8
-
1
(27)
=
[
s
4
-
1
]
[
s
4
+
1
]
=
[
s
4
-
1
]
[
s
4
+
1
]
(28)
=
[
(
s
2
-
1
)
(
s
2
+
1
)
]
[
(
s
2
-
j
)
(
s
2
+
j
)
]
=
[
(
s
2
-
1
)
(
s
2
+
1
)
]
[
(
s
2
-
j
)
(
s
2
+
j
)
]
(29)
=
[
(
s
-
1
)
(
s
+
1
)
(
s
-
j
)
(
s
+
j
)
]
[
(
s
-
a
)
(
s
+
a
)
(
s
-
j
a
)
(
s
+
j
a
)
]
=
[
(
s
-
1
)
(
s
+
1
)
(
s
-
j
)
(
s
+
j
)
]
[
(
s
-
a
)
(
s
+
a
)
(
s
-
j
a
)
(
s
+
j
a
)
]
(30)
where a2=ja2=j. Reducing X(s)X(s) by the first factoring gives two
third degree polynomials
X
(
s
)
=
x
0
+
x
1
s
+
x
2
s
2
+
.
.
.
+
x
7
s
7
X
(
s
)
=
x
0
+
x
1
s
+
x
2
s
2
+
.
.
.
+
x
7
s
7
(31)
gives the residue polynomials
X
1
(
s
)
=
(
(
X
(
s
)
)
)
(
s
4
-
1
)
=
(
x
0
+
x
4
)
+
(
x
1
+
x
5
)
s
+
(
x
2
+
x
6
)
s
2
+
(
x
3
+
x
7
)
s
3
X
1
(
s
)
=
(
(
X
(
s
)
)
)
(
s
4
-
1
)
=
(
x
0
+
x
4
)
+
(
x
1
+
x
5
)
s
+
(
x
2
+
x
6
)
s
2
+
(
x
3
+
x
7
)
s
3
(32)
X
2
(
s
)
=
(
(
X
(
s
)
)
)
(
s
4
+
1
)
=
(
x
0
-
x
4
)
+
(
x
1
-
x
5
)
s
+
(
x
2
-
x
6
)
s
2
+
(
x
3
-
x
7
)
s
3
X
2
(
s
)
=
(
(
X
(
s
)
)
)
(
s
4
+
1
)
=
(
x
0
-
x
4
)
+
(
x
1
-
x
5
)
s
+
(
x
2
-
x
6
)
s
2
+
(
x
3
-
x
7
)
s
3
(33)
Two more levels of reduction are carried out to finally give the
DFT. Close examination shows the resulting algorithm to be the
decimation-in-frequency radix-2 Cooley-Tukey FFT [7], [2].
Martens [3] has used this approach to derive an efficient
DFT algorithm.
Other algorithms and types of FFT can be developed using polynomial
representations and some are presented in the generalization in DFT and FFT: An Algebraic View.
-
Blahut, Richard E. (1985). Fast Algorithms for Digital Signal Processing. Reading, Mass.: Addison-Wesley.
-
Burrus, C. S. and Parks, T. W. (1985). DFT/FFT and Convolution Algorithms. New York: John Wiley & Sons.
-
Martens, J. B. (1984, August). Recursive Cyclotomic Factorization – A New Algorithm for Calculating the Discrete Fourier Transform. IEEE Trans. on ASSP, 32(4), 750–762.
-
McClellan, J. H. and Rader, C. M. (1979). Number Theory in Digital Signal Processing. Englewood Cliffs, NJ: Prentice-Hall.
-
Nussbaumer, H. J. (1981, 1982). Fast Fourier Transform and Convolution Algorithms. (Second). Heidelberg, Germany: Springer-Verlag.
-
Niven, Ivan and Zuckerman, H. S. (1980). An Introduction to the Theory of Numbers. (Fourth). New York: John Wiley & Sons.
-
Oppenheim, A. V. and Schafer, R. W. (1999). Discrete-Time Signal Processing. (Second). [Earlier editions in 1975 and 1989]. Englewood Cliffs, NJ: Prentice-Hall.
"The Fast Fourier Transform (FFT) is a landmark algorithm used in fields ranging from signal processing to high-performance computing. First popularized by two American scientists in 1965, the […]"