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