path: root/src/loris/FourierTransform.C

/*
 * This is the Loris C++ Class Library, implementing analysis, 
 * manipulation, and synthesis of digitized sounds using the Reassigned 
 * Bandwidth-Enhanced Additive Sound Model.
 *
 * Loris is Copyright (c) 1999-2010 by Kelly Fitz and Lippold Haken
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or
 * (at your option) any later version.
 *
 * This program 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 General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 *
 *
 * FourierTransform.C
 *
 * Implementation of class Loris::FourierTransform, providing a simplified
 * uniform interface to the FFTW library (www.fftw.org), version 2.1.3
 * or newer (including version 3), or to the General Purpose FFT package
 * by Takuya OOURA, http://momonga.t.u-tokyo.ac.jp/~ooura/fft.html if
 * FFTW is unavailable. 
 *
 * Kelly Fitz, 2 Jun 2006
 * loris@cerlsoundgroup.org
 *
 * http://www.cerlsoundgroup.org/Loris/
 *
 */

#if HAVE_CONFIG_H
	#include "config.h"
#endif

#include "FourierTransform.h"
#include "LorisExceptions.h"
#include "Notifier.h"

#include <cmath>
#include <complex>

#if defined(HAVE_M_PI) && (HAVE_M_PI)
	const double Pi = M_PI;
#else
	const double Pi = std::acos( -1.0 );
#endif

#if defined(HAVE_FFTW3_H) && HAVE_FFTW3_H
    #include <fftw3.h>
#elif defined(HAVE_FFTW_H) && HAVE_FFTW_H
    #include <fftw.h>
#endif

// ---------------------------------------------------------------------------
//	isPO2 - return true if N is a power of two
// ---------------------------------------------------------------------------
//  If out_expon is non-zero, return the exponent in that address.
//
static bool isPO2( unsigned int N, int * out_expon = 0 )
{
    unsigned int M = 1;
    int exp = 0;
    while ( M < N )
    {
        M *= 2;
        ++exp;
    }
    if ( 0 != out_expon && M == N )
    {
        *out_expon = exp;
    }
    return M == N;
}

// ===========================================================================
// No longer matters that fftw and this class use the same floating point
// data format. The insulating implementation class now does its job of
// really insulating clients completely from FFTW, by copying data between
// buffers of std::complex< double > and fftw_complex, rather than 
// relying on those two complex types to have the same memory layout.
// The overhead of copying is probably not significant, compared to 
// the expense of computing the transform. 
//
// about complex math functions for fftw_complex:
//
// These functions are all defined as templates in <complex>.
// Regrettably, they are all implemented using real() and 
// imag() _member_ functions of the template argument, T. 
// If they had instead been implemented in terms of the real()
// and imag() (template) free functions, then I could just specialize
// those two for the fftw complex data type, and the other template
// functions would work. Instead, I have to specialize _all_ of
// those functions that I want to use. I hope this was a learning 
// experience for someone... In the mean time, the alternative I 
// have is to take advantage of the fact that fftw_complex and 
// std::complex<double> have the same footprint, so I can just
// cast back and forth between the two types. Its icky, but it 
// works, and its a lot faster than converting, and more palatable
// than redefining all those operators.
//
// On the subject of brilliant designs, fftw_complex is defined as
// a typedef of an anonymous struct, as in typedef struct {...} fftw_complex,
// so I cannot forward-declare that type.
//
// In other good news, the planning structure got a slight name change
// in version 3, making it even more important to remove all traces of
// FFTW from the FourierTransform class definition.
//
// ===========================================================================

//	begin namespace
namespace Loris {

using std::complex;
using std::vector;

// --- private implementation class ---

// ---------------------------------------------------------------------------
//  FTimpl
//
// Insulating implementation class to insulate clients
// completely from everything about the interaction between
// Loris and FFTW. There is more copying of data between buffers,
// but this is not the expensive part of computing Fourier transforms
// and we don't have to do unsavory things that require std::complex
// and fftw_complex to have the same memory layout (we could even get
// by with single-precision floats in FFTW if necessary). Also got
// rid of lots of shared buffer stuff that just made the implementation
// lots more complicated than necessary. This one is simple, if not
// as memory efficient.
//

#if defined(HAVE_FFTW3_H) && HAVE_FFTW3_H

class FTimpl    //  FFTW version 3
{
private:

