/*
 * 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
 *
 *
 * Analyzer.C
 *
 * Implementation of class Loris::Analyzer.
 *
 * Kelly Fitz, 5 Dec 99
 * loris@cerlsoundgroup.org
 *
 * http://www.cerlsoundgroup.org/Loris/
 *
 */

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

#include "Analyzer.h"

#include "AssociateBandwidth.h"
#include "Breakpoint.h"
#include "BreakpointEnvelope.h"
#include "Envelope.h"
#include "F0Estimate.h"
#include "LorisExceptions.h"
#include "KaiserWindow.h"
#include "Notifier.h"
#include "Partial.h"
#include "PartialPtrs.h"
#include "ReassignedSpectrum.h"
#include "SpectralPeakSelector.h"
#include "PartialBuilder.h"

#include "phasefix.h"   //  for frequency/phase fixing at end of analysis



#include <algorithm>
#include <cmath>
#include <functional>   //  for std::plus
#include <memory>
#include <numeric>      //  for std::inner_product
#include <utility>
#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 {

// ---------------------------------------------------------------------------
//  helpers, used below
// ---------------------------------------------------------------------------
static double accumPeakSquaredAmps( double init, 
                                    const SpectralPeak & pk )
{
    return init + (pk.amplitude() * pk.amplitude());
}

template < class Pair >
static double compare2nd( const Pair & p1, const Pair & p2 )
{
    return p1.second < p2.second;
}


// ---------------------------------------------------------------------------
//  FundamentalBuilder::build
// ---------------------------------------------------------------------------
//
void FundamentalBuilder::build( const Peaks & peaks, double frameTime )
{
    amplitudes.clear();
    frequencies.clear();
    for ( Peaks::const_iterator spkpos = peaks.begin(); spkpos != peaks.end(); ++spkpos )
    {
        if ( spkpos->amplitude() > mAmpThresh &&
             spkpos->frequency() < mFreqThresh )
        {
            amplitudes.push_back( spkpos->amplitude() );
            frequencies.push_back( spkpos->frequency() );
        }
    }
    if ( ! amplitudes.empty() )
    {
        const double fmin = mFminEnv->valueAt( frameTime );
        const double fmax = mFmaxEnv->valueAt( frameTime );
        
        //  estimate f0
        F0Estimate est( amplitudes, frequencies, fmin, fmax, 0.1 );
        
        if ( est.confidence() >= mMinConfidence &&
             est.frequency() > fmin && est.frequency() < fmax  )
        {
            // notifier << "f0 is " << est.frequency << endl;
            //  add breakpoint to fundamental envelope
            mEnvelope.insert( frameTime, est.frequency() );
        }
    }
    
}


// ---------------------------------------------------------------------------
//  AmpEnvBuilder::build
// ---------------------------------------------------------------------------
//
void AmpEnvBuilder::build( const Peaks & peaks, double frameTime )
{
    double x = std::accumulate( peaks.begin(), peaks.end(), 0.0, accumPeakSquaredAmps );
    mEnvelope.insert( frameTime, std::sqrt( x ) );
}


// ---------------------------------------------------------------------------
//  Analyzer constructor - frequency resolution only
// ---------------------------------------------------------------------------
//! Construct a new Analyzer configured with the given  
//! frequency resolution (minimum instantaneous frequency   
//! difference between Partials). All other Analyzer parameters     
//! are computed from the specified frequency resolution.   
//! 
//! \param resolutionHz is the frequency resolution in Hz.
//
Analyzer::Analyzer( double resolutionHz )
{
    configure( resolutionHz, 2.0 * resolutionHz );
}

// ---------------------------------------------------------------------------
//  Analyzer constructor
// ---------------------------------------------------------------------------
//! Construct a new Analyzer configured with the given  
//! frequency resolution (minimum instantaneous frequency   
//! difference between Partials) and analysis window width
//! (main lobe, zero-to-zero). All other Analyzer parameters    
//! are computed from the specified resolution and window width.    
//! 
//! \param resolutionHz is the frequency resolution in Hz.
//! \param windowWidthHz is the main lobe width of the Kaiser
//! analysis window in Hz.
//
Analyzer::Analyzer( double resolutionHz, double windowWidthHz )
{
    configure( resolutionHz, windowWidthHz );
}

// ---------------------------------------------------------------------------
//  Analyzer constructor
// ---------------------------------------------------------------------------
//! Construct a new Analyzer configured with the given time-varying
//! frequency resolution (minimum instantaneous frequency   
//! difference between Partials) and analysis window width
//! (main lobe, zero-to-zero). All other Analyzer parameters    
//! are computed from the specified resolution and window width.    
//! 
//! \param resolutionHz is the frequency resolution in Hz.
//! \param windowWidthHz is the main lobe width of the Kaiser
//! analysis window in Hz.
//
Analyzer::Analyzer( const Envelope & resolutionEnv, double windowWidthHz )
{
    configure( resolutionEnv, windowWidthHz );
}

