A fast Fourier transform (FFT) is an algorithm that computes the discrete Fourier transform (DFT) of a sequence, or its inverse (IDFT). Fourier analysis converts a signal from its original domain (often time or space) to a representation in the frequency domain and vice versa. The DFT is obtained by decomposing a sequence of values into components of different frequencies.^{[1]} This operation is useful in many fields, but computing it directly from the definition is often too slow to be practical. An FFT rapidly computes such transformations by factorizing the DFT matrix into a product of sparse (mostly zero) factors.^{[2]} As a result, it manages to reduce the complexity of computing the DFT from , which arises if one simply applies the definition of DFT, to , where is the data size. The difference in speed can be enormous, especially for long data sets where N may be in the thousands or millions. In the presence of roundoff error, many FFT algorithms are much more accurate than evaluating the DFT definition directly. There are many different FFT algorithms based on a wide range of published theories, from simple complexnumber arithmetic to group theory and number theory.
Fast Fourier transforms are widely used for applications in engineering, science, and mathematics. The basic ideas were popularized in 1965, but some algorithms had been derived as early as 1805.^{[1]} In 1994, Gilbert Strang described the FFT as "the most important numerical algorithm of our lifetime",^{[3]}^{[4]} and it was included in Top 10 Algorithms of 20th Century by the IEEE magazine Computing in Science & Engineering.^{[5]}
The bestknown FFT algorithms depend upon the factorization of N, but there are FFTs with O(N log N) complexity for all N, even for prime N. Many FFT algorithms only depend on the fact that is an Nth primitive root of unity, and thus can be applied to analogous transforms over any finite field, such as numbertheoretic transforms. Since the inverse DFT is the same as the DFT, but with the opposite sign in the exponent and a 1/N factor, any FFT algorithm can easily be adapted for it.
YouTube Encyclopedic

1/5Views:2 433 174302 134100 402111 64617 853

✪ But what is the Fourier Transform? A visual introduction.

✪ The Fast Fourier Transform Algorithm

✪ 3. Divide & Conquer: FFT

✪ Fast Fourier Transform (FFT) in DTSP Discrete Time Signals Processing

