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 )
{
#if 0
/* sur: do not use real-only optimized FFT */
    double xi[1 << MAXLOGR];

    int nfft;

    if ( logm > MAXLOGR )
	{
		fprintf(stderr, "rfft size too big\n");
		exit(1);
	}

    nfft = logm_to_nfft[logm];

    if ( nfft )
    {
        //unsigned int i;
        //for ( i = 0; i < nfft; i++ )
        //{
        //    xi[i] = 0.0;
        //}
	    memset( xi, 0, nfft * sizeof( xi[0] ) );

	    fft( fft_tables, x, xi, logm );

	    memcpy( x + nfft / 2, xi, ( nfft / 2 ) * sizeof(x[0]) );
    }
    else
    {
        fprintf( stderr, "bad config for logm = %d\n", logm);
        exit( 1 );
    }

#else
/* 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 );
    }
#endif
}

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 );
    }
}

/* sur: Trying to use cfft from libfaad2 -- it does not work 'from scratch' */
//
//#include "cfft/common.h"
//
//void fft_initialize( FFT_Tables *fft_tables )
//{
//    memset( fft_tables->cfft, 0, sizeof( fft_tables->cfft ) );
//}
//void fft_terminate( FFT_Tables *fft_tables )
//{
//    unsigned int i;
//    for ( i = 0; i < sizeof( fft_tables->cfft ) / sizeof( fft_tables->cfft[0] ); i++ )
//    {
//        if ( fft_tables->cfft[i] )
//        {
//            cfftu( fft_tables->cfft[i] );
//            fft_tables->cfft[i] = NULL;
//        }
//    }
//}
//
//void rfft( FFT_Tables *fft_tables, double *x, int logm )
//{
//	double xi[1 << MAXLOGR];
//
//    int nfft;
//
//    if ( logm > MAXLOGR )
//	{
//		fprintf(stderr, "rfft size too big\n");
//		exit(1);
//	}
//
//    nfft = logm_to_nfft[logm];
//
//    if ( nfft )
//    {
//        unsigned int i;
//
//        for ( i = 0; i < nfft; i++ )
//        {
//            xi[i] = 0.0;
//        }
//	    //memset( xi, 0, nfft * sizeof( xi[0] ) );
//
//	    fft( fft_tables, x, xi, logm );
//
//	    memcpy( x + nfft / 2, xi, ( nfft / 2 ) * sizeof(x[0]) );
//
//#ifdef SUR_DEBUG_FFT
//        {
//            FILE* f = fopen( "fft.log", "at" );
//            
//            fprintf( f, "RFFT(%d)\n", nfft );
//            
//            for ( i = 0; i < nfft; i++ )
//            {
//                fprintf( f, ";%d;%g;%g\n", i, x[i], xi[i] );
//            }
//
//            fclose( f );
//        }
//#endif
//    }
//    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;
//
//    complex_t c[1 << MAXLOGM];
//
//    if ( logm > MAXLOGM )
//	{
//		fprintf(stderr, "fft size too big\n");
//		exit(1);
//	}
//
//    nfft = logm_to_nfft[logm];
//
//    if ( fft_tables->cfft[logm] == NULL )
//    {
//        if ( nfft )
//        {
//            fft_tables->cfft[logm] = cffti( nfft );
//        }
//        else
//        {
//	    	fprintf(stderr, "bad logm = %d\n", logm);
//            exit( 1 );
//        }
//    }
//
//    if ( fft_tables->cfft[logm] )
//    {
//        unsigned int i;
//        
//        for ( i = 0; i < nfft; i++ )
//        {
//            RE( c[i] ) = xr[i];    
//            IM( c[i] ) = xi[i];
//        }
//        
//        cfftf( fft_tables->cfft[logm], c );
//
//        for ( i = 0; i < nfft; i++ )
//        {
//            xr[i]   = RE( c[i] );
//            xi[i]   = IM( c[i] );
//        }
//
//#ifdef SUR_DEBUG_FFT
//        {
//            FILE* f = fopen( "fft.log", "at" );
//            
//            fprintf( f, "FFT(%d)\n", nfft );
//            
//            for ( i = 0; i < nfft; i++ )
//            {
//                fprintf( f, ";%d;%g;%g\n", i, xr[i], xi[i] );
//            }
//
//            fclose( f );
//        }
//#endif
//    }
//    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;
//
//    complex_t c[1 << MAXLOGM];
//
//    if ( logm > MAXLOGM )
//	{
//		fprintf(stderr, "fft size too big\n");
//		exit(1);
//	}
//
//    nfft = logm_to_nfft[logm];
//
//    if ( fft_tables->cfft[logm] == NULL )
//    {
//        if ( nfft )
//        {
//            fft_tables->cfft[logm] = cffti( nfft );
//        }
//        else
//        {
//	    	fprintf(stderr, "bad logm = %d\n", logm);
//            exit( 1 );
//        }
//    }
//    
//    if ( fft_tables->cfft[logm] )
//    {
//        unsigned int i;
//        
//        for ( i = 0; i < nfft; i++ )
//        {
//            RE( c[i] )  = xr[i];    
//            IM( c[i] )  = xi[i];
//        }
//        
//        cfftb( fft_tables->cfft[logm], c );
//
//        for ( i = 0; i < nfft; i++ )
//        {
//            xr[i]   = RE( c[i] );
//            xi[i]   = IM( c[i] );
//        }
//
//#ifdef SUR_DEBUG_FFT
//        {
//            FILE* f = fopen( "fft.log", "at" );
//            
//            fprintf( f, "FFTI(%d)\n", nfft );
//            
//            for ( i = 0; i < nfft; i++ )
//            {
//                fprintf( f, ";%d;%g;%g\n", i, xr[i], xi[i] );
//            }
//
//            fclose( f );
//        }
//#endif
//    }
//    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 */

/*
$Log: fft.c,v $
Revision 1.12  2005/02/02 07:49:55  sur
Added interface to kiss_fft library to implement FFT for 960 transform length.

Revision 1.11  2004/04/02 14:56:17  danchr
fix name clash w/ libavcodec: fft_init -> fft_initialize
bump version number to 1.24 beta

Revision 1.10  2003/11/16 05:02:51  stux
moved global tables from fft.c into hEncoder FFT_Tables. Add fft_init and fft_terminate, flowed through all necessary changes. This should remove at least one instance of a memory leak, and fix some thread-safety problems. Version update to 1.23.3

Revision 1.9  2003/09/07 16:48:01  knik
reduced arrays size

Revision 1.8  2002/11/23 17:32:54  knik
rfft: made xi a local variable

Revision 1.7  2002/08/21 16:52:25  knik
new simplier and faster fft routine and correct real fft
new real fft is just a complex fft wrapper so it is slower than optimal but
by surprise it seems to be at least as fast as the old buggy function

*/