Audacious  $Id:Doxyfile42802007-03-2104:39:00Znenolod$
fft.c
Go to the documentation of this file.
1 /*
2  * fft.c
3  * Copyright 2011 John Lindgren
4  *
5  * This file is part of Audacious.
6  *
7  * Audacious is free software: you can redistribute it and/or modify it under
8  * the terms of the GNU General Public License as published by the Free Software
9  * Foundation, version 2 or version 3 of the License.
10  *
11  * Audacious is distributed in the hope that it will be useful, but WITHOUT ANY
12  * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
13  * A PARTICULAR PURPOSE. See the GNU General Public License for more details.
14  *
15  * You should have received a copy of the GNU General Public License along with
16  * Audacious. If not, see <http://www.gnu.org/licenses/>.
17  *
18  * The Audacious team does not consider modular code linking to Audacious or
19  * using our public API to be a derived work.
20  */
21 
22 #include <complex.h>
23 #include <math.h>
24 
25 #include "config.h"
26 #include "fft.h"
27 
28 #ifndef HAVE_CEXPF
29 /* e^(a+bi) = (e^a)(cos(b)+sin(b)i) */
30 #define cexpf(x) (expf(crealf(x))*(cosf(cimagf(x))+sinf(cimagf(x))*I))
31 #warning Your C library does not have cexpf(). Please update it.
32 #endif
33 
34 #define N 512 /* size of the DFT */
35 #define LOGN 9 /* log N (base 2) */
36 
37 static float hamming[N]; /* hamming window, scaled to sum to 1 */
38 static int reversed[N]; /* bit-reversal table */
39 static float complex roots[N / 2]; /* N-th roots of unity */
40 static char generated = 0; /* set if tables have been generated */
41 
42 /* Reverse the order of the lowest LOGN bits in an integer. */
43 
44 static int bit_reverse (int x)
45 {
46  int y = 0;
47 
48  for (int n = LOGN; n --; )
49  {
50  y = (y << 1) | (x & 1);
51  x >>= 1;
52  }
53 
54  return y;
55 }
56 
57 /* Generate lookup tables. */
58 
59 static void generate_tables (void)
60 {
61  if (generated)
62  return;
63 
64  for (int n = 0; n < N; n ++)
65  hamming[n] = 1 - 0.85 * cosf (2 * M_PI * n / N);
66  for (int n = 0; n < N; n ++)
67  reversed[n] = bit_reverse (n);
68  for (int n = 0; n < N / 2; n ++)
69  roots[n] = cexpf (2 * M_PI * I * n / N);
70 
71  generated = 1;
72 }
73 
74 /* Perform the DFT using the Cooley-Tukey algorithm. At each step s, where
75  * s=1..log N (base 2), there are N/(2^s) groups of intertwined butterfly
76  * operations. Each group contains (2^s)/2 butterflies, and each butterfly has
77  * a span of (2^s)/2. The twiddle factors are nth roots of unity where n = 2^s. */
78 
79 static void do_fft (float complex a[N])
80 {
81  int half = 1; /* (2^s)/2 */
82  int inv = N / 2; /* N/(2^s) */
83 
84  /* loop through steps */
85  while (inv)
86  {
87  /* loop through groups */
88  for (int g = 0; g < N; g += half << 1)
89  {
90  /* loop through butterflies */
91  for (int b = 0, r = 0; b < half; b ++, r += inv)
92  {
93  float complex even = a[g + b];
94  float complex odd = roots[r] * a[g + half + b];
95  a[g + b] = even + odd;
96  a[g + half + b] = even - odd;
97  }
98  }
99 
100  half <<= 1;
101  inv >>= 1;
102  }
103 }
104 
105 /* Input is N=512 PCM samples.
106  * Output is intensity of frequencies from 1 to N/2=256. */
107 
108 void calc_freq (const float data[N], float freq[N / 2])
109 {
110  generate_tables ();
111 
112  /* input is filtered by a Hamming window */
113  /* input values are in bit-reversed order */
114  float complex a[N];
115  for (int n = 0; n < N; n ++)
116  a[reversed[n]] = data[n] * hamming[n];
117 
118  do_fft (a);
119 
120  /* output values are divided by N */
121  /* frequencies from 1 to N/2-1 are doubled */
122  for (int n = 0; n < N / 2 - 1; n ++)
123  freq[n] = 2 * cabsf (a[1 + n]) / N;
124 
125  /* frequency N/2 is not doubled */
126  freq[N / 2 - 1] = cabsf (a[N / 2]) / N;
127 }