A major application of the FFT is fast convolution or fast filtering
where the DFT of the signal is multiplied term-by-term by the DFT of
the impulse (helps to be doing finite impulse response (FIR) filtering)
and the time-domain output is obtained by taking the inverse DFT of
that product. What is less well-known is the DFT can be calculated
by convolution. There are several different approaches to this,
each with different application.
In this section a method quite different from the index
mapping or polynomial evaluation is developed. Rather than dealing
with the DFT directly, it is converted into a cyclic convolution
which must then be carried out by some efficient means. Those means
will be covered later, but here the conversion will be explained.
This method requires use of some number theory, which can be
found in an accessible form in Entry 1, Entry 12, Entry 13 and is easy
enough to verify on one's own. A good general reference on number
theory is Entry 14.
The DFT and cyclic convolution are defined by
C
(
k
)
=
∑
n
=
0
N
-
1
x
(
n
)
W
n
k
C
(
k
)
=
∑
n
=
0
N
-
1
x
(
n
)
W
n
k
(1)
y
(
k
)
=
∑
n
=
0
N
-
1
x
(
n
)
h
(
k
-
n
)
y
(
k
)
=
∑
n
=
0
N
-
1
x
(
n
)
h
(
k
-
n
)
(2)
For both, the indices are evaluated modulo NN. In order to convert
the DFT in Equation 1 into the cyclic convolution of
Equation 2, the nknk product must be changed to the k-nk-n
difference. With real numbers, this can be done with logarithms,
but it is more complicated when working in a finite set of integers
modulo NN. From number theory Entry 1, Entry 12, Entry 13, Entry 14, it can
be shown that if the modulus is a prime number, a base (called a
primitive root) exists such that a form of integer logarithm can be
defined. This is stated in the following way. If NN is a prime
number, a number rr called a primitive roots exists such that the
integer equation
n
=
(
(
r
m
)
)
N
n
=
(
(
r
m
)
)
N
(3)
creates a unique, one-to-one map of
the N-1N-1 member set m={0,...,N-2}m={0,...,N-2} and the N-1N-1 member
set n={1,...,N-1}n={1,...,N-1}. This is because the multiplicative group
of integers modulo a prime, pp, is isomorphic to the additive group
of integers modulo (p-1)(p-1) and is illustrated for N=5N=5 below.
Table of Integers n=((rm))n=((rm)) modulo 5, [* not defined]
| r |
m= |
0 |
1 |
2 |
3 |
4 |
5 |
6 |
7 |
| 1 |
|
1 |
1 |
1 |
1 |
1 |
1 |
1 |
1 |
| 2 |
|
1 |
2 |
4 |
3 |
1 |
2 |
4 |
3 |
| 3 |
|
1 |
3 |
4 |
2 |
1 |
3 |
4 |
2 |
| 4 |
|
1 |
4 |
1 |
4 |
1 |
4 |
1 |
4 |
| 5 |
|
* |
0 |
0 |
0 |
* |
0 |
0 |
0 |
| 6 |
|
1 |
1 |
1 |
1 |
1 |
1 |
1 |
1 |
Table 1 is an array of values of rmrm modulo NN and it is
easy to see that there are two primitive
roots, 2 and 3, and equation Equation 3 defines a permutation of
the integers nn from the integers mm (except for zero). Equation
Equation 3 and a primitive root (usually chosen to be the smallest
of those that exist) can be used to convert the DFT in Equation 1
to the convolution in Equation 2. Since Equation 3 cannot give a
zero, a new length-(N-1) data sequence is defined from x(n)x(n) by
removing the term with index zero. Let
n
=
r
-
m
n
=
r
-
m
(4)
and
k
=
r
s
k
=
r
s
(5)
where the term with the negative exponent (the inverse) is
defined as the integer that satisfies
(
(
r
-
m
r
m
)
)
N
=
1
(
(
r
-
m
r
m
)
)
N
=
1
(6)
If NN is a prime number, r-mr-m always exists. For
example, ((2-1))5=3((2-1))5=3. Equation Equation 1 now becomes
C
(
r
s
)
=
∑
m
=
0
N
-
2
x
(
r
-
m
)
W
r
-
m
r
s
+
x
(
0
)
,
C
(
r
s
)
=
∑
m
=
0
N
-
2
x
(
r
-
m
)
W
r
-
m
r
s
+
x
(
0
)
,
(7)
for s=0,1,..,N-2s=0,1,..,N-2, and
C
(
0
)
=
∑
n
=
0
N
-
1
x
(
n
)
C
(
0
)
=
∑
n
=
0
N
-
1
x
(
n
)
(8)
New functions are defined, which are simply a permutation in the
order of the original functions, as
x
'
(
m
)
=
x
(
r
-
m
)
,
C
'
(
s
)
=
C
(
r
s
)
,
W
'
(
n
)
=
W
r
n
x
'
(
m
)
=
x
(
r
-
m
)
,
C
'
(
s
)
=
C
(
r
s
)
,
W
'
(
n
)
=
W
r
n
(9)
Equation Equation 7 then becomes
C
'
(
s
)
=
∑
m
=
0
N
-
2
x
'
(
m
)
W
'
(
s
-
m
)
+
x
(
0
)
C
'
(
s
)
=
∑
m
=
0
N
-
2
x
'
(
m
)
W
'
(
s
-
m
)
+
x
(
0
)
(10)
which is cyclic convolution of length N-1 (plus x(0)x(0)) and is
denoted as
C
'
(
k
)
=
x
'
(
k
)
*
W
'
(
k
)
+
x
(
0
)
C
'
(
k
)
=
x
'
(
k
)
*
W
'
(
k
)
+
x
(
0
)
(11)
Applying this change of variables (use of logarithms) to the DFT
can best be illustrated from the matrix formulation of the DFT.
Equation Equation 1 is written for a length-5 DFT as
C
(
0
)
C
(
1
)
C
(
2
)
C
(
3
)
C
(
4
)
=
0
0
0
0
0
0
1
2
3
4
0
2
4
1
3
0
3
1
4
2
0
4
3
2
1
x
(
0
)
x
(
1
)
x
(
2
)
x
(
3
)
x
(
4
)
C
(
0
)
C
(
1
)
C
(
2
)
C
(
3
)
C
(
4
)
=
0
0
0
0
0
0
1
2
3
4
0
2
4
1
3
0
3
1
4
2
0
4
3
2
1
x
(
0
)
x
(
1
)
x
(
2
)
x
(
3
)
x
(
4
)
(12)
where the square matrix should contain the terms of WnkWnk. For
clarity, only the exponents nknk are shown. Separating the x(0)x(0)
term, applying the mapping of Equation 9, and using the primitive
roots r=2r=2 (and r-1=3r-1=3) gives
C
(
1
)
C
(
2
)
C
(
4
)
C
(
3
)
=
1
3
4
2
2
1
3
4
4
2
1
3
3
4
2
1
x
(
1
)
x
(
3
)
x
(
4
)
x
(
2
)
+
x
(
0
)
x
(
0
)
x
(
0
)
x
(
0
)
C
(
1
)
C
(
2
)
C
(
4
)
C
(
3
)
=
1
3
4
2
2
1
3
4
4
2
1
3
3
4
2
1
x
(
1
)
x
(
3
)
x
(
4
)
x
(
2
)
+
x
(
0
)
x
(
0
)
x
(
0
)
x
(
0
)
(13)
and
C
(
0
)
=
x
(
0
)
+
x
(
1
)
+
x
(
2
)
+
x
(
3
)
+
x
(
4
)
C
(
0
)
=
x
(
0
)
+
x
(
1
)
+
x
(
2
)
+
x
(
3
)
+
x
(
4
)
(14)
which can be seen to be a reordering of the structure in
Equation 12. This is in the form of cyclic convolution as indicated
in Equation 10. Rader first showed this in 1968 Entry 12, stating
that a prime length-N DFT could be converted into a length-(N-1)
cyclic convolution of a permutation of the data with a permutation
of the W's. He also stated that a slightly more complicated
version of the same idea would work for a DFT with a length equal
to an odd prime to a power. The details of that theory can be found
in Entry 12, Entry 10.
Until 1976, this conversion approach received little
attention since it seemed to offer few advantages. It has
specialized applications in calculating the DFT if the cyclic
convolution is done by distributed arithmetic table look-up
Entry 5 or by use of number theoretic transforms
Entry 1, Entry 12, Entry 13. It and the Goertzel algorithm
Entry 16, Entry 3 are efficient when only a few DFT values need to be
calculated. It may also have advantages when used with pipelined or
vector hardware designed for fast inner products. One example is the
TMS320 signal processing microprocessor which is pipelined for inner
products. The general use of this scheme emerged when new fast
cyclic convolution algorithms were developed by Winograd
Entry 21.
The DFT of x(n)x(n) evaluates the Z-transform of x(n)x(n) on NN equally
spaced points on the unit circle in the zz plane. Using a nonlinear
change of variables, one can create a structure which is equivalent
to modulation and filtering x(n)x(n) by a “chirp" signal.
Entry 2, Entry 20, Entry 19, Entry 16, Entry 18, Entry 3.
The mathematical identity (k-n)2=k2-2kn+n2(k-n)2=k2-2kn+n2 gives
n
k
=
(
n
2
-
(
k
-
n
)
2
+
k
2
)
/
2
n
k
=
(
n
2
-
(
k
-
n
)
2
+
k
2
)
/
2
(15)
which substituted into the definition of the DFT in Multidimensional Index Mapping: Equation 1 gives
C
(
k
)
=
{
∑
n
=
0
N
-
1
[
x
(
n
)
W
n
2
/
2
]
W
-
(
k
-
n
)
2
/
2
}
W
k
2
/
2
C
(
k
)
=
{
∑
n
=
0
N
-
1
[
x
(
n
)
W
n
2
/
2
]
W
-
(
k
-
n
)
2
/
2
}
W
k
2
/
2
(16)
This equation can be interpreted as first multiplying (modulating) the data
x(n)x(n) by a chirp sequence (Wn2/2Wn2/2, then convolving (filtering) it, then
finally multiplying the filter output by the chirp sequence to give the DFT.
Define the chirp sequence or signal as h(n)=Wn2/2h(n)=Wn2/2 which is called
a chirp because the squared exponent gives a sinusoid with changing frequency.
Using this definition, Equation 16 becomes
C
(
n
)
=
{
[
x
(
n
)
h
(
n
)
]
*
h
-
1
}
h
(
n
)
C
(
n
)
=
{
[
x
(
n
)
h
(
n
)
]
*
h
-
1
}
h
(
n
)
(17)
We know that convolution can be carried out by multiplying the DFTs of the signals,
here we see that evaluation of the DFT can be carried out by convolution. Indeed,
the convolution represented by ** in Equation 17 can be carried out by DFTs (actually
FFTs) of a larger length. This allows a prime length DFT to be calculated by a
very efficient length-2M2M FFT. This becomes practical for large NN when a particular
non-composite (or NN with few factors) length is required.
As developed here, the chirp z-transform evaluates the z-transform at equally spaced
points on the unit circle. A slight modification allows evaluation on a spiral and
in segments Entry 19, Entry 16 and allows savings with only some input values are nonzero or
when only some output values are needed. The story of the development of this
transform is given in Entry 18.
Two Matlab programs to calculate an arbitrary length DFT using the chirp z-transform
is shown in Figure 1.
Goertzel's algorithm Entry 7, Entry 3, Entry 15 is another methods that
calculates the DFT by converting it into a digital filtering problem. The
method looks at the calculation of the DFT as the evaluation of a
polynomial on the unit circle in the complex plane. This evaluation is
done by Horner's method which is implemented recursively by an IIR filter.
The polynomial whose values on the unit circle are the DFT is a slightly
modified z-transform of x(n) given by
X
(
z
)
=
∑
n
=
0
N
-
1
x
(
n
)
z
-
n
X
(
z
)
=
∑
n
=
0
N
-
1
x
(
n
)
z
-
n
(18)
which for clarity in this development uses a positive exponent .
This is illustrated for a length-4 sequence as a third-order
polynomial by
X
(
z
)
=
x
(
3
)
z
3
+
x
(
2
)
z
2
+
x
(
1
)
z
+
x
(
0
)
X
(
z
)
=
x
(
3
)
z
3
+
x
(
2
)
z
2
+
x
(
1
)
z
+
x
(
0
)
(19)
The DFT is found by evaluating Equation 18 at z=Wkz=Wk, which
can be written as
C
(
k
)
=
X
(
z
)
|
z
=
W
k
=
D
F
T
{
x
(
n
)
}
C
(
k
)
=
X
(
z
)
|
z
=
W
k
=
D
F
T
{
x
(
n
)
}
(20)
where
W
=
e
-
j
2
π
/
N
W
=
e
-
j
2
π
/
N
(21)
The most efficient way of evaluating a general polynomial without any
pre-processing is by “Horner's rule" Entry 11 which is a nested
evaluation. This is illustrated for the polynomial in Equation 19 by
X
(
z
)
=
[
x
(
3
)
z
+
x
(
2
)
]
z
+
x
(
1
)
z
+
x
(
0
)
X
(
z
)
=
[
x
(
3
)
z
+
x
(
2
)
]
z
+
x
(
1
)
z
+
x
(
0
)
(22)
This nested sequence of operations can be written as a linear difference
equation in the form of
y
(
m
)
=
z
y
(
m
-
1
)
+
x
(
N
-
m
)
y
(
m
)
=
z
y
(
m
-
1
)
+
x
(
N
-
m
)
(23)
with initial condition y(0)=0y(0)=0, and the desired result being
the solution at m=Nm=N. The value of the polynomial is given by
X
(
z
)
=
y
(
N
)
.
X
(
z
)
=
y
(
N
)
.
(24)
Equation Equation 23 can be viewed as a first-order IIR filter with the
input being the data sequence in reverse order and the value of the
polynomial at zz being the filter output sampled at m=Nm=N. Applying
this to the DFT gives the Goertzel algorithm Entry 17, Entry 15 which is
y
(
m
)
=
W
k
y
(
m
-
1
)
+
x
(
N
-
m
)
y
(
m
)
=
W
k
y
(
m
-
1
)
+
x
(
N
-
m
)
(25)
with y(0)=0y(0)=0 and
C
(
k
)
=
y
(
N
)
C
(
k
)
=
y
(
N
)
(26)
where
C
(
k
)
=
∑
n
=
0
N
-
1
x
(
n
)
W
n
k
.
C
(
k
)
=
∑
n
=
0