Start adding Loris files

1 files changed, 712 insertions, 0 deletions

diff --git a/src/loris/Fundamental.C b/src/loris/Fundamental.C
new file mode 100644
index 0000000..cddb226
--- /dev/null
+++ b/src/loris/Fundamental.C
@@ -0,0 +1,712 @@
+/*
+ * 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
+ *
+ *
+ * Fundamental.C
+ *
+ * Definition of classes for computing an estimate of time-varying
+ * fundamental frequency from either a sequence of samples or a
+ * collection of Partials using a frequency domain maximum likelihood 
+ * algorithm adapted from Quatieri's speech signal processing textbook. 
+ *
+ * Kelly Fitz, 25 March 2008
+ * loris@cerlsoundgroup.org
+ *
+ * http://www.cerlsoundgroup.org/Loris/
+ *
+ */
+
+#if HAVE_CONFIG_H
+	#include "config.h"
+#endif
+
+#include "Fundamental.h"
+
+#include "LorisExceptions.h"
+#include "KaiserWindow.h"
+#include "LinearEnvelope.h"
+#include "Notifier.h"
+#include "PartialUtils.h"
+#include "ReassignedSpectrum.h"
+#include "SpectralPeakSelector.h"
+
+#include "F0Estimate.h" 
+
+#include <algorithm>
+#include <cmath>
+#include <vector>
+
+using namespace std;
+
+#if defined(HAVE_M_PI) && (HAVE_M_PI)
+	const double Pi = M_PI;
+#else
+	const double Pi = 3.14159265358979324;
+#endif
+
+//	begin namespace
+namespace Loris {
+
+
+
+#define VERIFY_ARG(func, test)											\
+	do {																\
+		if (!(test)) 													\
+			Throw( Loris::InvalidArgument, #func ": " #test  );			\
+	} while (false)
+
+
+// ---------------------------------------------------------------------------
+// ---------------------------------------------------------------------------
+//  FundamentalEstimator members
+// ---------------------------------------------------------------------------
+// ---------------------------------------------------------------------------
+
+
+//  -- lifecycle --
+
+// ---------------------------------------------------------------------------
+//  constructor (protected)
+// ---------------------------------------------------------------------------
+//! Construct a new estimator with specified precision and
+//! other parameters given default values.
+//!
+//! The specified precision is used to terminate the iterative
+//! estimation procedure. 
+//!
+//! \param precisionHz is the precision in Hz with which the 
+//! fundamental estimates will be made.
+
+FundamentalEstimator::FundamentalEstimator( double precisionHz ) :
+    m_precision( precisionHz ),
+    m_ampFloor( DefaultAmpFloor ),
+    m_ampRange( DefaultAmpRange ),
+    m_freqCeiling( DefaultFreqCeiling )
+{
+	VERIFY_ARG( FundamentalEstimator, precisionHz > 0 );
+}
+
+// ---------------------------------------------------------------------------
+//  destructor
+// ---------------------------------------------------------------------------
+    
+FundamentalEstimator::~FundamentalEstimator( void )
+{
+}
+
+//  -- spectral analysis parameter access --
+
+// ---------------------------------------------------------------------------
+//	ampFloor
+// ---------------------------------------------------------------------------
+//! Return the absolute amplitude threshold in (negative) dB, 
+//! below which spectral peaks will not be considered in the 
+//! estimation of the fundamental (default is 30 dB).            
+double 
+FundamentalEstimator::ampFloor( void ) const
+{ 
+    return m_ampFloor; 
+}
+
+// ---------------------------------------------------------------------------
+//	ampRange
+// ---------------------------------------------------------------------------
+//!	Return the amplitude range in dB, 
+//! relative to strongest peak in a frame, floating
+//! amplitude threshold below which spectral
+//! peaks will not be considered in the estimation of 
+//! the fundamental (default is 30 dB).			
+//
+double 
+FundamentalEstimator::ampRange( void ) const 
+{ 
+    return m_ampRange; 
+}
+
+// ---------------------------------------------------------------------------
+//	freqCeiling
+// ---------------------------------------------------------------------------
+//!	Return the frequency ceiling in Hz, the
+//! frequency threshold above which spectral
+//! peaks will not be considered in the estimation of 
+//! the fundamental (default is 10 kHz).			
+//
+double 
+FundamentalEstimator::freqCeiling( void ) const 
+{ 
+    return m_freqCeiling; 
+}
+
+// ---------------------------------------------------------------------------
+//	precision
+// ---------------------------------------------------------------------------
+//!	Return the precision of the estimate in Hz, the
+//! fundamental frequency will be estimated to 
+//! within this range (default is 0.1 Hz).
+//
+double 
+FundamentalEstimator::precision( void ) const 
+{ 
+    return m_precision; 
+}
+
+                       
+//  -- spectral analysis parameter mutation --
+
+// ---------------------------------------------------------------------------
+//	setAmpFloor
+// ---------------------------------------------------------------------------
+//! Set the absolute amplitude threshold in (negative) dB, 
+//! below which spectral peaks will not be considered in the 
+//! estimation of the fundamental (default is 30 dB).            
+//! 
+//! \param x is the new value of this parameter.            
+void 
+FundamentalEstimator::setAmpFloor( double x )
+{ 
+	VERIFY_ARG( setAmpFloor, x < 0 );
+    m_ampFloor = x; 
+}
+
+// ---------------------------------------------------------------------------
+//	setAmpRange
+// ---------------------------------------------------------------------------
+//!	Set the amplitude range in dB, 
+//! relative to strongest peak in a frame, floating
+//! amplitude threshold (negative) below which spectral
+//! peaks will not be considered in the estimation of 
+//! the fundamental (default is 30 dB).	
+//!	
+//!	\param x is the new value of this parameter. 				
+//
+void 
+FundamentalEstimator::setAmpRange( double x ) 
+{ 
+	VERIFY_ARG( setAmpRange, x > 0 );
+    m_ampRange = x; 
+}
+
+// ---------------------------------------------------------------------------
+//	setFreqCeiling
+// ---------------------------------------------------------------------------
+//!	Set the frequency ceiling in Hz, the
+//! frequency threshold above which spectral
+//! peaks will not be considered in the estimation of 
+//! the fundamental (default is 10 kHz). Must be
+//! greater than the lower bound.
+//!	
+//!	\param x is the new value of this parameter. 				
+//
+void 
+FundamentalEstimator::setFreqCeiling( double x ) 
+{ 
+    m_freqCeiling = x; 
+}
+
+// ---------------------------------------------------------------------------
+//	setPrecision
+// ---------------------------------------------------------------------------
+//!	Set the precision of the estimate in Hz, the
+//! fundamental frequency will be estimated to 
+//! within this range (default is 0.1 Hz).
+//!	
+//!	\param x is the new value of this parameter. 				
+//
+void 
+FundamentalEstimator::setPrecision( double x ) 
+{ 
+	VERIFY_ARG( setPrecision, x > 0 );
+    m_precision = x; 
+}
+
+
+// ---------------------------------------------------------------------------
+// ---------------------------------------------------------------------------
+//  FundamentalFromSamples members
+// ---------------------------------------------------------------------------
+// ---------------------------------------------------------------------------
+
+//  -- lifecycle --
+
+// ---------------------------------------------------------------------------
+//  constructor
+// ---------------------------------------------------------------------------
+//! Construct a new estimator configured with the given  
+//! analysis window width (main lobe, zero-to-zero). All other 
+//! spectrum analysis parameters are computed from the specified 
+//! window width. 
+//!
+//! The specified precision is used to terminate the iterative
+//! estimation procedure. If unspecified, the default value,
+//! DefaultPrecisionOver100 * 100 is used. 
+//! 
+//! \param windowWidthHz is the main lobe width of the Kaiser
+//! analysis window in Hz.
+//!
+//! \param precisionHz is the precision in Hz with which the 
+//! fundamental estimates will be made.
+
+FundamentalFromSamples::FundamentalFromSamples( double winWidthHz, 
+                                                double precisionHz ) :
+    m_cacheSampleRate( 0 ),
+    m_windowWidth( winWidthHz ), 
+    FundamentalEstimator( precisionHz )
+{
+	VERIFY_ARG( FundamentalFromSamples, winWidthHz > 0 );
+}
+
+// ---------------------------------------------------------------------------
+//  destructor
+// ---------------------------------------------------------------------------
+    
+FundamentalFromSamples::~FundamentalFromSamples( void )
+{
+}
+
+//  -- fundamental frequency estimation --
+
+
+
+// ---------------------------------------------------------------------------
+//  buildEnvelope
+// ---------------------------------------------------------------------------
+//! Construct a linear envelope from fundamental frequency 
+//! estimates taken at the specified interval in seconds.
+//! 
+
+LinearEnvelope 
+FundamentalFromSamples::buildEnvelope( const double * sampsBeg, 
+                                       const double * sampsEnd, 
+                                       double sampleRate, 
+                                       double tbeg, double tend, 
+                                       double interval,
+                                       double lowerFreqBound, double upperFreqBound, 
+                                       double confidenceThreshold )
+{
+    //  sanity check
+    if ( tbeg > tend )
+    {
+        std::swap( tbeg, tend );
+    }
+
+    LinearEnvelope env;
+    
+    std::vector< double > amplitudes, frequencies;
+
+    double time = tbeg;
+    while ( time < tend )
+    {
+        collectFreqsAndAmps( sampsBeg, sampsEnd-sampsBeg, sampleRate,
+                             frequencies, amplitudes, time );
+        if ( ! amplitudes.empty() )
+        {
+            F0Estimate est( amplitudes, frequencies, lowerFreqBound, upperFreqBound, 
+                            m_precision );
+
+            if ( est.confidence() >= confidenceThreshold )
+            {   
+                env.insert( time, est.frequency() );
+            }
+        }
+        
+        time += interval;
+    }
+    
+    return env;            
+}                                       
+                             
+                             
+// ---------------------------------------------------------------------------
+//  estimateAt
+// ---------------------------------------------------------------------------
+//! Return an estimate of the fundamental frequency computed 
+//! at the specified time. 
+
+FundamentalFromSamples::value_type 
+FundamentalFromSamples::estimateAt( const double * sampsBeg, 
+                                    const double * sampsEnd, 
+                                    double sampleRate, 
+                                    double time,
+                                    double lowerFreqBound, double upperFreqBound )
+{
+    std::vector< double > amplitudes, frequencies;
+    
+    collectFreqsAndAmps( sampsBeg, sampsEnd-sampsBeg, sampleRate,
+                         frequencies, amplitudes, time );
+                         
+    F0Estimate est( amplitudes, frequencies, lowerFreqBound, upperFreqBound, m_precision );
+
+    return est;
+}                                    
+
+//  -- spectral analysis parameter access/mutation --
+
+// ---------------------------------------------------------------------------
+//  windowWidth
+// ---------------------------------------------------------------------------
+//! Return the frequency-domain main lobe width (in Hz) (measured between 
+//! zero-crossings) of the analysis window used in spectral
+//! analysis.               
+double FundamentalFromSamples::windowWidth( void ) const
+{
+    return m_windowWidth;
+}
+
+// ---------------------------------------------------------------------------
+//  setWindowWidth
+// ---------------------------------------------------------------------------
+//! Set the frequency-domain main lobe width (in Hz) (measured between 
+//! zero-crossings) of the analysis window used in spectral
+//! analysis.   
+//! 
+//! \param x is the new value of this parameter.            
+void FundamentalFromSamples::setWindowWidth( double x )
+{
+	VERIFY_ARG( setWindowWidth, x > 0 );
+    m_windowWidth = x; 
+}
+
+
+
+//  -- private auxiliary functions --
+
+// ---------------------------------------------------------------------------
+//  buildSpectrumAnalyzer
+// ---------------------------------------------------------------------------
+//! Construct the ReassignedSpectrum that will be used to perform
+//! spectral analysis from which peak frequencies and amplitudes 
+//! will be drawn. This construction is performed in a lazy fashion,
+//! and needs to be done again when certain of the parameters change.
+//!
+//! The spectrum analyzer cannot be constructed without knowledge of
+//! the sample rate, specified in Hz, which is needed to determine the
+//! parameters of the analysis window. (The sample rate is cached in
+//! this class in order that it be possible to determine whether the
+//! spectrum analyzer can be reused from one estimate to another.)
+//
+void 
+FundamentalFromSamples::buildSpectrumAnalyzer( double srate )
+{
+ 	//	configure the reassigned spectral analyzer, 
+    //	always use odd-length windows:
+    const double sidelobeLevel = - m_ampFloor; // amp floor is negative
+    double winshape = KaiserWindow::computeShape( sidelobeLevel );
+    long winlen = KaiserWindow::computeLength( m_windowWidth / srate, winshape );    
+    if ( 1 != (winlen % 2) ) 
+    {
+        ++winlen;
+    }
+    
+    std::vector< double > window( winlen );
+    KaiserWindow::buildWindow( window, winshape );
+    
+    std::vector< double > windowDeriv( winlen );
+    KaiserWindow::buildTimeDerivativeWindow( windowDeriv, winshape );
+   
+    m_spectrum.reset( new ReassignedSpectrum( window, windowDeriv ) );    
+    
+    //  remember the sample rate used to build this spectrum
+    //  analyzer:
+    m_cacheSampleRate = srate;
+}
+
+// ---------------------------------------------------------------------------
+//	sort_peaks_greater_amplitude
+// ---------------------------------------------------------------------------
+//	predicate used for sorting peaks in order of decreasing amplitude:
+static bool sort_peaks_greater_amplitude( const SpectralPeak & lhs, 
+										  const SpectralPeak & rhs )
+{ 
+	return lhs.amplitude() > rhs.amplitude(); 
+}
+
+// ---------------------------------------------------------------------------
+//  collectFreqsAndAmps
+// ---------------------------------------------------------------------------
+//! Perform spectral analysis on a sequence of samples, using
+//! an analysis window centered at the specified time in seconds.
+//! Collect the frequencies and amplitudes of the peaks and return
+//! them in the vectors provided. 
+//
+   
+void 
+FundamentalFromSamples::collectFreqsAndAmps( const double * samps,
+                                             unsigned long nsamps,
+                                             double sampleRate,
+                                             std::vector< double > & frequencies, 
+                                             std::vector< double > & amplitudes,
+                                             double time )
+{
+    amplitudes.clear();
+    frequencies.clear();
+
+    //  build the spectrum analyzer if necessary:
+    if ( m_cacheSampleRate != sampleRate ||
+         0 == m_spectrum.get() )
+    {
+        buildSpectrumAnalyzer( sampleRate );
+    }
+    
+    
+    //	configure the peak selection and partial formation policies:    
+    unsigned long winlen = m_spectrum->window().size();
+    const double maxTimeCorrection = 0.25 * winlen / sampleRate;   //  one-quarter the window width
+    SpectralPeakSelector selector( sampleRate, maxTimeCorrection );  
+ 	
+    
+     
+    //	compute reassigned spectrum:
+    //  sampsBegin is the position of the first sample to be transformed,
+    //	sampsEnd is the position after the last sample to be transformed.
+    //	(these computations work for odd length windows only)
+    unsigned long winMiddle = (unsigned long)( sampleRate * time );
+    unsigned long sampsBegin = std::max( long(winMiddle) - long(winlen / 2), 0L );
+    unsigned long sampsEnd = std::min( winMiddle + (winlen / 2) + 1, nsamps );
+    
+    if ( winMiddle < sampsEnd )
+    {
+        m_spectrum->transform( samps + sampsBegin, samps + winMiddle, samps + sampsEnd );
+                     
+        //	extract peaks from the spectrum, no fading:
+        Peaks peaks = selector.selectPeaks( *m_spectrum ); 
+        
+        if ( ! peaks.empty() )
+        {
+            //  sort the peaks in order of decreasing amplitude
+            //
+            //  (HEY is there any reason to do this, other than to find the largest?)
+            //std::sort( peaks.begin(), peaks.end(), sort_peaks_greater_amplitude );
+            Peaks::iterator maxpos = std::max_element( peaks.begin(), peaks.end(), sort_peaks_greater_amplitude );
+            
+            //  determine the floating amplitude threshold
+            const double thresh = 
+                std::max( std::pow( 10.0, - 0.05 * - m_ampFloor ), 
+                          std::pow( 10.0, - 0.05 * m_ampRange ) * maxpos->amplitude() );
+                        
+            //  collect amplitudes and frequencies and try to 
+            //  estimate the fundamental
+            for ( Peaks::const_iterator spkpos = peaks.begin(); spkpos != peaks.end(); ++spkpos )
+            {
+                if ( spkpos->amplitude() > thresh &&
+                     spkpos->frequency() < m_freqCeiling )
+                {
+                    amplitudes.push_back( spkpos->amplitude() );
+                    frequencies.push_back( spkpos->frequency() );
+                }
+            }
+        }
+    }
+}
+
+// ---------------------------------------------------------------------------
+// ---------------------------------------------------------------------------
+//  FundamentalFromPartials members
+// ---------------------------------------------------------------------------
+// ---------------------------------------------------------------------------
+
+// ---------------------------------------------------------------------------
+//  constructor
+// ---------------------------------------------------------------------------
+//! Construct a new estimator configured with the given precision.
+//! The specified precision is used to terminate the iterative
+//! estimation procedure. If unspecified, the default value,
+//! DefaultPrecisionOver100 * 100 is used. 
+//!
+//! \param precisionHz is the precision in Hz with which the 
+//! fundamental estimates will be made.
+
+FundamentalFromPartials::FundamentalFromPartials( double precisionHz ) :
+    FundamentalEstimator( precisionHz )
+{
+}
+
+// ---------------------------------------------------------------------------
+//  copy constructor
+// ---------------------------------------------------------------------------
+//! Construct a copy of an estimator. Nothing much to do since this class
+//! has no data members.
+//
+
+FundamentalFromPartials::FundamentalFromPartials( const FundamentalFromPartials & rhs ) :
+    FundamentalEstimator( rhs )
+{
+}
+
+// ---------------------------------------------------------------------------
+//  destructor
+// ---------------------------------------------------------------------------
+    
+FundamentalFromPartials::~FundamentalFromPartials( void )
+{
+}
+
+// ---------------------------------------------------------------------------
+//  assignment
+// ---------------------------------------------------------------------------
+//! Pass the assignment opertion up to the base class.
+//
+
+FundamentalFromPartials &
+FundamentalFromPartials::operator=( const FundamentalFromPartials & rhs )
+{
+    FundamentalEstimator::operator=( rhs );
+    
+    return *this;
+}
+
+//  -- fundamental frequency estimation --
+
+// ---------------------------------------------------------------------------
+//  buildEnvelope
+// ---------------------------------------------------------------------------
+//! Construct a linear envelope from fundamental frequency 
+//! estimates taken at the specified interval in seconds.
+//! 
+
+LinearEnvelope 
+FundamentalFromPartials::buildEnvelope( PartialList::const_iterator begin_partials, 
+                                        PartialList::const_iterator end_partials,
+                                        double tbeg, double tend, 
+                                        double interval,
+                                        double lowerFreqBound, double upperFreqBound, 
+                                        double confidenceThreshold )
+{
+    //  sanity check
+    if ( tbeg > tend )
+    {
+        std::swap( tbeg, tend );
+    }
+
+    LinearEnvelope env;
+    
+    std::vector< double > amplitudes, frequencies;
+
+    double time = tbeg;
+    while ( time < tend )
+    {
+        collectFreqsAndAmps( begin_partials, end_partials, frequencies, amplitudes, time );
+                  
+        if (! amplitudes.empty() )
+        {
+            F0Estimate est( amplitudes, frequencies, lowerFreqBound, upperFreqBound, 
+                            m_precision );
+        
+            if ( est.confidence() >= confidenceThreshold )
+            {   
+                env.insert( time, est.frequency() );
+            }
+        }
+        
+        time += interval;
+    }
+    
+    return env;            
+}                         
+
+// ---------------------------------------------------------------------------
+//  estimateAt
+// ---------------------------------------------------------------------------
+//! Return an estimate of the fundamental frequency computed 
+//! at the specified time. 
+
+FundamentalFromPartials::value_type 
+FundamentalFromPartials::estimateAt( PartialList::const_iterator begin_partials, 
+                                     PartialList::const_iterator end_partials,
+                                     double time,
+                                     double lowerFreqBound, double upperFreqBound )
+{
+    std::vector< double > amplitudes, frequencies;
+    
+    collectFreqsAndAmps( begin_partials, end_partials, frequencies, amplitudes, time );
+                         
+    F0Estimate est( amplitudes, frequencies, lowerFreqBound, upperFreqBound, m_precision );
+
+    return est;
+}  
+
+
+//  -- private auxiliary functions --
+
+// ---------------------------------------------------------------------------
+//  collectFreqsAndAmps
+// ---------------------------------------------------------------------------
+//! Perform spectral analysis on a sequence of samples, using
+//! an analysis window centered at the specified time in seconds.
+//! Collect the frequencies and amplitudes of the peaks and return
+//! them in the vectors provided. 
+//
+   
+void 
+FundamentalFromPartials::collectFreqsAndAmps( PartialList::const_iterator begin_partials, 
+                                              PartialList::const_iterator end_partials,
+                                              std::vector< double > & frequencies, 
+                                              std::vector< double > & amplitudes,
+                                              double time )
+{
+    amplitudes.clear();
+    frequencies.clear();
+    
+    if ( begin_partials != end_partials )
+    {            
+        //  determine the absolute amplitude threshold 
+        double thresh = std::pow( 10.0, - 0.05 * - m_ampFloor );
+        
+        double max_amp = 0;        
+        for ( PartialList::const_iterator it = begin_partials; it != end_partials; ++it )
+        {
+            //  compute the sinusoidal amplitude (without bandwidth energy)
+            double sine_amp = std::sqrt(1 - it->bandwidthAt( time )) * it->amplitudeAt( time );        
+            double freq = it->frequencyAt( time );
+            
+            if ( sine_amp > thresh &&
+                 freq < m_freqCeiling )
+            {
+                amplitudes.push_back( sine_amp );
+                frequencies.push_back( freq );
+            }
+            
+            max_amp = std::max( sine_amp, max_amp );                        
+        }
+        
+        //  remove quietest ones - this isn't very efficient, 
+        //  but it is much faster than making two passes (and 
+        //  computing two sequences of sinusoidal amplitudes).
+        thresh = std::pow( 10.0, - 0.05 * m_ampRange ) * max_amp;
+        vector< double >::size_type N = amplitudes.size();
+        vector< double >::size_type k = 0;
+        while ( k < N )
+        {
+            if ( amplitudes[k] < thresh )
+            {
+                amplitudes.erase( amplitudes.begin() + k );
+                frequencies.erase( frequencies.begin() + k );
+                --N;
+            }
+            else
+            {
+                ++k;
+            }
+        }
+    }
+    
+}
+
+}   //  end of namespace Loris