Chirp Z-transform

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

$f_{j}=\sum _{k=0}^{n-1}x_{k}e^{-{\frac {2\pi i}{n}}jk}\qquad j=0,\dots ,n-1.$

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

$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}}\qquad j=0,\dots ,n-1.$

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

$a_{k}=x_{k}e^{-{\frac {\pi i}{n}}k^{2}}$

$b_{k}=e^{{\frac {\pi i}{n}}k^{2}},$

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, and in fact they cannot be since the convolution is linear and not circular. Rather, a linear convolution can always be computed exactly from FFTs 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, albeit several times slower than the Cooley-Tukey algorithm for composite sizes.

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:

$f_{j}=\sum _{k=0}^{n-1}x_{k}e^{i\alpha jk}\qquad j=0,\dots ,m-1,$

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.]