annotate audacious/fft.c @ 1178:ae173b6e65a3 trunk

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