1 /*
2 * fft.c
3 * Copyright 2011 John Lindgren
4 *
5 * Redistribution and use in source and binary forms, with or without
6 * modification, are permitted provided that the following conditions are met:
7 *
8 * 1. Redistributions of source code must retain the above copyright notice,
9 * this list of conditions, and the following disclaimer.
10 *
11 * 2. Redistributions in binary form must reproduce the above copyright notice,
12 * this list of conditions, and the following disclaimer in the documentation
13 * provided with the distribution.
14 *
15 * This software is provided "as is" and without any warranty, express or
16 * implied. In no event shall the authors be liable for any damages arising from
17 * the use of this software.
18 */
19
20 #include "internal.h"
21
22 #include <complex>
23 #include <math.h>
24
25 #define TWO_PI 6.2831853f
26
27 #define N 512 /* size of the DFT */
28 #define LOGN 9 /* log N (base 2) */
29
30 typedef std::complex<float> Complex;
31
32 static float hamming[N]; /* hamming window, scaled to sum to 1 */
33 static int reversed[N]; /* bit-reversal table */
34 static Complex roots[N / 2]; /* N-th roots of unity */
35 static char generated = 0; /* set if tables have been generated */
36
37 /* Reverse the order of the lowest LOGN bits in an integer. */
38
bit_reverse(int x)39 static int bit_reverse(int x)
40 {
41 int y = 0;
42
43 for (int n = LOGN; n--;)
44 {
45 y = (y << 1) | (x & 1);
46 x >>= 1;
47 }
48
49 return y;
50 }
51
52 /* Generate lookup tables. */
53
generate_tables()54 static void generate_tables()
55 {
56 if (generated)
57 return;
58
59 for (int n = 0; n < N; n++)
60 hamming[n] = 1 - 0.85f * cosf(n * (TWO_PI / N));
61 for (int n = 0; n < N; n++)
62 reversed[n] = bit_reverse(n);
63 for (int n = 0; n < N / 2; n++)
64 roots[n] = exp(Complex(0, n * (TWO_PI / N)));
65
66 generated = 1;
67 }
68
69 /* Perform the DFT using the Cooley-Tukey algorithm. At each step s, where
70 * s=1..log N (base 2), there are N/(2^s) groups of intertwined butterfly
71 * operations. Each group contains (2^s)/2 butterflies, and each butterfly has
72 * a span of (2^s)/2. The twiddle factors are nth roots of unity where n = 2^s.
73 */
74
do_fft(Complex a[N])75 static void do_fft(Complex a[N])
76 {
77 int half = 1; /* (2^s)/2 */
78 int inv = N / 2; /* N/(2^s) */
79
80 /* loop through steps */
81 while (inv)
82 {
83 /* loop through groups */
84 for (int g = 0; g < N; g += half << 1)
85 {
86 /* loop through butterflies */
87 for (int b = 0, r = 0; b < half; b++, r += inv)
88 {
89 Complex even = a[g + b];
90 Complex odd = roots[r] * a[g + half + b];
91 a[g + b] = even + odd;
92 a[g + half + b] = even - odd;
93 }
94 }
95
96 half <<= 1;
97 inv >>= 1;
98 }
99 }
100
101 /* Input is N=512 PCM samples.
102 * Output is intensity of frequencies from 1 to N/2=256. */
103
calc_freq(const float data[N],float freq[N/2])104 void calc_freq(const float data[N], float freq[N / 2])
105 {
106 generate_tables();
107
108 /* input is filtered by a Hamming window */
109 /* input values are in bit-reversed order */
110 Complex a[N];
111 for (int n = 0; n < N; n++)
112 a[reversed[n]] = data[n] * hamming[n];
113
114 do_fft(a);
115
116 /* output values are divided by N */
117 /* frequencies from 1 to N/2-1 are doubled */
118 for (int n = 0; n < N / 2 - 1; n++)
119 freq[n] = 2 * abs(a[1 + n]) / N;
120
121 /* frequency N/2 is not doubled */
122 freq[N / 2 - 1] = abs(a[N / 2]) / N;
123 }
124