annotate audacious/fft.c @ 2155:840263ee3a5c trunk

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