Decimation in time

The radix-2 decimation-in-time algorithm rearranges the discrete Fourier transform (DFT) equation into two parts: a sum over the even-numbered discrete-time indices n=024…N-2 n 0 2 4 … N 2 and a sum over the odd-numbered indices n=135…N-1 n 1 3 5 … N 1 as in Equation 1:

The mathematical simplifications in Equation 1 reveal that all DFT frequency outputs Xk X k can be computed as the sum of the outputs of two length- N2 N 2 DFTs, of the even-indexed and odd-indexed discrete-time samples, respectively, where the odd-indexed short DFT is multiplied by a so-called twiddle factor term W N k =ⅇ-ⅈ2πkN W N k 2 k N . This is called a decimation in time because the time samples are rearranged in alternating groups, and a radix-2 algorithm because there are two groups. Figure 1 graphically illustrates this form of the DFT computation, where for convenience the frequency outputs of the length- N2 N 2 DFT of the even-indexed time samples are denoted Gk G k and those of the odd-indexed samples as Hk H k . Because of the periodicity with N2 N 2 frequency samples of these length- N2 N 2 DFTs, Gk G k and Hk H k can be used to compute two of the length-NN DFT frequencies, namely Xk X k and Xk+N2 X k N 2 , but with a different twiddle factor. This reuse of these short-length DFT outputs gives the FFT its computational savings.

Figure 1: Decimation in time of a length-NN DFT into two length-N2N2 DFTs followed by a combining stage.

Whereas direct computation of all NN DFT frequencies according to the DFT equation would require N2 N 2 complex multiplies and N2-N N 2 N complex additions (for complex-valued data), by reusing the results of the two short-length DFTs as illustrated in Figure 1, the computational cost is now

This simple reorganization and reuse has reduced the total computation by almost a factor of two over direct DFT computation!

Additional Simplification

A basic butterfly operation is shown in Figure 2, which requires only N2 N 2 twiddle-factor multiplies per stage. It is worthwhile to note that, after merging the twiddle factors to a single term on the lower branch, the remaining butterfly is actually a length-2 DFT! The theory of multi-dimensional index maps shows that this must be the case, and that FFTs of any factorable length may consist of successive stages of shorter-length FFTs with twiddle-factor multiplications in between.

Figure 2: Radix-2 DIT butterfly simplification: both operations produce the same outputs

Subfigure 2.1 Subfigure 2.2

Radix-2 decimation-in-time FFT

The same radix-2 decimation in time can be applied recursively to the two length N2 N 2 DFTs to save computation. When successively applied until the shorter and shorter DFTs reach length-2, the result is the radix-2 DIT FFT algorithm.

Figure 3: Radix-2 Decimation-in-Time FFT algorithm for a length-8 signal

The full radix-2 decimation-in-time decomposition illustrated in Figure 3 using the simplified butterflies involves M=log2N M 2 N stages, each with N2 N 2 butterflies per stage. Each butterfly requires 11 complex multiply and 22 adds per butterfly. The total cost of the algorithm is thus

This is a remarkable savings over direct computation of the DFT. For example, a length-1024 DFT would require 10485761048576 complex multiplications and 10475521047552 complex additions with direct computation, but only 51205120 complex multiplications and 1024010240 complex additions using the radix-2 FFT, a savings by a factor of 100 or more. The relative savings increase with longer FFT lengths, and are less for shorter lengths.

Modest additional reductions in computation can be achieved by noting that certain twiddle factors, namely Using special butterflies for W N 0 W N 0 , W N N 2 W N N 2 , W N N 4 W N N 4 , W N N 8 W N N 8 , W N 3 N 8 W N 3 N 8 , require no multiplications, or fewer real multiplies than other ones. By implementing special butterflies for these twiddle factors as discussed in FFT algorithm and programming tricks, the computational cost of the radix-2 decimation-in-time FFT can be reduced to

It is important to note that, if the input signal data are placed in bit-reversed order before beginning the FFT computations, the outputs of each butterfly throughout the computation can be placed in the same memory locations from which the inputs were fetched, resulting in an in-place algorithm that requires no extra memory to perform the FFT. Most FFT implementations are in-place, and overwrite the input data with the intermediate values and finally the output.

Example FFT Code

The following function, written in the C programming language, implements a radix-2 decimation-in-time FFT. It is designed for computing the DFT of complex-valued inputs to produce complex-valued outputs, with the real and imaginary parts of each number stored in separate double-precision floating-point arrays. It is an in-place algorithm, so the intermediate and final output values are stored in the same array as the input data, which is overwritten. After initializations, the program first bit-reverses the discrete-time samples, as is typical with a decimation-in-time algorithm (but see alternate FFT structures for DIT algorithms with other input orders), then computes the FFT in stages according to the above description.

Ihis FFT program uses a standard three-loop structure for the main FFT computation. The outer loop steps through the stages (each column in Figure 3); the middle loop steps through "flights" (butterflies with the same twiddle factor from each short-length DFT at each stage), and the inner loop steps through the individual butterflies. This ordering minimizes the number of fetches or computations of the twiddle-factor values. Since the bit-reverse of a bit-reversed index is the original index, bit-reversal can be performed fairly simply by swapping pairs of data.

/**********************************************************/
/* fft.c                                                  */
/* (c) Douglas L. Jones                                   */
/* University of Illinois at Urbana-Champaign             */
/* January 19, 1992                                       */
/*                                                        */
/*   fft: in-place radix-2 DIT DFT of a complex input     */
/*                                                        */
/*   input:                                               */
/* n: length of FFT: must be a power of two               */
/* m: n = 2**m                                            */
/*   input/output                                         */
/* x: double array of length n with real part of data     */
/* y: double array of length n with imag part of data     */
/*                                                        */
/*   Permission to copy and use this program is granted   */
/*   under a Creative Commons "Attribution" license       */
/*   http://creativecommons.org/licenses/by/1.0/          */
/**********************************************************/
fft(n,m,x,y)
int n,m;
double x[],y[];
{
int i,j,k,n1,n2;
double c,s,e,a,t1,t2;        
         
  
j = 0; /* bit-reverse */
n2 = n/2;
for (i=1; i < n - 1; i++)
{
  n1 = n2;
  while ( j >= n1 )
   {
    j = j - n1;
    n1 = n1/2;
   }
  j = j + n1;
               
  if (i < j)
   {
    t1 = x[i];
    x[i] = x[j];
    x[j] = t1;
    t1 = y[i];
    y[i] = y[j];
    y[j] = t1;
   }
}
                                       
                                           
n1 = 0; /* FFT */
n2 = 1;
                                             
for (i=0; i < m; i++)
{
  n1 = n2;
  n2 = n2 + n2;
  e = -6.283185307179586/n2;
  a = 0.0;
                                             
  for (j=0; j < n1; j++)
   {
    c = cos(a);
    s = sin(a);
    a = a + e;
                                            
    for (k=j; k < n; k=k+n2)
     {
      t1 = c*x[k+n1] - s*y[k+n1];
      t2 = s*x[k+n1] + c*y[k+n1];
      x[k+n1] = x[k] - t1;
      y[k+n1] = y[k] - t2;
      x[k] = x[k] + t1;
      y[k] = y[k] + t2;
     }
   }
}
                                      
return;
}