1 /*****************************************************************************
2 * fft.c: Iterative implementation of a FFT
3 *****************************************************************************
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., 51 Franklin Street, Fifth Floor, Boston MA 02110-1301, USA.
24 *****************************************************************************/
34 #define PI 3.14159265358979323846 /* pi */
38 /******************************************************************************
40 *****************************************************************************/
41 static void fft_prepare(const sound_sample *input, float * re, float * im,
42 const unsigned int *bitReverse);
43 static void fft_calculate(float * re, float * im,
44 const float *costable, const float *sintable );
45 static void fft_output(const float *re, const float *im, float *output);
46 static int reverseBits(unsigned int initial);
48 /*****************************************************************************
49 * These functions are the ones called externally
50 *****************************************************************************/
53 * Initialisation routine - sets up tables and space to work in.
54 * Returns a pointer to internal state, to be used when performing calls.
55 * On error, returns NULL.
56 * The pointer should be freed when it is finished with, by fft_close().
58 fft_state *visual_fft_init(void)
63 p_state = malloc( sizeof(*p_state) );
67 for(i = 0; i < FFT_BUFFER_SIZE; i++)
69 p_state->bitReverse[i] = reverseBits(i);
71 for(i = 0; i < FFT_BUFFER_SIZE / 2; i++)
73 float j = 2 * PI * i / FFT_BUFFER_SIZE;
74 p_state->costable[i] = cos(j);
75 p_state->sintable[i] = sin(j);
82 * Do all the steps of the FFT, taking as input sound data (as described in
83 * sound.h) and returning the intensities of each frequency as floats in the
84 * range 0 to ((FFT_BUFFER_SIZE / 2) * 32768) ^ 2
86 * The input array is assumed to have FFT_BUFFER_SIZE elements,
87 * and the output array is assumed to have (FFT_BUFFER_SIZE / 2 + 1) elements.
88 * state is a (non-NULL) pointer returned by visual_fft_init.
90 void fft_perform(const sound_sample *input, float *output, fft_state *state) {
91 /* Convert data from sound format to be ready for FFT */
92 fft_prepare(input, state->real, state->imag, state->bitReverse );
94 /* Do the actual FFT */
95 fft_calculate(state->real, state->imag, state->costable, state->sintable);
97 /* Convert the FFT output into intensities */
98 fft_output(state->real, state->imag, output);
104 void fft_close(fft_state *state) {
108 /*****************************************************************************
109 * These functions are called from the other ones
110 *****************************************************************************/
113 * Prepare data to perform an FFT on
115 static void fft_prepare( const sound_sample *input, float * re, float * im,
116 const unsigned int *bitReverse ) {
121 /* Get input, in reverse bit order */
122 for(i = 0; i < FFT_BUFFER_SIZE; i++)
124 *p_real++ = input[bitReverse[i]];
130 * Take result of an FFT and calculate the intensities of each frequency
131 * Note: only produces half as many data points as the input had.
133 static void fft_output(const float * re, const float * im, float *output)
135 float *p_output = output;
136 const float *p_real = re;
137 const float *p_imag = im;
138 float *p_end = output + FFT_BUFFER_SIZE / 2;
140 while(p_output <= p_end)
142 *p_output = (*p_real * *p_real) + (*p_imag * *p_imag);
143 p_output++; p_real++; p_imag++;
145 /* Do divisions to keep the constant and highest frequency terms in scale
146 * with the other terms. */
153 * Actually perform the FFT
155 static void fft_calculate(float * re, float * im, const float *costable, const float *sintable )
157 unsigned int i, j, k;
158 unsigned int exchanges;
159 float fact_real, fact_imag;
160 float tmp_real, tmp_imag;
161 unsigned int factfact;
163 /* Set up some variables to reduce calculation in the loops */
165 factfact = FFT_BUFFER_SIZE / 2;
167 /* Loop through the divide and conquer steps */
168 for(i = FFT_BUFFER_SIZE_LOG; i != 0; i--) {
169 /* In this step, we have 2 ^ (i - 1) exchange groups, each with
170 * 2 ^ (FFT_BUFFER_SIZE_LOG - i) exchanges
172 /* Loop through the exchanges in a group */
173 for(j = 0; j != exchanges; j++) {
174 /* Work out factor for this exchange
175 * factor ^ (exchanges) = -1
176 * So, real = cos(j * PI / exchanges),
177 * imag = sin(j * PI / exchanges)
179 fact_real = costable[j * factfact];
180 fact_imag = sintable[j * factfact];
182 /* Loop through all the exchange groups */
183 for(k = j; k < FFT_BUFFER_SIZE; k += exchanges << 1) {
184 int k1 = k + exchanges;
185 tmp_real = fact_real * re[k1] - fact_imag * im[k1];
186 tmp_imag = fact_real * im[k1] + fact_imag * re[k1];
187 re[k1] = re[k] - tmp_real;
188 im[k1] = im[k] - tmp_imag;
198 static int reverseBits(unsigned int initial)
200 unsigned int reversed = 0, loop;
201 for(loop = 0; loop < FFT_BUFFER_SIZE_LOG; loop++) {
203 reversed += (initial & 1);