Start adding Loris files

1 files changed, 641 insertions, 0 deletions

diff --git a/src/loris/FourierTransform.C b/src/loris/FourierTransform.C
new file mode 100644
index 0000000..2a29238
--- /dev/null
+++ b/src/loris/FourierTransform.C
@@ -0,0 +1,641 @@
+/*
+ * 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 );
+        }
+        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