Bluestein's FFT algorithm
Bluestein's FFT algorithm (1968), also called the
chirpz algorithm (1969), is a
fast Fourier transform (FFT) algorithm that computes the
discrete Fourier transform (DFT) of arbitrary sizes (including
prime sizes) by reexpressing 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 +
j^{2}/2 +
k^{2}/2, we thus obtain:
This summation is precisely a linear convolution of the two sequences
a_{k} and
b_{k} of length
n (
k = 0,...,
n1) defined by:

with the output of the convolution multiplied by
n phase factors
b_{j}^{*}. That is:
This convolution, in turn, can be performed with a pair of FFTs (plus the precomputed FFT of
b_{k}) 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 zeropadding 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
CooleyTukey algorithm in O(
n log
n) time. Thus, Bluestein's algorithm provides an O(
n log
n) way to compute primesize DFTs, albeit several times slower than the CooleyTukey algorithm for composite sizes.
The use of zeropadding for the convolution in Bluestein's algorithm deserves some additional comment. Suppose we zeropad to a length N ≥ 2n–1. This means that a_{k} is extended to an array A_{k} of length N, where A_{k} = a_{k} for 0 ≤ k < n and A_{k} = 0 otherwise—the usual meaning of "zeropadding". However, because of the b_{j–k} term in the linear convolution, both positive and negative values of k are required for b_{k} (noting that b_{–k} = b_{k}). The periodic boundaries implied by the DFT of the zeropadded array mean that –k is equivalent to N–k. Thus, b_{k} is extended to an array B_{k} of length N, where B_{0} = b_{0}, B_{k} = B_{N–k} = b_{k} for 0 < k < n, and B_{k} = 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.
In fact, Bluestein's algorithm can be used to compute more general transforms than the DFT, called chirp ztransforms (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 transferfunction 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, 218219 (1968).
 Lawrence R. Rabiner, Ronald W. Schafer, and Charles M. Rader, "The chirp ztransform algorithm and its application," Bell Syst. Tech. J. 48, 12491292 (1969).
 D. H. Bailey and P. N. Swarztrauber, "The fractional Fourier transform and applications," SIAM Review 33, 389404 (1991). [Note that this terminology for the chirp ztransform is nonstandard: a fractional Fourier transform conventionally refers to an entirely different, continuous transform.]