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

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

#include "KaiserWindow.h"
#include "LorisExceptions.h"
#include <cmath>

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

using namespace std;

//  begin namespace
namespace Loris {

//  prototypes for static helpers, defined below
static double zeroethOrderBessel( double x );
static double firstOrderBessel( double x );

// ---------------------------------------------------------------------------
//  buildWindow
// ---------------------------------------------------------------------------
//! Build a new Kaiser analysis window having the specified shaping
//! parameter. See Oppenheim and Schafer:  "Digital Signal Processing" 
//! (1975), p. 452 for further explanation of the Kaiser window. Also, 
//! see Kaiser and Schafer, 1980.
//!
//! \param      win is the vector that will store the window
//!             samples. The number of samples computed will be
//!             equal to the length of this vector. Any previous
//!             contents will be overwritten.
//! \param      shape is the Kaiser shaping parameter, controlling
//!             the sidelobe rejection level.
//
void
KaiserWindow::buildWindow( vector< double > & win, double shape )
{   
    //  Pre-compute the shared denominator in the Kaiser equation. 
    const double oneOverDenom = 1.0 / zeroethOrderBessel( shape );
  
    const unsigned int N = win.size() - 1;
    const double oneOverN = 1.0 / N;
    
    for ( unsigned int n = 0; n <= N; ++n )
    {
        const double K = (2.0 * n * oneOverN) - 1.0;
        const double arg = sqrt( 1.0 - (K * K) );
        
        win[n] = zeroethOrderBessel( shape * arg ) * oneOverDenom;
    }
}

// ---------------------------------------------------------------------------
//  createDerivativeWindow
// ---------------------------------------------------------------------------
//! Build a new time-derivative Kaiser analysis window having the
//! specified shaping parameter, for computing frequency reassignment.
//! The closed form for the time derivative can be obtained from the
//! property of modified Bessel functions that the derivative of the
//! zeroeth order function is equal to the first order function.
//!
//! \param      win is the vector that will store the window
//!             samples. The number of samples computed will be
//!             equal to the length of this vector. Any previous
//!             contents will be overwritten.
//! \param      shape is the Kaiser shaping parameter, controlling
//!             the sidelobe rejection level.
//
void
KaiserWindow::buildTimeDerivativeWindow( vector< double > & win, double shape )
{   
    //  Pre-compute the common factor that does not depend on n.
    const unsigned int N = win.size() - 1;
    const double oneOverN = 1.0 / N;
    
    const double commonFac = - 2.0 * shape / (N * zeroethOrderBessel( shape ) );
  
    //  w'[0] = w'[N] = 0
    win[0] = win[N] = 0.0;
    
    for ( unsigned int n = 1; n < N; ++n )
    {
        const double K = (2.0 * n * oneOverN) - 1.0;
        const double arg = sqrt( 1.0 - (K * K) );
        
        win[n] = commonFac * firstOrderBessel( shape * arg ) * K / arg;
    }
}


// ---------------------------------------------------------------------------
//  zeroethOrderBessel
// ---------------------------------------------------------------------------
//  Compute the zeroeth order modified Bessel function of the first kind 
//  at x using the series expansion, used to compute the Kasier window
//  function.
//
static double zeroethOrderBessel( double x )
{
    const double eps = 0.000001;
    
    //  initialize the series term for m=0 and the result
    double besselValue = 0;
    double term = 1;
    double m = 0;
    
    //  accumulate terms as long as they are significant
    while(term  > eps * besselValue)
    {
        besselValue += term;
        
        //  update the term
        ++m;
        term *= (x*x) / (4*m*m);
    }
    
    return besselValue;
}

// ---------------------------------------------------------------------------
//  firstOrderBessel
// ---------------------------------------------------------------------------
//  Compute the first order modified Bessel function of the first kind 
//  at x using the series expansion, used to compute the time derivative
//  of the Kasier window function for computing frequency reassignment.
//
static double firstOrderBessel( double x )
{
    const double eps = 0.000001;
    
    //  initialize the series term for m=0 and the result
    double besselValue = 0;
    double term = .5*x;
    double m = 0;
    
    //  accumulate terms as long as they are significant
    while(term  > eps * besselValue)
    {
        besselValue += term;
        
        //  update the term
        ++m;
        term *= (x*x) / (4*m*(m+1));
    }
    
    return besselValue;
}

// ---------------------------------------------------------------------------
//  computeShape
// ---------------------------------------------------------------------------
//  Compute the Kaiser window shaping parameter from the specified attenuation 
//  of side lobes. This algorithm is given in Kaiser an Schafer,1980 and is 
//  supposed to give better than 0.36% accuracy (Kaiser and Schafer 1980).
//
double
KaiserWindow::computeShape( double atten )
{
    if ( atten < 0. )
    {
        Throw( InvalidArgument, 
        	   "Kaiser window shape must be computed from positive (> 0dB)"
        	   " sidelobe attenuation. (received attenuation < 0)" );
	}
	
    double alpha;
    
    if ( atten > 60.0 )
    {
        alpha = 0.12438 * (atten + 6.3);
    }
    else if ( atten > 13.26 )
    {
        alpha = 0.76609L * ( pow((atten - 13.26), 0.4) ) + 
                            0.09834L * (atten - 13.26L);
    }
    else
    {
        //  can't have less than 13dB.
        alpha = 0.0;
    }
    
    return alpha;
}
// ---------------------------------------------------------------------------
//  computeLength
// ---------------------------------------------------------------------------
//  Compute the length (in samples) of the Kaiser window from the desired 
//  (approximate) main lobe width and the control parameter. Of course, since 
//  the window must be an integer number of samples in length, your actual 
//  lobal mileage may vary. This equation appears in Kaiser and Schafer 1980
//  (on the use of the I0 window class for spectral analysis) as Equation 9.
//
//  The main width of the main lobe must be normalized by the sample rate,
//  that is, it is a fraction of the sample rate.
//
unsigned long
KaiserWindow::computeLength( double width, double alpha )
{
    //double alpha = computeShape( atten );

    //  The last 0.5 is cheap rounding.
    //  But I think I don't need cheap rounding because the equation 
    //  from Kaiser and Schafer has a +1 that appears to be a cheap
    //  ceiling function.
    return long(1.0 + (2. * sqrt((Pi*Pi) + (alpha*alpha)) / (Pi * width)) /* + 0.5 */);
}

}   //  end of namespace Loris