// ---------------------------------------------------------------------------
//  Analyzer copy constructor
// ---------------------------------------------------------------------------
//! Construct  a new Analyzer having identical
//! parameter configuration to another Analyzer. 
//! The list of collected Partials is not copied.       
//! 
//! \param other is the Analyzer to copy.   
//
Analyzer::Analyzer( const Analyzer & other ) :
    m_freqResolutionEnv( other.m_freqResolutionEnv->clone() ),
    m_ampFloor( other.m_ampFloor ),
    m_windowWidth( other.m_windowWidth ),
    m_freqFloor( other.m_freqFloor ),
    m_freqDrift( other.m_freqDrift ),
    m_hopTime( other.m_hopTime ),
    m_cropTime( other.m_cropTime ),
    m_bwAssocParam( other.m_bwAssocParam ),
    m_sidelobeLevel( other.m_sidelobeLevel ),
    m_phaseCorrect( other.m_phaseCorrect ),
    m_partials( other.m_partials )
{
    m_f0Builder.reset( other.m_f0Builder->clone() );
    m_ampEnvBuilder.reset( other.m_ampEnvBuilder->clone() );
}

// ---------------------------------------------------------------------------
//  Analyzer assignment
// ---------------------------------------------------------------------------
//! Construct  a new Analyzer having identical
//! parameter configuration to another Analyzer. 
//! The list of collected Partials is not copied.       
//! 
//! \param rhs is the Analyzer to copy. 
//
Analyzer & 
Analyzer::operator=( const Analyzer & rhs )
{
    if ( this != & rhs ) 
    {
        m_freqResolutionEnv.reset( rhs.m_freqResolutionEnv->clone() );
        m_ampFloor = rhs.m_ampFloor;
        m_windowWidth = rhs.m_windowWidth;
        m_freqFloor = rhs.m_freqFloor;  
        m_freqDrift = rhs.m_freqDrift;
        m_hopTime = rhs.m_hopTime;
        m_cropTime = rhs.m_cropTime;
        m_bwAssocParam = rhs.m_bwAssocParam;
        m_sidelobeLevel = rhs.m_sidelobeLevel;
        m_phaseCorrect = rhs.m_phaseCorrect;
        m_partials = rhs.m_partials;

        m_f0Builder.reset( rhs.m_f0Builder->clone() );
        m_ampEnvBuilder.reset( rhs.m_ampEnvBuilder->clone() );
                
    }
    return *this;
}

// ---------------------------------------------------------------------------
//  Analyzer destructor
// ---------------------------------------------------------------------------
//! Destroy this Analyzer.
//
Analyzer::~Analyzer( void )
{
}

// -- configuration --

// ---------------------------------------------------------------------------
//  configure
// ---------------------------------------------------------------------------
//! Configure this Analyzer with the given frequency resolution 
//! (minimum instantaneous frequency difference between Partials, 
//! in Hz). All other Analyzer parameters are (re-)computed from the 
//! frequency resolution, including the window width, which is
//!	twice the resolution.      
//! 
//! \param resolutionHz is the frequency resolution in Hz.
//
void
Analyzer::configure( double resolutionHz )
{
	configure( resolutionHz, 2.0 * resolutionHz );
}

// ---------------------------------------------------------------------------
//  configure
// ---------------------------------------------------------------------------
//! Configure this Analyzer with the given frequency resolution 
//! (minimum instantaneous frequency difference between Partials)
//! and analysis window width (main lobe, zero-to-zero, in Hz). 
//! All other Analyzer parameters are (re-)computed from the 
//! frequency resolution and window width.      
//! 
//! \param resolutionHz is the frequency resolution in Hz.
//! \param windowWidthHz is the main lobe width of the Kaiser
//! analysis window in Hz.
//!     
//! There are three categories of analysis parameters:
//! - the resolution, and params that are usually related to (or
//! identical to) the resolution (frequency floor and drift)
//! - the window width and params that are usually related to (or
//! identical to) the window width (hop and crop times)
//! - independent parameters (bw region width and amp floor)
//
void
Analyzer::configure( double resolutionHz, double windowWidthHz )
{
    //  use specified resolution:
    setFreqResolution( resolutionHz );
    
    //  floor defaults to -90 dB:
    setAmpFloor( -90. );
    
    //  window width should generally be approximately 
    //  equal to, and never more than twice the 
    //  frequency resolution:
    setWindowWidth( windowWidthHz );
    
    //  the Kaiser window sidelobe level can be the same
    //  as the amplitude floor (except in positive dB):
    setSidelobeLevel( - m_ampFloor );
    
    //  for the minimum frequency, below which no data is kept,
    //  use the frequency resolution by default (this makes 
    //  Lip happy, and is always safe?) and allow the client 
    //  to change it to anything at all.
    setFreqFloor( resolutionHz );
    
    //  frequency drift in Hz is the maximum difference
    //  in frequency between consecutive Breakpoints in
    //  a Partial, by default, make it equal to one half
    //  the frequency resolution:
    setFreqDrift( .5 * resolutionHz );
    
    //  hop time (in seconds) is the inverse of the
    //  window width....really. Smith and Serra (1990) cite 
    //  Allen (1977) saying: a good choice of hop is the window 
    //  length divided by the main lobe width in frequency samples,
    //  which turns out to be just the inverse of the width.
    setHopTime( 1. / m_windowWidth );
    
    //  crop time (in seconds) is the maximum allowable time
    //  correction, beyond which a reassigned spectral component
    //  is considered unreliable, and not considered eligible for
    //  Breakpoint formation in extractPeaks(). By default, use
    //  the hop time (should it be half that?):
    setCropTime( m_hopTime );
    
    //  bandwidth association region width 
    //  defaults to 2 kHz, corresponding to 
    //  1 kHz region center spacing:
    storeResidueBandwidth();

	//  configure the envelope builders using default 
	//  parameters:
	buildFundamentalEnv( 0.99 * resolutionHz,
                         1.5 * resolutionHz );
    m_ampEnvBuilder.reset( new AmpEnvBuilder );
    
    //  enable phase-correct Partial construction:
    m_phaseCorrect = true;
}

