FFT C code

Colored By Color Scripter™
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
/*
   Modification of Paul Bourkes FFT code by Peter Cusack 
   to utilise the Microsoft complex type.
 
   This computes an in-place complex-to-complex FFT 
   x and y are the real and imaginary arrays of 2^m points.
   dir =  1 gives forward transform
   dir = -1 gives reverse transform 
*/
void FFT(int dir, long m, complex <double> x[])
{
   long i, i1, i2,j, k, l, l1, l2, n;
   complex <double> tx, t1, u, c;
 
   /*Calculate the number of points */
   n = 1;
   for(i = 0; i < m; i++) 
      n <<= 1;   
 
   /* Do the bit reversal */
   i2 = n >> 1;
   j = 0;
 
   for (i = 0; i < n-1 ; i++)
   {
      if (i < j)
         swap(x[i], x[j]);
 
      k = i2;
 
      while (k <= j) 
      {
         j -= k;
         k >>= 1;
      }
 
      j += k;
   }
 
   /* Compute the FFT */
   c.real(-1.0);
   c.imag(0.0);
   l2 = 1;
   for (l = 0; l < m; l++) 
   {
      l1 = l2;
      l2 <<= 1;
      u.real(1.0);
      u.imag(0.0);
 
      for (j = 0; j < l1; j++) 
      {
         for (i = j; i < n; i += l2) 
         {
            i1 = i + l1;
            t1 = u * x[i1];
            x[i1] = x[i] - t1; 
            x[i] += t1;
         }
 
         u = u * c;
      }
 
      c.imag(sqrt((1.0 - c.real()) / 2.0));
      if (dir == 1)
         c.imag(-c.imag());
      c.real(sqrt((1.0 + c.real()) / 2.0));
   }
 
   /* Scaling for forward transform */
   if (dir == 1) 
   {
      for (i = 0; i < n; i++)
         x[i] /= n;      
   }   
   return;
}

출처 - http://paulbourke.net/miscellaneous/dft/