Some days ago I read a post on the comp.dsp newsgroup and, if I understood the poster's words, it seemed that the poster would benefit from knowing how to compute the twiddle factors of a radix-2 fast Fourier transform (FFT).

Then, later it occurred to me that it might be useful for this blog's readers to be aware of algorithms for computing FFT twiddle factors. So,... what follows are two algorithms showing how to compute the individual twiddle factors of an N-point decimation-in-frequency (DIF) and an N-point decimation-in-time (DIT) FFT.

The vast majority of FFT applications use (notice how I used the word "use" instead of the clumsy word "utilize") standard, pre-written, FFT software routines. However, there are non-standard FFT applications (for example, specialized harmonic analysis, transmultiplexers, or perhaps using an FFT to implement a bank of filters) where only a subset of the full N-sample complex FFT results are required. Those oddball FFT applications, sometimes called "pruned FFTs", require computation of individual FFT twiddle factors, and that's the purpose of this blog.

(If, by chance, the computation of FFT twiddle factors is of no interest to you, you might just scroll down to the "A Little History of the FFT" part of this blog.)

Before we present the two twiddle factor computation algorithms, let's understand the configuration of a single "butterfly" operation used in our radix-2 FFTs. We've all seen the signal flow drawings of FFTs with their arrays of butterfly operations. There are various ways of implementing a butterfly operation, but my favorites are the efficient single-complex-multiply butterflies shown in Figure 1. A DIF butterfly is shown in Figure 1(a), while a DIT butterfly is shown in Figure 1(b). In Figure 1 the twiddle factors are shown as e^–j2πQ/N, where variable Q is merely an integer in the range of 0 ≤ Q ≤ (N/2)–1.

To simplify this blog's follow-on figures, we'll use Figures 1(c) and 1(d) to represent the DIF and DIT butterflies. As such, Figure 1(c) is equivalent to Figure 1(a), and Figure 1(d) is equivalent to Figure 1(b).

Computing DIF Twiddle Factors

Take a look at Figure 2 showing the butterfly operations for an 8-point radix-2 DIF FFT.

For the radix-2 DIF FFT using the Figures 1(c) and 1(d) butterflies,

• The N-point DIF FFT has log₂(N) stages, numbered P = 1, 2, ..., log₂(N).
• Each stage comprises N/2 butterflies.
• Not counting the –1 twiddle factors, the Pth stage has N/2^P unique twiddle factors, numbered k = 0, 1, 2, ..., N/2^P–1 as indicated by the bold numbers above the upward-pointing arrows at the bottom of Figure 2.

Given those characteristics, the kth unique twiddle factor phase angle for the Pth stage is computed using:

          kth DIF twiddle factor angle = k•2^P/2     (1)

where 0 ≤ k ≤ N/2^P–1. For example, for the second stage (P = 2) of an N = 8-point DIF FFT, the unique Q factors are:

          k = 0,   Q = 0•2^P/2 = 0•4/2 = 0
          k = 1,   Q = 1•2^P/2 = 1•4/2 = 2.

Computing DIT Twiddle Factors

Here's an algorithm for computing the individual twiddle factor angles of a radix-2 DIT FFT. Consider Figure 3 showing the butterfly signal flow of an 8-point DIT FFT.

For the DIT FFT using the Figures 1(c) and 1(d) butterflies,

• The N-point DIT FFT has log₂(N) stages, numbered P = 1, 2, ..., log₂(N).
• Each stage comprises N/2 butterflies.
• Not counting the –1 twiddle factors, the Pth stage has N/2 twiddle factors, numbered k = 0, 1, 2, ..., N/2–1 as indicated by the upward arrows at the bottom of Figure 3.

Given those characteristics, the kth DIT twiddle Q factor for the Pth stage is computed using:

          kth DIT twiddle factor Q = [⌊k2^P/N⌋]_bit-rev     (2)

where 0 ≤ k ≤ N/2–1. The ⌊q⌋ operation means the integer part of q. The [z]_bit-rev function represents the three-step operation of:

           [1] convert decimal integer z to a binary number represented by log₂(N)–1 binary bits,

           [2] perform bit reversal on the binary number as discussed in Section 4.5, and

           [3] convert the bit reversed number back to a decimal integer.

As an example of using Eq.(2), for the second stage (P = 2) of an N = 8-point DIT FFT, the k = 3 twiddle Q factor is:

          3rd twiddle factor Q = [⌊3•2²/8⌋]_bit-rev

                    = [⌊1.5⌋]_bit-rev = [1]_bit-rev = 2.

The above [1]_bit-rev operation is: take the decimal number 1 and represent it with log₂(N)–1 = 2 bits, i.e., as 01₂. Next, reverse those bits to a binary 10₂ and convert that binary number to our desired decimal result of 2.

A Little History of the FFT

The radix-2 FFT has a very interesting history. For example, one of the driving forces behind the development of the FFT was the United State's desire to detect nuclear explosions inside the Soviet Union in the early 1960s. Also, if it hadn't been for the influence of a patent attorney, the Cooley-Tukey radix-2 FFT algorithm might well have been known as the Sande-Tukey algorithm, named after Gordon Sande and John Tukey. (That's the same Gordon Sande that occasionally posts on the comp.dsp newsgroup.) For those and other interesting FFT historical facts, see the following web sites.

Cooley and Tukey, "On the Origin of the FFT Paper", http://www.garfield.library.upenn.edu/classics1993/A1993MJ84500001.pdf

Rockmore, "The FFT - An Algorithm the Whole Family Can Use", http://www.cs.dartmouth.edu/~rockmore/cse-fft.pdf

Cipra, "The FFT: Making Technology Fly", http://compmack.googlecode.com/svn/marcoshack/calc4/The_FFT_Making_Technology_Fly.pdf

Cooley, Lewis, and Welch, "Historical Notes on the Fast Fourier Transform", http://www.signallake.com/innovation/FFTHistoryJun67.pdf