✪ Fourier Transform  Part 1
Transcription
This right here is what we're going to build to, this video: A certain animated approach to thinking about a superimportant idea from math: The Fourier transform. For anyone unfamiliar with what that is, my #1 goal here is just for the video to be an introduction to that topic. But even for those of you who are already familiar with it, I still think that there's something fun and enriching about seeing what all of its components actually look like. The central example, to start, is gonna be the classic one: Decomposing frequencies from sound. But after that, I also really wanna show a glimpse of how this idea extends well beyond sound and frequency, and to many seemingly disparate areas of math, and even physics. Really, it is crazy just how ubiquitous this idea is. Let's dive in. This sound right here is a pure A. 440 beats per second. Meaning, if you were to measure the air pressure right next to your headphones, or your speaker, as a function of time, it would oscillate up and down around its usual equilibrium, in this wave. making 440 oscillations each second. A lowerpitched note, like a D, has the same structure, just fewer beats per second. And when both of them are played at once, what do you think the resulting pressure vs. time graph looks like? Well, at any point in time, this pressure difference is gonna be the sum of what it would be for each of those notes individually. Which, let's face it, is kind of a complicated thing to think about. At some points, the peaks match up with each other, resulting in a really high pressure. At other points, they tend to cancel out. And all in all, what you get is a waveish pressure vs. time graph, that is not a pure sine wave; it's something more complicated. And as you add in other notes, the wave gets more and more complicated. But right now, all it is is a combination of four pure frequencies. So it seems...needlessly complicated, given the low amount of information put into it. A microphone recording any sound just picks up on the air pressure at many different points in time. It only "sees" the final sum. So our central question is gonna be how you can take a signal like this, and decompose it into the pure frequencies that make it up. Pretty interesting, right? Adding up those signals really mixes them all together. So pulling them back apart...feels akin to unmixing multiple paint colors that have all been stirred up together. The general strategy is gonna be to build for ourselves a mathematical machine that treats signals with a given frequency... ..differently from how it treats other signals. To start, consider simply taking a pure signal say, with a lowly three beats per second, so that we can plot it easily. And let's limit ourselves to looking at a finite portion of this graph. In this case, the portion between zero seconds, and 4.5 seconds. The key idea, is gonna be to take this graph, and sort of wrap it up around a circle. Concretely, here's what I mean by that. Imagine a little rotating vector where each point in time its length is equal to the height of our graph for that time. So, high points of the graph correspond to a greater disance from the origin, and low points end up closer to the origin. And right now, I'm drawing it in such a way that moving forward two seconds in time corresponds to a single rotation around the circle. Our little vector drawing this wound up graph is rotating at half a cycle per second. So, this is important. There are two different frequencies at play here: There's the frequency of our signal, which goes up and down, three times per second. And then, separately, there's the frequency with which we're wrapping the graph around the circle. Which, at the moment, is half of a rotation per second. But we can adjust that second frequency however we want. Maybe we want to wrap it around faster... ..or maybe we go and wrap it around slower. And that choice of winding frequency determines what the wound up graph looks like. Some of the diagrams that come out of this can be pretty complicated; although, they are very pretty. But it's important to keep in mind that all that's happening here is that we're wrapping the signal around a circle. The vertical lines that I'm drawing up top, by the way, are just a way to keep track of the distance on the original graph that corresponds to a full rotation around the circle. So, lines spaced out by 1.5 seconds would mean it takes 1.5 seconds to make one full revolution. And at this point, we might have some sort of vague sense that something special will happen when the winding frequency matches the frequency of our signal: three beats per second. All the high points on the graph happen on the right side of the circle And all of the low points happen on the left. But how precisely can we take advantage of that in our attempt to build a frequencyunmixing machine? Well, imagine this graph is having some kind of mass to it, like a metal wire. This little dot is going to represent the center of mass of that wire. As we change the frequency, and the graph winds up differently, that center of mass kind of wobbles around a bit. And for most of the winding frequencies, the peaks and valleys are all spaced out around the circle in such a way that the center of mass stays pretty close to the origin. But! When the winding frequency is the same as the frequency of our signal, in this case, three cycles per second, all of the peaks are on the right, and all of the valleys are on the left.. ..so the center of mass is unusually far to the right. Here, to capture this, let's draw some kind of plot that keeps track of where that center of mass is for each winding frequency. Of course, the center of mass is a twodimensional thing, and requires two coordinates to fully keep track of, but for the moment, let's only keep track of the x coordinate. So, for a frequency of 0, when everything is bunched up on the right, this x coordinate is relatively high. And then, as you increase that winding frequency, and the graph balances out around the circle, the x coordinate of that center of mass goes closer to 0, and it just kind of wobbles around a bit. But then, at three beats per second, there's a spike as everything lines up to the right. This right here is the central construct, so let's sum up what we have so far: We have that original intensity vs. time graph, and then we have the wound up version of that in some twodimensional plane, and then, as a third thing, we have a plot for how the winding frequency influences the center of mass of that graph. And by the way, let's look back at those really low frequencies near 0. This big spike around 0 in our new frequency plot just corresponds to the fact that the whole cosine wave is shifted up. If I had chosen a signal oscillates around 0, dipping into negative values, then, as we play around with various winding frequences, this plot of the winding frequencies vs. center of mass would only have a spike at the value of three. But, negative values are a little bit weird and messy to think about especially for a first example, so let's just continue thinking in terms of the shiftedup graph. I just want you to understand that that spike around 0 only corresponds to the shift. Our main focus, as far as frequency decomposition is concerned, is that bump at three. This whole plot is what I'll call the "Almost Fourier Transform" of the original signal. There's a couple small distinctions between this and the actual Fourier transform, which I'll get to in a couple minutes, but already, you might be able to see how this machine lets us pick out the frequency of a signal. Just to play around with it a little bit more, take a different pure signal, let's say with a lower frequency of two beats per second, and do the same thing. Wind it around a circle, imagine different potential winding frequencies, and as you do that keep track of where the center of mass of that graph is, and then plot the x coordinate of that center of mass as you adjust the winding frequency. Just like before, we get a spike when the winding frequency is the same as the signal frequency, which in this case, is when it equals two cycles per second. But the real key point, the thing that makes this machine so delightful, is how it enables us to take a signal consisting of multiple frequencies, and pick out what they are. Imagine taking the two signals we just looked at: The wave with three beats per second, and the wave with two beats per second, and add them up. Like I said earlier, what you get is no longer a nice, pure cosine wave; it's something a little more complicated. But imagine throwing this into our windingfrequency machine... ..it is certainly the case that as you wrap this thing around, it looks a lot more complicated; you have this chaos (1) and chaos (2) and chaos (3) and chaos (4) and then WOOP! Things seem to line up really nicely at two cycles per second, and as you continue on it's more chaos (5) and more chaos (6) more chaos (7) chaos (8), chaos (9), chaos (10), WOOP! Things nicely align again at three cycles per second. And, like I said before, the wound up graph can look kind of busy and complicated, but all it is is the relatively simple idea of wrapping the graph around a circle. It's just a more complicated graph, and a pretty quick winding frequency. Now what's going on here with the two different spikes, is that if you were to take two signals, and then apply this AlmostFourier transform to each of them individually, and then add up the results, what you get is the same as if you first added up the signals, and then applied this AlmostFourier transorm. And the attentive viewers among you might wanna pause and ponder, and... ..convince yourself that what I just said is actually true. It's a pretty good test to verify for yourself that it's clear what exactly is being measured inside this winding machine. Now this property makes things really useful to us, because the transform of a pure frequency is close to 0 everywhere except for a spike around that frequency. So when you add together two pure frequencies, the transform graph just has these little peaks above the frequencies that went into it. So this little mathematical machine does exactly what we wanted. It pulls out the original frequencies from their jumbled up sums, unmixing the mixed bucket of paint. And before continuing into the full math that describes this operation, let's just get a quick glimpse of one context where this thing is useful: Sound editing. Let's say that you have some recording, and it's got an annoying high pitch that you'd like to filter out. Well, at first, your signal is coming in as a function of various intensities over time. Different voltages given to your speaker from one millisecond to the next. But we want to think of this in terms of frequencies, so, when you take the Fourier transform of that signal, the annoying high pitch is going to show up just as a spike at some high frequency. Filtering that out, by just smushing the spike down, what you'd be looking at is the Fourier transform of a sound that's just like your recording, only without that high frequency. Luckily, there's a notion of an inverse Fourier transform that tells you which signal would have produced this as its Fourier transform. I'll be talking about inverse much more fully in the next video, but long story short, applying the Fourier transform to the Fourier transform gives you back something close to the original function. Mm, kind of... this is... ..a little bit of a lie, but it's in the direction of the truth. And most of the reason that it's a lie is that I still have yet to tell you what the actual Fourier Transform is, since it's a little more complex than this xcoordinateofthecenterofmass idea. First off, bringing back this wound up graph, and looking at its center of mass, the x coordinate is really only half the story, right? I mean, this thing is in two dimensions, it's got a y coordinate as well. And, as is typical in math, whenever you're dealing with something twodimensional, it's elegant to think of it as the complex plane, where this center of mass is gonna be a complex number, that has both a real and an imaginary part. And the reason for talking in terms of complex numbers, rather than just saying, "It has two coordinates," is that complex numbers lend themselves to really nice descriptions of things that have to do with winding, and rotation. For example: Euler's formula famously tells us that if you take e to some number times i, you're gonna land on the point that you get if you were to walk that number of units around a circle with radius 1, counterclockwise starting on the right. So, imagine you wanted to describe rotating at a rate of one cycle per second. One thing that you could do is take the expression "e^2π*i*t," where t is the amount of time that has passed. Since, for a circle with radius 1, 2π describes the full length of its circumference. And... this is a little bit dizzying to look at, so maybe you wanna describe a different frequency... ..something lower and more reasonable... ..and for that, you would just multiply that time t in the exponent by the frequency, f. For example, if f was one tenth, then this vector makes one full turn every ten seconds, since the time t has to increase all the way to ten before the full exponent looks like 2πi. I have another video giving some intuition on why this is the behavior of e^x for imaginary inputs, if you're curious 😉, but for right now, we're just gonna take it as a given. Now why am I telling you this you this, you might ask. Well, it gives us a really nice way to write down the idea of winding up the graph into a single, tight little formula. First off, the convention in the context of Fourier transforms is to think about rotating in the clockwise direction, so let's go ahead and throw a negative sign up into that exponent. Now, take some function describing a signal intensity vs. time, like this pure cosine wave we had before, and call it g(t). If you multiply this exponential expression times g(t), it means that the rotating complex number is getting scaled up and down according to the value of this function. So you can think of this little rotating vector with its changing length as drawing the wound up graph. So think about it, this is awesome. This really small expression is a superelegant way to encapsulate the whole idea of winding a graph around a circle with a variable frequency f. And remember, that thing we want to do with this wound up graph is to track its center of mass. So think about what formula is going to capture that. Well, to approximate it at least, you might sample a whole bunch of times from the original signal, see where those points end up on the wound up graph, and then just take an average. That is, add them all together, as complex numbers, and then divide by the number of points that you've sampled. This will become more accurate if you sample more points which are closer together. And in the limit, rather than looking at the sum of a whole bunch of points divided by the number of points, you take an integral of this function, divided by the size of the time interval that we're looking at. Now the idea of integrating a complexvalued function might seem weird, and to anyone who's shaky with calculus, maybe even intimidating, but the underlying meaning here really doesn't require any calculus knowledge. The whole expression is just the center of mass of the wound up graph. So... Great! Stepbystep, we have built up this kind of complicated, but, let's face it, surprisingly small expression for the whole winding machine idea that I talked about. And now, there is only one final distinction to point out between this and the actual, honesttogoodness Fourier transform. Namely, just don't divide out by the time interval. The Fourier transform is just the integral part of this. What that means is that instead of looking at the center of mass, you would scale it up by some amount. If the portion of the original graph you were using spanned three seconds, you would multiply the center of mass by three. If it was spanning six seconds, you would multiply the center of mass by six. Physically, this has the effect that when a certain frequency persists for a long time, then the magnitude of the Fourier transform at that frequency is scaled up more and more. For example, what we're looking at right here is how when you have a pure frequency of two beats per second, and you wind it around the graph at two cycles per second, the center of mass stays in the same spot, right? It's just tracing out the same shape. But the longer that signal persists, the larger the value of the Fourier transform, at that frequency. For other frequencies, though, even if you just increase it by a bit, this is cancelled out by the fact that for longer time intervals you're giving the wound up graph more of a chance to balance itself around the circle. That is...a lot of different moving parts, so let's step back and summarize what we have so far. The Fourier transform of an intensity vs. time function, like g(t), is a new function, which doesn't have time as an input, but instead takes in a frequency, what I've been calling "the winding frequency." In terms of notation, by the way, the common convention is to call this new function "ghat," with a little circumflex on top of it. Now the output of this function is a complex number, some point in the 2D plane, that corresponds to the strength of a given frequency in the original signal. The plot that I've been graphing for the Fourier transform, is just the real component of that output, the xcoordinate But you could also graph the imaginary component separately, if you wanted a fuller description. And all of this is being encapsulated inside that formula that we built up. And out of context, you can imagine how seeing this formula would seem sort of daunting. But if you understand how exponentials correspond to rotation... ..how multiplying that by the function g(t) means drawing a wound up version of the graph, and how an integral of a complexvalued function can be interpreted in terms of a centerofmass idea, you can see how this whole thing carries with it a very rich, intuitive meaning. And, by the way, one quick small note before we can call this wrapped up. Even though in practice, with things like sound editing, you'll be integrating over a finite time interval, the theory of Fourier transforms is often phrased where the bounds of this integral are ∞ and ∞. Concretely, what that means is that you consider this expression for all possible finite time intervals, and you just ask, "What is its limit as that time interval grows to ∞?" And...man, oh man, there is so much more to say! So much, I don't wanna call it done here. This transform extends to corners of math well beyond the idea of extracting frequencies from signal. So, the next video I put out is gonna go through a couple of these, and that's really where things start getting interesting. So, stay subscribed for when that comes out, or an alternate option is to just binge a couple 3blue1brown videos so that the YouTube recommender is more inclined to show you new things that come out... ..really, the choice is yours! And to close things off, I have something pretty fun: A mathematical puzzler from this video's sponsor, Jane Street, who's looking to recruit more technical talent. So, let's say that you have a closed, bounded convex set C sitting in 3D space, and then let B be the boundary of that space, the surface of your complex blob. Now imagine taking every possible pair of points on that surface, and adding them up, doing a vector sum. Let's name this set of all possible sums D. Your task is to prove that D is also a convex set. So, Jane Street is a quantitative trading firm, and if you're the kind of person who enjoys math and solving puzzles like this, the team there really values intellectual curiosity. So, they might be interested in hiring you. And they're looking both for fulltime employees and interns. For my part, I can say that some people I've interacted with there just seem to love math, and sharing math, and when they're hiring they look less at a background in finance than they do at how you think, how you learn, and how you solve problems, hence the sponsorship of a 3blue1brown video. If you want the answer to that puzzler, or to learn more about what they do, or to apply for open positions, go to janestreet.com/3b1b/
Contents
History
The development of fast algorithms for DFT can be traced to Gauss's unpublished work in 1805 when he needed it to interpolate the orbit of asteroids Pallas and Juno from sample observations.^{[6]}^{[7]} His method was very similar to the one published in 1965 by Cooley and Tukey, who are generally credited for the invention of the modern generic FFT algorithm. While Gauss's work predated even Fourier's results in 1822, he did not analyze the computation time and eventually used other methods to achieve his goal.
Between 1805 and 1965, some versions of FFT were published by other authors. Frank Yates in 1932 published his version called interaction algorithm, which provided efficient computation of Hadamard and Walsh transforms.^{[8]} Yates' algorithm is still used in the field of statistical design and analysis of experiments. In 1942, G. C. Danielson and Cornelius Lanczos published their version to compute DFT for xray crystallography, a field where calculation of Fourier transforms presented a formidable bottleneck.^{[9]}^{[10]} While many methods in the past had focused on reducing the constant factor for computation by taking advantage of "symmetries", Danielson and Lanczos realized that one could use the "periodicity" and apply a "doubling trick" to get runtime.^{[11]}
James Cooley and John Tukey published a more general version of FFT in 1965 that is applicable when N is composite and not necessarily a power of 2.^{[12]} Tukey came up with the idea during a meeting of President Kennedy's Science Advisory Committee where a discussion topic involved detecting nuclear tests by the Soviet Union by setting up sensors to surround the country from outside. To analyze the output of these sensors, a fast Fourier transform algorithm would be needed. In discussion with Tukey, Richard Garwin recognized the general applicability of the algorithm not just to national security problems, but also to a wide range of problems including one of immediate interest to him, determining the periodicities of the spin orientations in a 3D crystal of Helium3.^{[13]} Garwin gave Tukey's idea to Cooley (both worked at IBM's Watson labs) for implementation.^{[14]} Cooley and Tukey published the paper in a relatively short time of six months.^{[15]} As Tukey did not work at IBM, the patentability of the idea was doubted and the algorithm went into the public domain, which, through the computing revolution of the next decade, made FFT one of the indispensable algorithms in digital signal processing.
Definition
Let x_{0}, ...., x_{N−1} be complex numbers. The DFT is defined by the formula
where is a primitive Nth root of 1.
Evaluating this definition directly requires operations: there are N outputs X_{k}, and each output requires a sum of N terms. An FFT is any method to compute the same results in operations. All known FFT algorithms require Θ operations, although there is no known proof that a lower complexity score is impossible.^{[16]}
To illustrate the savings of an FFT, consider the count of complex multiplications and additions for N=4096 data points. Evaluating the DFT's sums directly involves N^{2} complex multiplications and N(N−1) complex additions, of which operations can be saved by eliminating trivial operations such as multiplications by 1, leaving about 30 million operations. On the other hand, the radix2 Cooley–Tukey algorithm, for N a power of 2, can compute the same result with only (N/2)log_{2}(N) complex multiplications (again, ignoring simplifications of multiplications by 1 and similar) and N log_{2}(N) complex additions, in total about 30,000 operations  a thousand times less than with direct evaluation. In practice, actual performance on modern computers is usually dominated by factors other than the speed of arithmetic operations and the analysis is a complicated subject (see, e.g., Frigo & Johnson, 2005),^{[17]} but the overall improvement from to remains.
Algorithms
Cooley–Tukey algorithm
By far the most commonly used FFT is the Cooley–Tukey algorithm. This is a divide and conquer algorithm that recursively breaks down a DFT of any composite size N = N_{1}N_{2} into many smaller DFTs of sizes N_{1} and N_{2}, along with O(N) multiplications by complex roots of unity traditionally called twiddle factors (after Gentleman and Sande, 1966^{[18]}).
This method (and the general idea of an FFT) was popularized by a publication of Cooley and Tukey in 1965,^{[12]} but it was later discovered^{[1]} that those two authors had independently reinvented an algorithm known to Carl Friedrich Gauss around 1805^{[19]} (and subsequently rediscovered several times in limited forms).
The best known use of the Cooley–Tukey algorithm is to divide the transform into two pieces of size N/2 at each step, and is therefore limited to poweroftwo sizes, but any factorization can be used in general (as was known to both Gauss and Cooley/Tukey^{[1]}). These are called the radix2 and mixedradix cases, respectively (and other variants such as the splitradix FFT have their own names as well). Although the basic idea is recursive, most traditional implementations rearrange the algorithm to avoid explicit recursion. Also, because the Cooley–Tukey algorithm breaks the DFT into smaller DFTs, it can be combined arbitrarily with any other algorithm for the DFT, such as those described below.
Other FFT algorithms
There are other FFT algorithms distinct from Cooley–Tukey.
Cornelius Lanczos did pioneering work on the FFT and FFS (fast Fourier sampling method) with G. C. Danielson (1940).
For N = N_{1}N_{2} with coprime N_{1} and N_{2}, one can use the primefactor (Good–Thomas) algorithm (PFA), based on the Chinese remainder theorem, to factorize the DFT similarly to Cooley–Tukey but without the twiddle factors. The Rader–Brenner algorithm (1976)^{[20]} is a Cooley–Tukeylike factorization but with purely imaginary twiddle factors, reducing multiplications at the cost of increased additions and reduced numerical stability; it was later superseded by the splitradix variant of Cooley–Tukey (which achieves the same multiplication count but with fewer additions and without sacrificing accuracy). Algorithms that recursively factorize the DFT into smaller operations other than DFTs include the Bruun and QFT algorithms. (The Rader–Brenner^{[20]} and QFT algorithms were proposed for poweroftwo sizes, but it is possible that they could be adapted to general composite N. Bruun's algorithm applies to arbitrary even composite sizes.) Bruun's algorithm, in particular, is based on interpreting the FFT as a recursive factorization of the polynomial z^{N} − 1, here into realcoefficient polynomials of the form z^{M} − 1 and z^{2M} + az^{M} + 1.
Another polynomial viewpoint is exploited by the Winograd FFT algorithm,^{[21]}^{[22]} which factorizes z^{N} − 1 into cyclotomic polynomials—these often have coefficients of 1, 0, or −1, and therefore require few (if any) multiplications, so Winograd can be used to obtain minimalmultiplication FFTs and is often used to find efficient algorithms for small factors. Indeed, Winograd showed that the DFT can be computed with only O(N) irrational multiplications, leading to a proven achievable lower bound on the number of multiplications for poweroftwo sizes; unfortunately, this comes at the cost of many more additions, a tradeoff no longer favorable on modern processors with hardware multipliers. In particular, Winograd also makes use of the PFA as well as an algorithm by Rader for FFTs of prime sizes.
Rader's algorithm, exploiting the existence of a generator for the multiplicative group modulo prime N, expresses a DFT of prime size N as a cyclic convolution of (composite) size N−1, which can then be computed by a pair of ordinary FFTs via the convolution theorem (although Winograd uses other convolution methods). Another primesize FFT is due to L. I. Bluestein, and is sometimes called the chirpz algorithm; it also reexpresses a DFT as a convolution, but this time of the same size (which can be zeropadded to a power of two and evaluated by radix2 Cooley–Tukey FFTs, for example), via the identity
Hexagonal Fast Fourier Transform aims at computing an efficient FFT for the hexagonally sampled data by using a new addressing scheme for hexagonal grids, called Array Set Addressing (ASA).
FFT algorithms specialized for real and/or symmetric data
In many applications, the input data for the DFT are purely real, in which case the outputs satisfy the symmetry
and efficient FFT algorithms have been designed for this situation (see e.g. Sorensen, 1987).^{[23]}^{[24]} One approach consists of taking an ordinary algorithm (e.g. Cooley–Tukey) and removing the redundant parts of the computation, saving roughly a factor of two in time and memory. Alternatively, it is possible to express an evenlength realinput DFT as a complex DFT of half the length (whose real and imaginary parts are the even/odd elements of the original real data), followed by O(N) postprocessing operations.
It was once believed that realinput DFTs could be more efficiently computed by means of the discrete Hartley transform (DHT), but it was subsequently argued that a specialized realinput DFT algorithm (FFT) can typically be found that requires fewer operations than the corresponding DHT algorithm (FHT) for the same number of inputs.^{[23]} Bruun's algorithm (above) is another method that was initially proposed to take advantage of real inputs, but it has not proved popular.
There are further FFT specializations for the cases of real data that have even/odd symmetry, in which case one can gain another factor of (roughly) two in time and memory and the DFT becomes the discrete cosine/sine transform(s) (DCT/DST). Instead of directly modifying an FFT algorithm for these cases, DCTs/DSTs can also be computed via FFTs of real data combined with O(N) pre/post processing.
Computational issues
Bounds on complexity and operation counts
Unsolved problem in computer science: What is the lower bound on the complexity of fast Fourier transform algorithms? Can they be faster than ? (more unsolved problems in computer science)

A fundamental question of longstanding theoretical interest is to prove lower bounds on the complexity and exact operation counts of fast Fourier transforms, and many open problems remain. It is not even rigorously proved whether DFTs truly require Ω(N log N) (i.e., order N log N or greater) operations, even for the simple case of power of two sizes, although no algorithms with lower complexity are known. In particular, the count of arithmetic operations is usually the focus of such questions, although actual performance on modernday computers is determined by many other factors such as cache or CPU pipeline optimization.
Following pioneering work by Winograd (1978),^{[21]} a tight Θ(N) lower bound is known for the number of real multiplications required by an FFT. It can be shown that only irrational real multiplications are required to compute a DFT of poweroftwo length . Moreover, explicit algorithms that achieve this count are known (Heideman & Burrus, 1986;^{[25]} Duhamel, 1990^{[26]}). Unfortunately, these algorithms require too many additions to be practical, at least on modern computers with hardware multipliers (Duhamel, 1990;^{[26]} Frigo & Johnson, 2005).^{[17]}
A tight lower bound is not known on the number of required additions, although lower bounds have been proved under some restrictive assumptions on the algorithms. In 1973, Morgenstern^{[27]} proved an Ω(N log N) lower bound on the addition count for algorithms where the multiplicative constants have bounded magnitudes (which is true for most but not all FFT algorithms). This result, however, applies only to the unnormalized Fourier transform (which is a scaling of a unitary matrix by a factor of ), and does not explain why the Fourier matrix is harder to compute than any other unitary matrix (including the identity matrix) under the same scaling. Pan (1986)^{[28]} proved an Ω(N log N) lower bound assuming a bound on a measure of the FFT algorithm's "asynchronicity", but the generality of this assumption is unclear. For the case of poweroftwo N, Papadimitriou (1979)^{[29]} argued that the number of complexnumber additions achieved by Cooley–Tukey algorithms is optimal under certain assumptions on the graph of the algorithm (his assumptions imply, among other things, that no additive identities in the roots of unity are exploited). (This argument would imply that at least real additions are required, although this is not a tight bound because extra additions are required as part of complexnumber multiplications.) Thus far, no published FFT algorithm has achieved fewer than complexnumber additions (or their equivalent) for poweroftwo N.
A third problem is to minimize the total number of real multiplications and additions, sometimes called the "arithmetic complexity" (although in this context it is the exact count and not the asymptotic complexity that is being considered). Again, no tight lower bound has been proven. Since 1968, however, the lowest published count for poweroftwo N was long achieved by the splitradix FFT algorithm, which requires real multiplications and additions for N > 1. This was recently reduced to (Johnson and Frigo, 2007;^{[16]} Lundy and Van Buskirk, 2007^{[30]}). A slightly larger count (but still better than split radix for N≥256) was shown to be provably optimal for N≤512 under additional restrictions on the possible algorithms (splitradixlike flowgraphs with unitmodulus multiplicative factors), by reduction to a satisfiability modulo theories problem solvable by brute force (Haynal & Haynal, 2011).^{[31]}
Most of the attempts to lower or prove the complexity of FFT algorithms have focused on the ordinary complexdata case, because it is the simplest. However, complexdata FFTs are so closely related to algorithms for related problems such as realdata FFTs, discrete cosine transforms, discrete Hartley transforms, and so on, that any improvement in one of these would immediately lead to improvements in the others (Duhamel & Vetterli, 1990).^{[32]}
Approximations
All of the FFT algorithms discussed above compute the DFT exactly (i.e. neglecting floatingpoint errors). A few "FFT" algorithms have been proposed, however, that compute the DFT approximately, with an error that can be made arbitrarily small at the expense of increased computations. Such algorithms trade the approximation error for increased speed or other properties. For example, an approximate FFT algorithm by Edelman et al. (1999)^{[33]} achieves lower communication requirements for parallel computing with the help of a fast multipole method. A waveletbased approximate FFT by Guo and Burrus (1996)^{[34]} takes sparse inputs/outputs (time/frequency localization) into account more efficiently than is possible with an exact FFT. Another algorithm for approximate computation of a subset of the DFT outputs is due to Shentov et al. (1995).^{[35]} The Edelman algorithm works equally well for sparse and nonsparse data, since it is based on the compressibility (rank deficiency) of the Fourier matrix itself rather than the compressibility (sparsity) of the data. Conversely, if the data are sparse—that is, if only K out of N Fourier coefficients are nonzero—then the complexity can be reduced to O(K log(N)log(N/K)), and this has been demonstrated to lead to practical speedups compared to an ordinary FFT for N/K > 32 in a largeN example (N = 2^{22}) using a probabilistic approximate algorithm (which estimates the largest K coefficients to several decimal places).^{[36]}
Accuracy
Even the "exact" FFT algorithms have errors when finiteprecision floatingpoint arithmetic is used, but these errors are typically quite small; most FFT algorithms, e.g. Cooley–Tukey, have excellent numerical properties as a consequence of the pairwise summation structure of the algorithms. The upper bound on the relative error for the Cooley–Tukey algorithm is O(ε log N), compared to O(εN^{3/2}) for the naïve DFT formula,^{[18]} where ε is the machine floatingpoint relative precision. In fact, the root mean square (rms) errors are much better than these upper bounds, being only O(ε √log N) for Cooley–Tukey and O(ε √N) for the naïve DFT (Schatzman, 1996).^{[37]} These results, however, are very sensitive to the accuracy of the twiddle factors used in the FFT (i.e. the trigonometric function values), and it is not unusual for incautious FFT implementations to have much worse accuracy, e.g. if they use inaccurate trigonometric recurrence formulas. Some FFTs other than Cooley–Tukey, such as the Rader–Brenner algorithm, are intrinsically less stable.
In fixedpoint arithmetic, the finiteprecision errors accumulated by FFT algorithms are worse, with rms errors growing as O(√N) for the Cooley–Tukey algorithm (Welch, 1969).^{[38]} Moreover, even achieving this accuracy requires careful attention to scaling to minimize loss of precision, and fixedpoint FFT algorithms involve rescaling at each intermediate stage of decompositions like Cooley–Tukey.
To verify the correctness of an FFT implementation, rigorous guarantees can be obtained in O(N log N) time by a simple procedure checking the linearity, impulseresponse, and timeshift properties of the transform on random inputs (Ergün, 1995).^{[39]}
Multidimensional FFTs
As defined in the multidimensional DFT article, the multidimensional DFT
transforms an array x_{n} with a ddimensional vector of indices by a set of d nested summations (over for each j), where the division n/N, defined as , is performed elementwise. Equivalently, it is the composition of a sequence of d sets of onedimensional DFTs, performed along one dimension at a time (in any order).
This compositional viewpoint immediately provides the simplest and most common multidimensional DFT algorithm, known as the rowcolumn algorithm (after the twodimensional case, below). That is, one simply performs a sequence of d onedimensional FFTs (by any of the above algorithms): first you transform along the n_{1} dimension, then along the n_{2} dimension, and so on (or actually, any ordering works). This method is easily shown to have the usual O(N log N) complexity, where is the total number of data points transformed. In particular, there are N/N_{1} transforms of size N_{1}, etcetera, so the complexity of the sequence of FFTs is:
In two dimensions, the x_{k} can be viewed as an matrix, and this algorithm corresponds to first performing the FFT of all the rows (resp. columns), grouping the resulting transformed rows (resp. columns) together as another matrix, and then performing the FFT on each of the columns (resp. rows) of this second matrix, and similarly grouping the results into the final result matrix.
In more than two dimensions, it is often advantageous for cache locality to group the dimensions recursively. For example, a threedimensional FFT might first perform twodimensional FFTs of each planar "slice" for each fixed n_{1}, and then perform the onedimensional FFTs along the n_{1} direction. More generally, an asymptotically optimal cacheoblivious algorithm consists of recursively dividing the dimensions into two groups and that are transformed recursively (rounding if d is not even) (see Frigo and Johnson, 2005).^{[17]} Still, this remains a straightforward variation of the rowcolumn algorithm that ultimately requires only a onedimensional FFT algorithm as the base case, and still has O(N log N) complexity. Yet another variation is to perform matrix transpositions in between transforming subsequent dimensions, so that the transforms operate on contiguous data; this is especially important for outofcore and distributed memory situations where accessing noncontiguous data is extremely timeconsuming.
There are other multidimensional FFT algorithms that are distinct from the rowcolumn algorithm, although all of them have O(N log N) complexity. Perhaps the simplest nonrowcolumn FFT is the vectorradix FFT algorithm, which is a generalization of the ordinary Cooley–Tukey algorithm where one divides the transform dimensions by a vector of radices at each step. (This may also have cache benefits.) The simplest case of vectorradix is where all of the radices are equal (e.g. vectorradix2 divides all of the dimensions by two), but this is not necessary. Vector radix with only a single nonunit radix at a time, i.e. , is essentially a rowcolumn algorithm. Other, more complicated, methods include polynomial transform algorithms due to Nussbaumer (1977),^{[40]} which view the transform in terms of convolutions and polynomial products. See Duhamel and Vetterli (1990)^{[32]} for more information and references.
Other generalizations
An O(N^{5/2}log N) generalization to spherical harmonics on the sphere S^{2} with N^{2} nodes was described by Mohlenkamp,^{[41]} along with an algorithm conjectured (but not proven) to have O(N^{2} log^{2}(N)) complexity; Mohlenkamp also provides an implementation in the libftsh library.^{[42]} A sphericalharmonic algorithm with O(N^{2}log N) complexity is described by Rokhlin and Tygert.^{[43]}
The fast folding algorithm is analogous to the FFT, except that it operates on a series of binned waveforms rather than a series of real or complex scalar values. Rotation (which in the FFT is multiplication by a complex phasor) is a circular shift of the component waveform.
Various groups have also published "FFT" algorithms for nonequispaced data, as reviewed in Potts et al. (2001).^{[44]} Such algorithms do not strictly compute the DFT (which is only defined for equispaced data), but rather some approximation thereof (a nonuniform discrete Fourier transform, or NDFT, which itself is often computed only approximately). More generally there are various other methods of spectral estimation.
Applications
FFT's importance derives from the fact that in signal processing and image processing it has made working in frequency domain equally computationally feasible as working in temporal or spatial domain. Some of the important applications of FFT includes,^{[15]}^{[45]}
 Fast large integer and polynomial multiplication
 Efficient matrixvector multiplication for Toeplitz, circulant and other structured matrices
 Filtering algorithms (see overlapadd and overlapsave methods)
 Fast algorithms for discrete cosine or sine transforms (example, Fast DCT used for JPEG, MP3/MPEG encoding)
 Fast Chebyshev approximation
 Fast discrete Hartley transform
 Solving difference equations
 Computation of isotopic distributions.^{[46]}
Research areas
 Big FFTs: With the explosion of big data in fields such as astronomy, the need for 512k FFTs has arisen for certain interferometry calculations. The data collected by projects such as WMAP and LIGO require FFTs of tens of billions of points. As this size does not fit into main memory, so called outofcore FFTs are an active area of research.^{[47]}
 Approximate FFTs: For applications such as MRI, it is necessary to compute DFTs for nonuniformly spaced grid points and/or frequencies. Multipole based approaches can compute approximate quantities with factor of runtime increase.^{[48]}
 Group FFTs: The FFT may also be explained and interpreted using group representation theory that allows for further generalization. A function on any compact group, including non cyclic, has an expansion in terms of a basis of irreducible matrix elements. It remains active area of research to find efficient algorithm for performing this change of basis. Applications including efficient spherical harmonic expansion, analyzing certain Markov processes, robotics etc.^{[49]}
 Quantum FFTs: Shor's fast algorithm for integer factorization on a quantum computer has a subroutine to compute DFT of a binary vector. This is implemented as sequence of 1 or 2bit quantum gates now known as quantum FFT, which is effectively the Cooley–Tukey FFT realized as a particular factorization of the Fourier matrix. Extension to these ideas is currently being explored.
Language reference
Language  Command/Method  Prerequisites 

R  stats::fft(x)  None 
Octave/MATLAB  fft(x)  None 
Python  fft.fft(x)  numpy 
Mathematica  Fourier[x]  None 
Julia  fft(A [,dims])  FFTW 
See also
FFTrelated algorithms:
 Cooley–Tukey FFT algorithm
 Primefactor FFT algorithm
 Bruun's FFT algorithm
 Rader's FFT algorithm
 Bluestein's FFT algorithm
 Goertzel algorithm – Computes individual terms of discrete Fourier transform
FFT implementations:
 ALGLIB – C++ and C# library with real/complex FFT implementation.
 FFTW "Fastest Fourier Transform in the West" – C library for the discrete Fourier transform (DFT) in one or more dimensions.
 FFTS – The Fastest Fourier Transform in the South.
 FFTPACK – another Fortran FFT library (public domain)
 Math Kernel Library
 cuFFT  FFT for GPU accelerated CUDA
Other links:
 Overlap add/Overlap save – efficient convolution methods using FFT for long signals
 Odlyzko–Schönhage algorithm applies the FFT to finite Dirichlet series.
 Schönhage–Strassen algorithm  asymptotically fast multiplication algorithm for large integers
 Butterfly diagram – a diagram used to describe FFTs.
 Spectral music (involves application of FFT analysis to musical composition)
 Spectrum analyzer – any of several devices that perform an FFT
 Time series
 Fast Walsh–Hadamard transform
 Generalized distributive law
 Multidimensional transform
 Multidimensional discrete convolution
 DFT matrix
References
 ^ ^{a} ^{b} ^{c} ^{d} Heideman, Michael T.; Johnson, Don H.; Burrus, Charles Sidney (1984). "Gauss and the history of the fast Fourier transform" (PDF). IEEE ASSP Magazine. 1 (4): 14–21. doi:10.1109/MASSP.1984.1162257.
 ^ Van Loan, Charles (1992). Computational Frameworks for the Fast Fourier Transform. SIAM.
 ^ Strang, Gilbert (May–June 1994). "Wavelets". American Scientist. 82 (3): 250–255. JSTOR 29775194.
 ^ Kent, Ray D.; Read, Charles (2002). Acoustic Analysis of Speech. ISBN 0769301126. ISBN 9780769301129.
 ^ Dongarra, Jack; Sullivan, Francis (January 2000). "Guest Editors Introduction to the top 10 algorithms". Computing in Science Engineering. 2 (1): 22–23. doi:10.1109/MCISE.2000.814652. ISSN 15219615.
 ^ Gauss, Carl Friedrich (1866). "Theoria interpolationis methodo nova tractata" [Theory regarding a new method of interpolation]. Nachlass (Unpublished manuscript). Werke (in Latin and German). 3. Göttingen, Germany: Königlichen Gesellschaft der Wissenschaften zu Göttingen. pp. 265–303.
 ^ Heideman, Michael T.; Johnson, Don H.; Burrus, Charles Sidney (19850901). "Gauss and the history of the fast Fourier transform". Archive for History of Exact Sciences. 34 (3): 265–277. CiteSeerX 10.1.1.309.181. doi:10.1007/BF00348431. ISSN 00039519.
 ^ Yates, Frank (1937). "The design and analysis of factorial experiments". Technical Communication no. 35 of the Commonwealth Bureau of Soils.
 ^ Danielson, Gordon C.; Lanczos, Cornelius (1942). "Some improvements in practical Fourier analysis and their application to xray scattering from liquids". Journal of the Franklin Institute. 233 (4): 365–380. doi:10.1016/S00160032(42)907671.
 ^ Lanczos, Cornelius (1956). Applied Analysis. Prentice–Hall.
 ^ Cooley, James W.; Lewis, Peter A. W.; Welch, Peter D. (June 1967). "Historical notes on the fast Fourier transform". IEEE Transactions on Audio and Electroacoustics. 15 (2): 76–79. CiteSeerX 10.1.1.467.7209. doi:10.1109/TAU.1967.1161903. ISSN 00189278.
 ^ ^{a} ^{b} Cooley, James W.; Tukey, John W. (1965). "An algorithm for the machine calculation of complex Fourier series". Mathematics of Computation. 19 (90): 297–301. doi:10.1090/S00255718196501785861. ISSN 00255718.
 ^ Cooley, James W. (1987). The ReDiscovery of the Fast Fourier Transform Algorithm (PDF). Microchimica Acta. III. Vienna, Austria. pp. 33–45.
 ^ Garwin, Richard (June 1969). "The Fast Fourier Transform As an Example of the Difficulty in Gaining Wide Use for a New Technique" (PDF). IEEE Transactions on Audio and Electroacoustics. AU17 (2): 68–72.
 ^ ^{a} ^{b} Rockmore, Daniel N. (January 2000). "The FFT: an algorithm the whole family can use". Computing in Science Engineering. 2 (1): 60–64. CiteSeerX 10.1.1.17.228. doi:10.1109/5992.814659. ISSN 15219615.
 ^ ^{a} ^{b} Frigo, Matteo; Johnson, Steven G. (January 2007) [20061219]. "A Modified SplitRadix FFT With Fewer Arithmetic Operations". IEEE Transactions on Signal Processing. 55 (1): 111–119. CiteSeerX 10.1.1.582.5497. doi:10.1109/tsp.2006.882087.
 ^ ^{a} ^{b} ^{c} Frigo, Matteo; Johnson, Steven G. (2005). "The Design and Implementation of FFTW3" (PDF). Proceedings of the IEEE. 93 (2): 216–231. CiteSeerX 10.1.1.66.3097. doi:10.1109/jproc.2004.840301.
 ^ ^{a} ^{b} Gentleman, W. Morven; Sande, G. (1966). "Fast Fourier transforms—for fun and profit". Proceedings of the AFIPS. 29: 563–578. doi:10.1145/1464291.1464352.
 ^ Gauss, Carl Friedrich (1866) [1805]. Theoria interpolationis methodo nova tractata. Werke (in Latin and German). 3. Göttingen, Germany: Königliche Gesellschaft der Wissenschaften. pp. 265–327.
 ^ ^{a} ^{b} Brenner, Norman M.; Rader, Charles M. (1976). "A New Principle for Fast Fourier Transformation". IEEE Transactions on Acoustics, Speech, and Signal Processing. 24 (3): 264–266. doi:10.1109/TASSP.1976.1162805.
 ^ ^{a} ^{b} Winograd, Shmuel (1978). "On computing the discrete Fourier transform". Mathematics of Computation. 32 (141): 175–199. doi:10.1090/S00255718197804683064. JSTOR 2006266. PMC 430186. PMID 16592303.
 ^ Winograd, Shmuel (1979). "On the multiplicative complexity of the discrete Fourier transform". Advances in Mathematics. 32 (2): 83–117. doi:10.1016/00018708(79)900379.
 ^ ^{a} ^{b} Sorensen, Henrik V.; Jones, Douglas L.; Heideman, Michael T.; Burrus, Charles Sidney (1987). "Realvalued fast Fourier transform algorithms". IEEE Transactions on Acoustics, Speech, and Signal Processing. 35 (6): 849–863. CiteSeerX 10.1.1.205.4523. doi:10.1109/TASSP.1987.1165220.
 ^ Sorensen, Henrik V.; Jones, Douglas L.; Heideman, Michael T.; Burrus, Charles Sidney (1987). "Corrections to "Realvalued fast Fourier transform algorithms"". IEEE Transactions on Acoustics, Speech, and Signal Processing. 35 (9): 1353. doi:10.1109/TASSP.1987.1165284.
 ^ Heideman, Michael T.; Burrus, Charles Sidney (1986). "On the number of multiplications necessary to compute a length2^{n} DFT". IEEE Transactions on Acoustics, Speech, and Signal Processing. 34 (1): 91–95. doi:10.1109/TASSP.1986.1164785.
 ^ ^{a} ^{b} Duhamel, Pierre (1990). "Algorithms meeting the lower bounds on the multiplicative complexity of length2^{n} DFTs and their connection with practical algorithms". IEEE Transactions on Acoustics, Speech, and Signal Processing. 38 (9): 1504–1511. doi:10.1109/29.60070.
 ^ Morgenstern, Jacques (1973). "Note on a lower bound of the linear complexity of the fast Fourier transform". Journal of the ACM. 20 (2): 305–306. doi:10.1145/321752.321761.
 ^ Pan, Victor Ya. (19860102). "The tradeoff between the additive complexity and the asynchronicity of linear and bilinear algorithms". Information Processing Letters. 22 (1): 11–14. doi:10.1016/00200190(86)900359. Retrieved 20171031.
 ^ Papadimitriou, Christos H. (1979). "Optimality of the fast Fourier transform". Journal of the ACM. 26: 95–102. doi:10.1145/322108.322118.
 ^ Lundy, Thomas J.; Van Buskirk, James (2007). "A new matrix approach to real FFTs and convolutions of length 2^{k}". Computing. 80 (1): 23–45. doi:10.1007/s0060700702226.
 ^ Haynal, Steve; Haynal, Heidi (2011). "Generating and Searching Families of FFT Algorithms" (PDF). Journal on Satisfiability, Boolean Modeling and Computation. 7: 145–187. Archived from the original (PDF) on 20120426.
 ^ ^{a} ^{b} Duhamel, Pierre; Vetterli, Martin (1990). "Fast Fourier transforms: a tutorial review and a state of the art". Signal Processing. 19 (4): 259–299. doi:10.1016/01651684(90)90158U.
 ^ Edelman, Alan; McCorquodale, Peter; Toledo, Sivan (1999). "The Future Fast Fourier Transform?" (PDF). SIAM Journal on Scientific Computing. 20 (3): 1094–1114. CiteSeerX 10.1.1.54.9339. doi:10.1137/S1064827597316266.
 ^ Guo, Haitao; Burrus, Charles Sidney (1996). "Fast approximate Fourier transform via wavelets transform". Proceedings of SPIE. Wavelet Applications in Signal and Image Processing IV. 2825: 250–259. CiteSeerX 10.1.1.54.3984. doi:10.1117/12.255236.
 ^ Shentov, Ognjan V.; Mitra, Sanjit K.; Heute, Ulrich; Hossen, Abdul N. (1995). "Subband DFT. I. Definition, interpretations and extensions". Signal Processing. 41 (3): 261–277. doi:10.1016/01651684(94)001037.
 ^ Hassanieh, Haitham; Indyk, Piotr; Katabi, Dina; Price, Eric (January 2012). "Simple and Practical Algorithm for Sparse Fourier Transform" (PDF). ACMSIAM Symposium on Discrete Algorithms (SODA). (NB. See also the sFFT Web Page.)
 ^ Schatzman, James C. (1996). "Accuracy of the discrete Fourier transform and the fast Fourier transform". SIAM Journal on Scientific Computing. 17 (5): 1150–1166. CiteSeerX 10.1.1.495.9184. doi:10.1137/s1064827593247023.
 ^ Welch, Peter D. (1969). "A fixedpoint fast Fourier transform error analysis". IEEE Transactions on Audio and Electroacoustics. 17 (2): 151–157. doi:10.1109/TAU.1969.1162035.
 ^ Ergün, Funda (1995). Testing multivariate linear functions: Overcoming the generator bottleneck. Proceedings of the 27th ACM Symposium On the Theory of Computing. Kyoto, Japan. pp. 407–416. doi:10.1145/225058.225167. ISBN 9780897917186.
 ^ Nussbaumer, Henri J. (1977). "Digital filtering using polynomial transforms". Electronics Letters. 13 (13): 386–387. doi:10.1049/el:19770280.
 ^ Mohlenkamp, Martin J. (1999). "A Fast Transform for Spherical Harmonics" (PDF). Journal of Fourier Analysis and Applications. 5 (2–3): 159–184. CiteSeerX 10.1.1.135.9830. doi:10.1007/BF01261607. Retrieved 20180111.
 ^ libftsh library
 ^ Rokhlin, Vladimir; Tygert, Mark (2006). "Fast Algorithms for Spherical Harmonic Expansions" (PDF). SIAM Journal on Scientific Computing. 27 (6): 1903–1928. CiteSeerX 10.1.1.125.7415. doi:10.1137/050623073. Retrieved 20140918. [1]
 ^ Potts, Daniel; Steidl, Gabriele; Tasche, Manfred (2001). "Fast Fourier transforms for nonequispaced data: A tutorial" (PDF). In Benedetto, J. J.; Ferreira, P. (eds.). Modern Sampling Theory: Mathematics and Applications. Birkhäuser.
 ^ Chu, Eleanor; George, Alan (19991111) [19991111]. "Chapter 16". Inside the FFT Black Box: Serial and Parallel Fast Fourier Transform Algorithms. CRC Press. pp. 153–168. ISBN 9781420049961.
 ^ FernandezdeCossio Diaz, Jorge; FernandezdeCossio, Jorge (20120808). "Computation of Isotopic Peak CenterMass Distribution by Fourier Transform". Analytical Chemistry. 84 (16): 7052–7056. doi:10.1021/ac301296a. ISSN 00032700. PMID 22873736.
 ^ Cormen, Thomas H.; Nicol, David M. (1998). "Performing outofcore FFTs on parallel disk systems" (PDF). Parallel Computing. 24 (1): 5–20. CiteSeerX 10.1.1.44.8212. doi:10.1016/S01678191(97)001142.
 ^ Dutt, Alok; Rokhlin, Vladimir (19931101). "Fast Fourier Transforms for Nonequispaced Data". SIAM Journal on Scientific Computing. 14 (6): 1368–1393. doi:10.1137/0914081. ISSN 10648275.
 ^ Rockmore, Daniel N. (2004). Byrnes, Jim (ed.). Recent Progress and Applications in Group FFTs. Computational Noncommutative Algebra and Applications. NATO Science Series II: Mathematics, Physics and Chemistry. 136. Springer Netherlands. pp. 227–254. CiteSeerX 10.1.1.324.4700. doi:10.1007/1402023073_9. ISBN 9781402019821.
Further reading
 Brigham, E. Oran (2002). "The Fast Fourier Transform". New York, USA: PrenticeHall.
 Cormen, Thomas H.; Leiserson, Charles E.; Rivest, Ronald L.; Stein, Clifford (2001). "Chapter 30: Polynomials and the FFT". Introduction to Algorithms (2 ed.). MIT Press / McGrawHill. ISBN 0262032937. ISBN 9780262032933.
 Elliott, Douglas F.; Rao, K. Ramamohan (1982). Fast transforms: Algorithms, analyses, applications. New York, USA: Academic Press.
 Guo, Haitao; Sitton, Gary A.; Burrus, Charles Sidney (1994). The Quick Discrete Fourier Transform. Proceedings on the IEEE Conference on Acoustics, Speech, and Signal Processing (ICASSP). 3. pp. 445–448. doi:10.1109/ICASSP.1994.389994. ISBN 9780780317758.
 Johnson, Steven G.; Frigo, Matteo (2007). "A modified splitradix FFT with fewer arithmetic operations" (PDF). IEEE Transactions on Signal Processing. 55 (1): 111–119. CiteSeerX 10.1.1.582.5497. doi:10.1109/tsp.2006.882087.
 Press, William H.; Teukolsky, Saul A.; Vetterling, William T.; Flannery, Brian P. (2007). "Chapter 12. Fast Fourier Transform". Numerical Recipes: The Art of Scientific Computing (3 ed.). New York, USA: Cambridge University Press. ISBN 9780521880688.
 Singleton, Richard Collom (June 1969). "A Short Bibliography on the Fast Fourier Transform". Special Issue on Fast Fourier Transform. IEEE Transactions on Audio and Electroacoustics. AU17. IEEE Audio and Electroacoustics Group. pp. 166–169. Retrieved 20171031. (NB. Contains extensive bibliography.)
External links
 Fast Fourier Algorithm
 Fast Fourier Transforms, Connexions online book edited by Charles Sidney Burrus, with chapters by Charles Sidney Burrus, Ivan Selesnick, Markus Pueschel, Matteo Frigo, and Steven G. Johnson (2008).
 Links to FFT code and information online.
 National Taiwan University – FFT
 FFT programming in C++ — Cooley–Tukey algorithm.
 Online documentation, links, book, and code.
 Using FFT to construct aggregate probability distributions
 Sri Welaratna, "Thirty years of FFT analyzers", Sound and Vibration (January 1997, 30th anniversary issue). A historical review of hardware FFT devices.
 FFT Basics and Case Study Using MultiInstrument
 FFT Textbook notes, PPTs, Videos at Holistic Numerical Methods Institute.
 ALGLIB FFT Code GPL Licensed multilanguage (VBA, C++, Pascal, etc.) numerical analysis and data processing library.
 MIT's sFFT MIT Sparse FFT algorithm and implementation.
 VB6 FFT VB6 optimized library implementation with source code.
 Fast Fourier transform illustrated Demo examples and FFT calculator.
 Discrete Fourier Transform (Forward) A JavaScript implementation of FFTPACK code by Swarztrauber.
 Fourier Transforms of Discrete Signals (Microlink IT College)
 Interactive FFT Tutorial A visual interactive intro to Fourier transforms and FFT methods.