00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022 #include <complex.h>
00023 #include <math.h>
00024
00025 #include "config.h"
00026 #include "fft.h"
00027
00028 #ifndef HAVE_CEXPF
00029
00030 #define cexpf(x) (expf(crealf(x))*(cosf(cimagf(x))+sinf(cimagf(x))*I))
00031 #warning Your C library does not have cexpf(). Please update it.
00032 #endif
00033
00034 #define N 512
00035 #define LOGN 9
00036
00037 static float hamming[N];
00038 static int reversed[N];
00039 static float complex roots[N / 2];
00040 static char generated = 0;
00041
00042
00043
00044 static int bit_reverse (int x)
00045 {
00046 int y = 0;
00047
00048 for (int n = LOGN; n --; )
00049 {
00050 y = (y << 1) | (x & 1);
00051 x >>= 1;
00052 }
00053
00054 return y;
00055 }
00056
00057
00058
00059 static void generate_tables (void)
00060 {
00061 if (generated)
00062 return;
00063
00064 for (int n = 0; n < N; n ++)
00065 hamming[n] = 1 - 0.85 * cosf (2 * M_PI * n / N);
00066 for (int n = 0; n < N; n ++)
00067 reversed[n] = bit_reverse (n);
00068 for (int n = 0; n < N / 2; n ++)
00069 roots[n] = cexpf (2 * M_PI * I * n / N);
00070
00071 generated = 1;
00072 }
00073
00074
00075
00076
00077
00078
00079 static void do_fft (float complex a[N])
00080 {
00081 int half = 1;
00082 int inv = N / 2;
00083
00084
00085 while (inv)
00086 {
00087
00088 for (int g = 0; g < N; g += half << 1)
00089 {
00090
00091 for (int b = 0, r = 0; b < half; b ++, r += inv)
00092 {
00093 float complex even = a[g + b];
00094 float complex odd = roots[r] * a[g + half + b];
00095 a[g + b] = even + odd;
00096 a[g + half + b] = even - odd;
00097 }
00098 }
00099
00100 half <<= 1;
00101 inv >>= 1;
00102 }
00103 }
00104
00105
00106
00107
00108 void calc_freq (const float data[N], float freq[N / 2])
00109 {
00110 generate_tables ();
00111
00112
00113
00114 float complex a[N];
00115 for (int n = 0; n < N; n ++)
00116 a[reversed[n]] = data[n] * hamming[n];
00117
00118 do_fft (a);
00119
00120
00121
00122 for (int n = 0; n < N / 2 - 1; n ++)
00123 freq[n] = 2 * cabsf (a[1 + n]) / N;
00124
00125
00126 freq[N / 2 - 1] = cabsf (a[N / 2]) / N;
00127 }