DSP Tricks: Computing Fast Fourier Transform Twiddle Factors
Typical applications using an N-point radix-2 FFT accept N input time
samples, x(n), and compute N frequency-domain samples X(m).
However, there are non-standard FFT applications (for example, specialized harmonic analysis, or perhaps using an FFT to implement a bank of filters) where only a subset of the X(m) results are required.
Consider Figure 13"39 below showing a 16-point radix-2 decimation in time FFT, and assume we only need to compute the X(1), X(5), X(9) and X(13) output samples. In this case, rather than compute the entire FFT we only need perform the computations indicated by the heavy lines.
|Figure 13"39 16-point radix-2 FFT signal flow diagram.|
Reduced-computation FFTs are often called pruned FFTs. (for simplicity the butterflies in Figure 13"39 above only show the twiddle phase angle factors and not the entire twiddle phase angles.)
To implement pruned FFTs, we need to know the twiddle phase angles associated with each necessary butterfly computation as shown in Figure 13"40(a) below.
Here we give an interesting algorithm for computing the (2 pi x A1)/N and (2 pi x A2)/N twiddle phase angles for an arbitrary-size FFT. The algorithm draws upon the following characteristics of a radix-2 decimation-in-time FFT:
A general FFT butterfly is that shown in Figure 13-40(a).
The A1 and A2 angle factors are the integer values shown in Figure 13-39.
An N-point radix-2 FFT has M stages (shown at the top of Figure 13"39) where M = log2(N).
Each stage comprises N/2 butterflies.
The A1 phase angle factor of an arbitrary butterfly is
where S is the index of the M stages over the range 1 S M. Similarly, B serves as the index for the N/2 butterflies in each stage, where 1 B N/2.
|Figure 13"40 Two forms of a radix-2 butterfly.|
B = 1 means the top butterfly within a stage. The [q] operation is a function that returns the smallest integer q.
Brev[z] represents the three-step operation of: convert decimal integer z to a binary number represented by M"1 binary bits, perform bit reversal on the binary number, and convert the bit reversed number back to a decimal integer yielding angle factor A1. Angle factor A2, in Figure 13"40(a) above, is then computed using
The algorithm can also be used to compute the single twiddle angle factor of the optimized butterfly shown in Figure 13"40(b). Below is a code listing, in MATLAB, implementing our twiddle angle factor computational algorithm.
The variables in the code with start and stop in their names allow computation of a subset of the of all the twiddle angle factors in an N-point FFT.
For example, to compute the twiddle angle factors for the fifth and sixth butterflies in the third stage of a 32-point FFT, we can assign N = 32, Sstart = 3, Sstop = 3, Bstart = 5, and Bstop = 6, and run the code.
Used with the permission of the publisher, Prentice Hall, this on-going series of articles on Embedded.com is based on copyrighted material from "Understanding Digital Signal Processing, Second Edition" by Richard G. Lyons. The book can be purchased on line.
Richard Lyons is a consulting systems engineer and lecturer with Besser Associates. As a lecturer with Besser and an instructor for the University of California Santa Cruz Extension, Lyons has delivered digitasl signal processing seminars and training course at technical conferences as well at companies such as Motorola, Freescale, Lockheed Martin, Texas Instruments, Conexant, Northrop Grumman, Lucent, Nokia, Qualcomm, Honeywell, National Semiconductor, General Dynamics and Infinion.