Audacious  $Id:Doxyfile42802007-03-2104:39:00Znenolod$
fft.c
Go to the documentation of this file.
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 }