/*****************************************************************************
* fft.c: Iterative implementation of a FFT
*****************************************************************************
- * $Id: fft.c,v 1.2 2003/10/24 17:43:51 sam Exp $
+ * $Id: fft.c,v 1.3 2003/12/22 14:32:56 sam Exp $
*
* Mainly taken from XMMS's code
- *
+ *
* Authors: Richard Boulton <richard@tartarus.org>
* Ralph Loader <suckfish@ihug.co.nz>
*
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation; either version 2 of the License, or
* (at your option) any later version.
- *
+ *
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
#endif
/******************************************************************************
- * Local prototypes
+ * Local prototypes
*****************************************************************************/
static void fft_prepare(const sound_sample *input, float * re, float * im);
static void fft_calculate(float * re, float * im);
* On error, returns NULL.
* The pointer should be freed when it is finished with, by fft_close().
*/
-fft_state *visual_fft_init(void)
+fft_state *visual_fft_init(void)
{
fft_state *p_state;
unsigned int i;
p_state = (fft_state *) malloc (sizeof(fft_state));
- if(! p_state )
+ if(! p_state )
return NULL;
- for(i = 0; i < FFT_BUFFER_SIZE; i++)
+ for(i = 0; i < FFT_BUFFER_SIZE; i++)
{
- bitReverse[i] = reverseBits(i);
+ bitReverse[i] = reverseBits(i);
}
- for(i = 0; i < FFT_BUFFER_SIZE / 2; i++)
+ for(i = 0; i < FFT_BUFFER_SIZE / 2; i++)
{
- float j = 2 * PI * i / FFT_BUFFER_SIZE;
- costable[i] = cos(j);
- sintable[i] = sin(j);
+ float j = 2 * PI * i / FFT_BUFFER_SIZE;
+ costable[i] = cos(j);
+ sintable[i] = sin(j);
}
return p_state;
unsigned int i;
float *p_real = re;
float *p_imag = im;
-
+
/* Get input, in reverse bit order */
- for(i = 0; i < FFT_BUFFER_SIZE; i++)
+ for(i = 0; i < FFT_BUFFER_SIZE; i++)
{
- *p_real++ = input[bitReverse[i]];
- *p_imag++ = 0;
+ *p_real++ = input[bitReverse[i]];
+ *p_imag++ = 0;
}
}
const float *p_real = re;
const float *p_imag = im;
float *p_end = output + FFT_BUFFER_SIZE / 2;
-
- while(p_output <= p_end)
+
+ while(p_output <= p_end)
{
- *p_output = (*p_real * *p_real) + (*p_imag * *p_imag);
- p_output++; p_real++; p_imag++;
+ *p_output = (*p_real * *p_real) + (*p_imag * *p_imag);
+ p_output++; p_real++; p_imag++;
}
/* Do divisions to keep the constant and highest frequency terms in scale
* with the other terms. */
/*
* Actually perform the FFT
*/
-static void fft_calculate(float * re, float * im)
+static void fft_calculate(float * re, float * im)
{
unsigned int i, j, k;
unsigned int exchanges;
float fact_real, fact_imag;
float tmp_real, tmp_imag;
unsigned int factfact;
-
+
/* Set up some variables to reduce calculation in the loops */
exchanges = 1;
factfact = FFT_BUFFER_SIZE / 2;
/* Loop through the divide and conquer steps */
for(i = FFT_BUFFER_SIZE_LOG; i != 0; i--) {
- /* In this step, we have 2 ^ (i - 1) exchange groups, each with
- * 2 ^ (FFT_BUFFER_SIZE_LOG - i) exchanges
- */
- /* Loop through the exchanges in a group */
- for(j = 0; j != exchanges; j++) {
- /* Work out factor for this exchange
- * factor ^ (exchanges) = -1
- * So, real = cos(j * PI / exchanges),
- * imag = sin(j * PI / exchanges)
- */
- fact_real = costable[j * factfact];
- fact_imag = sintable[j * factfact];
-
- /* Loop through all the exchange groups */
- for(k = j; k < FFT_BUFFER_SIZE; k += exchanges << 1) {
- int k1 = k + exchanges;
- tmp_real = fact_real * re[k1] - fact_imag * im[k1];
- tmp_imag = fact_real * im[k1] + fact_imag * re[k1];
- re[k1] = re[k] - tmp_real;
- im[k1] = im[k] - tmp_imag;
- re[k] += tmp_real;
- im[k] += tmp_imag;
- }
- }
- exchanges <<= 1;
- factfact >>= 1;
+ /* In this step, we have 2 ^ (i - 1) exchange groups, each with
+ * 2 ^ (FFT_BUFFER_SIZE_LOG - i) exchanges
+ */
+ /* Loop through the exchanges in a group */
+ for(j = 0; j != exchanges; j++) {
+ /* Work out factor for this exchange
+ * factor ^ (exchanges) = -1
+ * So, real = cos(j * PI / exchanges),
+ * imag = sin(j * PI / exchanges)
+ */
+ fact_real = costable[j * factfact];
+ fact_imag = sintable[j * factfact];
+
+ /* Loop through all the exchange groups */
+ for(k = j; k < FFT_BUFFER_SIZE; k += exchanges << 1) {
+ int k1 = k + exchanges;
+ tmp_real = fact_real * re[k1] - fact_imag * im[k1];
+ tmp_imag = fact_real * im[k1] + fact_imag * re[k1];
+ re[k1] = re[k] - tmp_real;
+ im[k1] = im[k] - tmp_imag;
+ re[k] += tmp_real;
+ im[k] += tmp_imag;
+ }
+ }
+ exchanges <<= 1;
+ factfact >>= 1;
}
}
-static int reverseBits(unsigned int initial)
+static int reverseBits(unsigned int initial)
{
unsigned int reversed = 0, loop;
for(loop = 0; loop < FFT_BUFFER_SIZE_LOG; loop++) {
- reversed <<= 1;
- reversed += (initial & 1);
- initial >>= 1;
+ reversed <<= 1;
+ reversed += (initial & 1);
+ initial >>= 1;
}
return reversed;
}
projects / vlc / blobdiff