// ---------------------------------------------------------------------------
//  configure
// ---------------------------------------------------------------------------
//! Configure this Analyzer with the given time-varying frequency resolution 
//! (minimum instantaneous frequency difference between Partials)
//! and analysis window width (main lobe, zero-to-zero, in Hz). 
//! All other Analyzer parameters are (re-)computed from the 
//! frequency resolution and window width.      
//! 
//! \param resolutionEnv is the time-varying frequency resolution 
//!	in Hz.
//! \param windowWidthHz is the main lobe width of the Kaiser
//! analysis window in Hz.
//!     
//! There are three categories of analysis parameters:
//! - the resolution, and params that are usually related to (or
//! identical to) the resolution (frequency floor and drift)
//! - the window width and params that are usually related to (or
//! identical to) the window width (hop and crop times)
//! - independent parameters (bw region width and amp floor)
//
void
Analyzer::configure( const Envelope & resolutionEnv, double windowWidthHz )
{
    //  use specified resolution:
    setFreqResolution( resolutionEnv );
    
    //  floor defaults to -90 dB:
    setAmpFloor( -90. );
    
    //  window width should generally be approximately 
    //  equal to, and never more than twice the 
    //  frequency resolution:
    setWindowWidth( windowWidthHz );
    
    //  the Kaiser window sidelobe level can be the same
    //  as the amplitude floor (except in positive dB):
    setSidelobeLevel( - m_ampFloor );
    
    //  for the minimum frequency, below which no data is kept,
    //  use the frequency resolution by default (this makes 
    //  Lip happy, and is always safe?) and allow the client 
    //  to change it to anything at all.
    setFreqFloor( windowWidthHz * 0.5 );		//	!!!!!
    
    //  frequency drift in Hz is the maximum difference
    //  in frequency between consecutive Breakpoints in
    //  a Partial, by default, make it equal to one half
    //  the frequency resolution:
    setFreqDrift( windowWidthHz * 0.25 );		//	!!!!!
    
    //  hop time (in seconds) is the inverse of the
    //  window width....really. Smith and Serra (1990) cite 
    //  Allen (1977) saying: a good choice of hop is the window 
    //  length divided by the main lobe width in frequency samples,
    //  which turns out to be just the inverse of the width.
    setHopTime( 1. / m_windowWidth );
    
    //  crop time (in seconds) is the maximum allowable time
    //  correction, beyond which a reassigned spectral component
    //  is considered unreliable, and not considered eligible for
    //  Breakpoint formation in extractPeaks(). By default, use
    //  the hop time (should it be half that?):
    setCropTime( m_hopTime );
    
    //  bandwidth association region width 
    //  defaults to 2 kHz, corresponding to 
    //  1 kHz region center spacing:
    storeResidueBandwidth();

	//  configure the envelope builders using default 
	//  parameters:
	/*
	buildFundamentalEnv( *m_freqResolutionEnv * 0.99,
                         *m_freqResolutionEnv * 1.5 );
      
     */	//	!!!!!!!
	m_f0Builder.reset( 
        new FundamentalBuilder( *m_freqResolutionEnv * 0.99,
        						*m_freqResolutionEnv * 1.5,
        						-60., 8000. ) );
        
    m_ampEnvBuilder.reset( new AmpEnvBuilder );
    
    //  enable phase-correct Partial construction:
    m_phaseCorrect = true;
}

// -- analysis --
// ---------------------------------------------------------------------------
//  analyze
// ---------------------------------------------------------------------------
//! Analyze a vector of (mono) samples at the given sample rate         
//! (in Hz) and store the extracted Partials in the Analyzer's
//! PartialList (std::list of Partials). 
//! 
//! \param vec is a vector of floating point samples
//! \param srate is the sample rate of the samples in the vector 
//
void 
Analyzer::analyze( const std::vector<double> & vec, double srate )      
{ 
    BreakpointEnvelope reference( 1.0 );
    analyze( &(vec[0]),  &(vec[0]) + vec.size(), srate, reference ); 
}

// ---------------------------------------------------------------------------
//  analyze
// ---------------------------------------------------------------------------
//! Analyze a range of (mono) samples at the given sample rate      
//! (in Hz) and store the extracted Partials in the Analyzer's
//! PartialList (std::list of Partials). 
//! 
//! \param bufBegin is a pointer to a buffer of floating point samples
//! \param bufEnd is (one-past) the end of a buffer of floating point 
//! samples
//! \param srate is the sample rate of the samples in the buffer
//
void 
Analyzer::analyze( const double * bufBegin, const double * bufEnd, double srate )
{ 
    BreakpointEnvelope reference( 1.0 );
    analyze( bufBegin,  bufEnd, srate, reference ); 
}

// ---------------------------------------------------------------------------
//  analyze
// ---------------------------------------------------------------------------
//! Analyze a vector of (mono) samples at the given sample rate         
//! (in Hz) and store the extracted Partials in the Analyzer's
//! PartialList (std::list of Partials). Use the specified envelope
//! as a frequency reference for Partial tracking.
//!
//! \param vec is a vector of floating point samples
//! \param srate is the sample rate of the samples in the vector
//! \param reference is an Envelope having the approximate
//! frequency contour expected of the resulting Partials.
//
void 
Analyzer::analyze( const std::vector<double> & vec, double srate, 
                   const Envelope & reference )     
{ 
    analyze( &(vec[0]),  &(vec[0]) + vec.size(), srate, reference ); 
}


