1. technology and computing

IEEE TRANSBCTIONS ON AUDIO AND ELECTROACOUSTICS,
VOL. AU-15,NO.
2,
JUNE
1967
45
W h a t is the Fast Fourier Transform?
G-AE Subcommittee on Measurement Concepts
WILLIAM T. COCHRAN
JAMES W. COOLEY
DAVIb L. FAVIN, MEMBER, B E E
HOWARD D. HELMS, MEMBER, IEEE
REGINALD A. KAENEL, SENIOR MEMBER, IEEE
WILLIAM w. LANG, SENIOR MEMBER, E E E
GEORGE C. ML4L1NG, JE., ASSOCIATE B B ~ ~ ~ B EIEEE
R,
DAVIb E. NELSON, MEMBER, BEE,
CHARLES M. RADER, MEMBER, IEEE
PETER D. WELCH
Absfracf-The fastFouriertransform is a computational tool
which facilitates signal analysissuch as power spectrum analysis and
filter simulation by means of digital computers. It is a method for
efficiently computing the discrete Fourier transform of a series of
data samples (referred to
a s a time series). In this paper, the discrete
Fourier transfom of a time series is defined, some its properties
are discussed, the associated fast method (fast Fourier transform)
for computing this transform is derived, and same of &e compai.ational aspects the method arepresented. Examples are included to
demonstrate the concepts involved.
iblemappingoperationfortime
series. As the name
implies, it has mathematical properties that are entirely
analogous to those of the Fourier integral transform. In
particular, i t defines a spectrum of a time series; multiplicationof
thetransform
of twotimeseriescorresponds to convolving the time series.
If digitalanalysistechniques
aretobe
usedfor
analyzing a corltinuous waveform then i t is necessary
that the data be
sampled(usually at equallyspaced
INTRODUCTION
intervals of time) in order
to produce a time series of
discrete
samples
which
can
be fed into a digital comN ALGORITHM for the computation of Fourier
puter.
As
is
w
~
l
l
known
[
6
] , such a timeseriescomcoefficients which requires much less computapletely
represents
the
continuous
waveform, provided
tional effort than was required in the past was
this
waveform
is
frequency
band-limited
and
the
reported by Cooley and Tukey [ l ] in 1965. This method
samples
are
taken
a
t
a
rate
that
is
at
least
twice
the
is now widely known as the "fast Fourier transform,"
highest
frequency
present
in
the
waveform.
When
these
and has produced major changes in computational techsamples
are
eqeally
spaced
they
are
known
as
Nyquist
niques used in digital spectral analysis, filter simulation,
s a ~ ~ pIkt ~ . be shown that the DFT of such a time
andrefated fields. Thetechniquehasa
Iong and interesting history that has been summarized by Cooley, series is closely related to the Fourier transform of the
continuouswaveform fromwhich sampleshavebeen
Lewis, and Llielch in this issue[2].
takentoformthetime
series. This makes the DFT
The fast Fourier transform (FFT)
is a method for
particularly usefulforpower
spectrumanalysisand
efficiently computingthediscreteFouriertransform
(DFT)of a timeseries(discrete
data samples). The filter simulation on digital computers.
The fast Fourier transform (FFT), then,
a highly
efficiency of this method is such that solutions to many
D F T of a time
problems can now be obtained substantially more eco- efficient procedure for computing the
nomically than in the past. This is the
reason for the series. I t takes advantage of the fact that the calculation of the coefficients of t h e D F T can be carried out
very great current interest in this technique.
iteratively,
which results in a considerablesavingsof
The discrete Fourier transform (DFT) is a transform
computation
time, This manipulation is not intuitively
in its own right such as the Fourier integra,' traxwform
obvious,
perhaps
explainingwhythisapproachwas
or the Fourier series transform. I t is a pourerful reversoverlooked for such a longtime. Specifically,if the
timeseriesconsists
of N = 2" samples,thenabout
Manuscript received March 10,1967.
W. T. Cochran, D. L. Favip, and R. A. Kaenel are with Bell 2nN=2N*10g2 N arithmetic operations will be shown
Telephone Laboratories, Inc., Murray Hill, N. J.
to evaluate all N associated D F T COH. D. Helms iswith Bell Telephone Laboratories,IIK., Whippany, to be required
N.
efficients. I n comparison with the number of operations
J. W. Cooley and P. D. Welch are with the IBM Research Center, required for the calculationof t h e D F Tcoefficients with
Yorktown Heights, N. Y .
1%'. Lann and G. C. Malina are with the IBM CorDoration, straightforwardprocedures
( W ) , thisnumber
is so
Poughkeepsie, iiJ. Y .
small
when
N
is
large
as
t
o
completely
change
the
comC. M. Rader is with Lincoln Laboratory, Massachusetts Institute
of Technology, Lexington, Mass. (Operated with support from the
putationally economical approach to various problems.
U. S. Air Force.)
that for N=8192
D. E. Nelson is with the Electronics Division of the General Forexample, it has beenreported
Dynamics Corporation, Rochester, N. Y .
samples, the computations require about
fiveseconds
2 -
I
Authorized licensd use limted to: IE Xplore. Downlade on May 10,2 at 19:056 UTC from IE Xplore. Restricon aply.
IEEE
ELECTROACOUSTICS,
AND
TRANSBCTIOKS
AUDIOON
JUNE
46
for the evaluation ofall 8192 D F T coefficients on an
IBM 7094 computer. Conventional procedures take on
the order of half an hour.
T h e known applications wherea substantial reduction
incomputationtimehasbeen
achieved include: 1 )
computation of the power spectra and autocorrelation
functions of sampled data
2 ) simulation of filters
3j pattern recognition b y using a two-dimensional
form of the DFT; 4 ) computation of bispectra, crosscovariancefunctions,cepstraandrelatedfunctions;
and 5) decomposing of convolved functions.
THEDISCRETEFOURIER
TRANSFORM
(DFT)
Definition of the D F T and it:, Inverse
Since the FFT is an efficient method for computing
t h e D F T i t is appropriate to begin by discussing the
DFT and some of the properties that make it so useful
a transformation. The DFT is defined by1
N-1
A,
X k exp ( - 2 ~ r j r k I N ) r
M
0,
1 (1)
k=O
where A, is the rth coefficient of the DFT and X k denotes the kth sample of the time series which consists
of N samples and j = 47.
T h e Xk’s can be complex
numbers and the A,‘s are almost always complex. For
notationalconvenience ( 1 ) is often written as
1967
N-
exp (2~rj(n m ) r / N )
LIT, if n
m mod
r=O
0, otherwise
(6)
establishes that the right side
of (5) is in fact equal
t o Xk.
I t is useful to extend the range of definition of A, to
all integers(positiveandnegative).Withinthis
definition it follows that
A,
A.V+~
A,hr+,
7)
x1
XW+l
N,hT+l
(8)
Similarly,
Relationshipsbetweenthe DFT andtheFourierTransf o r m of a Continuous Waveform
An importantpropertythatmakestheDFT
so
eminently useful is the relationship between the
DFT
of a sequence of Nyquist samples and the Fourier transform of a continuous waveform, t h a t is represented by
theNyquist samples. T o recognize thisrelationship,
consider a frequency band-limited waveform
whose
Nyquist samples, X k , vanish outside the time interval
05tlNT
N-
A,
(Xk)W,k
0,
N
1
(2)
k=O
where
W
exp ( - 2 ~ r j l N ) .
(3)
Since the X k ’ s are often values of a function a t discrete
time points, the index
r issometimescalled
the “frequency” of the DFT. T h e D F T h a salso been called the
“discreteFouriertransform”orthe“discretetime,
finite range Fourier transform.”
There exists the usual inverse of the DFT and, because the form is very similar t o t h a t of the DFT, the
FFT may be used to compute it.
T h e inverse of ( 2 ) is
A,W4
(l/N)
Let the Fourier transform
of g(t) be Gcf). As is well
known [6], thistransformisexactly
specified at discretefrequenciesbythecomplexFourierseries
coefficients of gp(t). From this it follows:
NT
0,
NT
gp(t) .exp ( - 2 ~ r j n t / N T )
(l/NT)
N-1
XI
where T is thetimespacingbetweenthe
samples. A
periodic repetition of
can be constructed that has
identically the same Nyquist samples in the time interval 0
NT
0,
N
1.
(4)
r=O
N-
(l/NT)
Xk.exp ( - 2 r j n k T I N( T1 )1 )
This relationship is called the inverse discrete Fourier
k=O
transform (IDFT). Itis easy to show that this inversion
where nl < N / 2 due to the spectral bandwidth limitais valid by inserting (2) into (4)
tion implicitly assumed by the sampling theorem underN - 1 AT-1
lying the validityof Nyquist samples.
x1
(Xk/*W)W‘(k-l).
(5)
Comparing ( 1 1 ) and ( 1 ) i t is seen that they are exr=o k=O
actly the same except for a factor of N T and ( Y , n) are
Interchangingin (5) theorder of summingoverthe
both unbounded. T h a t is,
indices r and k, and using the orthogonality relation
N . A, D, for r n and T 1 second.
(12)
The definition of the DFTis not uniform in the literature.Some
authors use A , / N
the DFT coefficients, others use A,/dT, still
others use a positive exponent.
Authorized licensd use limted to: IE Xplore. Downlade on May 10,2 at 19:056 UTC from IE Xplore. Restricon aply.
The bounds specified for r and n require a correspondence which depends on (7)
G-AE SUBCOMMITTEE: T H E TRANSFORM
FAST FOURIER
47
ing between even and odd values
of N , a distinction
avoided by keeping r positive.
A waveform of the type considered by (9) is shown
where
in Fig. l ( e ) . I t is usually obtained as an approximation
of a frequency band-limited source waveform
[such as
n
r
for n
q N/2)
the one sketched in Fig.1(a)] by truncating the Nyquist
and
sample series of this waveform, and reconstructing the
continuouswaveformcorrespondingtothetruncated
n
N
Y for n
1, 2,
q
N / 2 (13) Nyquistsampleseries[Fig.1
(b), (d), and(e)].Notwithstanding
the
identity
of
the
Nyquist samplesof this
and
reconstructed
waveform
and
the
frequency
bandG(n/NT)
limited source waveform, these waveforms differ in the
D, 1V.A,/2
for n
N/2.
truncation interval [Fig. l(c) and (e)]. The difference
NT
is usually referred to as aliasing distortion; the mechanics of this distortion is most apparent in the frequency
Equations (13) and (14) give a direct relationship bedomain [Fig. l(c)-(e)]. I t can be made negligibly small
tween t h e D F T coefficients and the Fourier transform
by choosing a sufficiently largeproduct
at discrete frequencies for the waveform stipulated by
of the fre(9). A one-to-one correspondence could have been obquencybandwidth of the sourcewaveform
andthe
tained if the running variable r had been bounded by duration of thetruncationinterval
[6] (eg., N is
N / 2 . This, however, would have required distinguish- greater than ten).
I
*f
(Frequency)
T
(a) Frequency-band-limited source waveform.
(b)Nyquist samples of the frequencyband-limited
source waveform.
(c) Truncated source waveform.
(d) Truncated series of Nyquist samples of the source
waveform.
(e) Frequency-band-limited waveform whose Nyquist
samples are identical tothetruncated
series of
Nyquist samples of the source waveform.
(f) Periodic continuation of thetruncated
source
waveform.
(g) Periodic continuation of thetruncated series of
Nyquist samples of the source waveform.
(h) DFT coefficients interpreted as Fourier series coefficients producing complex waveform.
Fig. 1. Related waveforms and their corresponding spectra asdefined by the Fourier transforms (integral transforms for energy-limited waveforms; series transform
for periodic waveforms).
Authorized licensd use limted to: IE Xplore. Downlade on May 10,2 at 19:056 UTC from IE Xplore. Restricon aply.
IEEE TRANSACTIONS ON AUDIO AND ELECTROACOUSTICS,
JGNE
48
These aliasing distortions are carried over directly to
the discrete spectra of the periodically repeated waveforms [Fig. l(f) and (g)], and appear correspondingly
in the DFT of the truncated series of Nyquist samples
[Fig. l(h)]. It may be of interest to observe that the
waveform corresponding to the D F T coefficients interpretedasFourierseries
coefficients is complex[Fig.
1(h) 1.
1967
D F T of a time series more rapidly than do other algorithms available. Thepossibility of computing the DFT
bysuch a fastalgorithmmakestheDFTtechnique
important. A comparison of the cornputational savings
that may beachieved through use of the FFT is summarized in Table I for various computations that are frecomquently performed. I t is important to add that the
putationaleffortslistedrepresentcomparableupper
bounds; the actual efforts dependon the number N and
Some Useful Properties of the D F T
the programming ingenuity applied[7].
I t m a y b euseful to point out that the FFT not only
Anotherpropertythatmakesthe
DFT eminently
reduces
thecomputationtime;
i t alsosubstantially
useful is the convolution relationship.T h a t is, t h e I D F T
reduces
round-off
errors
associated
with these computaof the product of two DFTsis the periodic mean convotions.
In
fact,
both
computation
time
and
round-off
lution of the two time seriesof the DFTs. This relationerror
essentially
are
reduced
by
a
factor
of
(logz
N)/N
ship proves very useful when computing the filter outwhere
N
is
the
number
of
data
samples
in
the
time
p u ta s a result of aninputwaveform;it
becomes
especiallyeffective
whencomputedbythe
FFT. A series. Forexample, if N = 1024 21°, then N . log,
N = 10 240 [7], [9]. Conventional methods for computderivation of this property is given in Appendix
A.
ing (1) for N = 1024 would require an effort proportional
Other properties of the DFT are in agreement with
t o N 2 1 048 576, more than 50 times that required with
thecorrespondingproperties
of theFourierintegral
transform,perhapswithslight
modifications. For ex- the FFT.
T h e FFT is a clever computational technique of seample, the DFTof a time series circularly shifted by
h is
quentially
combining
progressively larger
weighted
t h e D F Tof the time series multiplied byWmrh.
Furthersums
of
data
samples
so
as
to
produce
the
DFT coeffimore, the DFTof the sum of two functions is the sum
of
cients
as
defined
b
y
(2).
The
technique
can
be intert h e D F T of thetwofunctions.Thesepropertiesare
preted
in
terms
of
combining
the
DFTs
of
the
individual
readily derived using the definition of the DFT. These
data
samples
such
that
the
occurrence
times
of these
and other properties have beencompiled by Gentleman
samples
are
taken
into
account
sequentially
and applied
and Sande [7].
to the DFTs of progressively larger mutually exclusive
subgroups of data samples, which are combined to ultiTHEFASTFOURIER
TRANSFORM
mately produce the D F T of the complete series of data
General Description of the F F T
samples. T h e explanation of the FFT algorithm adopted
As mentionedin theIntroduction,the
FFT is an in this paperis believed to be particularly descriptive
for
algorithm that makes possible the computation of the programming purposes.
TABLE I
COMPARISON OF T H E x U M B E R O F hIULTIPLICATIONS REQUIRED
USING “DIREcr’’ AND
Operation
Discrete Fourier Transform (DFT)
Filtering (Convolution)
Autocorrelation Functions
Two-Dimensional FourierTransform(Pattern Analysis:
Two-Dimensional Filtering
Authorized licensd use limted to: IE Xplore. Downlade on May 10,2 at 19:056 UTC from IE Xplore. Restricon aply.
FFT
METHODS
E R F A S TG-AE S U
TB
HC
EO M M I T T E E :
Conventional Forms of the F F T
which, using (16), may be written in the
following form:
Decimation in T i m e : T h e D F T [asper (2)] and its
inverse [see (4)] are of the same form so that a procedure, machine, or sub-routine capableof computing one
can be used for computing the other by simply exchanging the roles of X k and A,, andmakingappropriate
scale-factor and sign changes. The two basic forms of
the FFT, each with its several modifications, are thereforeequivalent.However,
i t isworthdistinguishing
between them and discussing them separately. Let
us
first consider the form used by Cooley and Tukey
[l]
which shall be called decimation in time. Reversing the
roles of A, and Xkgives the form called decimation in
frequency, which will be considered afterwards.
Suppose a time series having N samples [such as Xk
shown in Fig. 2 ( a ) ] is divided into two functions,Y k and
zk,each of which has only half as many points ( N / 2 ) .
T h e function Y k iscomposed
of the even-numbered
points ( X o ,X 2 , X g
and
is composed of the odd
numbered points (XI,
X p ,X5
These functions are
showninFig.
2 (b) and (c), and we maywritethem
formally as
A,.
yk
X2k
K=0,1,2,...
zk
(15)
1.
' 2
X2k+l
Since Y k and
are sequences of N / 2 points each, they
have discrete Fourier transformsdefined by
h=O
r=0,1,2;..
(N/2)--1
C,
N
(16)
' 2
zkexp (-4rjrK/;V)
T h e discreteFouriertransform
t h a t we want is A,,
whichwecanwrite
in terms of theodd-andevennumbered points
25
[2K
N
exp
11))
1
(17)
Zk exp ( - 4 n j r k / N )
(18)
2,
Y
N
or
(N/2)--1
A,
Y k
exp
4.-jrk/Ar)
k=O
(Ni2)-1
exp
2~jr/N)
k=O
Authorized licensd use limted to: IE Xplore. Downlade on May 10,2 at 19:056 UTC from IE Xplore. Restricon aply.
B,
exp
27rjr/N)C,
O l r < N / 2(. 1 9 )
For values of r greater than N / 2 , the discrete Fourier
transforms B , and C, repeatperiodicallythevalues
taken on when r N / 2 . Therefore, substituting r + N / 2
for r in (19), we obtain
Ar+N/2
B,
exp (-2=j
B,
r
$]/N)
C,
5Y
05Y
exp (-2n-jr/AT)C,
N/2
N / 2 .( 2 0 )
By using (3), ( 1 9 ) and ( 2 0 ) may be written as
A,
B,
4- W'C,
A,+N/2
B,
W'C,
5Y
05 r
(21)
0
hT/2.
(22)
From ( 2 1 ) and ( 2 2 ) , the first N / 2 and last N / 2 points of
the discrete Fourier transformof Xk (a sequence having
N samples) can be simply obtained from the DFT of Y k
and Z k , both sequences of N / 2 samples.
Assuming t h a t we have a methodwhichcomputes
discrete Fourier transforms ina time proportional to the
square of the number of samples, we can use this algorithm to compute the transformsof Yk and Zk,requiring
a time proportional to 2 ( N / 2 ) 2 , and use ( 2 1 ) and ( 2 2 ) to
find A , with additional N operations. This is illustrated
in the signal flow graph of Fig. 3. The points on the left
are the values of Xk (Le., Yk and Zk),and the points on
the right are the points of the discrete Fourier transform, A,. For simplicity, Fig. 3 is drawn for the case
where Xk is an eight-point function, and advantage is
taken of the fact that Wn
Wn-N'2,as per (3).
However, since Yk and zk are to be transformed, and
since we have shown that the computation of t h e D F T
of N samples can be reduced to computing the DFTsof
two sequences of N / 2 samples each, the computationof
B k [or C k > can be reduced to the computation
of sequences of N/4 samples. These reductions can be carried out as long as each function has a number of samples that is divisible by 2. Thus, if N = 2" we can make
such reductions, applying (IS), ( 2 1 ) , and ( 2 2 ) first for
N , then for N / 2 ,
and finally for a two-point function. The discreteFouriertransform
of a one-point
function is, of course, the sample itself. The successive
reduction of an eight-point discrete Fourier transform,
begun in Fig. 3, is continued in Figs. 4 and 5. In Fig. 5
the operation has been completely reduced to complex
multiplicationsandadditions.Fromthesignal
flow
graph there are 8 by 3 terminal nodes and 2 by 8 by 3
arrows, correspondingto 24 additions and48 multiplications. Half of the multiplications can be omitted since
the transmission indicated by the arrow is unity. Half
of
the remaining multiplications are also easily eliminated,
we shall see below. Thus, in general, N-logz N complex additions and, at most, 4N.Iog2 N complex multi-
IEEE TRANSACTIONS ON AUDIO AND ELECTROACOUSTICS, JUNE 1967
50
k
(b)
o
Fig. 2. Decomposition of a time series into two part-time
series, each of which consists of half the samples.
Fig. 4. Signal flow graph illustrating further reduction of
the DFT computation suggested by Fig. 3.
c3
Fig. 3. Signal flow graphillustrating
the reduction of endpoint
D F T to two DFTs of N / 2 points each, using decimation in time.
The signal flow graph may be unfamiliar to some readers. Basically it is composed of dots (or nodes) and arrows(transmissions). Each node represents a variable, and the arrows terminating at that node originate at the nodes whose variables contribute
to the value of the variable at that node. The contributions are
additive, and theweight of each contribution, if other than unity,
is indicated by the constant written close to the arrowhead of
the transmission. Thus, in this example, the quantity A7 a t t h e
bottom right node is equal t o B3fW7XCa. Operations other than Fig. 5. Signal flow graph illustrating the computation of the DFT
when the 0perations:involved are completely reduced to multipliaddition and constantmultiplication must be clearly indicated by
cations and additions
symbols other than or
Authorized licensd use limted to: IE Xplore. Downlade on May 10,2 at 19:056 UTC from IE Xplore. Restricon aply.
ER
51
FASTG-AE SUBCOMMITTEE:
THE
plications are required for computation of the discrete
Fourier transform of an N point sequence, where N is a
power of 2.
When N is not a power of 2 , but has a factor p , the
development of equationsanalogousto
(15) through
( 2 2 ) is possible byforming
p differentsequences,
Yk(i)
=Xpk+;,each having N / p samples. Each of these
sequences has aD F T Br(i),and the DFTof the sequence
Xkcan be computed from thep simpler DFTs with p N
complex multiplications and additions. That
is,
p-1
AT+rn(NlP)
&(i)WiI‘+”(”/P)I
m=0,1,2,...,p-l
9.=0,1,2,...
N
’ P
(23)
The computation of the DFTs can be further simplified
if N has additional prime factors.
Further information about thefast Fourier transform
can be extracted from Fig.
For example, if the input
sequence Xk is stored in computer memory in the order
xo, x4,
x 2 , X6,
XI, x5,x 3 X, ? ,
(24)
For this version of the algorithm, the initial shuffling
of the data sequence, X,, wasnecessaryfor
the“in
place”computation.This
shufflingis duetothe
repeated movement of odd-numbered members of a sequence to the end of the sequence during each stage of
the reduction, as shown in Figs. 3, 4, and 5. This shuffling has been called bit reversal2 because the samples are
stored in bit-reversed order; i.e.,
X4=X(100)2is stored
inposition (011)z = 3 , etc. Note that the initial data
shuffling can also be done (‘in place.”
Variations of Decimation in T i m e : If one so desires,
the signal flow graph shown in Fig.
5 can be manipulated to yield different forms of the decimation in time
version of thealgorithm.
If oneimagines
that in
Fig. 5 all the nodeson thesamehorizontal
level as
A1 areinterchangedwith
all the nodeson thesame
horizontal level as A 4 , and all the nodes on the level of
As are interchanged with the nodes on the level of Ag,
withthearrows carriedalong withthenodes, then one
obtains aflow graph like that
of Fig. 6.
For this rearrangement one need not shuffle the original data into the bit-reversed order, but the resulting
spectrum needs to be unshufled. An additional disadvantage might be that the powers of W needed in the
computation are in bit-reversed order. Cooley’s original
description of the algorithm [ l ] corresponds to the flow
graph of Fig. 6 .
A somewhat more complicated rearrangementof Fig.
5and
yields the signal flow graph of Fig. 7. For this case
both the input data and the resulting spectrum are in
“natural” order, and thecoefficients in the computation
are also used in a natural order. However, the computation may no longer be done “in place.” Therefore,
at
leastoneotherarray
of registersmustbeprovided.
This signal flow graph, and a procedure corresponding
to it, are due to Stockham[8].
Decimation in Frequency: Let us nowconsider
a
second, quite distinct,
form of the fast Fourier transform algorithm, decimation in frequency. This form was
foundindependentlybySande
and Cooley and
Stockham [SI. Let the time series X k have a DFT A,.
The series andthe DFT bothcontain N terms. As
before, we divide X k into two sequenceshaving N / 2
pointseach.However,thefirstsequence,
Yk,is now
composed of the first N / 2 points in X L ,and the second,
Z k , is composed of the last N / 2 points in Xk. Formally,
then
as inFig. 5, ihe computation of the discrete Fourier
transform may be done “in place,”that is, by writing all
intermediate results over the original data sequence,
writing the final answer over the intermediate results.
Thus, no storage is needed beyond t h a t required for the
original N complex numbers. T o see this, suppose that
eachnodecorrespondstotwomemoryregisters(the
quantities to be stored are complex). The eight nodes
farthest to the leftin Fig. 5 then represent the registers
containing the shuffled order input data. The first step
in the computation is to compute the contents of the
registers represented by the eight nodesjust to the right
of the input nodes. But each pair of input nodes affects
only the corresponding pairof nodes immediately to the
right, and if the computation deals with two nodes at a
time, the newly computed quantities may
be written
intotheregisters
fromwhich
theinputvalues
were
taken, since the input values are no longer needed for
further computation. The second step, computation of
the quantities associated with the next vertical arrayof
nodes to the right, also involves pairs of nodes although
these pairs are now two locations apart instead of one.
Yr, X k
This fact does not change the property of “in place”
computation, since each pair of nodes affects only the
pair of nodes immediately to the right. After a new pair
of results is computed, it may be storedin the registers
which held the old results that are nolonger needed. In
thecomputationforthefinalarray
ofnodes,corresponding to the values
of the DFT, the computation
This is a special case of digit reversal where the radix of the
involves pairs of nodes separated by four locations,but
address is 2; more general digit reversalsare available for transforms
the “in place” property still holds.
with other radices.
Authorized licensd use limted to: IE Xplore. Downlade on May 10,2 at 19:056 UTC from IE Xplore. Restricon aply.
IEEE TRANSSCTIONS
ON
52
Fig. 6 . Rearrangement of the flow graph of Fig. 5 illustrating the
D F T computation from naturally ordered
time
samples.
Fig. 7.
AUDIO AND ELECTROACOUSTICS, JUNE 1967
Rearrangement of the flow graph of Fig. 5 illustrating the
DFT computation without
bit
reversal.
Fig. 8. Signal flow graph illustrating the reduction of endpoint
DFT to two DFTs of
points each, using decimation in frequency.
Authorized licensd use limted to: IE Xplore. Downlade on May 10,2 at 19:056 UTC from IE Xplore. Restricon aply.
53
G-AE S U B C O M M I T T E E : T H E F A S T F O U R I E R T R A N S F O R M
TRANSFORM
AI
-w'=
Fig. 9.
Signal flow graph illustrating further reduction of
the DFT computationsuggested by Fig. 8.
Fig. 11. Rearrangement of the flow graph of Fig. 10 illustrating the
computation of the DFT to yield naturally ordered D F T coefficients.
Authorized licensd use limted to: IE Xplore. Downlade on May 10,2 at 19:056 UTC from IE Xplore. Restricon aply.
Fig. 10. Signal flow graph illustrating the computation of the DFT
when the operations involved are completely reduced to multiplications and additions.
Fig. 12.
Rearrangement of the flow graph of Fig. 10 illustrating
the DFT computation without bitreversal.
IEEE TRANSACTIONS
ELECTROACOUSTICS,
ON
AND
AUDIO
54
J U N E 1967
Examination of Fig.10gives
us muchinformation
about the method of decimation in frequency, and allows us to compare it with decimation in time. Both
methods requireN / 2 .log N complex additions, complex
subtractions, and complex multiplications. Both computations can be done in place. If the coefficients in the
computation are to be used in a “natural” rather than
“bit-reversed” order, as in Figs. 5 and 10, then the decimation in frequency method works on time samples in
(N/2)--1
unshuffled order and yields frequency samples in shufA,
Y k [exp (-rjr)]Zk] exp (-2?rjrL/N). (27)
fled (bit-reversed) order. Recall that Fig. 5 yielded the
k=O
opposite result.
Let us consider separately the even-numbered and oddWe are also able to rearrange the nodes in Fig. 10 to
numbered points of the transform. Let the even-numobtain the signal flow graph, Fig. 11, which works on
bered points be X, and the odd-numbered points be S,, shuffled time samples and yields naturally ordered frewhere
quency samples, but the coefficients are needed by the
computation in bit-reversedorder.
T h e geometry of
Rr A Zr
this signal flow graph is identical to the geometry of
05 r
N/2.
(28) Fig. 5, just as the geometry of Fig. 10 is identical to the
geometry of Fig. 6. T h e differences lie in the transmisSr
Asrtl
sions.
I t is this step that may called
be decifnation in frequency.
A somewhatmorecomplicatedrearrangement
of
Note that for computing the even-numbered spectrum
Fig. 10 (shown in Fig. 12) yields a signal flow graph that
points, (27) becomes simply
takes unshuffled samples of thetimeseriesandproduces a set of Fourier coefficients that are not in bit(N/2)-1
reversedorder. The computation cannot, however, be
R,
z k } e ( - 2 ~ j r k ) / ( ~ / 2 ) (29)
done “in place,” and
a t least one other arrayof registers
k=O
must be provided. The method is similar to that shown
which we recognize as the N / 2 point D F T of the func- in Fig. 7 for decimation in time. T h e forms of Figs. 5, 6,
tion Yk+Zk),
the sum of the first N / 2 and the lastN / 2 7, 10, 11, and 1 2 constitute a set of what we might call
timesamples.Similarly,fortheodd-numbered
spec- canonic forms of the fast Fourier transform. We may
trum points, (27) becomes
choose among these forms to find an algorithm with the
properties of “in place” computation, normally ordered
(W/Z)-l
input, normally ordered output,. or normally ordered
Sr
A2r+1
vk
exp ( - ~ j [ 2 ~ 11)
coefficients, but not all four at once. T o achieve“in
place” computation, we must deal with bit reversal, and
.exp ( - 2 r j [ 2 ~ I]R/N)
to eliminate bit reversal
we must give up “in place”
(N/Z)-l
computation.Thetwomethodsmost
effective when
ylc ~ , ] ~ ( - 2 r j i i ) : ~ ~ [ ( ~ ~ r r j r k ) / ( i \ ‘ /(30)
2)1
using homogeneous storage facilities are those providing
k=O
in right order the sine and cosine coefficients needed in
which we recognize as the N / 2 point D F T of the func- the computation. The other methods seem
less desirtion Y k Z,)exp 27rjjk/N).
ablesincetheyrequirewastefultables.Still,
all six
I t can be concluded from (29) and (30) t h a t t h e D F T methods have about equal usefulness, and the method
of an N-samplesequence, X k , may be determined as usedbest will dependontheproblem
a t hand.For
follows. For even-numbered transform points,it maybe example,themethodshowninFig.
10 maybeused
computed as anN / 2 point D F T of a simple combination to transform from the time to the frequency domain,
of the first N / 2 and last N / 2 samples of X k . For odd- and the method shown in Fig. 4 may be used for the
numberedtransformpoints,itmaybecomputed
as inverse transform. Any of the methods described above
another N / 2 point D F T of a different simple combina- may be used for the inverse discrete Fourier transform
if
tion of the first and lastN / 2 samples of X k . This is illus- the coefficients are replaced by their complexconjutrated in the signal flow graph of Fig. 8 for an eight- gates, and if the result of the computation is multiplied
point function. W has been defined in (3).
by 1/N.
As was the case with decimation in time, we can reT h e six forms mentioned are, in a sense, canonic, b u t
placeeach of theDFTsindicatedinFig.
8 bytwo
one could also employ a combination of decimation in
2-point DFTs, and each of the 2-point DFTs by two
time and decimation in frequency at different stages in
1-point transforms, these last being equivalency opera- thereductionprocess,yielding
a hybridsignal
flow
tions. These steps are indicatedin Figs. 9 and 10.
graph.
zk
Authorized licensd use limted to: IE Xplore. Downlade on May 10,2 at 19:056 UTC from IE Xplore. Restricon aply.
G-BE SUBCOMMITTEE: THE FAST FOURIER TRANSFORM
A Useful Computational Variation: I t may be worth
pointing out here how some programming simplicity is
realized when the factors p and p = N / p are relatively
prime. As described by Cooley, Lewis, and Welch, [2],
the “twiddle factor” Wi’ of (23) can be eliminated by
choosingsubsequences of the X k ’ s thataredifferent
thanthose usedbefore. T h e DFT computationsare
then conveniently performed in two stages.
1) Compute the q-point transforms
of each of the
55
is alsofrequency-band-limitedto
l/2T Hz and completely specified byitsNyquistsamples
spaced T
second apart
The convolution relationship facilitates computation of
this equation.
T o prove the convolution relationship, let the
D F T of
the XL’S be A , andcorrespondingly the D F T of the
Yh’s be B,. T h e I D F T of the product of A,. B , then
becomes [see
]
p sequences
2) Compute, then, the p-point transforms
(33)
of the
where
sequences B,(%),
(+)-Zs
perturbation term.
(36)
is meant to represent the rei.e., thesolution of ~ ( p ) ~ - ~ >If1 the first N / 2 samples of each of the two time series
(X,)and ( Y h ) are assumed to be identically zero, then
the
perturbation term of (36) is zero so t h a t t h e I D F of
T
CONCLUSION
the productof the two DFTs multiplied by N is equal to
T h e integral transform method has been one
of the the convolution product
of (35). Since it is always
foundations of analysis for many years because of the possible to select the time series to be convolved such
easewithwhichthetransformedexpressionsmaybe
that half of the samples are zero, the convolution relamanipulated, particularly in such diverse areasa s acous- tionship for the DFT can be used to compute the contic wave propagation. speech transmission, linear netvolution product [see (35)] of two time series.
work theory, transport phenomena, optics, and electroI t is useful to point out that
if A,=B,, aperiodic
magnetic theory. Many problems which are particularly autocorrelation function emerges.
amenabletosolutionbyintegraltransformmethods
REFERENCES
have not been attacked by this method in the past because of the high cost of obtaining numerical results this [ l ] J. W. Cpoley and J. W. Tukey, ILAnalgorithm for the machine
calculation of complex Fourier series,” Math. of Comput., vol. 19,
way.
pp. 297-301, April 1965.
The fast Fourier transform has certainly
modified the [2] J. W. Cooley, P. A. Mr. Lewis, and P. D. Welch, “Historical notes
the fast Fourier transform” thisissue, p. 76-79.
economics of solution by ‘transform methods. Some new [3] on
R.B. Blackman and J. W. Tukey, TheMeasurement of Power
applicationsarepresented
in thisspecialissue,and
Spectra. New York: Dover, 1959.
C. Bingham, M. D. Godfrey, and J. W. Tukey, “Modern techfurther interesting and profitable applications probably [4] niques
of power spectrum estimation,” this issue, p. 56-66.
will be found during the next few years.
[5] T. G. Stockham, “High speed convolution and correlation,” 1966
and the notation (p)pl
ciprocal of p, mod
(mod
SpringJoint ComputerConf., A F I P S Proc., vol. 28. Washington, D. C.: Spartan, 1966, pp. 229-233.
[6] W. T. Cochran, J. J. Downing, D. L. Davin, H. D. Helms, R. A.
Kaenel, W. W. Lang, and D. E. Nelson, “Burst measurements in
As is well known, if the filter impulse response is frethe frequency domain,” Proc. IEEE, vol. 54, pp. 830-841, June
quency-band-limitedto1/2THzand
is given by its
1966.
Nyquist samples Yh spaced T second apart, and further- [7] W. M. Gentleman and G. Sande, “Fast Fourier transforms-for
fun and profit,” 1966Fall Joint ComputerConf. A F I P S Proc.,
more, if the input waveform is alsofrequency-bandvol. 29. Washington, D. C.: Spartan, 1966, pp. 563-578.
Private communication.
limited t o 1/2T Hz andgiven by its Nyquist samples
x k
[9] J. W. Cooley, “Applications of the fast Fourier
transform method,’’
spaced T second apart, then the filter output waveform
Proc. of the I B M Scientific Computing Symp., June 1966.
APPENDIX
Authorized licensd use limted to: IE Xplore. Downlade on May 10,2 at 19:056 UTC from IE Xplore. Restricon aply.