A Guided Tour of The Fast Fourier Transform
A Guided Tour of The Fast Fourier Transform
A Guided Tour of The Fast Fourier Transform
Fourier transform
The fast Fourier transform algorithm can reduce the time
involved in finding a discrete Fourier transform from several
minutes to less than a second, and also can lower the
cost from several dollars to several cents
G. D. Bergland
For some time the Fourier transform has served as a
bridge between the time domain and the frequency
domain. It is now possible to go back and forth between waveform and spectrum with enough speed and
economy to create a whole new range of applications
for this classic mathematical device. This article is intended as a primer on the fast Fourier transform,
which has revolutionized the digital processing of
waveforms. The reader's attention is especially
directed to the IEEE Transactions on Audio and
Electroacoustics for June 1969, a special issue devoted to the fast Fourier transform.
This article is written as an introduction to the fast
Fourier transform. The need for an FFT primer is apparent from the barrage of questions asked by each new
person entering the field. Eventually, most of these questions are answered when the person gains an understanding of some relatively simple concept that is taken for
granted by all but the uninitiated. Here the basic concepts will be introduced by the use of specific examples.
The discussion is centered around these questions:
1. What is the fast Fourier transform?
2. What can it do?
3. What are the pitfalls in using it?
4. How has it been implemented?
Representative references are cited for each topic
covered so that the reader can conveniently interrupt this
fast guided tour for a more detailed study.
What is the fast Fourier transform?
The Fourier transform has long been used for characterizing linear systems and for identifying the frequency components making up a continuous waveform.
However, when the waveform is sampled, or the system
is to be analyzed on a digital computer, it is the finite, disIEEE
x(t) =
X(J)ei2tft dl
(2)
1N-i
X(j) = - E x(k)ei2j /N
N k 0
(3)
x(k) = E X(I)ei2Xjk/N
(4)
N-i1
j=0
for]j = 0, 1, * * , N - 1; k = 0, 1, * * , N - 1. Both
41
1
X(j) =-N E x(k) WN-k
(5)
k=o
x(k)
N-1
E X(J) WNk
(6)
x(k A T)
An example of a real-valued time series and its associated DFT is shown in Fig. 1. The time series x(kAT)
is assumed to be periodic in the time domain of period T
seconds, and the set of Fourier coefficients X(jfo) is assumed to be periodic over the sample frequency fi. Only
one complete period of each function is shown.
The fundamental frequency Jo and the sample period
AT do not appear explicitly in Eqs. (5) and (6), but each
j should still be interpreted as a harmonic number and
each k still refers to a sample period number. That is, the
0true frequency is the product of j and fo and the true
time is the product of k and AT.
When the x(k) series is real, the real part of X(j) is
symmetric about the folding frequency ff (where f =
and the imaginary part is antisymmetric. Since X(j)
has
interpreted as being periodic, these symmetries
ha been
enitrrtda
engproi,teesmere
are equivalent to saying that the real part of X(j) is an
even function, and that the imaginary part of X(j) is an
odd function. This also means that the Fourier coefficients between N/2 and N 1 can be viewed as the
/"negative frequency" harmonics between - N/2 and -1.
Likewise, the last half of the time series can be interpreted
as negative time (that is, as occurring before t 0= ).
Derivation of the Cooley-Tukey FFT algorithm. A derivation of the Cooley-Tukey FFT algorithm8 for evaluating Eq. (6) is given in this section for the example of
N = 8. This derivation is also appropriate to the forward
transform, since Eq. (5) can be rewritten in the form9
if/2)
I
___ I / .
t=O
k =0O
AT
1
2 AT
2
3 AT
3
* *
* * *
NAT
N
X(jf0)
-.JX
=
gX(X)
f=0
fo
2fo
j =0
fo = 1/T
fs
NFo
3fo... ff...
3
N/2
x[IEx(k)*WNik]
(7)
~~
where the asterisk refers to the complex conjugate operation. Alternatively, the FFT algorithm used for comput-
1/AT
FIGURE 2. A flow diagram of the Cooley-Tukey FFT algorithm for performing an eight-point transform.
AO(OOO)
Al(000)
0
A2(000)
AO(001)
Al(001)
0414
A2(001)
AO(010)
A0(01 1)
A1(010)
0I4
A2(010)
Al(011)
AO(101)
AO(100)
A2(011)
Al(100)
Al(101)
04
A2(100)
A2(101)
Al(110)
Al(111)
67
A2(110)
A2(111)
:=l
no
42
A0(1 11)
A0(1 10)
A(00
A(01
^300
^301
X(OOO)
X(OO1)
X(O10)
X(O1
A(0)
(O)1)
A(1)
X()
A(10
A(1)
X101)
X(1 11)
N-I
X(j) =E A(k) Wk
(8)
k 0
X(j2,jl,jO)
1
= Z
1
E
1
Z
A(k2,k,ko)W(324J1a+)o)(k24k2+k^)
ko=O k1=O
k2= 0
Noting that
Wm'+
w -
(10)
Wn, we have
W(j24+i12+io)k24W(i24+12+jo)ki2W(i24+j)2+io);o (11)
W
If we look at these terms individually, it is clear that
they can be written in the form
w(i24+ji2+io)124 = [W8(i22+ii)k2] Wiok24
(12)
-L[W8i2kli w(i12+io)ki
(13)
W(i24+i12+-o)ko= W(j24+i12+Jo)ko
(14)
[e-27/8]8= e2r=
-;Q2Jl,JO)
ko=O
E
ko =0 A2(joJ,]ko)
X(jjljo)
A3(jOJ2)
Wj4+1
i)k
(324i1aJo)o
(19)
(20)
What
can it do?
The operations usually associated with the FFT are:
(1) computing a spectrogram (a display of the shortterm power spectrum as a function of time); (2) the convolution of two time series to perform digital filtering;
and (3) the correlation of two time series. Although all of
these operations can be performed without the FFT,
its computational savings have significantly increased the
interest in performing these operations digitally.
Spectrograms. The diagram in Fig. 4 represents a
method of obtaining estimates of the power spectrum of a
time signal through the use of the fast Fourier transform.
(15)
Ed
=.....d
.Eki=O
A3(io,J2) =
FIGURE 3. The number of operations required for computing the discrete Fourier transform using the FFT algo-
I
Wjok24 W(ji2+jo) ki2w(Y4+Yi2+Yo)kg
EA(k2,ki,ko)
1024
k2=
A i(jo,ki,ko)
A2(jo,l,ko)I
Aa(joj,j.j2)
A3
(16)
Ai(]o,ki,ko) = E A(k2,k1k0)WJok24
(17)
(18)
12=0
ki=0
.0~
~~
Direct calculation
of DFT
0~
-/
X
512
256_
DFT via FFT
128 _/
64\
64 128 256
512
N
1024
43
In this case the square of the magnitude of the set of complex Fourier coefficients (that is, the periodogram) is used
to estimate the power spectrum of the original signal.
A snapshot of the spectrum of the signal can always be
computed from the last T seconds of data. By taking a
series of these snapshots, estimates of the power spectrum
can be displayed as a function of time as shown in Fig. 5.
When the audio range is displayed, this is usually called a
sound spectrogram or sonogram.
XR(jfOI
XI (ifo)
'X(ifo) 2
O
12
N--r=r_ --r-
-r iT
FIGURE 4. The power spectrum of a real function computed by taking the sum of the squares of the real and
imaginary components of the discrete Fourier transform
Fourier coefficients.
| t | | | /l
V | K1
- l -,| i1 '-,
l
i l l
|E
||
/
1
;
| i
I NE
r(k)
N(k)
(21)
(22)
where C() is the DFT of c(k), G(j) is the DFT of g(k), and
t~~~~~c()
g(r)h(k T)
g(ET)W i )(N
hE 6) W _j
(24)
| | l |!
| |
| |
~~~~IV1N- N- I
Ik
N-
N IrI 11
1111gIl|rg(r)h(f)
W+j(/
[N
E -t-^r)](25)
I 1111 I I l i 11 li * l g _lV=li
FIGURE 6. The response of a linear system to a driving
function g(k) expressed as the convolution of the input
sgnal with the impulse response of the system.
Digital filter
g(k)
h(r)
c(k) =
Time
44
h>
N-1
:cik)
g(r) h(k - .)
= g*h
IEEE spectrum JULY 1969
i=O
(26)
= k - r, for
otherwise
(27)
C(j)
G(j).H*(j)
(28)
ceding s
welobtant sult
1
N-
g(r)h(T- k)
c(k) = N
-E
=o
(29)
tion involves.
In Fig. 7, the results of the DFT are treated as a corrupted estimate of the CFT. The example considers a
cosine-wave input, but the results can be extended to any
function expressed as a sum of sine and cosine waves.
Line (a) of Fig. 7 represents the input signal s(t) and its
continuous Fourier transform S(f). Since s(t) is shown to
be a real cosine waveform, the continuous Fourier transform consists of two impulse functions that are symmetric
/\J
\t
1V
t\
(a)
W(f)
w(t)
(b)
T/2
-T/2
S(f) W(f)
s(t) . w(t)
(C)
0\
'
C(f)
.4
A
s(k)
1 11t
lSr
*
llSl 21
1
''~j~
-2tem
u-N termsA
(f)
c(t)
(d)
(e)
AA
A
S(j)
* .t
, *j
,T
as*s^.*l,X
45
of leakage, aliasing, and picket-fence effects are associated with variations from this ideal condition.
Since both s(k) and S(j) are periodic, with a period N,
any set of N adjacent terms can be used to characterize
either set. The magnitude of S(j) is the same for samples
of si(k) that are symmetric about the origin as it is for the
DFT of a set of samples starting at the origin. The
change in the phase of S(j) for different time origins can
be determined by application of the DFT shifting
theorem.9
f(t)
1'
46
-f
(k) = 0
(30)
.rgnlst
The
additional Fourier coefficients that result are interleve wit th
AvAs shown
wit tin Fig.iginal
set.
13, the ripple in the power spectrum
has been reduced fro
60 percent to 20
approximately 60
reduce from apomately
per to 2020
harben
larger
smaller
percent.ethe pe e
or
than
zeros.
In practice the picket-fence problem is not as great as
this discussion implies. In many cases the signal being
processed will not be a pure sinusoid but will be broad
enough to fill several of the original spectral windows.
Moreover, the use of any data window other than the
rectangular (or boxcar) data window discussed here,
1_
0J
L
9N
N-
10
(picket-fence effect).
Independent filters
1.0
0.6
0___l_V_V_V_
0
1
2
3
V_V_V_V_\_\_
4
5
6
Harmonic number
jPower
1.0 F
0.
>f
response
\f/\/I\lI\/I\I\/I\/I\
_F IV Y51V V V
1
0
Yg Vg
f
47
;;,
1.0
0o9 A A NAAAAMAAAAAA
1.0
if
>
0.8
f
FIGURE 13. The reduction of the picket-fence effect,
brought about by computing redundant overlapping sets of
Fourier coefficients.
>
g(t
0 < k <N
N < k < 2N
= g(k)
'(k) = 0
(31)
Filtered by
h/t
Past 4
> Past
Future
Present
Present
g(T)
h(t-r)~~ht-7
- FutuJre
Past
Future|
Past
present
Present
IEEE spectrum JULY 1969
k(k) = g(k)
=
(32)
< kk<2N(32
< 2N
O0.
0 < k<N
(33)
2N
<2N
= 0N < k <
A
A AA
/\/\/\ ~~~~~~~(r)~
h(t)
g(t)
L1
-.
-NM
FIGURE 17. A method for convolving a finite impulse response with an infinite time function by performing a series
samples
Digital filtering
A
g
g(t)
g(t)
l
1
2 ll
4
3
Al'6
8
7
htr)
Ksamples
--
~ ~-
2N terms-
H- 2T seconds ---
49
H(f)
+f
-f
Software. The number of variations of the FFT algorithm appears to be directly proportional to the number
of people usingit. 3 4-53Most of these algorithms are based
on either the Cooley-Tukey or the Sande-Tukey algorithm," but are formulated to exploit properties of the
series being analyzed or properties of the computer being
implies
h(t)
'Nvs
Jr-\ /
\ , -'
FIGURE 19. The (sin x)/x impulse response of a filter implied when a rectangular frequency response is specified.
FIGURE 20. Block diagrams illustrating the use of the FFT
for simulating a radar signal processing system.
used.
Equations (17) through (20) represent the CooleyTukey formulation of the FFT algorithm. By separating
the components of ] instead of the components of k in
Eq. (11), the Sande-Tukey equations would have resulted.
In either case, the recursive equations can be written in
the form of a two-point transform followed by a rereferencing (or twiddling9) operation. The resulting
Aj(jo,k,,ko)
A2(jo,jl,ko)
{E
J -i
(35)
Dechirp
s(t)
Taylor
~~~~filter
weighting
f(t)
w(t)
O
r()ko=
x(t)s(t)
50
s(t) Y
~~~~Dechirp
~~~filter
Taylor
weightin
r(t) = [l - x(t)]
A2(]o,l]ko) W2i2jo}
A3(jo,ji,js)
(36)
(37)
>Cr
A3(jOjl,j2)
r(t)
(38)
_k = k2(r,r3) + k1(r:3) + ko
* , r, - 1;]', k1 = 0, 1, * * , r2 -1;
where]jO, k2 = 0, 1,
X(j) =
X(j)
The author would like to thank D. E. Wilson for the loan of the
title;
Kaenel and W. W. Lang for their M. Rader, and
and B.R.P. A.
Bogert, W. T. Hartwell, H. D. Helms, C.encouragement;
P. T. Rux for suggesting several improvements, which were incorporated in this article.
i N-1
x(k) W
N k=O
(39)
REFERENCES
NXo x(k)W-jk+[(k2-k2+j2ij2)/21
(40)
tice-Hall, 1966.
X(j) = W
Q)=
(41)
X(j) Q
= -
E g(k)h(]-k)
(42)
trans-
Paper
RC-1743, Feb. 1967.
8. Cooley, J. W., and Tukey,
cross-cepstrum and saphe-cracking," Time Series Analysis, Murray Rosenblatt, ed. New York: Wiley, 1963, pp. 201-243.
14. Noll, A. Michael, "Short-time spectrum and 'cepstrum' tech-
niques for vocal-pitch detection," J. Acoust. Soc. Am., vol. 36, pp.
296-302, 1964.
15. Oppenheim, A. V., Schafer, R. W., and Stockham, T. G., Jr.,
"sNonlinear filtering of multiplied and convolved signals,"
Proc. IEEE, vol. 56, pp. 1264-1291, Aug. 1968. (Reprinted in
IEEE
Trans. Audio and Electroacoustics, vol. AU-16, pp. 437-465,
Sept. 1968.)
16. Parzen, -E., "Statistical spectral analysis (single channel case)
in 1968," Tech. Rep. 11, ONR Contract Nonr-225(80)(NR-042234), Stanford University, Dept. of Statistics, Stanford, Calif.,
17. Richards, P. I., "Computing reliable power spectra," IEEE
51
IS. Tukey. J. W., "An introduction to the measurement of spectra," in Probabiliti anid Statistics, Ulf Greniander, ed. New York:
Wiley, 1959, pp. 300-330.
19. Tukey, J. W., "An initroduction to thle calculations of niumerical spectrumn analysis," in Spectral Anati-sis of Time Series, Beriiard Harris, ed. New York: Wiley, 1967, pp. 25-46.
20. Welch, P. D., "A direct digital method of power spectrum estimationi," IBM J. Res. Develop., vol. 5, pp. 141-156, 1961.
21. Welch, P. D., "Thie use of the faist Fourier transform for the
estimiation of power spectra: a methodi based on timne averaging
over shiort, modified periodogratms," IEEE Traits. Aud(io anid
Flectroacouistics, vol. AU-I5, pp. 70-73. Junie 1967.
Use of the fiast Foutrier trantsformi int couvoluttioo, correlatio,,, dligital
filterinig, etc.
22. Cooley, J. W., "Applications of thec fast Fourier trantsform
miethod." Proc. IBM Scienitific Conipuptintg S.t'nip. oni Digital Si,nutlatio,, of Co,,ttittuons Srvstemns, Thomas J. Wattsoni Researchi Center,
Yorktowni Heights, N.Y., 1966.
23. Cooley, J. W., Lewis, P. A. W., anid Welch., P. D., "Applicationi of the faist Fourier trainsformi to comiputation of Fourier integratls, Fourier series, anid conivolution initegrals." IEEE Tranis.
Auedio anid Electroaconstics, vol. AU-IS. pp. 79-84, Junie 1967.
24. Helmis, H. D.. "Fast Fourier tranisform miethod of computing
difrerenice equations and simulating tilters." IFEE Trants. Auldio
anid Electroacou(stics, vol. A U-IS5, pp. 85-90, June 1967.
25. Rabiner, L. R., Schiafer, R. W., anid Rader, C. M.. "Thie chiirp
Z-tranisforni algorithmi antd its applications.'" Bell SyVstemi Tech.
J., vol. 48. pP). 1249-1292, May-June 1969.
26. Sanide, G., "Oni ani alternativec miethiod of calculatinig covarianice funictions," unipublished techniical niote, Prineceton Uniiversity,
Priniceton, N.J., 1965.
27. Sinigleton., R. C., "Algorithim 345. ani Algol conivolution p)roceduire based oni the fakst Fourier tranisformi," Cominuniii. Assoc.
.
Contipit. Mach.. vol. 12. Mar. 1969.
28. Stockhami, T. G., "High-speed conivolutioni anid correlation."
1966 S5printg Jointt Conmputer Coutf.. AFIPS Proc., vol. 28. Wastiinigton, D.C.: Spatrtant Books. pp. 229-233.
The picket-ftnuie e/i'ct
29. Hartwell. W. T., "An alternate approach to thec use of discrete
Fourier tranisformis." to be publishied.
44. Gentleman, W. M., anid Sande, G., "Fast Fourier transforms -ror funi anid profit." 1966 Fall Jointt Comiputer Conf.,
AFIPS Proc., vol. 29. Washington. D.C.: Spartani Books.
45. Good, I. J., "The interactioni algorithm and pracetical Fourier
series," J. Royv. Statist. Soc., vol. 20, series B, pp. 361-372, 1958;
addendum, vol. 22, pp. 372-375, 1960.
46. Pease, M. C., "Ani adaption of the fast Fourier tranisform for
parallel processing," J. Assoc. Comiput. Mach., vol. 15, pp. 252264. Apr. 1968.
47. Rader, C. M., "Discrete Fourier transformns wheni the number
of data samples is prime," Proc. IEEE, vol. 56, pp. 1107-1108,
June 1968.
48. Sanide. G., "Arbitrary radix one-dimensionial fast Fourier
transformn subroutinies," University of Chicago. Ill., 1968.
49. Sinigleton., R. C.. "On computing the fast Fourier trainsform." CTomninii. Assoc. Comtpuit. Maclh., vol. 10, pp. 647-654,
Oct. 1967.
50. Singleton. R. C.. "Algorithml 338. Algol procedures for the
fitst Fourier tranisformi." Com,in,mi. Assoc. Comnput. Mach., vol. II,
pp. 647-654, Nov. 1968.
SI1. Singletoni, R. C., "Algorithim 339, an Algol procedure for the
faist Fourier tranisform with arbitrary factors," Comm,uuu. Assoc.
Comput. Maclh., vol. I 1, pp. 776-779, Nov. 1968.
52. Sinigletoii, R. C... "Algorithmn 345. an Algol conivolutionl procedure based oni the fast Fourier tranisformn," Commntii. Assoc.
Comtpuit. Mfach., vol. 12, Mar. 1969.
53. Yavne. R., "Ani econiomiical mnethiod for catlculatinig thec discrete Fourier tranisform," 1968 Fatll Joinit Computer Con!f, 1FfIPS
Proc.. vol. 33. Washinigton., D.C.: Spartani Books. pp. 115-125.
Fatouirrnsrmhdae
Berglanid,
hardware
tin-
plementattionis
-at survev," IEEE Trails. Audtio and Electroacontstics, vol. AU-17, Junie 1969.
Berglanid,
"Fatst
hairdware
56.
G. D.,
Fourier transform
implementations---an overview,"' IEEE Trants. Aufdio anid Electroacontstics, vol. AU-17. Junte 1969.
57. McCullough, R. B.. "A real-time digital spectrumi anialyzer,"
Frequenc~r-doinain filter (lesigitStaniford Electronics Laboratories Sci. Rept. 23, Stanford UniiFrequenci-domnninfilter rlesign
~~~versity. Calif., Nov. 1967.
30. Helmis. H. D., "Fast Fourier transform miethod for compuiting
M. C., Ill,
Goldberg. J.,
study of a
58.
dilferencec equattionis anid simulatinig filters,"~IEEE Tratits. Autdio
spciluroedgtlom trfro--inFuiranlys,
aiidElecroaoitsies,vol AU-5, l). 5-90 Juie
No. 989. Advanced Research Projects Agency, Washington,
31. Helms, H. D., "'Nonirecursive digitail filters: design miethiods
D.C., Maty 1967.
for achieving specificattionis oni frequenicy response." IEEE Trunis.
in real
processor to
59. Shively, R. R., "A
Aiti(iaii Elctrttettsies vl. U-1, p. 36-32, ept 198.
time," IEEE Tranis. Comnputers, vol. C-17, pp. 485-491, May 1968.
32. Kaiser, J. F., "Digital filters," in Sriste,ip Anialrsis bhi Digital
Fourier
Bell
60. Smith, R. A., "A
TlpoeLbrtre,Ic,Wipiy
.. processor,"
Comiputter, F. F. Kuo anid J. F, Kaiser, eds. New York: Wiley,
97
Peatse,
167.Order
anid
"Featsibility
digitatl
faist
genierate spectrat
tranisform
1966.TleloiLaoztrs,Ic.WhpiiyN..197
34. Berglanid, G. D., '"Tlie fast Fourier tratisformn recursive equations for arbitrary lenigthi records," Maith. Comiput., vol. 21, pp.
236--238, Apr. 1967.
35. Berglanid, G. D., "~A fast Fourier trantsform algorithmn usinig
base 8 iterattionis.' Math. Compupit., vol. 22, pp. 275-279, Apr. 1968.
36. Berglanid, G. D., "A fatst Fourier transform algorithim for
real-valued series." Comntu,iit. Assoc. Conitpit. Muich., vol. 11, pp.
703-710, Oct. 1968.
37. Berglanid. G. D., and Wilsoin, D. E., "Ani FFT algorithmt for a
globail, highly pairallel processor," IEEE Traits. Audcio amid Electroacoustics, vol. AU-17. Junie 1969).
38. Bluestein., L. I., "'A liniear filtering approach to the computationi of the discrete Fourier trainsforml," 1968 NEREM Record,
pP). 218-219.
39. Brennier, N. M., "~Three Fortrani programiis thatt performi thle
Coolcv-Tukey Fourier trainsformi," Tech. Note 1967-2, Lincoln
Laiborattory, M.I.T., Lexingtoni, Mass., July 1967.
40. Cooley, J. W., anid Tukey, J. W.. "An algorithm for the
machinie calculation of complex Fourier series," Math. C'omtput.,
vol. 19, pp. 297-301, Apr. 1965.
41. Cooley, J. W., "Hairmionic anialysis complex Fourier series,"
SHARE Doe. 3425, Feb. 7, 1966.
42. Coolcy, J. W., "Complex finiite Fourier transform subroutinec,"
SHARE Doc. 3465, Sept. 8, 1966.
43. Danielsoni, G. C.. anid Laniezos. C.. "Some
in
practical Fourier anialysis and their applicationi to X-ray scattering
from liquids," J. Franiklini liust., vol. 233, pp. 365-380, 435-452.
improve!ments
52
Other
61. Gilbert, S. M., Privatte commiunicationi, Bell Teleplione Laboratories, Inc., Whiippaniy. N.J.
62. Klauder, J. R., Price. A. C.. Darlingtoni. S., antd Albersheim,
W. J., "Thie theory and designi of chiirp radars.,' Bell Systemn Tech.
J. vol. 39, pp. 745-808. Julv 1960.
G. D.
grees from Iowa State University in 1962, 1964, and 1966, respectively. While at Iowa State he was a teaching assistant
inteEcrcaEgneigDprmntndaesrh
asistanletrincth Engineering DexperimentSation. Frome1964
tos1966
he thelaEnationaleSienxeienFoundation.
traimeeship
t 96h edaNtoa cec onaintanehp
1966
he
In
joined the Digital Systems Laboratory at Bell
Telephone Laboratories, Inc., Whippany, N.J., where he
data-processing techniques. He
is a member of Sigma Xi, Phi
Kappa Phi, Eta Kappa Nu, and