ref: 31a27a38271fdff180f07c621fde650cb31977bd
dir: /libfaac/fft.c/
/*
* FAAC - Freeware Advanced Audio Coder
* $Id: fft.c,v 1.12 2005/02/02 07:49:55 sur Exp $
* Copyright (C) 2002 Krzysztof Nikiel
*
* This library is free software; you can redistribute it and/or
* modify it under the terms of the GNU Lesser General Public
* License as published by the Free Software Foundation; either
* version 2.1 of the License, or (at your option) any later version.
*
* This library is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
* Lesser General Public License for more details.
* You should have received a copy of the GNU Lesser General Public
* License along with this library; if not, write to the Free Software
* Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
*
*/
#include <math.h>
#include <stdlib.h>
#include <stdio.h>
#include "fft.h"
#include "util.h"
#define MAXLOGM 9
#define MAXLOGR 8
#if defined DRM && !defined DRM_1024
#include "kiss_fft/kiss_fft.h"
#include "kiss_fft/kiss_fftr.h"
static const int logm_to_nfft[] =
{
/* 0 1 2 3 */
0, 0, 0, 0,
/* 4 5 6 7 */
0, 0, 60, 0,
/* 8 9 */
240, 480
};
void fft_initialize( FFT_Tables *fft_tables )
{
memset( fft_tables->cfg, 0, sizeof( fft_tables->cfg ) );
}
void fft_terminate( FFT_Tables *fft_tables )
{
unsigned int i;
for ( i = 0; i < sizeof( fft_tables->cfg ) / sizeof( fft_tables->cfg[0] ); i++ )
{
if ( fft_tables->cfg[i][0] )
{
free( fft_tables->cfg[i][0] );
fft_tables->cfg[i][0] = NULL;
}
if ( fft_tables->cfg[i][1] )
{
free( fft_tables->cfg[i][1] );
fft_tables->cfg[i][1] = NULL;
}
}
}
void rfft( FFT_Tables *fft_tables, double *x, int logm )
{
/* sur: use real-only optimized FFT */
int nfft = 0;
kiss_fft_scalar fin[1 << MAXLOGR];
kiss_fft_cpx fout[1 << MAXLOGR];
if ( logm > MAXLOGR )
{
fprintf(stderr, "fft size too big\n");
exit(1);
}
nfft = logm_to_nfft[logm];
if ( fft_tables->cfg[logm][0] == NULL )
{
if ( nfft )
{
fft_tables->cfg[logm][0] = kiss_fftr_alloc( nfft, 0, NULL, NULL );
}
else
{
fprintf(stderr, "bad logm = %d\n", logm);
exit( 1 );
}
}
if ( fft_tables->cfg[logm][0] )
{
unsigned int i;
for ( i = 0; i < nfft; i++ )
{
fin[i] = x[i];
}
kiss_fftr( (kiss_fftr_cfg)fft_tables->cfg[logm][0], fin, fout );
for ( i = 0; i < nfft / 2; i++ )
{
x[i] = fout[i].r;
x[i + nfft / 2] = fout[i].i;
}
}
else
{
fprintf( stderr, "bad config for logm = %d\n", logm);
exit( 1 );
}
}
void fft( FFT_Tables *fft_tables, double *xr, double *xi, int logm )
{
int nfft = 0;
kiss_fft_cpx fin[1 << MAXLOGM];
kiss_fft_cpx fout[1 << MAXLOGM];
if ( logm > MAXLOGM )
{
fprintf(stderr, "fft size too big\n");
exit(1);
}
nfft = logm_to_nfft[logm];
if ( fft_tables->cfg[logm][0] == NULL )
{
if ( nfft )
{
fft_tables->cfg[logm][0] = kiss_fft_alloc( nfft, 0, NULL, NULL );
}
else
{
fprintf(stderr, "bad logm = %d\n", logm);
exit( 1 );
}
}
if ( fft_tables->cfg[logm][0] )
{
unsigned int i;
for ( i = 0; i < nfft; i++ )
{
fin[i].r = xr[i];
fin[i].i = xi[i];
}
kiss_fft( (kiss_fft_cfg)fft_tables->cfg[logm][0], fin, fout );
for ( i = 0; i < nfft; i++ )
{
xr[i] = fout[i].r;
xi[i] = fout[i].i;
}
}
else
{
fprintf( stderr, "bad config for logm = %d\n", logm);
exit( 1 );
}
}
void ffti( FFT_Tables *fft_tables, double *xr, double *xi, int logm )
{
int nfft = 0;
kiss_fft_cpx fin[1 << MAXLOGM];
kiss_fft_cpx fout[1 << MAXLOGM];
if ( logm > MAXLOGM )
{
fprintf(stderr, "fft size too big\n");
exit(1);
}
nfft = logm_to_nfft[logm];
if ( fft_tables->cfg[logm][1] == NULL )
{
if ( nfft )
{
fft_tables->cfg[logm][1] = kiss_fft_alloc( nfft, 1, NULL, NULL );
}
else
{
fprintf(stderr, "bad logm = %d\n", logm);
exit( 1 );
}
}
if ( fft_tables->cfg[logm][1] )
{
unsigned int i;
double fac = 1.0 / (double)nfft;
for ( i = 0; i < nfft; i++ )
{
fin[i].r = xr[i];
fin[i].i = xi[i];
}
kiss_fft( (kiss_fft_cfg)fft_tables->cfg[logm][1], fin, fout );
for ( i = 0; i < nfft; i++ )
{
xr[i] = fout[i].r * fac;
xi[i] = fout[i].i * fac;
}
}
else
{
fprintf( stderr, "bad config for logm = %d\n", logm);
exit( 1 );
}
}
#else /* !defined DRM || defined DRM_1024 */
void fft_initialize( FFT_Tables *fft_tables )
{
int i;
fft_tables->costbl = AllocMemory( (MAXLOGM+1) * sizeof( fft_tables->costbl[0] ) );
fft_tables->negsintbl = AllocMemory( (MAXLOGM+1) * sizeof( fft_tables->negsintbl[0] ) );
fft_tables->reordertbl = AllocMemory( (MAXLOGM+1) * sizeof( fft_tables->reordertbl[0] ) );
for( i = 0; i< MAXLOGM+1; i++ )
{
fft_tables->costbl[i] = NULL;
fft_tables->negsintbl[i] = NULL;
fft_tables->reordertbl[i] = NULL;
}
}
void fft_terminate( FFT_Tables *fft_tables )
{
int i;
for( i = 0; i< MAXLOGM+1; i++ )
{
if( fft_tables->costbl[i] != NULL )
FreeMemory( fft_tables->costbl[i] );
if( fft_tables->negsintbl[i] != NULL )
FreeMemory( fft_tables->negsintbl[i] );
if( fft_tables->reordertbl[i] != NULL )
FreeMemory( fft_tables->reordertbl[i] );
}
FreeMemory( fft_tables->costbl );
FreeMemory( fft_tables->negsintbl );
FreeMemory( fft_tables->reordertbl );
fft_tables->costbl = NULL;
fft_tables->negsintbl = NULL;
fft_tables->reordertbl = NULL;
}
static void reorder( FFT_Tables *fft_tables, double *x, int logm)
{
int i;
int size = 1 << logm;
unsigned short *r; //size
if ( fft_tables->reordertbl[logm] == NULL ) // create bit reversing table
{
fft_tables->reordertbl[logm] = AllocMemory(size * sizeof(*(fft_tables->reordertbl[0])));
for (i = 0; i < size; i++)
{
int reversed = 0;
int b0;
int tmp = i;
for (b0 = 0; b0 < logm; b0++)
{
reversed = (reversed << 1) | (tmp & 1);
tmp >>= 1;
}
fft_tables->reordertbl[logm][i] = reversed;
}
}
r = fft_tables->reordertbl[logm];
for (i = 0; i < size; i++)
{
int j = r[i];
double tmp;
if (j <= i)
continue;
tmp = x[i];
x[i] = x[j];
x[j] = tmp;
}
}
static void fft_proc(
double *xr,
double *xi,
fftfloat *refac,
fftfloat *imfac,
int size)
{
int step, shift, pos;
int exp, estep;
estep = size;
for (step = 1; step < size; step *= 2)
{
int x1;
int x2 = 0;
estep >>= 1;
for (pos = 0; pos < size; pos += (2 * step))
{
x1 = x2;
x2 += step;
exp = 0;
for (shift = 0; shift < step; shift++)
{
double v2r, v2i;
v2r = xr[x2] * refac[exp] - xi[x2] * imfac[exp];
v2i = xr[x2] * imfac[exp] + xi[x2] * refac[exp];
xr[x2] = xr[x1] - v2r;
xr[x1] += v2r;
xi[x2] = xi[x1] - v2i;
xi[x1] += v2i;
exp += estep;
x1++;
x2++;
}
}
}
}
static void check_tables( FFT_Tables *fft_tables, int logm)
{
if( fft_tables->costbl[logm] == NULL )
{
int i;
int size = 1 << logm;
if( fft_tables->negsintbl[logm] != NULL )
FreeMemory( fft_tables->negsintbl[logm] );
fft_tables->costbl[logm] = AllocMemory((size / 2) * sizeof(*(fft_tables->costbl[0])));
fft_tables->negsintbl[logm] = AllocMemory((size / 2) * sizeof(*(fft_tables->negsintbl[0])));
for (i = 0; i < (size >> 1); i++)
{
double theta = 2.0 * M_PI * ((double) i) / (double) size;
fft_tables->costbl[logm][i] = cos(theta);
fft_tables->negsintbl[logm][i] = -sin(theta);
}
}
}
void fft( FFT_Tables *fft_tables, double *xr, double *xi, int logm)
{
if (logm > MAXLOGM)
{
fprintf(stderr, "fft size too big\n");
exit(1);
}
if (logm < 1)
{
//printf("logm < 1\n");
return;
}
check_tables( fft_tables, logm);
reorder( fft_tables, xr, logm);
reorder( fft_tables, xi, logm);
fft_proc( xr, xi, fft_tables->costbl[logm], fft_tables->negsintbl[logm], 1 << logm );
}
void rfft( FFT_Tables *fft_tables, double *x, int logm)
{
double xi[1 << MAXLOGR];
if (logm > MAXLOGR)
{
fprintf(stderr, "rfft size too big\n");
exit(1);
}
memset(xi, 0, (1 << logm) * sizeof(xi[0]));
fft( fft_tables, x, xi, logm);
memcpy(x + (1 << (logm - 1)), xi, (1 << (logm - 1)) * sizeof(*x));
}
void ffti( FFT_Tables *fft_tables, double *xr, double *xi, int logm)
{
int i, size;
double fac;
double *xrp, *xip;
fft( fft_tables, xi, xr, logm);
size = 1 << logm;
fac = 1.0 / size;
xrp = xr;
xip = xi;
for (i = 0; i < size; i++)
{
*xrp++ *= fac;
*xip++ *= fac;
}
}
#endif /* defined DRM && !defined DRM_1024 */