1 /*
2 LICENSE
3 -------
4 Copyright 2005-2013 Nullsoft, Inc.
5 All rights reserved.
6
7 Redistribution and use in source and binary forms, with or without modification,
8 are permitted provided that the following conditions are met:
9
10 * Redistributions of source code must retain the above copyright notice,
11 this list of conditions and the following disclaimer.
12
13 * Redistributions in binary form must reproduce the above copyright notice,
14 this list of conditions and the following disclaimer in the documentation
15 and/or other materials provided with the distribution.
16
17 * Neither the name of Nullsoft nor the names of its contributors may be used to
18 endorse or promote products derived from this software without specific prior written permission.
19
20 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR
21 IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND
22 FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR
23 CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
24 DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
25 DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER
26 IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
27 OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
28 */
29
30 #include <math.h>
31 #include <memory.h>
32 #include "fft.h"
33
34 #define PI 3.141592653589793238462643383279502884197169399f
35
36 #define SafeDeleteArray(x) { if (x) { delete [] x; x = 0; } }
37
38 /*****************************************************************************/
39
FFT()40 FFT::FFT()
41 {
42 NFREQ = 0;
43
44 envelope = 0;
45 equalize = 0;
46 bitrevtable = 0;
47 cossintable = 0;
48 temp1 = 0;
49 temp2 = 0;
50 }
51
52 /*****************************************************************************/
53
~FFT()54 FFT::~FFT()
55 {
56 CleanUp();
57 }
58
59 /*****************************************************************************/
60
Init(int samples_in,int samples_out,int bEqualize,float envelope_power)61 void FFT::Init(int samples_in, int samples_out, int bEqualize, float envelope_power)
62 {
63 // samples_in: # of waveform samples you'll feed into the FFT
64 // samples_out: # of frequency samples you want out; MUST BE A POWER OF 2.
65 // bEqualize: set to 1 if you want the magnitude of the basses and trebles
66 // to be roughly equalized; 0 to leave them untouched.
67 // envelope_power: set to -1 to disable the envelope; otherwise, specify
68 // the envelope power you want here. See InitEnvelopeTable for more info.
69
70 CleanUp();
71
72 m_samples_in = samples_in;
73 NFREQ = samples_out*2;
74
75 InitBitRevTable();
76 InitCosSinTable();
77 if (envelope_power > 0)
78 InitEnvelopeTable(envelope_power);
79 if (bEqualize)
80 InitEqualizeTable();
81 temp1 = new float[NFREQ];
82 temp2 = new float[NFREQ];
83 }
84
85 /*****************************************************************************/
86
CleanUp()87 void FFT::CleanUp()
88 {
89 SafeDeleteArray(envelope);
90 SafeDeleteArray(equalize);
91 SafeDeleteArray(bitrevtable);
92 SafeDeleteArray(cossintable);
93 SafeDeleteArray(temp1);
94 SafeDeleteArray(temp2);
95 }
96
97 /*****************************************************************************/
98
InitEqualizeTable()99 void FFT::InitEqualizeTable()
100 {
101 int i;
102 float scaling = -0.02f;
103 float inv_half_nfreq = 1.0f/(float)(NFREQ/2);
104
105 equalize = new float[NFREQ/2];
106
107 for (i=0; i<NFREQ/2; i++)
108 equalize[i] = scaling * logf( (float)(NFREQ/2-i)*inv_half_nfreq );
109 }
110
111 /*****************************************************************************/
112
InitEnvelopeTable(float power)113 void FFT::InitEnvelopeTable(float power)
114 {
115 // this precomputation is for multiplying the waveform sample
116 // by a bell-curve-shaped envelope, so we don't see the ugly
117 // frequency response (oscillations) of a square filter.
118
119 // a power of 1.0 will compute the FFT with exactly one convolution.
120 // a power of 2.0 is like doing it twice; the resulting frequency
121 // output will be smoothed out and the peaks will be "fatter".
122 // a power of 0.5 is closer to not using an envelope, and you start
123 // to get the frequency response of the square filter as 'power'
124 // approaches zero; the peaks get tighter and more precise, but
125 // you also see small oscillations around their bases.
126
127 int i;
128 float mult = 1.0f/(float)m_samples_in * 6.2831853f;
129
130 envelope = new float[m_samples_in];
131
132 if (power==1.0f)
133 for (i=0; i<m_samples_in; i++)
134 envelope[i] = 0.5f + 0.5f*sinf(i*mult - 1.5707963268f);
135 else
136 for (i=0; i<m_samples_in; i++)
137 envelope[i] = powf(0.5f + 0.5f*sinf(i*mult - 1.5707963268f), power);
138 }
139
140 /*****************************************************************************/
141
InitBitRevTable()142 void FFT::InitBitRevTable()
143 {
144 int i,j,m,temp;
145 bitrevtable = new int[NFREQ];
146
147 for (i=0; i<NFREQ; i++)
148 bitrevtable[i] = i;
149
150 for (i=0,j=0; i < NFREQ; i++)
151 {
152 if (j > i)
153 {
154 temp = bitrevtable[i];
155 bitrevtable[i] = bitrevtable[j];
156 bitrevtable[j] = temp;
157 }
158
159 m = NFREQ >> 1;
160
161 while (m >= 1 && j >= m)
162 {
163 j -= m;
164 m >>= 1;
165 }
166
167 j += m;
168 }
169 }
170
171 /*****************************************************************************/
172
InitCosSinTable()173 void FFT::InitCosSinTable()
174 {
175
176 int i,dftsize,tabsize;
177 float theta;
178
179 dftsize = 2;
180 tabsize = 0;
181 while (dftsize <= NFREQ)
182 {
183 tabsize++;
184 dftsize <<= 1;
185 }
186
187 cossintable = new float[tabsize][2];
188
189 dftsize = 2;
190 i = 0;
191 while (dftsize <= NFREQ)
192 {
193 theta = (float)(-2.0f*PI/(float)dftsize);
194 cossintable[i][0] = (float)cosf(theta);
195 cossintable[i][1] = (float)sinf(theta);
196 i++;
197 dftsize <<= 1;
198 }
199 }
200
201 /*****************************************************************************/
202
time_to_frequency_domain(float * in_wavedata,float * out_spectraldata)203 void FFT::time_to_frequency_domain(float *in_wavedata, float *out_spectraldata)
204 {
205 // Converts time-domain samples from in_wavedata[]
206 // into frequency-domain samples in out_spectraldata[].
207 // The array lengths are the two parameters to Init().
208
209 // The last sample of the output data will represent the frequency
210 // that is 1/4th of the input sampling rate. For example,
211 // if the input wave data is sampled at 44,100 Hz, then the last
212 // sample of the spectral data output will represent the frequency
213 // 11,025 Hz. The first sample will be 0 Hz; the frequencies of
214 // the rest of the samples vary linearly in between.
215 // Note that since human hearing is limited to the range 200 - 20,000
216 // Hz. 200 is a low bass hum; 20,000 is an ear-piercing high shriek.
217 // Each time the frequency doubles, that sounds like going up an octave.
218 // That means that the difference between 200 and 300 Hz is FAR more
219 // than the difference between 5000 and 5100, for example!
220 // So, when trying to analyze bass, you'll want to look at (probably)
221 // the 200-800 Hz range; whereas for treble, you'll want the 1,400 -
222 // 11,025 Hz range.
223 // If you want to get 3 bands, try it this way:
224 // a) 11,025 / 200 = 55.125
225 // b) to get the number of octaves between 200 and 11,025 Hz, solve for n:
226 // 2^n = 55.125
227 // n = log 55.125 / log 2
228 // n = 5.785
229 // c) so each band should represent 5.785/3 = 1.928 octaves; the ranges are:
230 // 1) 200 - 200*2^1.928 or 200 - 761 Hz
231 // 2) 200*2^1.928 - 200*2^(1.928*2) or 761 - 2897 Hz
232 // 3) 200*2^(1.928*2) - 200*2^(1.928*3) or 2897 - 11025 Hz
233
234 // A simple sine-wave-based envelope is convolved with the waveform
235 // data before doing the FFT, to emeliorate the bad frequency response
236 // of a square (i.e. nonexistent) filter.
237
238 // You might want to slightly damp (blur) the input if your signal isn't
239 // of a very high quality, to reduce high-frequency noise that would
240 // otherwise show up in the output.
241
242 int j, m, i, dftsize, hdftsize, t;
243 float wr, wi, wpi, wpr, wtemp, tempr, tempi;
244
245 if (!bitrevtable) return;
246 //if (!envelope) return;
247 //if (!equalize) return;
248 if (!temp1) return;
249 if (!temp2) return;
250 if (!cossintable) return;
251
252 // 1. set up input to the fft
253 if (envelope)
254 {
255 for (i=0; i<NFREQ; i++)
256 {
257 int idx = bitrevtable[i];
258 if (idx < m_samples_in)
259 temp1[i] = in_wavedata[idx] * envelope[idx];
260 else
261 temp1[i] = 0;
262 }
263 }
264 else
265 {
266 for (i=0; i<NFREQ; i++)
267 {
268 int idx = bitrevtable[i];
269 if (idx < m_samples_in)
270 temp1[i] = in_wavedata[idx];// * envelope[idx];
271 else
272 temp1[i] = 0;
273 }
274 }
275 memset(temp2, 0, sizeof(float)*NFREQ);
276
277 // 2. perform FFT
278 float *real = temp1;
279 float *imag = temp2;
280 dftsize = 2;
281 t = 0;
282 while (dftsize <= NFREQ)
283 {
284 wpr = cossintable[t][0];
285 wpi = cossintable[t][1];
286 wr = 1.0f;
287 wi = 0.0f;
288 hdftsize = dftsize >> 1;
289
290 for (m = 0; m < hdftsize; m+=1)
291 {
292 for (i = m; i < NFREQ; i+=dftsize)
293 {
294 j = i + hdftsize;
295 tempr = wr*real[j] - wi*imag[j];
296 tempi = wr*imag[j] + wi*real[j];
297 real[j] = real[i] - tempr;
298 imag[j] = imag[i] - tempi;
299 real[i] += tempr;
300 imag[i] += tempi;
301 }
302
303 wr = (wtemp=wr)*wpr - wi*wpi;
304 wi = wi*wpr + wtemp*wpi;
305 }
306
307 dftsize <<= 1;
308 t++;
309 }
310
311 // 3. take the magnitude & equalize it (on a log10 scale) for output
312 if (equalize)
313 for (i=0; i<NFREQ/2; i++)
314 out_spectraldata[i] = equalize[i] * sqrtf(temp1[i]*temp1[i] + temp2[i]*temp2[i]);
315 else
316 for (i=0; i<NFREQ/2; i++)
317 out_spectraldata[i] = sqrtf(temp1[i]*temp1[i] + temp2[i]*temp2[i]);
318 }
319
320 /*****************************************************************************/