1 /*
2 * FAAC - Freeware Advanced Audio Coder
3 * $Id: fft.c,v 1.12 2005/02/02 07:49:55 sur Exp $
4 * Copyright (C) 2002 Krzysztof Nikiel
5 *
6 * This library is free software; you can redistribute it and/or
7 * modify it under the terms of the GNU Lesser General Public
8 * License as published by the Free Software Foundation; either
9 * version 2.1 of the License, or (at your option) any later version.
10 *
11 * This library is distributed in the hope that it will be useful,
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
14 * Lesser General Public License for more details.
15
16 * You should have received a copy of the GNU Lesser General Public
17 * License along with this library; if not, write to the Free Software
18 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
19 *
20 */
21
22 #include <math.h>
23 #include <stdlib.h>
24 #include <stdio.h>
25
26 #include "fft.h"
27 #include "util.h"
28
29 #define MAXLOGM 9
30 #define MAXLOGR 8
31
32 #if defined DRM && !defined DRM_1024
33
34 #include "kiss_fft/kiss_fft.h"
35 #include "kiss_fft/kiss_fftr.h"
36
37 static const int logm_to_nfft[] =
38 {
39 /* 0 1 2 3 */
40 0, 0, 0, 0,
41 /* 4 5 6 7 */
42 0, 0, 60, 0,
43 /* 8 9 */
44 240, 480
45 };
46
fft_initialize(FFT_Tables * fft_tables)47 void fft_initialize( FFT_Tables *fft_tables )
48 {
49 memset( fft_tables->cfg, 0, sizeof( fft_tables->cfg ) );
50 }
fft_terminate(FFT_Tables * fft_tables)51 void fft_terminate( FFT_Tables *fft_tables )
52 {
53 unsigned int i;
54 for ( i = 0; i < sizeof( fft_tables->cfg ) / sizeof( fft_tables->cfg[0] ); i++ )
55 {
56 if ( fft_tables->cfg[i][0] )
57 {
58 free( fft_tables->cfg[i][0] );
59 fft_tables->cfg[i][0] = NULL;
60 }
61 if ( fft_tables->cfg[i][1] )
62 {
63 free( fft_tables->cfg[i][1] );
64 fft_tables->cfg[i][1] = NULL;
65 }
66 }
67 }
68
rfft(FFT_Tables * fft_tables,double * x,int logm)69 void rfft( FFT_Tables *fft_tables, double *x, int logm )
70 {
71 /* sur: use real-only optimized FFT */
72
73 int nfft = 0;
74
75 kiss_fft_scalar fin[1 << MAXLOGR];
76 kiss_fft_cpx fout[1 << MAXLOGR];
77
78 if ( logm > MAXLOGR )
79 {
80 fprintf(stderr, "fft size too big\n");
81 exit(1);
82 }
83
84 nfft = logm_to_nfft[logm];
85
86 if ( fft_tables->cfg[logm][0] == NULL )
87 {
88 if ( nfft )
89 {
90 fft_tables->cfg[logm][0] = kiss_fftr_alloc( nfft, 0, NULL, NULL );
91 }
92 else
93 {
94 fprintf(stderr, "bad logm = %d\n", logm);
95 exit( 1 );
96 }
97 }
98
99 if ( fft_tables->cfg[logm][0] )
100 {
101 unsigned int i;
102
103 for ( i = 0; i < nfft; i++ )
104 {
105 fin[i] = x[i];
106 }
107
108 kiss_fftr( (kiss_fftr_cfg)fft_tables->cfg[logm][0], fin, fout );
109
110 for ( i = 0; i < nfft / 2; i++ )
111 {
112 x[i] = fout[i].r;
113 x[i + nfft / 2] = fout[i].i;
114 }
115 }
116 else
117 {
118 fprintf( stderr, "bad config for logm = %d\n", logm);
119 exit( 1 );
120 }
121 }
122
fft(FFT_Tables * fft_tables,double * xr,double * xi,int logm)123 void fft( FFT_Tables *fft_tables, double *xr, double *xi, int logm )
124 {
125 int nfft = 0;
126
127 kiss_fft_cpx fin[1 << MAXLOGM];
128 kiss_fft_cpx fout[1 << MAXLOGM];
129
130 if ( logm > MAXLOGM )
131 {
132 fprintf(stderr, "fft size too big\n");
133 exit(1);
134 }
135
136 nfft = logm_to_nfft[logm];
137
138 if ( fft_tables->cfg[logm][0] == NULL )
139 {
140 if ( nfft )
141 {
142 fft_tables->cfg[logm][0] = kiss_fft_alloc( nfft, 0, NULL, NULL );
143 }
144 else
145 {
146 fprintf(stderr, "bad logm = %d\n", logm);
147 exit( 1 );
148 }
149 }
150
151 if ( fft_tables->cfg[logm][0] )
152 {
153 unsigned int i;
154
155 for ( i = 0; i < nfft; i++ )
156 {
157 fin[i].r = xr[i];
158 fin[i].i = xi[i];
159 }
160
161 kiss_fft( (kiss_fft_cfg)fft_tables->cfg[logm][0], fin, fout );
162
163 for ( i = 0; i < nfft; i++ )
164 {
165 xr[i] = fout[i].r;
166 xi[i] = fout[i].i;
167 }
168 }
169 else
170 {
171 fprintf( stderr, "bad config for logm = %d\n", logm);
172 exit( 1 );
173 }
174 }
175
ffti(FFT_Tables * fft_tables,double * xr,double * xi,int logm)176 void ffti( FFT_Tables *fft_tables, double *xr, double *xi, int logm )
177 {
178 int nfft = 0;
179
180 kiss_fft_cpx fin[1 << MAXLOGM];
181 kiss_fft_cpx fout[1 << MAXLOGM];
182
183 if ( logm > MAXLOGM )
184 {
185 fprintf(stderr, "fft size too big\n");
186 exit(1);
187 }
188
189 nfft = logm_to_nfft[logm];
190
191 if ( fft_tables->cfg[logm][1] == NULL )
192 {
193 if ( nfft )
194 {
195 fft_tables->cfg[logm][1] = kiss_fft_alloc( nfft, 1, NULL, NULL );
196 }
197 else
198 {
199 fprintf(stderr, "bad logm = %d\n", logm);
200 exit( 1 );
201 }
202 }
203
204 if ( fft_tables->cfg[logm][1] )
205 {
206 unsigned int i;
207 double fac = 1.0 / (double)nfft;
208
209 for ( i = 0; i < nfft; i++ )
210 {
211 fin[i].r = xr[i];
212 fin[i].i = xi[i];
213 }
214
215 kiss_fft( (kiss_fft_cfg)fft_tables->cfg[logm][1], fin, fout );
216
217 for ( i = 0; i < nfft; i++ )
218 {
219 xr[i] = fout[i].r * fac;
220 xi[i] = fout[i].i * fac;
221 }
222 }
223 else
224 {
225 fprintf( stderr, "bad config for logm = %d\n", logm);
226 exit( 1 );
227 }
228 }
229
230 #else /* !defined DRM || defined DRM_1024 */
231
fft_initialize(FFT_Tables * fft_tables)232 void fft_initialize( FFT_Tables *fft_tables )
233 {
234 int i;
235 fft_tables->costbl = AllocMemory( (MAXLOGM+1) * sizeof( fft_tables->costbl[0] ) );
236 fft_tables->negsintbl = AllocMemory( (MAXLOGM+1) * sizeof( fft_tables->negsintbl[0] ) );
237 fft_tables->reordertbl = AllocMemory( (MAXLOGM+1) * sizeof( fft_tables->reordertbl[0] ) );
238
239 for( i = 0; i< MAXLOGM+1; i++ )
240 {
241 fft_tables->costbl[i] = NULL;
242 fft_tables->negsintbl[i] = NULL;
243 fft_tables->reordertbl[i] = NULL;
244 }
245 }
246
fft_terminate(FFT_Tables * fft_tables)247 void fft_terminate( FFT_Tables *fft_tables )
248 {
249 int i;
250
251 for( i = 0; i< MAXLOGM+1; i++ )
252 {
253 if( fft_tables->costbl[i] != NULL )
254 FreeMemory( fft_tables->costbl[i] );
255
256 if( fft_tables->negsintbl[i] != NULL )
257 FreeMemory( fft_tables->negsintbl[i] );
258
259 if( fft_tables->reordertbl[i] != NULL )
260 FreeMemory( fft_tables->reordertbl[i] );
261 }
262
263 FreeMemory( fft_tables->costbl );
264 FreeMemory( fft_tables->negsintbl );
265 FreeMemory( fft_tables->reordertbl );
266
267 fft_tables->costbl = NULL;
268 fft_tables->negsintbl = NULL;
269 fft_tables->reordertbl = NULL;
270 }
271
reorder(FFT_Tables * fft_tables,double * x,int logm)272 static void reorder( FFT_Tables *fft_tables, double *x, int logm)
273 {
274 int i;
275 int size = 1 << logm;
276 unsigned short *r; //size
277
278
279 if ( fft_tables->reordertbl[logm] == NULL ) // create bit reversing table
280 {
281 fft_tables->reordertbl[logm] = AllocMemory(size * sizeof(*(fft_tables->reordertbl[0])));
282
283 for (i = 0; i < size; i++)
284 {
285 int reversed = 0;
286 int b0;
287 int tmp = i;
288
289 for (b0 = 0; b0 < logm; b0++)
290 {
291 reversed = (reversed << 1) | (tmp & 1);
292 tmp >>= 1;
293 }
294 fft_tables->reordertbl[logm][i] = reversed;
295 }
296 }
297
298 r = fft_tables->reordertbl[logm];
299
300 for (i = 0; i < size; i++)
301 {
302 int j = r[i];
303 double tmp;
304
305 if (j <= i)
306 continue;
307
308 tmp = x[i];
309 x[i] = x[j];
310 x[j] = tmp;
311 }
312 }
313
fft_proc(double * xr,double * xi,fftfloat * refac,fftfloat * imfac,int size)314 static void fft_proc(
315 double *xr,
316 double *xi,
317 fftfloat *refac,
318 fftfloat *imfac,
319 int size)
320 {
321 int step, shift, pos;
322 int exp, estep;
323
324 estep = size;
325 for (step = 1; step < size; step *= 2)
326 {
327 int x1;
328 int x2 = 0;
329 estep >>= 1;
330 for (pos = 0; pos < size; pos += (2 * step))
331 {
332 x1 = x2;
333 x2 += step;
334 exp = 0;
335 for (shift = 0; shift < step; shift++)
336 {
337 double v2r, v2i;
338
339 v2r = xr[x2] * refac[exp] - xi[x2] * imfac[exp];
340 v2i = xr[x2] * imfac[exp] + xi[x2] * refac[exp];
341
342 xr[x2] = xr[x1] - v2r;
343 xr[x1] += v2r;
344
345 xi[x2] = xi[x1] - v2i;
346
347 xi[x1] += v2i;
348
349 exp += estep;
350
351 x1++;
352 x2++;
353 }
354 }
355 }
356 }
357
check_tables(FFT_Tables * fft_tables,int logm)358 static void check_tables( FFT_Tables *fft_tables, int logm)
359 {
360 if( fft_tables->costbl[logm] == NULL )
361 {
362 int i;
363 int size = 1 << logm;
364
365 if( fft_tables->negsintbl[logm] != NULL )
366 FreeMemory( fft_tables->negsintbl[logm] );
367
368 fft_tables->costbl[logm] = AllocMemory((size / 2) * sizeof(*(fft_tables->costbl[0])));
369 fft_tables->negsintbl[logm] = AllocMemory((size / 2) * sizeof(*(fft_tables->negsintbl[0])));
370
371 for (i = 0; i < (size >> 1); i++)
372 {
373 double theta = 2.0 * M_PI * ((double) i) / (double) size;
374 fft_tables->costbl[logm][i] = cos(theta);
375 fft_tables->negsintbl[logm][i] = -sin(theta);
376 }
377 }
378 }
379
fft(FFT_Tables * fft_tables,double * xr,double * xi,int logm)380 void fft( FFT_Tables *fft_tables, double *xr, double *xi, int logm)
381 {
382 if (logm > MAXLOGM)
383 {
384 fprintf(stderr, "fft size too big\n");
385 exit(1);
386 }
387
388 if (logm < 1)
389 {
390 //printf("logm < 1\n");
391 return;
392 }
393
394 check_tables( fft_tables, logm);
395
396 reorder( fft_tables, xr, logm);
397 reorder( fft_tables, xi, logm);
398
399 fft_proc( xr, xi, fft_tables->costbl[logm], fft_tables->negsintbl[logm], 1 << logm );
400 }
401
rfft(FFT_Tables * fft_tables,double * x,int logm)402 void rfft( FFT_Tables *fft_tables, double *x, int logm)
403 {
404 double xi[1 << MAXLOGR];
405
406 if (logm > MAXLOGR)
407 {
408 fprintf(stderr, "rfft size too big\n");
409 exit(1);
410 }
411
412 memset(xi, 0, (1 << logm) * sizeof(xi[0]));
413
414 fft( fft_tables, x, xi, logm);
415
416 memcpy(x + (1 << (logm - 1)), xi, (1 << (logm - 1)) * sizeof(*x));
417 }
418
ffti(FFT_Tables * fft_tables,double * xr,double * xi,int logm)419 void ffti( FFT_Tables *fft_tables, double *xr, double *xi, int logm)
420 {
421 int i, size;
422 double fac;
423 double *xrp, *xip;
424
425 fft( fft_tables, xi, xr, logm);
426
427 size = 1 << logm;
428 fac = 1.0 / size;
429 xrp = xr;
430 xip = xi;
431
432 for (i = 0; i < size; i++)
433 {
434 *xrp++ *= fac;
435 *xip++ *= fac;
436 }
437 }
438
439 #endif /* defined DRM && !defined DRM_1024 */
440