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