Bluestein's FFT algorithm
Bluestein's FFT algorithm (1968), also called the
chirp-z algorithm (1969), is a
fast Fourier transform (FFT) algorithm that computes the
discrete Fourier transform (DFT) of arbitrary sizes (including
prime sizes) by re-expressing the DFT as a linear
convolution. (The other algorithm for FFTs of prime sizes,
Rader's algorithm, also works by rewriting the DFT as a convolution.)
Recall that the DFT is defined by the formula
If we replace the product
jk in the exponent by the identity
jk = -(
j-
k)
2/2 +
j2/2 +
k2/2, we thus obtain:
This summation is precisely a linear convolution of the two sequences
ak and
bk of length
n (
k = 0,...,
n-1) defined by:
-
with the output of the convolution multiplied by
n phase factors
bj*. That is:
This convolution, in turn, can be performed with a pair of FFTs (plus the pre-computed FFT of
bk) via the
convolution theorem. The key point is that these FFTs are not of the same length
n: such a linear convolution can be computed exactly from FFTs only by zero-padding it to a length greater than or equal to 2
n–1. In particular, one can pad to a power of two or some other highly composite size, for which the FFT can be efficiently performed by e.g. the
Cooley-Tukey algorithm in O(
n log
n) time. Thus, Bluestein's algorithm provides an O(
n log
n) way to compute prime-size DFTs, albeit several times slower than the Cooley-Tukey algorithm for composite sizes.
The use of zero-padding for the convolution in Bluestein's algorithm deserves some additional comment. Suppose we zero-pad to a length N ≥ 2n–1. This means that ak is extended to an array Ak of length N, where Ak = ak for 0 ≤ k < n and Ak = 0 otherwise—the usual meaning of "zero-padding". However, because of the bj–k term in the linear convolution, both positive and negative values of k are required for bk (noting that b–k = bk). The periodic boundaries implied by the DFT of the zero-padded array mean that –k is equivalent to N–k. Thus, bk is extended to an array Bk of length N, where B0 = b0, Bk = BN–k = bk for 0 < k < n, and Bk = 0 otherwise. A and B are then FFTed, multiplied pointwise, and inverse FFTed to obtain the linear convolution of a and b, according to the usual convolution theorem.
Chirp z-Transforms
In fact, Bluestein's algorithm can be used to compute more general transforms than the DFT, called chirp z-transforms (Rabiner, 1969); this is any transform of the form:
for an arbitrary complex number α, and for differing numbers n and m of inputs and outputs. Given Bluestein's algorithm, such a transform can be used, for example, to obtain a more finely spaced interpolation of some portion of the spectrum (although the frequency resolution is still limited by the total sampling time), enhance arbitrary poles in transfer-function analyses, etcetera.
References:
- Leo I. Bluestein, "A linear filtering approach to the computation of the discrete Fourier transform," Northeast Electronics Research and Engineering Meeting Record 10, 218-219 (1968).
- 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).
- D. H. Bailey and P. N. Swarztrauber, "The fractional Fourier transform and applications," SIAM Review 33, 389-404 (1991). [Note that this terminology for the chirp z-transform is nonstandard: a fractional Fourier transform conventionally refers to an entirely different, continuous transform.]