Audacious $Id:Doxyfile42802007-03-2104:39:00Znenolod$
|
00001 /* 00002 * fft.c 00003 * Copyright 2011 John Lindgren 00004 * 00005 * This file is part of Audacious. 00006 * 00007 * Audacious is free software: you can redistribute it and/or modify it under 00008 * the terms of the GNU General Public License as published by the Free Software 00009 * Foundation, version 2 or version 3 of the License. 00010 * 00011 * Audacious is distributed in the hope that it will be useful, but WITHOUT ANY 00012 * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR 00013 * A PARTICULAR PURPOSE. See the GNU General Public License for more details. 00014 * 00015 * You should have received a copy of the GNU General Public License along with 00016 * Audacious. If not, see <http://www.gnu.org/licenses/>. 00017 * 00018 * The Audacious team does not consider modular code linking to Audacious or 00019 * using our public API to be a derived work. 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 /* e^(a+bi) = (e^a)(cos(b)+sin(b)i) */ 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 /* size of the DFT */ 00035 #define LOGN 9 /* log N (base 2) */ 00036 00037 static float hamming[N]; /* hamming window, scaled to sum to 1 */ 00038 static int reversed[N]; /* bit-reversal table */ 00039 static float complex roots[N / 2]; /* N-th roots of unity */ 00040 static char generated = 0; /* set if tables have been generated */ 00041 00042 /* Reverse the order of the lowest LOGN bits in an integer. */ 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 /* Generate lookup tables. */ 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 /* Perform the DFT using the Cooley-Tukey algorithm. At each step s, where 00075 * s=1..log N (base 2), there are N/(2^s) groups of intertwined butterfly 00076 * operations. Each group contains (2^s)/2 butterflies, and each butterfly has 00077 * a span of (2^s)/2. The twiddle factors are nth roots of unity where n = 2^s. */ 00078 00079 static void do_fft (float complex a[N]) 00080 { 00081 int half = 1; /* (2^s)/2 */ 00082 int inv = N / 2; /* N/(2^s) */ 00083 00084 /* loop through steps */ 00085 while (inv) 00086 { 00087 /* loop through groups */ 00088 for (int g = 0; g < N; g += half << 1) 00089 { 00090 /* loop through butterflies */ 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 /* Input is N=512 PCM samples. 00106 * Output is intensity of frequencies from 1 to N/2=256. */ 00107 00108 void calc_freq (const float data[N], float freq[N / 2]) 00109 { 00110 generate_tables (); 00111 00112 /* input is filtered by a Hamming window */ 00113 /* input values are in bit-reversed order */ 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 /* output values are divided by N */ 00121 /* frequencies from 1 to N/2-1 are doubled */ 00122 for (int n = 0; n < N / 2 - 1; n ++) 00123 freq[n] = 2 * cabsf (a[1 + n]) / N; 00124 00125 /* frequency N/2 is not doubled */ 00126 freq[N / 2 - 1] = cabsf (a[N / 2]) / N; 00127 }