// ---------------------------------------------------------------------------
//  analyze
// ---------------------------------------------------------------------------
//! Analyze a range of (mono) samples at the given sample rate      
//! (in Hz) and store the extracted Partials in the Analyzer's
//! PartialList (std::list of Partials). Use the specified envelope
//! as a frequency reference for Partial tracking.
//! 
//! \param bufBegin is a pointer to a buffer of floating point samples
//! \param bufEnd is (one-past) the end of a buffer of floating point 
//! samples
//! \param srate is the sample rate of the samples in the buffer
//! \param reference is an Envelope having the approximate
//! frequency contour expected of the resulting Partials.
//
void 
Analyzer::analyze( const double * bufBegin, const double * bufEnd, double srate,
                   const Envelope & reference )
{ 
    //  configure the reassigned spectral analyzer, 
    //  always use odd-length windows:

    //  Kaiser window
    double winshape = KaiserWindow::computeShape( sidelobeLevel() );
    long winlen = KaiserWindow::computeLength( windowWidth() / srate, winshape );    
    if (! (winlen % 2)) 
    {
        ++winlen;
    }
    debugger << "Using Kaiser window of length " << winlen << endl;
    
    std::vector< double > window( winlen );
    KaiserWindow::buildWindow( window, winshape );
    
    std::vector< double > windowDeriv( winlen );
    KaiserWindow::buildTimeDerivativeWindow( windowDeriv, winshape );
       
    ReassignedSpectrum spectrum( window, windowDeriv );   
    
    //  configure the peak selection and partial formation policies:
    SpectralPeakSelector selector( srate, m_cropTime );
    PartialBuilder builder( m_freqDrift, reference );
    
    //  configure bw association policy, unless
    //  bandwidth association is disabled:
    std::auto_ptr< AssociateBandwidth > bwAssociator;
    if( m_bwAssocParam > 0 )
    {
        debugger << "Using bandwidth association regions of width " 
                 << bwRegionWidth() << " Hz" << endl;
        bwAssociator.reset( new AssociateBandwidth( bwRegionWidth(), srate ) );
    }
    else
    {
        debugger << "Bandwidth association disabled" << endl;
    }

    //  reset envelope builders:
    m_ampEnvBuilder->reset();
    m_f0Builder->reset();
    
    m_partials.clear();
        
    try 
    { 
        const double * winMiddle = bufBegin; 

        //  loop over short-time analysis frames:
        while ( winMiddle < bufEnd )
        {
            //  compute the time of this analysis frame:
            const double currentFrameTime = long(winMiddle - bufBegin) / srate;
            
            //  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)
            const double * sampsBegin = std::max( winMiddle - (winlen / 2), bufBegin );
            const double * sampsEnd = std::min( winMiddle + (winlen / 2) + 1, bufEnd );
            spectrum.transform( sampsBegin, winMiddle, sampsEnd );
            
             
            //  extract peaks from the spectrum, and thin
            Peaks peaks = selector.selectPeaks( spectrum, m_freqFloor ); 
			Peaks::iterator rejected = thinPeaks( peaks, currentFrameTime );

            //	fix the stored bandwidth values
            //	KLUDGE: need to do this before the bandwidth
            //	associator tries to do its job, because the mixed
            //	derivative is temporarily stored in the Breakpoint 
            //	bandwidth!!! FIX!!!!
            fixBandwidth( peaks );
            
            if ( m_bwAssocParam > 0 )
            {
                bwAssociator->associateBandwidth( peaks.begin(), rejected, peaks.end() );
            }
            
            //  remove rejected Breakpoints (needed above to 
            //  compute bandwidth envelopes):
            peaks.erase( rejected, peaks.end() );
            
            //  estimate the amplitude in this frame:
            m_ampEnvBuilder->build( peaks, currentFrameTime );
                        
            //  collect amplitudes and frequencies and try to 
            //  estimate the fundamental
            m_f0Builder->build( peaks, currentFrameTime );          

            //  form Partials from the extracted Breakpoints:
            builder.buildPartials( peaks, currentFrameTime );
            
            //  slide the analysis window:
            winMiddle += long( m_hopTime * srate ); //  hop in samples, truncated

        }   //  end of loop over short-time frames
        
        //  unwarp the Partial frequency envelopes:
        builder.finishBuilding( m_partials );
        
        //  fix the frequencies and phases to be consistent.
        if ( m_phaseCorrect )
        {
            fixFrequency( m_partials.begin(), m_partials.end() );
        }
        
        
        //  for debugging:
        /*
        if ( ! m_ampEnv.empty() )
        {
            LinearEnvelope::iterator peakpos = 
                std::max_element( m_ampEnv.begin(), m_ampEnv.end(), 
                                  compare2nd<LinearEnvelope::iterator::value_type> );
            notifier << "Analyzer found amp peak at time : " << peakpos->first
                     << " value: " << peakpos->second << endl;
        }
        */
    }
    catch ( Exception & ex ) 
    {
        ex.append( "analysis failed." );
        throw;
    }
}

// -- parameter access --

// ---------------------------------------------------------------------------
//  ampFloor
// ---------------------------------------------------------------------------
//! Return the amplitude floor (lowest detected spectral amplitude),            
//! in (negative) dB, for this Analyzer.                
//
double 
Analyzer::ampFloor( void ) const 
{ 
    return m_ampFloor; 
}

