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  * 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 <complex.h>
21 #include <math.h>
22 
23 #include "config.h"
24 #include "fft.h"
25 
26 #ifndef HAVE_CEXPF
27 /* e^(a+bi) = (e^a)(cos(b)+sin(b)i) */
28 #define cexpf(x) (expf(crealf(x))*(cosf(cimagf(x))+sinf(cimagf(x))*I))
29 #warning Your C library does not have cexpf(). Please update it.
30 #endif
31 
32 #define N 512 /* size of the DFT */
33 #define LOGN 9 /* log N (base 2) */
34 
35 static float hamming[N]; /* hamming window, scaled to sum to 1 */
36 static int reversed[N]; /* bit-reversal table */
37 static float complex roots[N / 2]; /* N-th roots of unity */
38 static char generated = 0; /* set if tables have been generated */
39 
40 /* Reverse the order of the lowest LOGN bits in an integer. */
41 
42 static int bit_reverse (int x)
43 {
44  int y = 0;
45 
46  for (int n = LOGN; n --; )
47  {
48  y = (y << 1) | (x & 1);
49  x >>= 1;
50  }
51 
52  return y;
53 }
54 
55 /* Generate lookup tables. */
56 
57 static void generate_tables (void)
58 {
59  if (generated)
60  return;
61 
62  for (int n = 0; n < N; n ++)
63  hamming[n] = 1 - 0.85 * cosf (2 * M_PI * n / N);
64  for (int n = 0; n < N; n ++)
65  reversed[n] = bit_reverse (n);
66  for (int n = 0; n < N / 2; n ++)
67  roots[n] = cexpf (2 * M_PI * I * n / N);
68 
69  generated = 1;
70 }
71 
72 /* Perform the DFT using the Cooley-Tukey algorithm. At each step s, where
73  * s=1..log N (base 2), there are N/(2^s) groups of intertwined butterfly
74  * operations. Each group contains (2^s)/2 butterflies, and each butterfly has
75  * a span of (2^s)/2. The twiddle factors are nth roots of unity where n = 2^s. */
76 
77 static void do_fft (float complex a[N])
78 {
79  int half = 1; /* (2^s)/2 */
80  int inv = N / 2; /* N/(2^s) */
81 
82  /* loop through steps */
83  while (inv)
84  {
85  /* loop through groups */
86  for (int g = 0; g < N; g += half << 1)
87  {
88  /* loop through butterflies */
89  for (int b = 0, r = 0; b < half; b ++, r += inv)
90  {
91  float complex even = a[g + b];
92  float complex odd = roots[r] * a[g + half + b];
93  a[g + b] = even + odd;
94  a[g + half + b] = even - odd;
95  }
96  }
97 
98  half <<= 1;
99  inv >>= 1;
100  }
101 }
102 
103 /* Input is N=512 PCM samples.
104  * Output is intensity of frequencies from 1 to N/2=256. */
105 
106 void calc_freq (const float data[N], float freq[N / 2])
107 {
108  generate_tables ();
109 
110  /* input is filtered by a Hamming window */
111  /* input values are in bit-reversed order */
112  float complex a[N];
113  for (int n = 0; n < N; n ++)
114  a[reversed[n]] = data[n] * hamming[n];
115 
116  do_fft (a);
117 
118  /* output values are divided by N */
119  /* frequencies from 1 to N/2-1 are doubled */
120  for (int n = 0; n < N / 2 - 1; n ++)
121  freq[n] = 2 * cabsf (a[1 + n]) / N;
122 
123  /* frequency N/2 is not doubled */
124  freq[N / 2 - 1] = cabsf (a[N / 2]) / N;
125 }