/* * 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 "util.h"#include "fft.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 */
/*
* 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 "util.h"
#include "fft.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 */