// ---------------------------------------------------------------------------
//  cropTime
// ---------------------------------------------------------------------------
//! Return the crop time (maximum temporal displacement of a time-
//! frequency data point from the time-domain center of the analysis
//! window, beyond which data points are considered "unreliable")
//! for this Analyzer.
//
double 
Analyzer::cropTime( void ) const 
{ 
    // debugger << "Analyzer::cropTime() is a deprecated member, and will be removed in a future Loris release." << endl;
    return m_cropTime; 
}

// ---------------------------------------------------------------------------
//  freqDrift
// ---------------------------------------------------------------------------
//! Return the maximum allowable frequency difference 
//! consecutive Breakpoints in a Partial envelope for this Analyzer.                
//
double 
Analyzer::freqDrift( void ) const 
{ 
    return m_freqDrift;
}

// ---------------------------------------------------------------------------
//  freqFloor
// ---------------------------------------------------------------------------
//! Return the frequency floor (minimum instantaneous Partial               
//! frequency), in Hz, for this Analyzer.               
//
double 
Analyzer::freqFloor( void ) const 
{ 
    return m_freqFloor; 
}

// ---------------------------------------------------------------------------
//  freqResolution
// ---------------------------------------------------------------------------
//! Return the frequency resolution (minimum instantaneous frequency        
//! difference between Partials) for this Analyzer at the specified
//! time in seconds. If no time is specified, then the initial resolution
//!	(at 0 seconds) is returned.
//! 
//! \param time is the time in seconds at which to evaluate the 
//!		   frequency resolution
//
double 
Analyzer::freqResolution( double time /* = 0.0 */ ) const 
{ 
    return m_freqResolutionEnv->valueAt( time ); 
}

// ---------------------------------------------------------------------------
//  hopTime
// ---------------------------------------------------------------------------
//! Return the hop time (which corresponds approximately to the 
//! average density of Partial envelope Breakpoint data) for this 
//! Analyzer.
//
double 
Analyzer::hopTime( void ) const 
{ 
    return m_hopTime; 
}

// ---------------------------------------------------------------------------
//  sidelobeLevel
// ---------------------------------------------------------------------------
//! Return the sidelobe attenutation level for the Kaiser analysis window in
//! positive dB. Larger numbers (e.g. 90) give very good sidelobe 
//! rejection but cause the window to be longer in time. Smaller numbers 
//! (like 60) raise the level of the sidelobes, increasing the likelihood
//! of frequency-domain interference, but allow the window to be shorter
//! in time.
//
double 
Analyzer::sidelobeLevel( void ) const 
{ 
    return m_sidelobeLevel; 
}

// ---------------------------------------------------------------------------
//  windowWidth
// ---------------------------------------------------------------------------
//! Return the frequency-domain main lobe width (measured between 
//! zero-crossings) of the analysis window used by this Analyzer.               
//
double 
Analyzer::windowWidth( void ) const 
{ 
    return m_windowWidth; 
}

// ---------------------------------------------------------------------------
//  phaseCorrect
// ---------------------------------------------------------------------------
//! Return true if the phases and frequencies of the constructed
//! partials should be modified to be consistent at the end of the
//! analysis, and false otherwise. (Default is true.)
//!
//! \param  TF is a flag indicating whether or not to construct
//!         phase-corrected Partials
bool 
Analyzer::phaseCorrect( void ) const
{
    return m_phaseCorrect;
}

// -- parameter mutation --

