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