Mercurial > audlegacy
annotate src/audacious/fft.c @ 3121:3b6d316f8b09 trunk
GPL3 relicensing.
author | William Pitcock <nenolod@atheme-project.org> |
---|---|
date | Fri, 20 Jul 2007 08:59:47 -0500 |
parents | 3149d4b1a9a9 |
children | f1c756f39e6c |
rev | line source |
---|---|
2313 | 1 /* Audacious - Cross-platform multimedia player |
2 * Copyright (C) 2005-2007 Audacious development team | |
3 * | |
4 * Copyright (C) 1999 Richard Boulton <richard@tartarus.org> | |
5 * Convolution stuff by Ralph Loader <suckfish@ihug.co.nz> | |
6 * | |
7 * This program is free software; you can redistribute it and/or modify | |
8 * it under the terms of the GNU General Public License as published by | |
3121
3b6d316f8b09
GPL3 relicensing.
William Pitcock <nenolod@atheme-project.org>
parents:
2313
diff
changeset
|
9 * the Free Software Foundation; under version 3 of the License. |
2313 | 10 * |
11 * This program is distributed in the hope that it will be useful, | |
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of | |
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the | |
14 * GNU General Public License for more details. | |
15 * | |
16 * You should have received a copy of the GNU General Public License | |
3121
3b6d316f8b09
GPL3 relicensing.
William Pitcock <nenolod@atheme-project.org>
parents:
2313
diff
changeset
|
17 * along with this program. If not, see <http://www.gnu.org/licenses>. |
2313 | 18 */ |
19 | |
20 /* fft.c: iterative implementation of a FFT */ | |
21 | |
22 /* | |
23 * TODO | |
24 * Remove compiling in of FFT_BUFFER_SIZE? (Might slow things down, but would | |
25 * be nice to be able to change size at runtime.) | |
26 * Finish making / checking thread-safety. | |
27 * More optimisations. | |
28 */ | |
29 | |
30 #ifdef HAVE_CONFIG_H | |
31 # include "config.h" | |
32 #endif | |
33 | |
34 #include "fft.h" | |
35 | |
36 #include <glib.h> | |
37 #include <stdlib.h> | |
38 #include <math.h> | |
39 #ifndef PI | |
40 #ifdef M_PI | |
41 #define PI M_PI | |
42 #else | |
43 #define PI 3.14159265358979323846 /* pi */ | |
44 #endif | |
45 #endif | |
46 | |
47 /* ########### */ | |
48 /* # Structs # */ | |
49 /* ########### */ | |
50 | |
51 struct _struct_fft_state { | |
52 /* Temporary data stores to perform FFT in. */ | |
53 float real[FFT_BUFFER_SIZE]; | |
54 float imag[FFT_BUFFER_SIZE]; | |
55 }; | |
56 | |
57 /* ############################# */ | |
58 /* # Local function prototypes # */ | |
59 /* ############################# */ | |
60 | |
61 static void fft_prepare(const sound_sample * input, float *re, float *im); | |
62 static void fft_calculate(float *re, float *im); | |
63 static void fft_output(const float *re, const float *im, float *output); | |
64 static int reverseBits(unsigned int initial); | |
65 | |
66 /* #################### */ | |
67 /* # Global variables # */ | |
68 /* #################### */ | |
69 | |
70 /* Table to speed up bit reverse copy */ | |
71 static unsigned int bitReverse[FFT_BUFFER_SIZE]; | |
72 | |
73 /* The next two tables could be made to use less space in memory, since they | |
74 * overlap hugely, but hey. */ | |
75 static float sintable[FFT_BUFFER_SIZE / 2]; | |
76 static float costable[FFT_BUFFER_SIZE / 2]; | |
77 | |
78 /* ############################## */ | |
79 /* # Externally called routines # */ | |
80 /* ############################## */ | |
81 | |
82 /* --------- */ | |
83 /* FFT stuff */ | |
84 /* --------- */ | |
85 | |
86 /* | |
87 * Initialisation routine - sets up tables and space to work in. | |
88 * Returns a pointer to internal state, to be used when performing calls. | |
89 * On error, returns NULL. | |
90 * The pointer should be freed when it is finished with, by fft_close(). | |
91 */ | |
92 fft_state * | |
93 fft_init(void) | |
94 { | |
95 fft_state *state; | |
96 unsigned int i; | |
97 | |
98 state = (fft_state *) g_malloc(sizeof(fft_state)); | |
99 if (!state) | |
100 return NULL; | |
101 | |
102 for (i = 0; i < FFT_BUFFER_SIZE; i++) { | |
103 bitReverse[i] = reverseBits(i); | |
104 } | |
105 for (i = 0; i < FFT_BUFFER_SIZE / 2; i++) { | |
106 float j = 2 * PI * i / FFT_BUFFER_SIZE; | |
107 costable[i] = cos(j); | |
108 sintable[i] = sin(j); | |
109 } | |
110 | |
111 return state; | |
112 } | |
113 | |
114 /* | |
115 * Do all the steps of the FFT, taking as input sound data (as described in | |
116 * sound.h) and returning the intensities of each frequency as floats in the | |
117 * range 0 to ((FFT_BUFFER_SIZE / 2) * 32768) ^ 2 | |
118 * | |
119 * FIXME - the above range assumes no frequencies present have an amplitude | |
120 * larger than that of the sample variation. But this is false: we could have | |
121 * a wave such that its maximums are always between samples, and it's just | |
122 * inside the representable range at the places samples get taken. | |
123 * Question: what _is_ the maximum value possible. Twice that value? Root | |
124 * two times that value? Hmmm. Think it depends on the frequency, too. | |
125 * | |
126 * The input array is assumed to have FFT_BUFFER_SIZE elements, | |
127 * and the output array is assumed to have (FFT_BUFFER_SIZE / 2 + 1) elements. | |
128 * state is a (non-NULL) pointer returned by fft_init. | |
129 */ | |
130 void | |
131 fft_perform(const sound_sample * input, float *output, fft_state * state) | |
132 { | |
133 /* Convert data from sound format to be ready for FFT */ | |
134 fft_prepare(input, state->real, state->imag); | |
135 | |
136 /* Do the actual FFT */ | |
137 fft_calculate(state->real, state->imag); | |
138 | |
139 /* Convert the FFT output into intensities */ | |
140 fft_output(state->real, state->imag, output); | |
141 } | |
142 | |
143 /* | |
144 * Free the state. | |
145 */ | |
146 void | |
147 fft_close(fft_state * state) | |
148 { | |
149 if (state) | |
150 free(state); | |
151 } | |
152 | |
153 /* ########################### */ | |
154 /* # Locally called routines # */ | |
155 /* ########################### */ | |
156 | |
157 /* | |
158 * Prepare data to perform an FFT on | |
159 */ | |
160 static void | |
161 fft_prepare(const sound_sample * input, float *re, float *im) | |
162 { | |
163 unsigned int i; | |
164 float *realptr = re; | |
165 float *imagptr = im; | |
166 | |
167 /* Get input, in reverse bit order */ | |
168 for (i = 0; i < FFT_BUFFER_SIZE; i++) { | |
169 *realptr++ = input[bitReverse[i]]; | |
170 *imagptr++ = 0; | |
171 } | |
172 } | |
173 | |
174 /* | |
175 * Take result of an FFT and calculate the intensities of each frequency | |
176 * Note: only produces half as many data points as the input had. | |
177 * This is roughly a consequence of the Nyquist sampling theorm thingy. | |
178 * (FIXME - make this comment better, and helpful.) | |
179 * | |
180 * The two divisions by 4 are also a consequence of this: the contributions | |
181 * returned for each frequency are split into two parts, one at i in the | |
182 * table, and the other at FFT_BUFFER_SIZE - i, except for i = 0 and | |
183 * FFT_BUFFER_SIZE which would otherwise get float (and then 4* when squared) | |
184 * the contributions. | |
185 */ | |
186 static void | |
187 fft_output(const float *re, const float *im, float *output) | |
188 { | |
189 float *outputptr = output; | |
190 const float *realptr = re; | |
191 const float *imagptr = im; | |
192 float *endptr = output + FFT_BUFFER_SIZE / 2; | |
193 | |
194 #ifdef DEBUG | |
195 unsigned int i, j; | |
196 #endif | |
197 | |
198 while (outputptr <= endptr) { | |
199 *outputptr = (*realptr * *realptr) + (*imagptr * *imagptr); | |
200 outputptr++; | |
201 realptr++; | |
202 imagptr++; | |
203 } | |
204 /* Do divisions to keep the constant and highest frequency terms in scale | |
205 * with the other terms. */ | |
206 *output /= 4; | |
207 *endptr /= 4; | |
208 | |
209 #ifdef DEBUG | |
210 printf("Recalculated input:\n"); | |
211 for (i = 0; i < FFT_BUFFER_SIZE; i++) { | |
212 float val_real = 0; | |
213 float val_imag = 0; | |
214 for (j = 0; j < FFT_BUFFER_SIZE; j++) { | |
215 float fact_real = cos(-2 * j * i * PI / FFT_BUFFER_SIZE); | |
216 float fact_imag = sin(-2 * j * i * PI / FFT_BUFFER_SIZE); | |
217 val_real += fact_real * re[j] - fact_imag * im[j]; | |
218 val_imag += fact_real * im[j] + fact_imag * re[j]; | |
219 } | |
220 printf("%5d = %8f + i * %8f\n", i, | |
221 val_real / FFT_BUFFER_SIZE, val_imag / FFT_BUFFER_SIZE); | |
222 } | |
223 printf("\n"); | |
224 #endif | |
225 } | |
226 | |
227 /* | |
228 * Actually perform the FFT | |
229 */ | |
230 static void | |
231 fft_calculate(float *re, float *im) | |
232 { | |
233 unsigned int i, j, k; | |
234 unsigned int exchanges; | |
235 float fact_real, fact_imag; | |
236 float tmp_real, tmp_imag; | |
237 unsigned int factfact; | |
238 | |
239 /* Set up some variables to reduce calculation in the loops */ | |
240 exchanges = 1; | |
241 factfact = FFT_BUFFER_SIZE / 2; | |
242 | |
243 /* Loop through the divide and conquer steps */ | |
244 for (i = FFT_BUFFER_SIZE_LOG; i != 0; i--) { | |
245 /* In this step, we have 2 ^ (i - 1) exchange groups, each with | |
246 * 2 ^ (FFT_BUFFER_SIZE_LOG - i) exchanges | |
247 */ | |
248 /* Loop through the exchanges in a group */ | |
249 for (j = 0; j != exchanges; j++) { | |
250 /* Work out factor for this exchange | |
251 * factor ^ (exchanges) = -1 | |
252 * So, real = cos(j * PI / exchanges), | |
253 * imag = sin(j * PI / exchanges) | |
254 */ | |
255 fact_real = costable[j * factfact]; | |
256 fact_imag = sintable[j * factfact]; | |
257 | |
258 /* Loop through all the exchange groups */ | |
259 for (k = j; k < FFT_BUFFER_SIZE; k += exchanges << 1) { | |
260 int k1 = k + exchanges; | |
261 /* newval[k] := val[k] + factor * val[k1] | |
262 * newval[k1] := val[k] - factor * val[k1] | |
263 **/ | |
264 #ifdef DEBUG | |
265 printf("%d %d %d\n", i, j, k); | |
266 printf("Exchange %d with %d\n", k, k1); | |
267 printf("Factor %9f + i * %8f\n", fact_real, fact_imag); | |
268 #endif | |
269 /* FIXME - potential scope for more optimization here? */ | |
270 tmp_real = fact_real * re[k1] - fact_imag * im[k1]; | |
271 tmp_imag = fact_real * im[k1] + fact_imag * re[k1]; | |
272 re[k1] = re[k] - tmp_real; | |
273 im[k1] = im[k] - tmp_imag; | |
274 re[k] += tmp_real; | |
275 im[k] += tmp_imag; | |
276 #ifdef DEBUG | |
277 for (k1 = 0; k1 < FFT_BUFFER_SIZE; k1++) { | |
278 printf("%5d = %8f + i * %8f\n", k1, real[k1], imag[k1]); | |
279 } | |
280 #endif | |
281 } | |
282 } | |
283 exchanges <<= 1; | |
284 factfact >>= 1; | |
285 } | |
286 } | |
287 | |
288 static int | |
289 reverseBits(unsigned int initial) | |
290 { | |
291 unsigned int reversed = 0, loop; | |
292 for (loop = 0; loop < FFT_BUFFER_SIZE_LOG; loop++) { | |
293 reversed <<= 1; | |
294 reversed += (initial & 1); | |
295 initial >>= 1; | |
296 } | |
297 return reversed; | |
298 } |