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., 51 Franklin Street, Fifth Floor, Boston MA 02110-1301, USA.
  24  *****************************************************************************/
  25
  26 #include <stdlib.h>
  27 #include "fft.h"
  28
  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                         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);
  47
  48 /*****************************************************************************
  49  * These functions are the ones called externally
  50  *****************************************************************************/
  51
  52 /*
  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().
  57  */
  58 fft_state *visual_fft_init(void)
  59 {
  60     fft_state *p_state;
  61     unsigned int i;
  62
  63     p_state = malloc( sizeof(*p_state) );
  64     if(! p_state )
  65         return NULL;
  66
  67     for(i = 0; i < FFT_BUFFER_SIZE; i++)
  68     {
  69         p_state->bitReverse[i] = reverseBits(i);
  70     }
  71     for(i = 0; i < FFT_BUFFER_SIZE / 2; i++)
  72     {
  73         float j = 2 * PI * i / FFT_BUFFER_SIZE;
  74         p_state->costable[i] = cos(j);
  75         p_state->sintable[i] = sin(j);
  76     }
  77
  78     return p_state;
  79 }
  80
  81 /*
  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
  85  *
  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.
  89  */
  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 );
  93
  94     /* Do the actual FFT */
  95     fft_calculate(state->real, state->imag, state->costable, state->sintable);
  96
  97     /* Convert the FFT output into intensities */
  98     fft_output(state->real, state->imag, output);
  99 }
 100
 101 /*
 102  * Free the state.
 103  */
 104 void fft_close(fft_state *state) {
 105     free( state );
 106 }
 107
 108 /*****************************************************************************
 109  * These functions are called from the other ones
 110  *****************************************************************************/
 111
 112 /*
 113  * Prepare data to perform an FFT on
 114  */
 115 static void fft_prepare( const sound_sample *input, float * re, float * im,
 116                          const unsigned int *bitReverse ) {
 117     unsigned int i;
 118     float *p_real = re;
 119     float *p_imag = im;
 120
 121     /* Get input, in reverse bit order */
 122     for(i = 0; i < FFT_BUFFER_SIZE; i++)
 123     {
 124         *p_real++ = input[bitReverse[i]];
 125         *p_imag++ = 0;
 126     }
 127 }
 128
 129 /*
 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.
 132  */
 133 static void fft_output(const float * re, const float * im, float *output)
 134 {
 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;
 139
 140     while(p_output <= p_end)
 141     {
 142         *p_output = (*p_real * *p_real) + (*p_imag * *p_imag);
 143         p_output++; p_real++; p_imag++;
 144     }
 145     /* Do divisions to keep the constant and highest frequency terms in scale
 146      * with the other terms. */
 147     *output /= 4;
 148     *p_end /= 4;
 149 }
 150
 151
 152 /*
 153  * Actually perform the FFT
 154  */
 155 static void fft_calculate(float * re, float * im, const float *costable, const float *sintable )
 156 {
 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;
 162
 163     /* Set up some variables to reduce calculation in the loops */
 164     exchanges = 1;
 165     factfact = FFT_BUFFER_SIZE / 2;
 166
 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
 171          */
 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)
 178              */
 179             fact_real = costable[j * factfact];
 180             fact_imag = sintable[j * factfact];
 181
 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;
 189                 re[k]  += tmp_real;
 190                 im[k]  += tmp_imag;
 191             }
 192         }
 193         exchanges <<= 1;
 194         factfact >>= 1;
 195     }
 196 }
 197
 198 static int reverseBits(unsigned int initial)
 199 {
 200     unsigned int reversed = 0, loop;
 201     for(loop = 0; loop < FFT_BUFFER_SIZE_LOG; loop++) {
 202         reversed <<= 1;
 203         reversed += (initial & 1);
 204         initial >>= 1;
 205     }
 206     return reversed;
 207 }