• Main Page
  • Related Pages
  • Classes
  • Files
  • File List
  • File Members

fft.c

Go to the documentation of this file.
00001 /* Audacious - Cross-platform multimedia player
00002  * Copyright (C) 2005-2007  Audacious development team
00003  *
00004  * Copyright (C) 1999 Richard Boulton <richard@tartarus.org>
00005  * Convolution stuff by Ralph Loader <suckfish@ihug.co.nz>
00006  *
00007  *  This program is free software; you can redistribute it and/or modify
00008  *  it under the terms of the GNU General Public License as published by
00009  *  the Free Software Foundation; under version 3 of the License.
00010  *
00011  *  This program is distributed in the hope that it will be useful,
00012  *  but WITHOUT ANY WARRANTY; without even the implied warranty of
00013  *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00014  *  GNU General Public License for more details.
00015  *
00016  *  You should have received a copy of the GNU General Public License
00017  *  along with this program.  If not, see <http://www.gnu.org/licenses>.
00018  *
00019  *  The Audacious team does not consider modular code linking to
00020  *  Audacious or using our public API to be a derived work.
00021  */
00022 
00023 /* fft.c: iterative implementation of a FFT */
00024 
00025 /*
00026  * TODO
00027  * Remove compiling in of FFT_BUFFER_SIZE?  (Might slow things down, but would
00028  * be nice to be able to change size at runtime.)
00029  * Finish making / checking thread-safety.
00030  * More optimisations.
00031  */
00032 
00033 #ifdef HAVE_CONFIG_H
00034 #  include "config.h"
00035 #endif
00036 
00037 #include "fft.h"
00038 
00039 #include <glib.h>
00040 #include <stdlib.h>
00041 #include <math.h>
00042 #ifndef PI
00043 #ifdef M_PI
00044 #define PI M_PI
00045 #else
00046 #define PI            3.14159265358979323846    /* pi */
00047 #endif
00048 #endif
00049 
00050 /* ########### */
00051 /* # Structs # */
00052 /* ########### */
00053 
00054 struct _struct_fft_state {
00055     /* Temporary data stores to perform FFT in. */
00056     float real[FFT_BUFFER_SIZE];
00057     float imag[FFT_BUFFER_SIZE];
00058 };
00059 
00060 /* ############################# */
00061 /* # Local function prototypes # */
00062 /* ############################# */
00063 
00064 static void fft_prepare(const sound_sample * input, float *re, float *im);
00065 static void fft_calculate(float *re, float *im);
00066 static void fft_output(const float *re, const float *im, float *output);
00067 static int reverseBits(unsigned int initial);
00068 
00069 /* #################### */
00070 /* # Global variables # */
00071 /* #################### */
00072 
00073 /* Table to speed up bit reverse copy */
00074 static unsigned int bitReverse[FFT_BUFFER_SIZE];
00075 
00076 /* The next two tables could be made to use less space in memory, since they
00077  * overlap hugely, but hey. */
00078 static float sintable[FFT_BUFFER_SIZE / 2];
00079 static float costable[FFT_BUFFER_SIZE / 2];
00080 
00081 /* ############################## */
00082 /* # Externally called routines # */
00083 /* ############################## */
00084 
00085 /* --------- */
00086 /* FFT stuff */
00087 /* --------- */
00088 
00089 /*
00090  * Initialisation routine - sets up tables and space to work in.
00091  * Returns a pointer to internal state, to be used when performing calls.
00092  * On error, returns NULL.
00093  * The pointer should be freed when it is finished with, by fft_close().
00094  */
00095 fft_state *
00096 fft_init(void)
00097 {
00098     fft_state *state;
00099     unsigned int i;
00100 
00101     state = (fft_state *) g_malloc(sizeof(fft_state));
00102     if (!state)
00103         return NULL;
00104 
00105     for (i = 0; i < FFT_BUFFER_SIZE; i++) {
00106         bitReverse[i] = reverseBits(i);
00107     }
00108     for (i = 0; i < FFT_BUFFER_SIZE / 2; i++) {
00109         float j = 2 * PI * i / FFT_BUFFER_SIZE;
00110         costable[i] = cos(j);
00111         sintable[i] = sin(j);
00112     }
00113 
00114     return state;
00115 }
00116 
00117 /*
00118  * Do all the steps of the FFT, taking as input sound data (as described in
00119  * sound.h) and returning the intensities of each frequency as floats in the
00120  * range 0 to ((FFT_BUFFER_SIZE / 2) * 32768) ^ 2
00121  *
00122  * FIXME - the above range assumes no frequencies present have an amplitude
00123  * larger than that of the sample variation.  But this is false: we could have
00124  * a wave such that its maximums are always between samples, and it's just
00125  * inside the representable range at the places samples get taken.
00126  * Question: what _is_ the maximum value possible.  Twice that value?  Root
00127  * two times that value?  Hmmm.  Think it depends on the frequency, too.
00128  *
00129  * The input array is assumed to have FFT_BUFFER_SIZE elements,
00130  * and the output array is assumed to have (FFT_BUFFER_SIZE / 2 + 1) elements.
00131  * state is a (non-NULL) pointer returned by fft_init.
00132  */
00133 void
00134 fft_perform(const sound_sample * input, float *output, fft_state * state)
00135 {
00136     /* Convert data from sound format to be ready for FFT */
00137     fft_prepare(input, state->real, state->imag);
00138 
00139     /* Do the actual FFT */
00140     fft_calculate(state->real, state->imag);
00141 
00142     /* Convert the FFT output into intensities */
00143     fft_output(state->real, state->imag, output);
00144 }
00145 
00146 /*
00147  * Free the state.
00148  */
00149 void
00150 fft_close(fft_state * state)
00151 {
00152     if (state)
00153         free(state);
00154 }
00155 
00156 /* ########################### */
00157 /* # Locally called routines # */
00158 /* ########################### */
00159 
00160 /*
00161  * Prepare data to perform an FFT on
00162  */
00163 static void
00164 fft_prepare(const sound_sample * input, float *re, float *im)
00165 {
00166     unsigned int i;
00167     float *realptr = re;
00168     float *imagptr = im;
00169 
00170     /* Get input, in reverse bit order */
00171     for (i = 0; i < FFT_BUFFER_SIZE; i++) {
00172         *realptr++ = input[bitReverse[i]];
00173         *imagptr++ = 0;
00174     }
00175 }
00176 
00177 /*
00178  * Take result of an FFT and calculate the intensities of each frequency
00179  * Note: only produces half as many data points as the input had.
00180  * This is roughly a consequence of the Nyquist sampling theorm thingy.
00181  * (FIXME - make this comment better, and helpful.)
00182  * 
00183  * The two divisions by 4 are also a consequence of this: the contributions
00184  * returned for each frequency are split into two parts, one at i in the
00185  * table, and the other at FFT_BUFFER_SIZE - i, except for i = 0 and
00186  * FFT_BUFFER_SIZE which would otherwise get float (and then 4* when squared)
00187  * the contributions.
00188  */
00189 static void
00190 fft_output(const float *re, const float *im, float *output)
00191 {
00192     float *outputptr = output;
00193     const float *realptr = re;
00194     const float *imagptr = im;
00195     float *endptr = output + FFT_BUFFER_SIZE / 2;
00196 
00197     while (outputptr <= endptr) {
00198         *outputptr = (*realptr * *realptr) + (*imagptr * *imagptr);
00199         outputptr++;
00200         realptr++;
00201         imagptr++;
00202     }
00203     /* Do divisions to keep the constant and highest frequency terms in scale
00204      * with the other terms. */
00205     *output /= 4;
00206     *endptr /= 4;
00207 }
00208 
00209 /*
00210  * Actually perform the FFT
00211  */
00212 static void
00213 fft_calculate(float *re, float *im)
00214 {
00215     unsigned int i, j, k;
00216     unsigned int exchanges;
00217     float fact_real, fact_imag;
00218     float tmp_real, tmp_imag;
00219     unsigned int factfact;
00220 
00221     /* Set up some variables to reduce calculation in the loops */
00222     exchanges = 1;
00223     factfact = FFT_BUFFER_SIZE / 2;
00224 
00225     /* Loop through the divide and conquer steps */
00226     for (i = FFT_BUFFER_SIZE_LOG; i != 0; i--) {
00227         /* In this step, we have 2 ^ (i - 1) exchange groups, each with
00228          * 2 ^ (FFT_BUFFER_SIZE_LOG - i) exchanges
00229          */
00230         /* Loop through the exchanges in a group */
00231         for (j = 0; j != exchanges; j++) {
00232             /* Work out factor for this exchange
00233              * factor ^ (exchanges) = -1
00234              * So, real = cos(j * PI / exchanges),
00235              *     imag = sin(j * PI / exchanges)
00236              */
00237             fact_real = costable[j * factfact];
00238             fact_imag = sintable[j * factfact];
00239 
00240             /* Loop through all the exchange groups */
00241             for (k = j; k < FFT_BUFFER_SIZE; k += exchanges << 1) {
00242                 int k1 = k + exchanges;
00243                 /* newval[k]  := val[k] + factor * val[k1]
00244                  * newval[k1] := val[k] - factor * val[k1]
00245                  **/
00246                 /* FIXME - potential scope for more optimization here? */
00247                 tmp_real = fact_real * re[k1] - fact_imag * im[k1];
00248                 tmp_imag = fact_real * im[k1] + fact_imag * re[k1];
00249                 re[k1] = re[k] - tmp_real;
00250                 im[k1] = im[k] - tmp_imag;
00251                 re[k] += tmp_real;
00252                 im[k] += tmp_imag;
00253             }
00254         }
00255         exchanges <<= 1;
00256         factfact >>= 1;
00257     }
00258 }
00259 
00260 static int
00261 reverseBits(unsigned int initial)
00262 {
00263     unsigned int reversed = 0, loop;
00264     for (loop = 0; loop < FFT_BUFFER_SIZE_LOG; loop++) {
00265         reversed <<= 1;
00266         reversed += (initial & 1);
00267         initial >>= 1;
00268     }
00269     return reversed;
00270 }

Generated on Wed Apr 6 2011 for Audacious by  doxygen 1.7.1