The Goertzel Algorithm and the Chirp Transform

Massachusetts Institute of Technology Department of Electrical Engineering and Computer Science 6.341: Discrete-Time Signal Processing OpenCourseWare 2006 Lecture 20 The Goertzel Algorithm and the Chirp Transform Reading: Sections 9.0 - 9.2 and 9.6 in Oppenheim, Schafer & Buck (OSB). In the previous lecture we discussed a well-known class of algorithms for computing the DFT efficiently. While the fast Fourier transform’s various incarnations have gained considerable popularity, careful selection of an appropriate algorithm for computing the DFT in practice need not be limited to choosing between these so-called “fast” implementations. We’ll there­ fore focus in this lecture on two other techniques for computing the DFT: the Goertzel algorithm and the chirp transform. Before going further, think about the following question: Given a signal x[n] with corresponding 8-point DFT X[k], what is the most efficient way, in terms of total multiplications required, to compute X[5] from x[n]? To start, note that using a radix-2 FFT algorithm isn’t the best choice; judicious application of the canonical DFT equation is enough to beat it. Computing X[5] requires roughly 8 log 2 8 = 24 complex multiplications when employing a radix-2 FFT algorithm, since X[k] for k = 0, 1, . . . , 7 must first be computed. However, calculating X[5] = 7n=0 x[n]e−j(2�/8)5n using OSB Equation 8.11 requires only 8 complex multiplications and is strikingly simple compared to the radix-2 FFT algorithms, depicted in OSB Figure 9.10. This example illustrates the point that while the FFT algorithms may be useful for a wide class of problems, they are by no means the most efficient choice in all situations. 1 The Goertzel Algorithm We’ll now discuss the Goertzel Algorithm, an efficient method (in terms of multiplications) for computing X[k] for a given k. The derivation of the algorithm, which is developed in OSB Section 9.2, begins by noting that the DFT can be formulated in terms of a convolution. X[k] = = N −1 � n=0 N −1 � x[n]WNnk , −k(N −n) x[n]WN n=0 1 2� WN = e−j N nk , WN−kN = 1 = N −1 � = � −k(N −r) x[r]WN r=0 = �� � x[n] � WN−nk � n=N � � ��� � −nk x[n] � WN u[n] � n=N Specifically, processing a signal x[n] through an LTI filter with impulse response h[n] = W N−nk u[n] and evaluating the result, yk [n], at n = N will give the corresponding N -point DFT coefficient X[k] = yk [N ]. This LTI filtering process is illustrated below. � � x[n] −� h[n] = WN−nk u[n] −� yk [n] = x[n] � WN−nk u[n] Representing the filter by its z-transform, we have x[n] −� � m=0 WN−mk z −m −� yk [n], and since H(z) = � � WN−mk z −m m=0 = � 1 − WN−k z −1 � 1 = = = WN−mk z −m − WN−k z −1 m=0 � −mk −m z m=0 WN − 1− 1 −(m+1)k −(m+1) � z m=0 WN −k −1 WN z WN0 z 0 − WN−k z −1 1 , 1 − WN−k z −1 the filtering operation can be equivalently performed by the system x[n] −� 1 1 − WN−k z −1 −� yk [n], with the associated flow graph depicted in OSB Figure 9.10. Noting that H(z) = 1 1 − WN−k z −1 2 = 1 − WNk z −1 1 k −1 1 − WN z 1 − WN−k z −1 = 1 − WNk z −1 � � , 1 − 2 cos 2�k z −1 + z −2 N the filtering operation can also be equivalently performed by x[n] −� 1 − WNk z −1 � −1 � z + z −2 1 − 2 cos 2�k N −� yk [n], with the associated flow graph depicted in OSB Figure 9.20. Since we’re only evaluating the output of this filter at y k [N ], the multiplier −WNk need only be used at time n = N . With this in mind, the algorithm requires N real multiplications and a single complex multiplication to compute X[k] for a given k. This compares favorably to N complex multiplications as required by the canonical DFT equation and approximately N log 2 N complex multiplications as required by the radix-2 FFT algorithms, when computing X[k] for a given k. A final note about the Goertzel algorithm: since it is not restricted to −1 nk computing N n=0 x[n]WN for integer values of k, the algorithm has the convenient property that it can efficiently compute X(ej� ) for arbitrary choice of �. 2 The Chirp Transform Algorithm The chirp transform algorithm, which is derived in detail in OSB Subsection 9.6.2, computes X(ej� ) for uniformly-spaced samples of � anywhere along the unit circle, as depicted in OSB Figure 9.25. The algorithm therefore computes X(ej(�0 +k��) ) = N −1 n=0 x[n]e−j(�0 +k��)n , k = 0, 1, . . . , M − 1 , or equivalently, X(ej�k ) = N −1 � x[n]e−j�k n , k = 0, 1, . . . , M − 1, n=0 where �k = �0 + k��, k = 0, 1, . . . , M − 1. Defining W = e−j�� , this becomes X(ej�k ) = = N −1 � n=0 N −1 � e−j�0 n W nk x[n]e−j�0 n W n 2 /2 Wk 2 /2 W −(k−n) 2 /2 n=0 = Wn 2 /2 �� e−j�0 n W n 3 2 /2 � � 2 x[n] � W −n /2 . By defining g[n] = e−j�0 n W n 2 /2 x[n], the above equation is equivalently represented as � � 2 2 X(ej�k ) = W n /2 g[n] � W −n /2 , which leads to the block diagram depicted in OSB Figure 9.26. This system isn’t realizable as-is, since it is neither causal nor stable. If, however, we restrict ourselves to operating on a finite-length causal signal x[n], the system can be implemented, since we’re also only interested its output y[n] over a given, finite time interval. There are several possible causal realizations of this, and in each case the implementation relies on appropriate bookkeeping. There may also be some question as to how implementing X(e j�k ) as in OSB Figure 9.26 could be practically useful. Note however that for a given ��, the analysis parameter � 0 appears only at the first modulator block. As long as the frequency sample spacing �� is fixed, the system can therefore be easily adjusted to compute for frequency samples beginning at any point on the unit circle. Another advantage of this technique is that it computes DTFT samples using a fixed convolution block, a property which is convenient when using discrete-time analog hard­ ware, such as charge-coupled devices (CCDs). 4