Add a growl (http://growl.info/) notification plugin (most of code from msn.c)
[vlc.git] / 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 "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 }