git.sesse.net Git - vlc/blob - modules/visualization/visual/fft.c

   1 /*****************************************************************************
   2  * fft.c: Iterative implementation of a FFT
   3  *****************************************************************************
   4  * $Id$
   5  *
   6  * Mainly taken from XMMS's code
   7  *
   8  * Authors: Richard Boulton <richard@tartarus.org>
   9  *          Ralph Loader <suckfish@ihug.co.nz>
  10  *
  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.
  15  *
  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.
  20  *
  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  *****************************************************************************/
  25
  26 #include "fft.h"
  27
  28 #include <stdlib.h>
  29 #include <math.h>
  30 #ifndef PI
  31  #ifdef M_PI
  32   #define PI M_PI
  33  #else
  34   #define PI            3.14159265358979323846  /* pi */
  35  #endif
  36 #endif
  37
  38 /******************************************************************************
  39  * Local prototypes
  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);
  45
  46
  47 /* Table to speed up bit reverse copy */
  48 static unsigned int bitReverse[FFT_BUFFER_SIZE];
  49
  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];
  54
  55 /*****************************************************************************
  56  * These functions are the ones called externally
  57  *****************************************************************************/
  58
  59 /*
  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().
  64  */
  65 fft_state *visual_fft_init(void)
  66 {
  67     fft_state *p_state;
  68     unsigned int i;
  69
  70     p_state = (fft_state *) malloc (sizeof(fft_state));
  71     if(! p_state )
  72         return NULL;
  73
  74     for(i = 0; i < FFT_BUFFER_SIZE; i++)
  75     {
  76         bitReverse[i] = reverseBits(i);
  77     }
  78     for(i = 0; i < FFT_BUFFER_SIZE / 2; i++)
  79     {
  80         float j = 2 * PI * i / FFT_BUFFER_SIZE;
  81         costable[i] = cos(j);
  82         sintable[i] = sin(j);
  83     }
  84
  85     return p_state;
  86 }
  87
  88 /*
  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
  92  *
  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.
  96  */
  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);
 100
 101     /* Do the actual FFT */
 102     fft_calculate(state->real, state->imag);
 103
 104     /* Convert the FFT output into intensities */
 105     fft_output(state->real, state->imag, output);
 106 }
 107
 108 /*
 109  * Free the state.
 110  */
 111 void fft_close(fft_state *state) {
 112     if(state) free(state);
 113 }
 114
 115 /*****************************************************************************
 116  * These functions are called from the other ones
 117  *****************************************************************************/
 118
 119 /*
 120  * Prepare data to perform an FFT on
 121  */
 122 static void fft_prepare(const sound_sample *input, float * re, float * im) {
 123     unsigned int i;
 124     float *p_real = re;
 125     float *p_imag = im;
 126
 127     /* Get input, in reverse bit order */
 128     for(i = 0; i < FFT_BUFFER_SIZE; i++)
 129     {
 130         *p_real++ = input[bitReverse[i]];
 131         *p_imag++ = 0;
 132     }
 133 }
 134
 135 /*
 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.
 138  */
 139 static void fft_output(const float * re, const float * im, float *output)
 140 {
 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;
 145
 146     while(p_output <= p_end)
 147     {
 148         *p_output = (*p_real * *p_real) + (*p_imag * *p_imag);
 149         p_output++; p_real++; p_imag++;
 150     }
 151     /* Do divisions to keep the constant and highest frequency terms in scale
 152      * with the other terms. */
 153     *output /= 4;
 154     *p_end /= 4;
 155 }
 156
 157
 158 /*
 159  * Actually perform the FFT
 160  */
 161 static void fft_calculate(float * re, float * im)
 162 {
 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;
 168
 169     /* Set up some variables to reduce calculation in the loops */
 170     exchanges = 1;
 171     factfact = FFT_BUFFER_SIZE / 2;
 172
 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
 177          */
 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)
 184              */
 185             fact_real = costable[j * factfact];
 186             fact_imag = sintable[j * factfact];
 187
 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;
 195                 re[k]  += tmp_real;
 196                 im[k]  += tmp_imag;
 197             }
 198         }
 199         exchanges <<= 1;
 200         factfact >>= 1;
 201     }
 202 }
 203
 204 static int reverseBits(unsigned int initial)
 205 {
 206     unsigned int reversed = 0, loop;
 207     for(loop = 0; loop < FFT_BUFFER_SIZE_LOG; loop++) {
 208         reversed <<= 1;
 209         reversed += (initial & 1);
 210         initial >>= 1;
 211     }
 212     return reversed;
 213 }