	fftw_plan plan;
	FourierTransform::size_type N;
	fftw_complex * ftIn;   
	fftw_complex * ftOut;

public:
   
	// Construct an implementation instance:
	// allocate an input buffer, and an output buffer
	// and make a plan.
	FTimpl( FourierTransform::size_type sz ) : 
	  plan( 0 ), N( sz ), ftIn( 0 ), ftOut( 0 ) 
	{      
		// allocate buffers:
		ftIn = (fftw_complex *)fftw_malloc( sizeof( fftw_complex ) * N );
		ftOut = (fftw_complex *)fftw_malloc( sizeof( fftw_complex ) * N );
		if ( 0 == ftIn || 0 == ftOut )
		{
			fftw_free( ftIn );
			fftw_free( ftOut );
			throw RuntimeError( "cannot allocate Fourier transform buffers" );
		}
	  
		//	create a plan:
		plan = fftw_plan_dft_1d( N, ftIn, ftOut, FFTW_FORWARD, FFTW_ESTIMATE );

		//	verify:
		if ( 0 == plan )
		{
			Throw( RuntimeError, "FourierTransform could not make a (fftw) plan." );
		}
	}
   
	// Destroy the implementation instance:
	// dump the plan.
	~FTimpl( void )
	{
		if ( 0 != plan )
		{
            fftw_destroy_plan( plan );
		}         
		
		fftw_free( ftIn );
		fftw_free( ftOut );
	}
	
	// Copy complex< double >'s from a buffer into ftIn, 
	// the buffer must be as long as ftIn.
	void loadInput( const complex< double > * bufPtr )
	{
		for ( FourierTransform::size_type k = 0; k < N; ++k )
		{
			ftIn[ k ][0] = bufPtr->real();  //  real part
			ftIn[ k ][1] = bufPtr->imag();  //  imaginary part
			++bufPtr;
		}
	}
   
	// Copy complex< double >'s from ftOut into a buffer,
	// which must be as long as ftOut.
	void copyOutput( complex< double > * bufPtr ) const
	{
		for ( FourierTransform::size_type k = 0; k < N; ++k )
		{
			*bufPtr = complex< double >( ftOut[ k ][0], ftOut[ k ][1] );
			++bufPtr;
		}
	}
    
    // Compute a forward transform.
    void forward( void )
    {
        fftw_execute( plan );
    }
    
}; // end of class FTimpl for FFTW version 3

#elif defined(HAVE_FFTW_H) && HAVE_FFTW_H

//	"die hook" for FFTW, which otherwise try to write to a
//	non-existent console.
static void fftw_die_Loris( const char * s )
{
	using namespace std; // exit might be hidden in there
	
	notifier << "The FFTW library used by Loris has encountered a fatal error: " << s << endl;
	exit(EXIT_FAILURE);
}

class FTimpl    //  FFTW version 2
{
private:

	fftw_plan plan;
	FourierTransform::size_type N;
	fftw_complex * ftIn;   
	fftw_complex * ftOut;
   
public:

	// Construct an implementation instance:
	// allocate an input buffer, and an output buffer
	// and make a plan.
	FTimpl( FourierTransform::size_type sz ) : 
	  plan( 0 ), N( sz ), ftIn( 0 ), ftOut( 0 ) 
	{      
		// allocate buffers:
		ftIn = (fftw_complex *)fftw_malloc( sizeof( fftw_complex ) * N );
		ftOut = (fftw_complex *)fftw_malloc( sizeof( fftw_complex ) * N );
		if ( 0 == ftIn || 0 == ftOut )
		{
			fftw_free( ftIn );
			fftw_free( ftOut );
			Throw( RuntimeError, "cannot allocate Fourier transform buffers" );
		}
	  
		//	create a plan:
		plan = fftw_create_plan_specific( N, FFTW_FORWARD, FFTW_ESTIMATE,
                                          ftIn, 1, ftOut, 1 );

		//	verify:
		if ( 0 == plan )
		{
			Throw( RuntimeError, "FourierTransform could not make a (fftw) plan." );
		}

        //	FFTW calls fprintf a lot, which may be a problem in
        //	non-console-enabled applications. Catch fftw_die()
        //	calls by routing the error message to our own Notifier
        //	and exiting, using the function defined above.
        //
        //	(version 2 only)
        fftw_die_hook = fftw_die_Loris;
	}
   
	// Destroy the implementation instance:
	// dump the plan.
	~FTimpl( void )
	{
		if ( 0 != plan )
		{
            fftw_destroy_plan( plan );
		}         
		
		fftw_free( ftIn );
		fftw_free( ftOut );
	}
	
	// Copy complex< double >'s from a buffer into ftIn, 
	// the buffer must be as long as ftIn.
	void loadInput( const complex< double > * bufPtr )
	{
		for ( FourierTransform::size_type k = 0; k < N; ++k )
		{
			c_re( ftIn[ k ] ) = bufPtr->real();
			c_im( ftIn[ k ] ) = bufPtr->imag();
			++bufPtr;
		}
	}
   
	// Copy complex< double >'s from ftOut into a buffer,
	// which must be as long as ftOut.
	void copyOutput( complex< double > * bufPtr ) const
	{
		for ( FourierTransform::size_type k = 0; k < N; ++k )
		{
			*bufPtr = complex< double >( c_re( ftOut[ k ] ), c_im( ftOut[ k ] ) );
			++bufPtr;
		}
	}
    
    // Compute a forward transform.
    void forward( void )
    {
        fftw_one( plan, ftIn, ftOut );	
    }
    
}; // end of class FTimpl for FFTW version 2

#else

#define SORRY_NO_FFTW  1

//  function prototype, definition in fftsg.c
/* extern "C" void cdft(int, int, double *, int *, double *); */

//  function prototype, definition below
static void slowDFT( double * in, double * out, int N );

//  Uses General Purpose FFT (Fast Fourier/Cosine/Sine Transform) Package
//  by Takuya OOURA, http://momonga.t.u-tokyo.ac.jp/~ooura/fft.html defined
//  in fftsg.c.
//
//  In the event that the size is not a power of two, uses a (very) slow
//  direct DFT computation, defined below. In this case, the workspace
//  array is not used, and the twiddle factor array is used to store the
//  transform result.

class FTimpl    //  platform-neutral stand-alone implementation
{
private:

	double * mTxInOut;      //	input/output buffer for in-place transform                                
	double * mTwiddle;      //	storage for twiddle factors
	int * mWorkspace;		//	workspace storage

	FourierTransform::size_type N;
    
    bool mIsPO2;
   
public:

	// Construct an implementation instance:
	// allocate buffers and workspace, and
	// initialize the twiddle factors.
	FTimpl( FourierTransform::size_type sz ) : 
	  mTxInOut( 0 ), mTwiddle( 0 ), mWorkspace( 0 ), N( sz ), mIsPO2( isPO2( sz ) )
	{      
        mTxInOut = new double[ 2*N ]; 	
            //	input/output buffer for in-place transform
            
        if ( mIsPO2 )
        {    
            mTwiddle = new double[ N/2 ]; 		
                //	storage for twiddle factors
                
            mWorkspace = new int[ 2*int( std::sqrt((double)N) + 0.5 ) ];		
                //	workspace 
                
            if ( 0 == mTxInOut || 0 == mTwiddle || 0 == mWorkspace )
            {
                delete [] mTxInOut;
                delete [] mTwiddle;
                delete [] mWorkspace;
                Throw( RuntimeError, "FourierTransform: could not initialize tranform" );
            }

            mWorkspace[0] = 0;  // first time only, triggers setup    
                                // no need to do it now, it will happen the 
                                // first time a transform is computed
            // cdft_double( 2*N, -1, mTxInOut, mWorkspace, mTwiddle );
        }
        else
        {
            mTwiddle = new double[ 2*N ]; 	
                //	use for result in slowDFT 
                
            if ( 0 == mTxInOut || 0 == mTwiddle )
            {
                delete [] mTxInOut;
                delete [] mTwiddle;
                Throw( RuntimeError, "FourierTransform: could not initialize tranform" );
            }
        }
	}
   
