If we replace the product jk in the exponent by the identity jk = -(j-k)²/2 + j²/2 + k²/2, we thus obtain:

<math> f_j = e^{-\frac{\pi i}{n} j^2 } \sum_{k=0}^{n-1} \left( x_k e^{-\frac{\pi i}{n} k^2 } \right) e^{\frac{\pi i}{n} (j-k)^2 }

This summation is precisely a convolution of the two sequences a_k and b_k of length n (k = 0,...,n-1) defined by:

<math>a_k = x_k e^{-\frac{\pi i}{n} k^2 }</math>

<math>b_k = e^{\frac{\pi i}{n} k^2 },</math>

with the output of the convolution multiplied by n phase factors b_j^*.

This convolution, in turn, can be performed with a pair of FFTs (plus the pre-computed FFT of b_k) via the convolution theorem. Although this may seem circular, the key point is that these FFTs need not be of the same length n. Rather, a convolution can always be computed exactly by zero-padding it to any size greater than or equal to 2n-1. In particular, one can pad to a power of two or some other highly composite size, for which the FFT can be efficiently computed 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.

In fact, Bluestein's algorithm can be used to compute more general transforms than the DFT: any transform of the form:

<math> f_j = \sum_{k=0}^{n-1} x_k e^{i \alpha jk }

for an arbitrary complex number α. This is called a chirp-z transform or, for real α, a fractional Fourier transform. This 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).

• L. I. Bluestein, Northeast Electronics Res. and Eng. Meeting Record 10, 218-219 (1968).
• D. H. Bailey and P. N. Swarztrauber, "The fractional Fourier transform and applications," SIAM Review 33, 389-404 (1991).