fft.c

   1 #define _GNU_SOURCE
   2 #include <stdio.h>
   3 #include <math.h>
   4 #include <complex.h>
   5
   6 static int rev(int i, int n)
   7 {
   8   int j = 0;
   9   n--;
  10   while (n)
  11     {
  12       j = (j << 1) | (i & 1);
  13       i >>= 1;
  14       n >>= 1;
  15     }
  16   return j;
  17 }
  18
  19 static complex omega(int k, int n)
  20 {
  21   return cos(2*M_PI*k/n) + I*sin(2*M_PI*k/n);
  22 }
  23
  24 static void fft(double *x, double *y, int n)
  25 {
  26   complex omegas[n];
  27   for (int i=0; i<n; i++)
  28     omegas[i] = omega(i, n);
  29
  30   for (int i=0; i<n; i++)
  31     y[rev(i,n)] = x[i];
  32
  33   for (int i=1; i<n; i*=2)
  34     for (int j=0; j<n; j+=2*i)
  35       for (int k=0; k<i; k++)
  36         {
  37           complex a = y[j+k];
  38           complex b = y[j+k+i];
  39           complex o = omegas[(n/(2*i))*k % n];
  40           y[j+k] = a + b*o;
  41           y[j+k+i] = a - b*o;
  42         }
  43 }
  44
  45 int main(void)
  46 {
  47   double x[8] = { 0, 1, 0, 1, 0, 1, 0, 1 };
  48   double y[8];
  49   fft(x, y, 8);
  50   for (int i=0; i<8; i++)
  51     printf("%f + %fi\n", creal(y[i]), cimag(y[i]));
  52   return 0;
  53 }