The Cooley–Tukey algorithm, named after
J. W. Cooley and
John Tukey, is the most common
fast Fourier transform (FFT) algorithm. It re-expresses the
discrete Fourier transform
In mathematics, the discrete Fourier transform (DFT) converts a finite sequence of equally-spaced samples of a function into a same-length sequence of equally-spaced samples of the discrete-time Fourier transform (DTFT), which is a comple ...
(DFT) of an arbitrary
composite
Composite or compositing may refer to:
Materials
* Composite material, a material that is made from several different substances
** Metal matrix composite, composed of metal and other parts
** Cermet, a composite of ceramic and metallic materials
...
size
in terms of ''N''
1 smaller DFTs of sizes ''N''
2,
recursively
Recursion (adjective: ''recursive'') occurs when a thing is defined in terms of itself or of its type. Recursion is used in a variety of disciplines ranging from linguistics to logic. The most common application of recursion is in mathematics ...
, to reduce the computation time to O(''N'' log ''N'') for highly composite ''N'' (
smooth number
In number theory, an ''n''-smooth (or ''n''-friable) number is an integer whose prime factors are all less than or equal to ''n''. For example, a 7-smooth number is a number whose every prime factor is at most 7, so 49 = 72 and 15750 = 2 × 32 ...
s). Because of the algorithm's importance, specific variants and implementation styles have become known by their own names, as described below.
Because the Cooley–Tukey algorithm breaks the DFT into smaller DFTs, it can be combined arbitrarily with any other algorithm for the DFT. For example,
Rader's or
Bluestein's algorithm can be used to handle large prime factors that cannot be decomposed by Cooley–Tukey, or the
prime-factor algorithm can be exploited for greater efficiency in separating out
relatively prime
In mathematics, two integers and are coprime, relatively prime or mutually prime if the only positive integer that is a divisor of both of them is 1. Consequently, any prime number that divides does not divide , and vice versa. This is equivale ...
factors.
The algorithm, along with its recursive application, was invented by
Carl Friedrich Gauss
Johann Carl Friedrich Gauss (; german: Gauß ; la, Carolus Fridericus Gauss; 30 April 177723 February 1855) was a German mathematician and physicist who made significant contributions to many fields in mathematics and science. Sometimes refer ...
. Cooley and Tukey independently rediscovered and popularized it 160 years later.
History
This algorithm, including its recursive application, was invented around 1805 by
Carl Friedrich Gauss
Johann Carl Friedrich Gauss (; german: Gauß ; la, Carolus Fridericus Gauss; 30 April 177723 February 1855) was a German mathematician and physicist who made significant contributions to many fields in mathematics and science. Sometimes refer ...
, who used it to interpolate the trajectories of the
asteroids
Pallas
Pallas may refer to:
Astronomy
* 2 Pallas asteroid
** Pallas family, a group of asteroids that includes 2 Pallas
* Pallas (crater), a crater on Earth's moon
Mythology
* Pallas (Giant), a son of Uranus and Gaia, killed and flayed by Athena
* Pa ...
and
Juno
Juno commonly refers to:
*Juno (mythology), the Roman goddess of marriage and queen of the gods
*Juno (film), ''Juno'' (film), 2007
Juno may also refer to:
Arts, entertainment and media Fictional characters
*Juno, in the film ''Jenny, Juno''
*Ju ...
, but his work was not widely recognized (being published only posthumously and in
neo-Latin
New Latin (also called Neo-Latin or Modern Latin) is the revival of Literary Latin used in original, scholarly, and scientific works since about 1500. Modern scholarly and technical nomenclature, such as in zoological and botanical taxonomy ...
).
[Heideman, M. T., D. H. Johnson, and C. S. Burrus,]
Gauss and the history of the fast Fourier transform
" IEEE ASSP Magazine, 1, (4), 14–21 (1984) Gauss did not analyze the asymptotic computational time, however. Various limited forms were also rediscovered several times throughout the 19th and early 20th centuries.
[ FFTs became popular after James Cooley of IBM and John Tukey of ]Princeton
Princeton University is a private research university in Princeton, New Jersey. Founded in 1746 in Elizabeth as the College of New Jersey, Princeton is the fourth-oldest institution of higher education in the United States and one of the ni ...
published a paper in 1965 reinventing the algorithm and describing how to perform it conveniently on a computer.
Tukey reportedly came up with the idea during a meeting of President Kennedy's Science Advisory Committee discussing ways to detect nuclear-weapon tests in the Soviet Union
The Soviet Union,. officially the Union of Soviet Socialist Republics. (USSR),. was a List of former transcontinental countries#Since 1700, transcontinental country that spanned much of Eurasia from 1922 to 1991. A flagship communist state, ...
by employing seismometers located outside the country. These sensors would generate seismological time series. However, analysis of this data would require fast algorithms for computing DFTs due to the number of sensors and length of time. This task was critical for the ratification of the proposed nuclear test ban so that any violations could be detected without need to visit Soviet facilities. Another participant at that meeting, Richard Garwin
Richard Lawrence Garwin (born April 19, 1928) is an American physicist, best known as the author of the first hydrogen bomb design.
In 1978, Garwin was elected a member of the National Academy of Engineering for contributing to the application ...
of IBM, recognized the potential of the method and put Tukey in touch with Cooley. However, Garwin made sure that Cooley did not know the original purpose. Instead, Cooley was told that this was needed to determine periodicities of the spin orientations in a 3-D crystal of helium-3. Cooley and Tukey subsequently published their joint paper, and wide adoption quickly followed due to the simultaneous development of Analog-to-digital converter
In electronics, an analog-to-digital converter (ADC, A/D, or A-to-D) is a system that converts an analog signal, such as a sound picked up by a microphone or light entering a digital camera, into a digital signal. An ADC may also provide ...
s capable of sampling at rates up to 300 kHz.
The fact that Gauss had described the same algorithm (albeit without analyzing its asymptotic cost) was not realized until several years after Cooley and Tukey's 1965 paper.[ Their paper cited as inspiration only the work by I. J. Good on what is now called the ]prime-factor FFT algorithm
The prime-factor algorithm (PFA), also called the Good–Thomas algorithm (1958/1963), is a fast Fourier transform (FFT) algorithm that re-expresses the discrete Fourier transform (DFT) of a size ''N'' = ''N''1''N''2 as a two-dimensional ''N''1× ...
(PFA);[ although Good's algorithm was initially thought to be equivalent to the Cooley–Tukey algorithm, it was quickly realized that PFA is a quite different algorithm (working only for sizes that have ]relatively prime
In mathematics, two integers and are coprime, relatively prime or mutually prime if the only positive integer that is a divisor of both of them is 1. Consequently, any prime number that divides does not divide , and vice versa. This is equivale ...
factors and relying on the Chinese Remainder Theorem, unlike the support for any composite size in Cooley–Tukey).
The radix-2 DIT case
A radix-2 decimation-in-time (DIT) FFT is the simplest and most common form of the Cooley–Tukey algorithm, although highly optimized Cooley–Tukey implementations typically use other forms of the algorithm as described below. Radix-2 DIT divides a DFT of size ''N'' into two interleaved DFTs (hence the name "radix-2") of size ''N''/2 with each recursive stage.
The discrete Fourier transform (DFT) is defined by the formula:
:
where is an integer ranging from 0 to .
Radix-2 DIT first computes the DFTs of the even-indexed inputs
and of the odd-indexed inputs , and then combines those two results to produce the DFT of the whole sequence. This idea can then be performed recursively
Recursion (adjective: ''recursive'') occurs when a thing is defined in terms of itself or of its type. Recursion is used in a variety of disciplines ranging from linguistics to logic. The most common application of recursion is in mathematics ...
to reduce the overall runtime to O(''N'' log ''N''). This simplified form assumes that ''N'' is a power of two
A power of two is a number of the form where is an integer, that is, the result of exponentiation with number two as the base and integer as the exponent.
In a context where only integers are considered, is restricted to non-negativ ...
; since the number of sample points ''N'' can usually be chosen freely by the application (e.g. by changing the sample rate or window, zero-padding, etcetera), this is often not an important restriction.
The radix-2 DIT algorithm rearranges the DFT of the function into two parts: a sum over the even-numbered indices and a sum over the odd-numbered indices :
:
One can factor a common multiplier out of the second sum, as shown in the equation below. It is then clear that the two sums are the DFT of the even-indexed part and the DFT of odd-indexed part of the function . Denote the DFT of the ''E''ven-indexed inputs by and the DFT of the ''O''dd-indexed inputs by and we obtain:
:
Note that the equalities hold for
but the crux is that and are calculated in this way for only.
Thanks to the periodicity of the complex exponential, is also obtained from and :
:
We can rewrite and as:
:
This result, expressing the DFT of length ''N'' recursively in terms of two DFTs of size ''N''/2, is the core of the radix-2 DIT fast Fourier transform. The algorithm gains its speed by re-using the results of intermediate computations to compute multiple DFT outputs. Note that final outputs are obtained by a +/− combination of and , which is simply a size-2 DFT (sometimes called a butterfly
Butterflies are insects in the macrolepidopteran clade Rhopalocera from the order Lepidoptera, which also includes moths. Adult butterflies have large, often brightly coloured wings, and conspicuous, fluttering flight. The group comprise ...
in this context); when this is generalized to larger radices below, the size-2 DFT is replaced by a larger DFT (which itself can be evaluated with an FFT).
This process is an example of the general technique of divide and conquer algorithm
In computer science, divide and conquer is an algorithm design paradigm. A divide-and-conquer algorithm recursively breaks down a problem into two or more sub-problems of the same or related type, until these become simple enough to be solved direc ...
s; in many conventional implementations, however, the explicit recursion is avoided, and instead one traverses the computational tree in breadth-first
Breadth-first search (BFS) is an algorithm for searching a tree data structure for a node that satisfies a given property. It starts at the tree root and explores all nodes at the present depth prior to moving on to the nodes at the next de ...
fashion.
The above re-expression of a size-''N'' DFT as two size-''N''/2 DFTs is sometimes called the Danielson
Danielson is an American rock band from Clarksboro, New Jersey, that plays indie pop gospel music. The group consists of frontman Daniel Smith and a number of various artists with whom he collaborates. Smith has also released solo work as Br ...
–Lanczos
__NOTOC__
Cornelius (Cornel) Lanczos ( hu, Lánczos Kornél, ; born as Kornél Lőwy, until 1906: ''Löwy (Lőwy) Kornél''; February 2, 1893 – June 25, 1974) was a Hungarian-American and later Hungarian-Irish mathematician and physicist. Acco ...
lemma, since the identity was noted by those two authors in 1942 (influenced by Runge's 1903 work[). They applied their lemma in a "backwards" recursive fashion, repeatedly ''doubling'' the DFT size until the transform spectrum converged (although they apparently didn't realize the ]linearithmic
In computer science, the time complexity is the computational complexity that describes the amount of computer time it takes to run an algorithm. Time complexity is commonly estimated by counting the number of elementary operations performed by ...
.e., order ''N'' log ''N''asymptotic complexity they had achieved). The Danielson–Lanczos work predated widespread availability of mechanical or electronic computers and required manual calculation (possibly with mechanical aids such as adding machine
An adding machine is a class of mechanical calculator, usually specialized for bookkeeping calculations.
In the United States, the earliest adding machines were usually built to read in dollars and cents. Adding machines were ubiquitous of ...
s); they reported a computation time of 140 minutes for a size-64 DFT operating on real inputs to 3–5 significant digits. Cooley and Tukey's 1965 paper reported a running time of 0.02 minutes for a size-2048 complex DFT on an IBM 7094 (probably in 36-bit single precision
Single-precision floating-point format (sometimes called FP32 or float32) is a computer number format, usually occupying 32 bits in computer memory; it represents a wide dynamic range of numeric values by using a floating radix point.
A floatin ...
, ~8 digits).[ Rescaling the time by the number of operations, this corresponds roughly to a speedup factor of around 800,000. (To put the time for the hand calculation in perspective, 140 minutes for size 64 corresponds to an average of at most 16 seconds per floating-point operation, around 20% of which are multiplications.)
]
Pseudocode
In pseudocode, the below procedure could be written:
''X''0,...,''N''−1 ← ditfft2(''x'', ''N'', ''s''): ''DFT of (x''0, ''xs'', ''x''2''s'', ..., ''x''(''N''-1)''s''):
if ''N'' = 1 then
''X''0 ← ''x''0 ''trivial size-1 DFT base case''
else
''X''0,...,''N''/2−1 ← ditfft2(''x'', ''N''/2, 2''s'') ''DFT of (x''0, ''x''2''s'', ''x''4''s'', ..., ''x''(''N''-2)''s'')
''XN''/2,...,''N''−1 ← ditfft2(''x''+s, ''N''/2, 2''s'') ''DFT of (xs'', ''xs''+2''s'', ''xs''+4''s'', ..., ''x''(''N''-1)''s'')
for ''k'' = 0 to ''N''/2−1 do ''combine DFTs of two halves into full DFT:''
p ← ''Xk''
q ← exp(−2π''i''/''N'' ''k'') ''Xk''+''N''/2
''Xk'' ← p + q
''Xk''+''N''/2 ← p − q
end for
end if
Here, ditfft2
(''x'',''N'',1), computes ''X''=DFT(''x'') out-of-place by a radix-2 DIT FFT, where ''N'' is an integer power of 2 and ''s''=1 is the stride
Stride or STRIDE may refer to:
Computing
* STRIDE (security), spoofing, tampering, repudiation, information disclosure, denial of service, elevation of privilege
* Stride (software), a successor to the cloud-based HipChat, a corporate cloud-based ...
of the input ''x'' array
An array is a systematic arrangement of similar objects, usually in rows and columns.
Things called an array include:
{{TOC right
Music
* In twelve-tone and serial composition, the presentation of simultaneous twelve-tone sets such that the ...
. ''x''+''s'' denotes the array starting with ''xs''.
(The results are in the correct order in ''X'' and no further bit-reversal permutation is required; the often-mentioned necessity of a separate bit-reversal stage only arises for certain in-place algorithms, as described below.)
High-performance FFT implementations make many modifications to the implementation of such an algorithm compared to this simple pseudocode. For example, one can use a larger base case than ''N''=1 to amortize the overhead of recursion, the twiddle factor
A twiddle factor, in fast Fourier transform (FFT) algorithms, is any of the trigonometric constant coefficients that are multiplied by the data in the course of the algorithm. This term was apparently coined by Gentleman & Sande in 1966, and has ...
s