Bluestein's FFT algorithm
The Equation
The discrete fourier transform is defined as
\[X_k = \sum_{n=0}^{N-1} x_n e^{\frac{-i 2 \pi n k}{N}} \]using the identity \(nk = -(k-n)^2/2 + n^2/2 + k^2/2\) this becomes
\[X_k = e^{-\frac{\pi i}{N}k^2}\sum_{n=0}^{N-1} (x_n e^{-\frac{\pi i}{N}n^2})e^{\frac{\pi i}{N}(k-n)^2} = b(k)^* \sum_{n=0}^{N-1} a(n)b(k-n) \]which expresses the transform as a convolution. The symmetry in b
\[b(-k) = b(k)\]allows the simple transformation into a circular convolution of size M
\[(a * b_M)_n = \sum_{m=0}^{N-1} a_m \cdot b_{n-m \textrm{ mod } M} = F^{-1}(A \cdot B)_n\]where a and b are the same for 0 - N-1, b(M-k) = b(-k) = b(k), for 0 < k < N, 0 otherwise, and a(k) = 0 for k > N-1. For convenience M is calculated as the next largest power of 2 using this function
unsigned int nlpo2(register unsigned int x) {
x |= (x >> 1);
x |= (x >> 2);
x |= (x >> 4);
x |= (x >> 8);
x |= (x >> 16);
return(x+1);
}
The order of operations for the transform is: padding, DFT, pointwise multiplication of spectra, iDFT, pointwise multiplication with b(n)*. The DFT, pointwise mult and iDFT are wrapped up in another function fft_convolve
. The output is still N values, so after the iDFT when we have the result of the convolution, we only care about the first N values and ignore the rest. We'll use the Cooley-Tukey radix 2 to perform the DFTs.
void fft_bluestein(float *inx, float *iny, unsigned int N) {
unsigned int n, M = 2*N - 1; /* padded length */
M = nlpo2(M); /* next largest power of 2 */
float *ax, *ay, *bx, *by, theta, t;
/* allocate padded arrays */
ax = malloc(sizeof(float)*M);
ay = malloc(sizeof(float)*M);
bx = malloc(sizeof(float)*M);
by = malloc(sizeof(float)*M);
if(!(ax && ay && bx && by)) {
free(ax);
free(ay);
free(bx);
free(by);
return;
}
/* fill padded arrays */
for(n = 0; n < N; n++) {
theta = (M_PI * n / N) * n;
bx[n] = cos(theta);
by[n] = sin(theta);
ax[n] = inx[n] * bx[n] + iny[n] * by[n];
ay[n] = iny[n] * bx[n] - inx[n] * by[n];
}
for(n = 1; n < N; n++) { /* do the symmetry */
bx[M-n] = bx[n];
by[M-n] = by[n];
}
/* zero pad */
memset(ax + N, 0, (M-N)*sizeof(float));
memset(ay + N, 0, (M-N)*sizeof(float));
memset(bx + N, 0, (M-N-N+1)*sizeof(float));
memset(by + N, 0, (M-N-N+1)*sizeof(float));
/* convolve */
fft_convolve(ax, ay, bx, by, M);
/* pointwise multiplication with b(n) */
for(n = 0; n < N; n++) {
theta = (M_PI * n / N) * n;
inx[n] = ax[n] * cos(theta) - ay[n] * sin(theta);
iny[n] = ay[n] * cos(theta) + ax[n] * sin(theta);
}
free(ax);
free(ay);
free(bx);
free(by);
return;
}
We could have stored b in a separate array so we didn't have to recalculate it in the last step. A downside to this algorithm is the need to allocate memory for the computation. You could preallocate it, but that would only save time if you were performing multiple DFTs all with 2N - 1 less than or equal to a certain power of 2. The algorithm gives output in the correct order, but we had to bit reverse twice inside... also the error seems to grow faster, ~5e-5 for 1000, 1e-4 for 2000, 2e-4 for 4000 and 4e-4 for 8000. The implementation probably needs work, the last pointwise mult isn't as it should be, but it gives the right answer
Exercises / Questions
- The error growth of this algorithm, depending on your version of libm, may be linear because of the large values being fed to
sin
andcos
. How would you give these functions smaller values while still getting the correct outputs?
References
- http://en.wikipedia.org/wiki/Bluestein%27s_FFT_algorithm - general description of algorithm.
- http://aggregate.org/MAGIC/ - for next greatest power of 2 function, and power of 2 test.
- Lawrence R. Rabiner, Ronald W. Schafer, and Charles M. Rader, The chirp z-transform algorithm and its application, Bell Syst. Tech. J. 48, 1249-1292 (1969). Also published in: Rabiner, Shafer, and Rader, "The chirp z-transform algorithm," IEEE Trans. Audio Electroacoustics 17 (2), 86–92 (1969). - more advanced than the wikipedia version, excellent description of steps on page 1256 (8 in the pdf) with graphs showing each step on the next page.
- Uses javascript SyntaxHighlighter for syntax highlighting.
- Uses javascript MathJax for LateX rendering.
To the extent possible under law,
Craig Hughes
has waived all copyright and related or neighboring rights to
the text and source code in this example/tutorial.