Performance/bin/FFT.c

0001 #include <stdio.h>
0002 #include <stdlib.h>
0003 #include <math.h>
0004
0005 #include "FFT.h"
0006
0007 #define PI 3.1415926535897932
0008
0009 /*-----------------------------------------------------------------------*/
0010
0011 static int int_log2(int n);
0012
0013 double FFT_num_flops(int N) {
0014   double Nd = (double)N;
0015   double logN = (double)int_log2(N);
0016
0017   return (5.0 * Nd - 2) * logN + 2 * (Nd + 1);
0018 }
0019
0020 static int int_log2(int n) {
0021   int k = 1;
0022   int log = 0;
0023   for (/*k=1*/; k < n; k *= 2, log++)
0024     ;
0025   if (n != (1 << log)) {
0026     printf("FFT: Data length is not a power of 2!: %d ", n);
0027     exit(1);
0028   }
0029   return log;
0030 }
0031
0032 static void FFT_transform_internal(int N, double *data, int direction) {
0033   int n = N / 2;
0034   int bit = 0;
0035   int logn;
0036   int dual = 1;
0037
0038   if (n == 1)
0039     return; /* Identity operation! */
0040   logn = int_log2(n);
0041
0042   if (N == 0)
0043     return;
0044
0045   /* bit reverse the input data for decimation in time algorithm */
0046   FFT_bitreverse(N, data);
0047
0048   /* apply fft recursion */
0049   /* this loop executed int_log2(N) times */
0050   for (bit = 0; bit < logn; bit++, dual *= 2) {
0051     double w_real = 1.0;
0052     double w_imag = 0.0;
0053     int a;
0054     int b;
0055
0056     double theta = 2.0 * direction * PI / (2.0 * (double)dual);
0057     double s = sin(theta);
0058     double t = sin(theta / 2.0);
0059     double s2 = 2.0 * t * t;
0060
0061     for (a = 0, b = 0; b < n; b += 2 * dual) {
0062       int i = 2 * b;
0063       int j = 2 * (b + dual);
0064
0065       double wd_real = data[j];
0066       double wd_imag = data[j + 1];
0067
0068       data[j] = data[i] - wd_real;
0069       data[j + 1] = data[i + 1] - wd_imag;
0070       data[i] += wd_real;
0071       data[i + 1] += wd_imag;
0072     }
0073
0074     /* a = 1 .. (dual-1) */
0075     for (a = 1; a < dual; a++) {
0076       /* trignometric recurrence for w-> exp(i theta) w */
0077       {
0078         double tmp_real = w_real - s * w_imag - s2 * w_real;
0079         double tmp_imag = w_imag + s * w_real - s2 * w_imag;
0080         w_real = tmp_real;
0081         w_imag = tmp_imag;
0082       }
0083       for (b = 0; b < n; b += 2 * dual) {
0084         int i = 2 * (b + a);
0085         int j = 2 * (b + a + dual);
0086
0087         double z1_real = data[j];
0088         double z1_imag = data[j + 1];
0089
0090         double wd_real = w_real * z1_real - w_imag * z1_imag;
0091         double wd_imag = w_real * z1_imag + w_imag * z1_real;
0092
0093         data[j] = data[i] - wd_real;
0094         data[j + 1] = data[i + 1] - wd_imag;
0095         data[i] += wd_real;
0096         data[i + 1] += wd_imag;
0097       }
0098     }
0099   }
0100 }
0101
0102 void FFT_bitreverse(int N, double *data) {
0103   /* This is the Goldrader bit-reversal algorithm */
0104   int n = N / 2;
0105   int nm1 = n - 1;
0106   int i = 0;
0107   int j = 0;
0108   for (; i < nm1; i++) {
0109     /*int ii = 2*i; */
0110     int ii = i << 1;
0111
0112     /*int jj = 2*j; */
0113     int jj = j << 1;
0114
0115     /* int k = n / 2 ; */
0116     int k = n >> 1;
0117
0118     if (i < j) {
0119       double tmp_real = data[ii];
0120       double tmp_imag = data[ii + 1];
0121       data[ii] = data[jj];
0122       data[ii + 1] = data[jj + 1];
0123       data[jj] = tmp_real;
0124       data[jj + 1] = tmp_imag;
0125     }
0126
0127     while (k <= j) {
0128       /*j = j - k ; */
0129       j -= k;
0130
0131       /*k = k / 2 ;  */
0132       k >>= 1;
0133     }
0134     j += k;
0135   }
0136 }
0137
0138 void FFT_transform(int N, double *data) { FFT_transform_internal(N, data, -1); }
0139
0140 void FFT_inverse(int N, double *data) {
0141   int n = N / 2;
0142   double norm = 0.0;
0143   int i = 0;
0144   FFT_transform_internal(N, data, +1);
0145
0146   /* Normalize */
0147
0148   norm = 1 / ((double)n);
0149   for (i = 0; i < N; i++)
0150     data[i] *= norm;
0151 }