/* * 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 #include #include 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