SSU ES 400 - Fast Fourier Transform (FFT)

Unformatted text preview:

http://www.relisoft.com/science/Physics/fft.html Fast Fourier Transform (FFT) The naive implementation of the N-point digital Fourier transform involves calculating the scalar product of the sample buffer (treated as an N-dimensional vector) with N separate basis vectors. Since each scalar product involves N multiplications and N additions, the total time is proportional to N2 (in other words, it's an O(N2) algorithm). However, it turns out that by cleverly re-arranging these operations, one can optimize the algorithm down to O(N log(N)), which for large N makes a huge difference. The optimized version of the algorithm is called the fast Fourier transform, or the FFT. Let's do some back of the envelope calculations. Suppose that we want to do a real-time Fourier transform of one channel of CD-quality sound. That's 44k samples per second. Suppose also that we have a 1k buffer that is being re-filled with data 44 times per second. To generate a 1000-point Fourier transform we would have to perform 2 million floating-point operations (1M multiplications and 1M additions). To keep up with incoming buffers, we would need at least 88M flops (floating-point operations per second). Now, if you are lucky enough to have a 100 Mflop machine, that might be fine, but consider that you'd be left with very little processing power to spare. Using the FFT, on the other hand, we would perform on the order of 2*1000*log2(1000) operations per buffer, which is more like 20,000. Which requires 880k flops--less than 1 Mflop! A hundred-fold speedup. The standard strategy to speed up an algorithm is to divide and conquer. We have to find some way to group the terms in the equation V[k] = Σn=0..N-1 WNkn v[n] Let's see what happens when we separate odd ns from even ns (from now on, let's assume that N is even): V[k] = Σn even WNkn v[n] + Σn odd WNkn v[n] = Σr=0..N/2-1 WNk(2r) v[2r] + Σr=0..N/2-1 WNk(2r+1) v[2r+1] = Σr=0..N/2-1 WNk(2r) v[2r] + Σr=0..N/2-1 WNk(2r) WNk v[2r+1] = Σr=0..N/2-1 WNk(2r) v[2r] + WNk Σr=0..N/2-1 WNk(2r) v[2r+1] = (Σr=0..N/2-1 WN/2kr v[2r]) + WNk (Σr=0..N/2-1 WN/2kr v[2r+1]) where we have used one crucial identity: WNk(2r) = e-2πi*2kr/N = e-2πi*kr/(N/2) = WN/2krNotice an interesting thing: the two sums are nothing else but N/2-point Fourier transforms of, respectively, the even subset and the odd subset of samples. Terms with k greater or equal N/2 can be reduced using another identity: WN/2m+N/2 = WN/2mWN/2N/2 = WN/2m which is true because Wmm = e-2πi = cos(-2π) + i sin(-2π)= 1. If we start with N that is a power of 2, we can apply this subdivision recursively until we get down to 2-point transforms. We can also go backwards, starting with the 2-point transform: Note W10*k =W20*k V[k] = W20*k v[0] + W21*k v[1], k=0,1 The two components are: V[0] = W20 v[0] + W20 v[1] = v[0] + W20 v[1] V[1] = W20 v[0] + W21 v[1] = v[0] + W21 v[1] We can represent the two equations for the components of the 2-point transform graphically using the, so called, butterfly Fig. Butterfly calculation Furthermore, using the divide and conquer strategy, a 4-point transform can be reduced to two 2-point transforms: one for even elements, one for odd elements. The odd one will be multiplied by W4k. Diagrammatically, this can be represented as two levels of butterflies. Notice that using the identity WN/2n = WN2n, we can always express all the multipliers as powers of the same WN (in this case we choose N=4).Fig. Diagrammatical representation of the 4-point Fourier transform calculation I encourage the reader to derive the analogous diagrammatical representation for N=8. What will become obvious is that all the butterflies have similar form: Fig. Generic butterfly graph This graph can be further simplified using this identity: WNs+N/2 = WNs WNN/2 = -WNs which is true because WNN/2 = e-2πi(N/2)/N = e-πi = cos(-π) + isin(-π) = -1 Here's the simplified butterfly:Fig. Simplified generic butterfly Using this result, we can now simplify our 4-point diagram. Fig. 4-point FFT calculation This diagram is the essence of the FFT algorithm. The main trick is that you don't calculate each component of the Fourier transform separately. That would involve unnecessary repetition of a substantial number of calculations. Instead, you do your calculations in stages. At each stage you start with N (in general complex) numbers and "butterfly" them to obtain a new set of N complex numbers. Those numbers, in turn, become the input for the next stage. The calculation of a 4-point FFT involves two stages. The input of the first stage are the 4 original samples. The output of the second stage are the 4 components of the Fourier transform. Notice that each stage involves N/2 complex multiplications (or N real multiplications), N/2 sign inversions (multiplication by -1), and N complex additions. So each stage can be done in O(N) time. The number of stages islog2N (which, since N is a power of 2, is the exponent m in N = 2m). Altogether, the FFT requires on the order of O(N logN) calculations. Moreover, the calculations can be done in-place, using a single buffer of N complex numbers. The trick is to initialize this buffer with appropriately scrambled samples. For N=4, the order of samples is v[0], v[2], v[1], v[3]. In general, according to our basic identity, we first divide the samples into two groups, even ones and odd ones. Applying this division recursively, we split these groups of samples into two groups each by selecting every other sample. For instance, the group (0, 2, 4, 6, 8, 10, ... 2N-2) will be split into (0, 4, 8, ...) and (2, 6, 10, ...), and so on. If you write these numbers in binary notation, you'll see that the first split (odd/even) is done according to the lowest bit; the second split is done according to the second lowest bit, and so on. So if we start with the sequence of, say, 8 consecutive binary numbers: 000, 001, 010, 011, 100, 101, 110, 111 we will first scramble them like this: [even] (000, 010, 100, 110), [odd] (001, 011, 101, 111) then we'll scramble the groups: ((000, 100), (010, 110)), (001, 101), (011, 111)) which gives the result: 000, 100, 010, 110, 001, 101, 011, 111 This is equivalent to sorting the numbers in bit-reversed order--if you reverse bits in each number (for instance, 110 becomes 011, and so on), you'll get a set of consecutive numbers. So this is how the FFT algorithm works (more precisely, this is the decimation-in-time in-place FFT


View Full Document

SSU ES 400 - Fast Fourier Transform (FFT)

Download Fast Fourier Transform (FFT)
Our administrator received your request to download this document. We will send you the file to your email shortly.
Loading Unlocking...
Login

Join to view Fast Fourier Transform (FFT) and access 3M+ class-specific study document.

or
We will never post anything without your permission.
Don't have an account?
Sign Up

Join to view Fast Fourier Transform (FFT) 2 2 and access 3M+ class-specific study document.

or

By creating an account you agree to our Privacy Policy and Terms Of Use

Already a member?