Page 283 -
P. 283

Chapter 7 ■ Image Restoration   257


                               of the loop) and partly mathematical techniques. There are a number of
                               very good references on the FFT [Bracewell, 1965; Brigham, 1974], but these
                               deal rather rigorously with the subject. Here, the code in Figure 7.3 will be
                               successively improved until it implements the basic FFT method.
                                 The first optimization involves moving the exponential calculation (cexp) to
                               a position outside of the inner loop. This is done by precomputing all the N
                               possible products:
                                              n − 1                 n − 1
                                                      (−2πj/n)(wk)
                                       F(w) =    f(k) e       f(k) =   pre[wk mod n] f(k)      (EQ 7.6)
                                              k = 0                 k = 0
                                 This reduces the strength of the operation within the loop from a com-
                               plex exponentiation to a complex product. For the transform computed in
                               Figure 7.2a–b, for example, the program (Figure 7.3) requires about 5 times the
                               CPU time as does the program slow2.c (Figure 7.4), which uses precomputed
                               exponentials.



                                void slowft(float *x, COMPLEX *y, int n)  /* Compute all Y values */
                                {
                                  COMPLEX tmp,z1, z2, z3, z4, pre[1024];    for (m = 0; m<=n; m++)
                                  int m, k, i, p;                         {
                                                                              cmplx (x[0], 0.0, &(y[m]));
                                /* Constant factor -2 pi */                   for (k=1; k<=n-1; k++)
                                  cmplx (0.0, atan(1.0)/n * -8.0, &z1);        {
                                  cexp (z1, &tmp);                               p = (k*m % n);
                                                                                 cmplx (x[k], 0.0, &z3);
                                /* Precompute most of the exponential */           cmult (z3, pre[p], &z4);
                                  cmplx (1.0, 0.0, &z1);/*z1=1.0; */             csum (y[m], z4 &z2);
                                  for (i=0; i<n; i++)                            y[m].real = z2.real;
                                  {                                              y[m].imag = z2.imag;
                                    cmplx(z1.real, z1.imag, &(pre[i]));        }
                                    cmult (z1, tmp, &z3);                    }
                                    cmplx (z3.real, z3.imag, &z1);      }
                                  }

                               Figure 7.4: The program slow2, obtained from slow by precomputing the complex
                               exponential entries.

                                 The next step involves the mathematical observation that the even-numbered
                               elements can be computed separately from the odd ones. In cases where n
                               is even, this will reduce the number of multiplications by half. The even
                               coefficients are found using:
                                                              n/2 − 1
                                                               
   (−2πj/n)(ak)
                                                      F(2a) =      e       Sum [k]
                                                               k = 0
                                                                                               (EQ 7.7)
                                                              (f(k) + f(k + m))
                                                    Sum[k] =
                                                                     2
   278   279   280   281   282   283   284   285   286   287   288