	// Destroy the implementation instance:
	~FTimpl( void )
	{
        delete [] mTxInOut;
        delete [] mTwiddle;
        delete [] mWorkspace;
	}
	
	// Copy complex< double >'s from a buffer into ftIn, 
	// the buffer must be as long as ftIn.
	void loadInput( const complex< double > * bufPtr )
	{
		for ( FourierTransform::size_type k = 0; k < N; ++k )
		{
			mTxInOut[ 2*k ] = bufPtr->real();
			mTxInOut[ 2*k+1 ] = bufPtr->imag();
			++bufPtr;
		}
	}
   
	//  Copy complex< double >'s from ftOut into a buffer,
	//  which must be as long as ftOut. Result might be
    //  stored in the twiddle factor array if this is not
    //  power of two length DFT.
	void copyOutput( complex< double > * bufPtr ) const
	{
        double * result = mTxInOut;
        if ( !mIsPO2 )
        {
            result = mTwiddle;
        }

		for ( FourierTransform::size_type k = 0; k < N; ++k )
		{
			*bufPtr = complex< double >( result[ 2*k ], result[ 2*k+1 ] );
			++bufPtr;
		}
	}
    
    // Compute a forward transform.
    void forward( void )
    {        
        if ( mIsPO2 )
        {
            /* cdft( 2*N, -1, mTxInOut, mWorkspace, mTwiddle ); */
            slowDFT( mTxInOut, mTwiddle, N );
        }
        else
        {
            slowDFT( mTxInOut, mTwiddle, N );
        }
    }
    
}; // end of class platform-neutral stand-alone FTimpl 

#endif

// --- FourierTransform members ---

// ---------------------------------------------------------------------------
//	FourierTransform constructor
// ---------------------------------------------------------------------------
//! Initialize a new FourierTransform of the specified size.
//!
//! \param  len is the length of the transform in samples (the
//!         number of samples in the transform)
//! \throw  RuntimeError if the necessary buffers cannot be 
//!         allocated, or there is an error configuring FFTW.
//
FourierTransform::FourierTransform( size_type len ) :
	_buffer( len ),
	_impl( new FTimpl( len ) )
{
	//	zero:
	std::fill( _buffer.begin(), _buffer.end(), 0. );
}

// ---------------------------------------------------------------------------
//	FourierTransform copy constructor
// ---------------------------------------------------------------------------
//! Initialize a new FourierTransform that is a copy of another,
//! having the same size and the same buffer contents.
//!
//! \param  rhs is the instance to copy
//! \throw  RuntimeError if the necessary buffers cannot be 
//!         allocated, or there is an error configuring FFTW.
//
FourierTransform::FourierTransform( const FourierTransform & rhs ) :
	_buffer( rhs._buffer ),
	_impl( new FTimpl( rhs._buffer.size() ) ) // not copied
{
}

// ---------------------------------------------------------------------------
//	FourierTransform destructor
// ---------------------------------------------------------------------------
//! Free the resources associated with this FourierTransform.
//
FourierTransform::~FourierTransform( void )
{	
   delete _impl;
}

// ---------------------------------------------------------------------------
//	FourierTransform assignment operator
// ---------------------------------------------------------------------------
//! Make this FourierTransform a copy of another, having
//! the same size and buffer contents.
//!
//! \param  rhs is the instance to copy
//! \return a refernce to this instance
//! \throw  RuntimeError if the necessary buffers cannot be 
//!         allocated, or there is an error configuring FFTW.
//
FourierTransform &
FourierTransform::operator=( const FourierTransform & rhs )
{
   if ( this != &rhs )
   {
      _buffer = rhs._buffer;
      
      // The implementation instance is not assigned, 
      // but a new one is created.
      delete _impl;
      _impl = 0;
      _impl = new FTimpl( _buffer.size() );
   }
   
   return *this;
}

// ---------------------------------------------------------------------------
//	size
// ---------------------------------------------------------------------------
//! Return the length of the transform (in samples).
//! 
//! \return the length of the transform in samples.
FourierTransform::size_type 
FourierTransform::size( void ) const 
{ 
   return _buffer.size(); 
}
	
// ---------------------------------------------------------------------------
//	transform
// ---------------------------------------------------------------------------
//! Compute the Fourier transform of the samples stored in the 
//! transform buffer. The samples stored in the transform buffer
//! (accessed by index or by iterator) are replaced by the 
//! transformed samples, in-place. 
//
void
FourierTransform::transform( void )
{
    // copy data into the transform input buffer:
    _impl->loadInput( &_buffer.front() );
    
    //	crunch:	
    _impl->forward();
    
    // copy the data out of the transform output buffer:
    _impl->copyOutput( &_buffer.front() );
}


// --- slow non-power-of-two DFT implementation ---

#if defined(SORRY_NO_FFTW) 

// ---------------------------------------------------------------------------
//	slowDFT
// ---------------------------------------------------------------------------
//  Non-power-of-two DFT implementation. in and out cannot be the same,
//  and each is 2*N long, storing interleaved complex numbers.
//  This version is only used when FFTW is unavailable.
//
static
void slowDFT( double * in, double * out, int N )
{
#if 1 
    // slow DFT 
    // This version of the direct DFT computation is tolerably
    // speedy, even for rather long transforms (like 8k).
    // There is only one expensive call to std::polar, and
    // twiddle factors are updated by multiplying. This
    // causes some loss in numerical accuracy, worse for
    // longer transforms, but even for 10k long transforms,
    // accuracy is within hundredths of a percent in my experiments.
    
    const std::complex< double > eminj2pioN = 
        std::polar( 1.0, -2.0 * Pi / N );
              
    std::complex< double > Wn = 1;
    for ( int n = 0; n < N; ++n )
    {
        std::complex< double > Wkn = 1;
        std::complex< double > Xkn = 0;
        for ( int k = 0; k < N; ++k )
        {
            Xkn += std::complex< double >( in[ 2*k ], in[ 2*k+1 ] ) * Wkn;
            Wkn *= Wn;
        }
        
        out[ 2*n ] = Xkn.real();
        out[ 2*n+1 ] = Xkn.imag();
        Wn *= eminj2pioN;
    }

#else 
    // very, very slow
    // This version of the direct DFT computation is slightly
    // more accurate, increasingly so for longer transforms, 
    // but it is so much slower (4-5x) that the small improvement
    // in accuracy is probably not worth the extra wait. For
    // short transforms, the accuracy of both algorithms is
    // very high, and for long transforms, this algorithm is
    // too slow.
    //
    // Both algorithms are retained here, so that this 
    // tradeoff can be re-evaluated if necessary.

    for ( int n = 0; n < N; ++n )
    {
        std::complex< double > Xkn = 0;
        for ( int k = 0; k < N; ++k )
        {
            std::complex< double > Wkn = std::polar( 1.0, -2.0 * Pi * k * n / N );
            Xkn += std::complex< double >( in[ 2*k ], in[ 2*k+1 ] ) * Wkn;
        }
        
        out[ 2*n ] = Xkn.real();
        out[ 2*n+1 ] = Xkn.imag();
    }
    
#endif    
}

#endif  //  defined(SORRY_NO_FFTW)

}	//	end of namespace Loris