cannam@0: /* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
cannam@0: 
cannam@0: /*
cannam@0:     QM DSP Library
cannam@0: 
cannam@0:     Centre for Digital Music, Queen Mary, University of London.
cannam@0:     This file is based on Don Cross's public domain FFT implementation.
cannam@0: */
cannam@0: 
cannam@0: #include "FFT.h"
cannam@55: 
cannam@55: #include "maths/MathUtilities.h"
cannam@55: 
cannam@0: #include <cmath>
cannam@0: 
cannam@55: #include <iostream>
cannam@55: 
cannam@65: //#define USE_BUILTIN_FFT 1
cannam@0: 
cannam@64: #ifdef USE_BUILTIN_FFT
cannam@64: 
cannam@64: FFT::FFT(unsigned int n) :
cannam@64:     m_n(n),
cannam@64:     m_private(0)
cannam@0: {
cannam@64:     if( !MathUtilities::isPowerOfTwo(m_n) )
cannam@64:     {
cannam@64:         std::cerr << "ERROR: FFT: Non-power-of-two FFT size "
cannam@64:                   << m_n << " not supported in this implementation"
cannam@64:                   << std::endl;
cannam@64: 	return;
cannam@64:     }
cannam@0: }
cannam@0: 
cannam@0: FFT::~FFT()
cannam@0: {
cannam@0: 
cannam@0: }
cannam@0: 
cannam@64: FFTReal::FFTReal(unsigned int n) :
cannam@64:     m_n(n),
cannam@64:     m_private(0)
cannam@0: {
cannam@64:     m_private = new FFT(m_n);
cannam@64: }
cannam@0: 
cannam@64: FFTReal::~FFTReal()
cannam@64: {
cannam@64:     delete (FFT *)m_private;
cannam@64: }
cannam@64: 
cannam@64: void
cannam@64: FFTReal::process(bool inverse,
cannam@64:                  const double *realIn,
cannam@64:                  double *realOut, double *imagOut)
cannam@64: {
cannam@64:     ((FFT *)m_private)->process(inverse, realIn, 0, realOut, imagOut);
cannam@64: }
cannam@64: 
cannam@64: static unsigned int numberOfBitsNeeded(unsigned int p_nSamples)
cannam@64: {	
cannam@64:     int i;
cannam@64: 
cannam@64:     if( p_nSamples < 2 )
cannam@64:     {
cannam@64: 	return 0;
cannam@64:     }
cannam@64: 
cannam@64:     for ( i=0; ; i++ )
cannam@64:     {
cannam@64: 	if( p_nSamples & (1 << i) ) return i;
cannam@64:     }
cannam@64: }
cannam@64: 
cannam@64: static unsigned int reverseBits(unsigned int p_nIndex, unsigned int p_nBits)
cannam@64: {
cannam@64:     unsigned int i, rev;
cannam@64: 
cannam@64:     for(i=rev=0; i < p_nBits; i++)
cannam@64:     {
cannam@64: 	rev = (rev << 1) | (p_nIndex & 1);
cannam@64: 	p_nIndex >>= 1;
cannam@64:     }
cannam@64: 
cannam@64:     return rev;
cannam@64: }
cannam@64: 
cannam@64: void
cannam@64: FFT::process(bool p_bInverseTransform,
cannam@64:              const double *p_lpRealIn, const double *p_lpImagIn,
cannam@64:              double *p_lpRealOut, double *p_lpImagOut)
cannam@64: {
cannam@66:     if (!p_lpRealIn || !p_lpRealOut || !p_lpImagOut) return;
cannam@66: 
cannam@66: //    std::cerr << "FFT::process(" << m_n << "," << p_bInverseTransform << ")" << std::endl;
cannam@0: 
cannam@0:     unsigned int NumBits;
cannam@0:     unsigned int i, j, k, n;
cannam@0:     unsigned int BlockSize, BlockEnd;
cannam@0: 
cannam@0:     double angle_numerator = 2.0 * M_PI;
cannam@0:     double tr, ti;
cannam@0: 
cannam@64:     if( !MathUtilities::isPowerOfTwo(m_n) )
cannam@0:     {
cannam@55:         std::cerr << "ERROR: FFT::process: Non-power-of-two FFT size "
cannam@64:                   << m_n << " not supported in this implementation"
cannam@55:                   << std::endl;
cannam@0: 	return;
cannam@0:     }
cannam@0: 
cannam@0:     if( p_bInverseTransform ) angle_numerator = -angle_numerator;
cannam@0: 
cannam@64:     NumBits = numberOfBitsNeeded ( m_n );
cannam@0: 
cannam@0: 
cannam@64:     for( i=0; i < m_n; i++ )
cannam@0:     {
cannam@0: 	j = reverseBits ( i, NumBits );
cannam@0: 	p_lpRealOut[j] = p_lpRealIn[i];
cannam@0: 	p_lpImagOut[j] = (p_lpImagIn == 0) ? 0.0 : p_lpImagIn[i];
cannam@0:     }
cannam@0: 
cannam@0: 
cannam@0:     BlockEnd = 1;
cannam@64:     for( BlockSize = 2; BlockSize <= m_n; BlockSize <<= 1 )
cannam@0:     {
cannam@0: 	double delta_angle = angle_numerator / (double)BlockSize;
cannam@0: 	double sm2 = -sin ( -2 * delta_angle );
cannam@0: 	double sm1 = -sin ( -delta_angle );
cannam@0: 	double cm2 = cos ( -2 * delta_angle );
cannam@0: 	double cm1 = cos ( -delta_angle );
cannam@0: 	double w = 2 * cm1;
cannam@0: 	double ar[3], ai[3];
cannam@0: 
cannam@64: 	for( i=0; i < m_n; i += BlockSize )
cannam@0: 	{
cannam@0: 
cannam@0: 	    ar[2] = cm2;
cannam@0: 	    ar[1] = cm1;
cannam@0: 
cannam@0: 	    ai[2] = sm2;
cannam@0: 	    ai[1] = sm1;
cannam@0: 
cannam@0: 	    for ( j=i, n=0; n < BlockEnd; j++, n++ )
cannam@0: 	    {
cannam@0: 
cannam@0: 		ar[0] = w*ar[1] - ar[2];
cannam@0: 		ar[2] = ar[1];
cannam@0: 		ar[1] = ar[0];
cannam@0: 
cannam@0: 		ai[0] = w*ai[1] - ai[2];
cannam@0: 		ai[2] = ai[1];
cannam@0: 		ai[1] = ai[0];
cannam@0: 
cannam@0: 		k = j + BlockEnd;
cannam@0: 		tr = ar[0]*p_lpRealOut[k] - ai[0]*p_lpImagOut[k];
cannam@0: 		ti = ar[0]*p_lpImagOut[k] + ai[0]*p_lpRealOut[k];
cannam@0: 
cannam@0: 		p_lpRealOut[k] = p_lpRealOut[j] - tr;
cannam@0: 		p_lpImagOut[k] = p_lpImagOut[j] - ti;
cannam@0: 
cannam@0: 		p_lpRealOut[j] += tr;
cannam@0: 		p_lpImagOut[j] += ti;
cannam@0: 
cannam@0: 	    }
cannam@0: 	}
cannam@0: 
cannam@0: 	BlockEnd = BlockSize;
cannam@0: 
cannam@0:     }
cannam@0: 
cannam@0: 
cannam@0:     if( p_bInverseTransform )
cannam@0:     {
cannam@64: 	double denom = (double)m_n;
cannam@0: 
cannam@64: 	for ( i=0; i < m_n; i++ )
cannam@0: 	{
cannam@0: 	    p_lpRealOut[i] /= denom;
cannam@0: 	    p_lpImagOut[i] /= denom;
cannam@0: 	}
cannam@0:     }
cannam@0: }
cannam@0: 
cannam@64: #else
cannam@0: 
cannam@64: #include "kissfft/kiss_fft.h"
cannam@64: #include "kissfft/kiss_fftr.h"
cannam@0: 
cannam@65: struct KissFFTRec {
cannam@65:     kiss_fft_cfg forward;
cannam@65:     kiss_fft_cfg inverse;
cannam@65:     kiss_fft_cpx *in;
cannam@65:     kiss_fft_cpx *out;
cannam@65: };
cannam@65: 
cannam@65: FFT::FFT(unsigned int n) :
cannam@65:     m_n(n),
cannam@65:     m_private(0)
cannam@65: {
cannam@65:     KissFFTRec *rec = new KissFFTRec;
cannam@65:     rec->forward = kiss_fft_alloc(m_n, 0, 0, 0);
cannam@65:     rec->inverse = kiss_fft_alloc(m_n, 1, 0, 0);
cannam@65:     rec->in = new kiss_fft_cpx[m_n];
cannam@65:     rec->out = new kiss_fft_cpx[m_n];
cannam@65:     m_private = rec;
cannam@65: }
cannam@65: 
cannam@65: FFT::~FFT()
cannam@65: {
cannam@65:     KissFFTRec *rec = (KissFFTRec *)m_private;
cannam@65:     kiss_fft_free(rec->forward);
cannam@65:     kiss_fft_free(rec->inverse);
cannam@65:     delete[] rec->in;
cannam@65:     delete[] rec->out;
cannam@65: }
cannam@65: 
cannam@65: void
cannam@65: FFT::process(bool inverse,
cannam@65:              const double *rin, const double *iin,
cannam@65:              double *rout, double *iout)
cannam@65: {
cannam@65:     KissFFTRec *rec = (KissFFTRec *)m_private;
cannam@65:     for (int i = 0; i < m_n; ++i) {
cannam@65:         rec->in[i].r = rin[i];
cannam@65:     }
cannam@65:     if (iin) {
cannam@65:         for (int i = 0; i < m_n; ++i) {
cannam@65:             rec->in[i].i = iin[i];
cannam@65:         }
cannam@65:     } else {
cannam@65:         for (int i = 0; i < m_n; ++i) {
cannam@65:             rec->in[i].i = 0.0;
cannam@65:         }
cannam@65:     }
cannam@65:     if (inverse) {
cannam@65:         kiss_fft(rec->inverse, rec->in, rec->out);
cannam@65:     } else {
cannam@65:         kiss_fft(rec->forward, rec->in, rec->out);
cannam@65:     }
cannam@65:     for (int i = 0; i < m_n; ++i) {
cannam@65:         rout[i] = rec->out[i].r;
cannam@65:         iout[i] = rec->out[i].i;
cannam@65:     }
cannam@65: }
cannam@65: 
cannam@65: struct KissFFTRealRec {
cannam@65:     kiss_fftr_cfg forward;
cannam@65:     kiss_fftr_cfg inverse;
cannam@65:     kiss_fft_cpx *out;
cannam@65: };
cannam@65: 
cannam@65: FFTReal::FFTReal(unsigned int n) :
cannam@65:     m_n(n),
cannam@65:     m_private(0)
cannam@65: {
cannam@65:     KissFFTRealRec *rec = new KissFFTRealRec;
cannam@65:     rec->forward = kiss_fftr_alloc(m_n, 0, 0, 0);
cannam@65:     rec->inverse = kiss_fftr_alloc(m_n, 1, 0, 0);
cannam@65:     rec->out = new kiss_fft_cpx[m_n];
cannam@65:     m_private = rec;
cannam@65: }
cannam@65: 
cannam@65: FFTReal::~FFTReal()
cannam@65: {
cannam@65:     KissFFTRealRec *rec = (KissFFTRealRec *)m_private;
cannam@65:     kiss_fftr_free(rec->forward);
cannam@65:     kiss_fftr_free(rec->inverse);
cannam@65:     delete[] rec->out;
cannam@65: }
cannam@65: 
cannam@65: void
cannam@65: FFTReal::process(bool inverse,
cannam@65:                  const double *rin,
cannam@65:                  double *rout, double *iout)
cannam@65: {
cannam@65:     KissFFTRealRec *rec = (KissFFTRealRec *)m_private;
cannam@65:     if (inverse) {
cannam@65:         kiss_fftr(rec->inverse, rin, rec->out);
cannam@73:         for (int i = 0; i < m_n; ++i) {
cannam@73:             rout[i] = rec->out[i].r;
cannam@73:             iout[i] = rec->out[i].i;
cannam@73:         }
cannam@65:     } else {
cannam@65:         kiss_fftr(rec->forward, rin, rec->out);
cannam@73:         rout[0] = rec->out[0].r;
cannam@73:         iout[0] = rec->out[0].i;
cannam@73:         for (int i = 1; i < m_n/2; ++i) {
cannam@73:             rout[m_n-i] = rout[i] = rec->out[i].r;
cannam@73:         }
cannam@73:         for (int i = 1; i < m_n/2; ++i) {
cannam@73:             iout[i] = rec->out[i].i;
cannam@73:             iout[m_n-i] = -iout[i];
cannam@73:         }
cannam@73:         rout[m_n/2] = rec->out[m_n/2].r;
cannam@73:         iout[m_n/2] = rec->out[m_n/2].i;
cannam@65:     }
cannam@65: }
cannam@65: 
cannam@64: #endif