#define VERIFY_ARG(func, test)                                          \
    do {                                                                \
        if (!(test))                                                    \
            Throw( Loris::InvalidArgument, #func ": " #test  );         \
    } while (false)


// ---------------------------------------------------------------------------
//  setAmpFloor
// ---------------------------------------------------------------------------
//! Set the amplitude floor (lowest detected spectral amplitude), in            
//! (negative) dB, for this Analyzer. 
//! 
//! \param x is the new value of this parameter.                
//
void 
Analyzer::setAmpFloor( double x ) 
{ 
    VERIFY_ARG( setAmpFloor, x < 0 );
    m_ampFloor = x; 
}


// ---------------------------------------------------------------------------
//  setCropTime
// ---------------------------------------------------------------------------
//! Set the crop time (maximum temporal displacement of a time-
//! frequency data point from the time-domain center of the analysis
//! window, beyond which data points are considered "unreliable")
//! for this Analyzer.
//! 
//! \param x is the new value of this parameter.
//
void 
Analyzer::setCropTime( double x ) 
{ 
    VERIFY_ARG( setCropTime, x > 0 );
   // debugger << "Analyzer::setCropTime() is a deprecated member, and will be removed in a future Loris release." << endl;
    m_cropTime = x; 
}

// ---------------------------------------------------------------------------
//  setFreqDrift
// ---------------------------------------------------------------------------
//! Set the maximum allowable frequency difference between                  
//! consecutive Breakpoints in a Partial envelope for this Analyzer.                
//! 
//! \param x is the new value of this parameter.            
//
void 
Analyzer::setFreqDrift( double x ) 
{ 
    VERIFY_ARG( setFreqDrift, x > 0 );
    m_freqDrift = x; 
}

// ---------------------------------------------------------------------------
//  setFreqFloor
// ---------------------------------------------------------------------------
//! Set the frequency floor (minimum instantaneous Partial                  
//! frequency), in Hz, for this Analyzer.
//! 
//! \param x is the new value of this parameter.                    
//
void 
Analyzer::setFreqFloor( double x ) 
{ 
    VERIFY_ARG( setFreqFloor, x >= 0 );
    m_freqFloor = x; 
}

// ---------------------------------------------------------------------------
//  setFreqResolution (constant)
// ---------------------------------------------------------------------------
//! Set the frequency resolution (minimum instantaneous frequency       
//! difference between Partials) for this Analyzer. (Does not cause     
//! other parameters to be recomputed.)                                     
//! 
//! \param x is the new value of this parameter.                                        
//
void 
Analyzer::setFreqResolution( double x ) 
{ 
    VERIFY_ARG( setFreqResolution, x > 0 );
    m_freqResolutionEnv.reset( new LinearEnvelope( x ) ); 
}

// ---------------------------------------------------------------------------
//  setFreqResolution (envelope)
// ---------------------------------------------------------------------------
//! Set the time-varying frequency resolution (minimum instantaneous frequency       
//! difference between Partials) for this Analyzer. (Does not cause     
//! other parameters to be recomputed.)                                     
//! 
//! \param e is the envelope to copy for this parameter.                                        
//
void 
Analyzer::setFreqResolution( const Envelope & e ) 
{ 
	//	No mechanism exists to verify that the envelope never
	//	drops below zero, this can only be checked at analysis-time.
    // VERIFY_ARG( setFreqResolution, x > 0 );
    m_freqResolutionEnv.reset( e.clone() ); 
}

// ---------------------------------------------------------------------------
//  setSidelobeLevel
// ---------------------------------------------------------------------------
//! Set the sidelobe attenutation level for the Kaiser analysis window in
//! positive dB. Higher numbers (e.g. 90) give very good sidelobe 
//! rejection but cause the window to be longer in time. Lower 
//! numbers raise the level of the sidelobes, increasing the likelihood
//! of frequency-domain interference, but allow the window to be shorter
//! in time.
//! 
//! \param x is the new value of this parameter.    
//
void 
Analyzer::setSidelobeLevel( double x ) 
{ 
    VERIFY_ARG( setSidelobeLevel, x > 0 );
    m_sidelobeLevel = x; 
}

// ---------------------------------------------------------------------------
//  setHopTime
// ---------------------------------------------------------------------------
//! Set the hop time (which corresponds approximately to the average
//! density of Partial envelope Breakpoint data) for this Analyzer.
//! 
//! \param x is the new value of this parameter.
//
void 
Analyzer::setHopTime( double x ) 
{ 
    VERIFY_ARG( setHopTime, x > 0 );
    m_hopTime = x; 
}

// ---------------------------------------------------------------------------
//  setWindowWidth
// ---------------------------------------------------------------------------
//! Set the frequency-domain main lobe width (measured between 
//! zero-crossings) of the analysis window used by this Analyzer.   
//! 
//! \param x is the new value of this parameter.            
//
void 
Analyzer::setWindowWidth( double x ) 
{ 
    VERIFY_ARG( setWindowWidth, x > 0 );
    m_windowWidth = x; 
}

// ---------------------------------------------------------------------------
//  setPhaseCorrect
// ---------------------------------------------------------------------------
//! Indicate whether the phases and frequencies of the constructed
//! partials should be modified to be consistent at the end of the
//! analysis. (Default is true.)
//!
//! \param  TF is a flag indicating whether or not to construct
//!         phase-corrected Partials
void 
Analyzer::setPhaseCorrect( bool TF )
{
    m_phaseCorrect = TF;
}

//  -- bandwidth envelope specification --


// ---------------------------------------------------------------------------
//  storeResidueBandwidth
// ---------------------------------------------------------------------------
//!	Construct Partial bandwidth envelopes during analysis
//!	by associating residual energy in the spectrum (after
//! peak extraction) with the selected spectral peaks that
//!	are used to construct Partials. 
//!	
//!	\param regionWidth is the width (in Hz) of the bandwidth 
//!	association regions used by this process, must be positive.
//!	If unspecified, a default value is used.
//
void 
Analyzer::storeResidueBandwidth( double regionWidth ) 
{ 
    VERIFY_ARG( storeResidueBandwidth, regionWidth > 0 );
    m_bwAssocParam = regionWidth; 
}   

// ---------------------------------------------------------------------------
//  storeConvergenceBandwidth
// ---------------------------------------------------------------------------
//!	Construct Partial bandwidth envelopes during analysis
//!	by storing the mixed derivative of short-time phase, 
//!	scaled and shifted so that a value of 0 corresponds
//!	to a pure sinusoid, and a value of 1 corresponds to a
//! bandwidth-enhanced sinusoid with maximal energy spread
//! (minimum sinusoidal convergence).
//!
//!	\param tolerance is the amount of range over which the 
//!	mixed derivative indicator should be allowed to drift away 
//!	from a pure sinusoid before saturating. This range is mapped
//!	to bandwidth values on the range [0,1]. Must be positive and 
//!	not greater than 1. If unspecified, a default value is used.
//
void 
Analyzer::storeConvergenceBandwidth( double tolerance ) 
{ 
	if ( 1.0 < tolerance )
	{
		//	notify and scale, in Loris 1.5, tolerance was
		//	specified as a percent
		notifier << "Analyzer::storeConvergenceBandwidth, conergence tolerance "
					"should be positive and less than 1.0, scaling by 1/100" << endl;
		tolerance *= 0.01;
	}
	
    VERIFY_ARG( storeConvergenceBandwidth, 
    			(tolerance > 0) && (tolerance <= 1) );
    			
	//	store a negative value so that it can be 
	//	identified when used:
    m_bwAssocParam = -tolerance; 
}   

// ---------------------------------------------------------------------------
//  storeNoBandwidth
// ---------------------------------------------------------------------------
//!	Disable bandwidth envelope construction. Bandwidth 
//!	will be zero for all Breakpoints in all Partials.
//
void 
Analyzer::storeNoBandwidth( void ) 
{ 
    m_bwAssocParam = 0; 
}   

// ---------------------------------------------------------------------------
//!	Return true if this Analyzer is configured to compute
//! bandwidth envelopes using the spectral residue after
//! peaks have been identified, and false otherwise.
// ---------------------------------------------------------------------------
bool 
Analyzer::bandwidthIsResidue( void ) const
{ 
    return m_bwAssocParam > 0.; 
}

// ---------------------------------------------------------------------------
//!	Return true if this Analyzer is configured to compute
//!	bandwidth envelopes using the mixed derivative convergence
//!	indicator, and false otherwise.
// ---------------------------------------------------------------------------
bool 
Analyzer::bandwidthIsConvergence( void ) const
{ 
    return m_bwAssocParam < 0.; 
}


// ---------------------------------------------------------------------------
//! Return the width (in Hz) of the Bandwidth Association regions
//! used by this Analyzer, only if the spectral residue method is
//!	used to compute bandwidth envelopes. Return zero if the mixed
//! derivative method is used, or if no bandwidth is computed.
// ---------------------------------------------------------------------------
double 
Analyzer::bwRegionWidth( void ) const
{
	if ( m_bwAssocParam > 0 )
	{
		return m_bwAssocParam;
	}
	return 0;
}

// ---------------------------------------------------------------------------
//!	Return the mixed derivative convergence tolerance (percent)
//!	only if the convergence indicator is used to compute
//!	bandwidth envelopes. Return zero if the spectral residue
//!	method is used or if no bandwidth is computed.
// ---------------------------------------------------------------------------
double 
Analyzer::bwConvergenceTolerance( void ) const
{
	if ( m_bwAssocParam < 0 )
	{
		return - m_bwAssocParam;
	}
	return 0;
}


// -- PartialList access --

// ---------------------------------------------------------------------------
//  partials
// ---------------------------------------------------------------------------
//! Return a mutable reference to this Analyzer's list of 
//! analyzed Partials. 
//
PartialList & 
Analyzer::partials( void ) 
{ 
    return m_partials; 
}

// ---------------------------------------------------------------------------
//  partials
// ---------------------------------------------------------------------------
//! Return an immutable (const) reference to this Analyzer's 
//! list of analyzed Partials. 
//
const PartialList & 
Analyzer::partials( void ) const
{ 
    return m_partials; 
}

// ---------------------------------------------------------------------------
//  buildFundamentalEnv
// ---------------------------------------------------------------------------
//! Specify parameters for constructing a fundamental frequency 
//! envelope for the analyzed sound during analysis. The fundamental 
//! frequency estimate can be accessed by fundamentalEnv() after the 
//! analysis is complete. 
//!
//! \param  fmin is the lower bound on the fundamental frequency estimate
//! \param  fmax is the upper bound on the fundamental frequency estimate
//! \param  threshDb is the lower bound on the amplitude of a spectral peak
//!         that will constribute to the fundamental frequency estimate (very
//!         low amplitude peaks tend to have less reliable frequency estimates).
//!         Default is -60 dB.
//! \param  threshHz is the upper bound on the frequency of a spectral
//!         peak that will constribute to the fundamental frequency estimate.
//!         Default is 8 kHz.
//
void Analyzer::buildFundamentalEnv( double fmin, double fmax, 
                                    double threshDb, double threshHz )
{
    m_f0Builder.reset( 
        new FundamentalBuilder( fmin, fmax, threshDb, threshHz ) );
}

// ---------------------------------------------------------------------------
//  fundamentalEnv
// ---------------------------------------------------------------------------
//! Return the fundamental frequency estimate envelope constructed
//! during the most recent analysis performed by this Analyzer.
//! Will be empty unless buildFundamentalEnv was invoked to enable the
//! construction of this envelope during analysis.
//
const LinearEnvelope &
Analyzer::fundamentalEnv( void ) const
{   
    return m_f0Builder->envelope(); 
}


// ---------------------------------------------------------------------------
//  ampEnv
// ---------------------------------------------------------------------------
//! Return the overall amplitude estimate envelope constructed
//! during the most recent analysis performed by this Analyzer.
//! Will be empty unless buildAmpEnv was invoked to enable the
//! construction of this envelope during analysis.
//
const LinearEnvelope & 
Analyzer::ampEnv( void ) const
{ 
    return m_ampEnvBuilder->envelope(); 
}

// -- private helpers --

// ---------------------------------------------------------------------------
//	can_mask
// ---------------------------------------------------------------------------
//	functor used for identying peaks that are too close
//	in frequency to another louder peak:
struct can_mask
{
	//	masking occurs if any (louder) peak falls
	//	in the frequency range delimited by fmin and fmax:
	bool operator()( const SpectralPeak & v )  const
	{ 
		return	( v.frequency() > _fmin ) && 
				( v.frequency() < _fmax ); 
	}
		
	//	constructor:
	can_mask( double x, double y ) : 
		_fmin( x ), _fmax( y ) 
		{ if (x>y) std::swap(x,y); }
		
	//	bounds:
private:
	double _fmin, _fmax;
};

// ---------------------------------------------------------------------------
//	negative_time
// ---------------------------------------------------------------------------
//	functor used to identify peaks that have reassigned times 
//	before 0:
struct negative_time
{
	//	negative times occur when the reassigned time
	// 	plus the current frame time is less than 0:
	bool operator()( const Peaks::value_type & v )  const
	{ 
		return 0 > ( v.time() + _frameTime );
	}
		
	//	constructor:
	negative_time( double t ) : 
		_frameTime( t ) {}
		
	//	bounds:
private:
	double _frameTime;
	
};


// ---------------------------------------------------------------------------
//	thinPeaks (HELPER)
// ---------------------------------------------------------------------------
//	Reject peaks that are too quiet (low amplitude). Peaks that are retained,
//	but are quiet enough to be in the specified fadeRange should be faded.
//	Peaks having negative times are also rejected.
//
//	This is exactly the same as the basic peak selection strategy, there
//	is no tracking here.
//	
//	Rejected peaks are placed at the end of the peak collection.
//	Return the first position in the collection containing a rejected peak,
//	or the end of the collection if no peaks are rejected.
//
//	This used to be part of SpectralPeakSelector, but it really had no place
//	there. It _should_ remove the rejected peaks, but for now, those are needed
//	by the bandwidth association strategy.
//
Peaks::iterator 
Analyzer::thinPeaks( Peaks & peaks, double frameTime  )
{
	const double ampFloordB = m_ampFloor;

	//  fade quiet peaks out over 10 dB:
	const double fadeRangedB = 10.0;

	//	compute absolute magnitude thresholds:
	const double threshold = std::pow( 10., 0.05 * ampFloordB );
	const double beginFade = std::pow( 10., 0.05 * (ampFloordB+fadeRangedB) );

	//	louder peaks are preferred, so consider them 
	//	in order of louder magnitude:
	std::sort( peaks.begin(), peaks.end(), SpectralPeak::sort_greater_amplitude );
	
    //  negative times are not real, but still might represent
    //  a noisy part of the spectrum...
	Peaks::iterator bogusTimes = 
		std::remove_if( peaks.begin(), peaks.end(), negative_time( frameTime ) );
		
	//	...get rid of them anyway
    peaks.erase( bogusTimes, peaks.end() );
    bogusTimes = peaks.end();
        
    
	Peaks::iterator it = peaks.begin();
	Peaks::iterator beginRejected = it;

    const double freqResolution = 
    	std::max( m_freqResolutionEnv->valueAt( frameTime ), 0.0 ); 
    
    
	while ( it != peaks.end() ) 
	{
		SpectralPeak & pk = *it;
		
		//	keep this peak if it is loud enough and not
		//	 too near in frequency to a louder one:
		double lower = pk.frequency() - freqResolution;
		double upper = pk.frequency() + freqResolution;
		if ( pk.amplitude() > threshold &&
			 beginRejected == std::find_if( peaks.begin(), beginRejected, can_mask(lower, upper) ) )
		{
			//	this peak is a keeper, fade its
			//	amplitude if it is too quiet:
			if ( pk.amplitude() < beginFade )
			{
				double alpha = (beginFade - pk.amplitude())/(beginFade - threshold);
				pk.setAmplitude( pk.amplitude() * (1. - alpha) );
			}
			
			//	keep retained peaks at the front of the collection:
			if ( it != beginRejected )
			{
				std::swap( *it, *beginRejected );
			}
			++beginRejected;
		}
		++it;
	}
	
	// debugger << "thinPeaks retained " << std::distance( peaks.begin(), beginRejected ) << endl;

	//  remove rejected Breakpoints:
	//peaks.erase( beginRejected, peaks.end() );
	
	return beginRejected;
}

// ---------------------------------------------------------------------------
//	fixBandwidth (HELPER)
// ---------------------------------------------------------------------------
//	Fix the bandwidth value stored in the specified Peaks. 
//	This function is invoked if the spectral residue method is
//	not used to compute bandwidth (that method overwrites the
//	bandwidth already). If the convergence method is used to 
//	compute bandwidth, the appropriate scaling is applied
//	to the stored mixed phase derivative. Otherwise, the
//	Peak bandwidth is set to zero.
//
//  The convergence value is on the range [0,1], 0 for a sinusoid, 
//  and 1 for an impulse. If convergence tolerance is specified (as
//  a negative value in m_bwAssocParam), it should be positive and 
//  less than 1, and specifies the convergence value that is to 
//  correspond to bandwidth equal to 1.0. This is achieved by scaling
//  the convergence by the inverse of the tolerance, and saturating
//  at 1.0.
void Analyzer::fixBandwidth( Peaks & peaks )
{
	
	if ( m_bwAssocParam < 0 )
	{
		double scale = 1.0 / (- m_bwAssocParam);	
			// m_bwAssocParam stores negative tolerance
	
		for ( Peaks::iterator it = peaks.begin(); it != peaks.end(); ++it )
		{
            SpectralPeak & pk = *it;
			pk.setBandwidth( std::min( 1.0, scale * pk.bandwidth() ) );
		}
	}
	else if ( m_bwAssocParam == 0 )
	{
		for ( Peaks::iterator it = peaks.begin(); it != peaks.end(); ++it )
		{
            SpectralPeak & pk = *it;
			pk.setBandwidth( 0 );
		}
	}
}

}   //  end of namespace Loris