Audacious
$Id:Doxyfile42802007-03-2104:39:00Znenolod$
|
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 }