1 /*****************************************************************************
2 * fft.c: Iterative implementation of a FFT
3 *****************************************************************************
4 * $Id: fft.c,v 1.3 2003/12/22 14:32:56 sam Exp $
6 * Mainly taken from XMMS's code
8 * Authors: Richard Boulton <richard@tartarus.org>
9 * Ralph Loader <suckfish@ihug.co.nz>
11 * This program is free software; you can redistribute it and/or modify
12 * it under the terms of the GNU General Public License as published by
13 * the Free Software Foundation; either version 2 of the License, or
14 * (at your option) any later version.
16 * This program is distributed in the hope that it will be useful,
17 * but WITHOUT ANY WARRANTY; without even the implied warranty of
18 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
19 * GNU General Public License for more details.
21 * You should have received a copy of the GNU General Public License
22 * along with this program; if not, write to the Free Software
23 * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111, USA.
24 *****************************************************************************/
34 #define PI 3.14159265358979323846 /* pi */
38 /******************************************************************************
40 *****************************************************************************/
41 static void fft_prepare(const sound_sample *input, float * re, float * im);
42 static void fft_calculate(float * re, float * im);
43 static void fft_output(const float *re, const float *im, float *output);
44 static int reverseBits(unsigned int initial);
47 /* Table to speed up bit reverse copy */
48 static unsigned int bitReverse[FFT_BUFFER_SIZE];
50 /* The next two tables could be made to use less space in memory, since they
51 * overlap hugely, but hey. */
52 static float sintable[FFT_BUFFER_SIZE / 2];
53 static float costable[FFT_BUFFER_SIZE / 2];
55 /*****************************************************************************
56 * These functions are the ones called externally
57 *****************************************************************************/
60 * Initialisation routine - sets up tables and space to work in.
61 * Returns a pointer to internal state, to be used when performing calls.
62 * On error, returns NULL.
63 * The pointer should be freed when it is finished with, by fft_close().
65 fft_state *visual_fft_init(void)
70 p_state = (fft_state *) malloc (sizeof(fft_state));
74 for(i = 0; i < FFT_BUFFER_SIZE; i++)
76 bitReverse[i] = reverseBits(i);
78 for(i = 0; i < FFT_BUFFER_SIZE / 2; i++)
80 float j = 2 * PI * i / FFT_BUFFER_SIZE;
89 * Do all the steps of the FFT, taking as input sound data (as described in
90 * sound.h) and returning the intensities of each frequency as floats in the
91 * range 0 to ((FFT_BUFFER_SIZE / 2) * 32768) ^ 2
93 * The input array is assumed to have FFT_BUFFER_SIZE elements,
94 * and the output array is assumed to have (FFT_BUFFER_SIZE / 2 + 1) elements.
95 * state is a (non-NULL) pointer returned by visual_fft_init.
97 void fft_perform(const sound_sample *input, float *output, fft_state *state) {
98 /* Convert data from sound format to be ready for FFT */
99 fft_prepare(input, state->real, state->imag);
101 /* Do the actual FFT */
102 fft_calculate(state->real, state->imag);
104 /* Convert the FFT output into intensities */
105 fft_output(state->real, state->imag, output);
111 void fft_close(fft_state *state) {
112 if(state) free(state);
115 /*****************************************************************************
116 * These functions are called from the other ones
117 *****************************************************************************/
120 * Prepare data to perform an FFT on
122 static void fft_prepare(const sound_sample *input, float * re, float * im) {
127 /* Get input, in reverse bit order */
128 for(i = 0; i < FFT_BUFFER_SIZE; i++)
130 *p_real++ = input[bitReverse[i]];
136 * Take result of an FFT and calculate the intensities of each frequency
137 * Note: only produces half as many data points as the input had.
139 static void fft_output(const float * re, const float * im, float *output)
141 float *p_output = output;
142 const float *p_real = re;
143 const float *p_imag = im;
144 float *p_end = output + FFT_BUFFER_SIZE / 2;
146 while(p_output <= p_end)
148 *p_output = (*p_real * *p_real) + (*p_imag * *p_imag);
149 p_output++; p_real++; p_imag++;
151 /* Do divisions to keep the constant and highest frequency terms in scale
152 * with the other terms. */
159 * Actually perform the FFT
161 static void fft_calculate(float * re, float * im)
163 unsigned int i, j, k;
164 unsigned int exchanges;
165 float fact_real, fact_imag;
166 float tmp_real, tmp_imag;
167 unsigned int factfact;
169 /* Set up some variables to reduce calculation in the loops */
171 factfact = FFT_BUFFER_SIZE / 2;
173 /* Loop through the divide and conquer steps */
174 for(i = FFT_BUFFER_SIZE_LOG; i != 0; i--) {
175 /* In this step, we have 2 ^ (i - 1) exchange groups, each with
176 * 2 ^ (FFT_BUFFER_SIZE_LOG - i) exchanges
178 /* Loop through the exchanges in a group */
179 for(j = 0; j != exchanges; j++) {
180 /* Work out factor for this exchange
181 * factor ^ (exchanges) = -1
182 * So, real = cos(j * PI / exchanges),
183 * imag = sin(j * PI / exchanges)
185 fact_real = costable[j * factfact];
186 fact_imag = sintable[j * factfact];
188 /* Loop through all the exchange groups */
189 for(k = j; k < FFT_BUFFER_SIZE; k += exchanges << 1) {
190 int k1 = k + exchanges;
191 tmp_real = fact_real * re[k1] - fact_imag * im[k1];
192 tmp_imag = fact_real * im[k1] + fact_imag * re[k1];
193 re[k1] = re[k] - tmp_real;
194 im[k1] = im[k] - tmp_imag;
204 static int reverseBits(unsigned int initial)
206 unsigned int reversed = 0, loop;
207 for(loop = 0; loop < FFT_BUFFER_SIZE_LOG; loop++) {
209 reversed += (initial & 1);