/*
 * 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
 *
 *
 * Oscillator.C
 *
 * Implementation of class Loris::Oscillator, a Bandwidth-Enhanced Oscillator.
 *
 * Kelly Fitz, 31 Aug 1999
 * loris@cerlsoundgroup.org
 *
 * http://www.cerlsoundgroup.org/Loris/
 *
 */
 
#if HAVE_CONFIG_H
    #include "config.h"
#endif

#include "Oscillator.h"

#include "Filter.h"
#include "Partial.h"
#include "Notifier.h"

#include <cmath>
#include <vector>

#if defined(HAVE_M_PI) && (HAVE_M_PI)
    const double Pi = M_PI;
#else
    const double Pi = 3.14159265358979324;
#endif
const double TwoPi = 2*Pi;

//  begin namespace
namespace Loris {


// ---------------------------------------------------------------------------
//  Oscillator construction
// ---------------------------------------------------------------------------
//  Initialize stochastic modulators and state variables.
//
Oscillator::Oscillator( void ) :
    m_modulator( 1.0 /* seed */ ),
    m_filter( prototype_filter() ),
    m_instfrequency( 0 ),
    m_instamplitude( 0 ),
    m_instbandwidth( 0 ),
    m_determphase( 0 )
{
}

// ---------------------------------------------------------------------------
//  resetEnvelopes
// ---------------------------------------------------------------------------
//  Reset the instantaneous envelope parameters 
//  (frequency, amplitude, bandwidth, and phase).
//  The sample rate is needed to convert the 
//  Breakpoint frequency (Hz) to radians per sample.
//
void 
Oscillator::resetEnvelopes( const Breakpoint & bp, double srate )
{
    //  Remember that the oscillator only knows about 
    //  radian frequency! Convert!
    m_instfrequency = bp.frequency() * TwoPi / srate;
    m_instamplitude = bp.amplitude();
    m_instbandwidth = bp.bandwidth();
    m_determphase = bp.phase();
    
    //  clamp bandwidth:
    if ( m_instbandwidth > 1. )
    {
        debugger << "clamping bandwidth at 1." << endl;
        m_instbandwidth = 1.;
    }
    else if ( m_instbandwidth < 0. )
    { 
        debugger << "clamping bandwidth at 0." << endl;
        m_instbandwidth = 0.;
    }

    //  don't alias:
    if ( m_instfrequency > Pi )
    { 
        debugger << "fading out aliasing Partial" << endl;
        m_instamplitude = 0.;
    }
    
    //  Reset the fitler state too.
    m_filter.clear();
    
}

// ---------------------------------------------------------------------------
//  m2pi
// ---------------------------------------------------------------------------
//  O'Donnell's phase wrapping function.
//
static inline double m2pi( double x )
{
    using namespace std; // floor should be in std
    #define ROUND(x) (floor(.5 + (x)))
    return x + ( TwoPi * ROUND(-x/TwoPi) );
}

// ---------------------------------------------------------------------------
//  setPhase
// ---------------------------------------------------------------------------
//  Reset the phase of the Oscillator to the specified
//  value, and clear the accumulated phase modulation. (?)
//  Or not.
//  This is done when the amplitude of a Partial goes to 
//  zero, so that onsets are preserved in distilled
//  and collated Partials.
//
void 
Oscillator::setPhase( double ph )
{
    m_determphase = m2pi(ph);
}

// ---------------------------------------------------------------------------
//  oscillate
// ---------------------------------------------------------------------------
//  Accumulate bandwidth-enhanced sinusoidal samples modulating the 
//  oscillator state from its current values of radian frequency,
//  amplitude, and bandwidth to the specified target values, into
//  the specified half-open range of doubles.
//
//  The caller must ensure that the range is valid. Target parameters
//  are bounds-checked. 
//
void
Oscillator::oscillate( double * begin, double * end,
                       const Breakpoint & bp, double srate )
{
    double targetFreq = bp.frequency() * TwoPi / srate;     //  radians per sample
    double targetAmp = bp.amplitude(); 
    double targetBw = bp.bandwidth();
    
    //  clamp bandwidth:
    if ( targetBw > 1. )
    {
        debugger << "clamping bandwidth at 1." << endl;
        targetBw = 1.;
    }
    else if ( targetBw < 0. )
    { 
        debugger << "clamping bandwidth at 0." << endl;
        targetBw = 0.;
    }
        
    //  don't alias:
    if ( targetFreq > Pi )  //  radian Nyquist rate
    {
        debugger << "fading out Partial above Nyquist rate" << endl;
        targetAmp = 0.;
    }

    //  compute trajectories:
    const double dTime = 1. / (end - begin);
    const double dFreqOver2 = 0.5 * (targetFreq - m_instfrequency) * dTime;
    	//	split frequency update in two steps, update phase using average
    	//	frequency, after adding only half the frequency step
    	
    const double dAmp = (targetAmp - m_instamplitude)  * dTime;
    const double dBw = (targetBw - m_instbandwidth)  * dTime;

    //  Use temporary local variables for speed.
    //  Probably not worth it when I am computing square roots 
    //  and cosines...
    double ph = m_determphase;
    double f = m_instfrequency;
    double a = m_instamplitude;
    double bw = m_instbandwidth;
    
    //	Also use a more efficient sample loop when the bandwidth is zero.
    if ( 0 < bw || 0 < dBw )
    {
		double am, nz;
		for ( double * putItHere = begin; putItHere != end; ++putItHere )
		{
			//  use math functions in namespace std:
			using namespace std;
	
			//  compute amplitude modulation due to bandwidth:
			//
			//  This will give the right amplitude modulation when scaled
			//  by the Partial amplitude:
			//
			//  carrier amp: sqrt( 1. - bandwidth ) * amp
			//  modulation index: sqrt( 2. * bandwidth ) * amp
			//
			nz = m_filter.apply( m_modulator.sample() );
			am = sqrt( 1. - bw ) + ( nz * sqrt( 2. * bw ) );  
					
			//  compute a sample and add it into the buffer:
			*putItHere += am * a * cos( ph );
				
			//  update the instantaneous oscillator state:
			f += dFreqOver2;
			ph += f;   //  frequency is radians per sample
			f += dFreqOver2;
			a += dAmp;
			bw += dBw;
			if (bw < 0.)
			{
				bw = 0.;
			}				
		}   // end of sample computation loop
	}
	else
	{
		for ( double * putItHere = begin; putItHere != end; ++putItHere )
		{
			//  use math functions in namespace std:
			using namespace std;
	
			//	no modulation when there is no bandwidth
			
			//  compute a sample and add it into the buffer:
			*putItHere += a * cos( ph );
				
			//  update the instantaneous oscillator state:
			f += dFreqOver2;
			ph += f;   //  frequency is radians per sample
			f += dFreqOver2;
			a += dAmp;
		}   // end of sample computation loop
	
	}
	
	
    //	copy out of the local variables?
    //	no need because we are assigning to the target
    //	values below:
    /*
    m_instfrequency = f;
    m_instamplitude = a;
    m_instbandwidth = bw;
    */
    
    //  wrap phase to prevent eventual loss of precision at
    //  high oscillation frequencies:
    //  (Doesn't really matter much exactly how we wrap it, 
    //  as long as it brings the phase nearer to zero.)
    m_determphase = m2pi( ph );
    
    //  set the state variables to their target values,
    //  just in case they didn't arrive exactly (overshooting
    //  amplitude or, especially, bandwidth, could be bad, and
    //  it does happen):
    m_instfrequency = targetFreq;
    m_instamplitude = targetAmp;
    m_instbandwidth = targetBw;
}

// ---------------------------------------------------------------------------
//  protoype filter (static member)
// ---------------------------------------------------------------------------
//  Static local function for obtaining a prototype Filter
//  to use in Oscillator construction. Eventually, allow
//  external (client) specification of the Filter prototype.
//
const Filter & 
Oscillator::prototype_filter( void )
{
    //  Chebychev order 3, cutoff 500, ripple -1.
    //
    //  Coefficients obtained from http://www.cs.york.ac.uk/~fisher/mkfilter/
    //  Digital filter designed by mkfilter/mkshape/gencode   A.J. Fisher
    //
    static const double Gain = 4.663939184e+04;
    static const double ExtraScaling = 6.;
    static const double MaCoefs[] = { 1., 3., 3., 1. }; 
    static const double ArCoefs[] = { 1., -2.9258684252, 2.8580608586, -0.9320209046 };

    static const Filter proto( MaCoefs, MaCoefs + 4, ArCoefs, ArCoefs + 4, ExtraScaling/Gain );
    return proto;
}



}   //  end of namespace Loris