// file: $isip_ifc/class/algo/Window/win_06.cc
// version: $Id: win_06.cc 10561 2006-04-21 04:12:06Z srinivas $
//

// isip include files
//
#include "Window.h"
#include "Bessel.h"
#include "Chebyshev.h"

// method: setConstants
//
// arguments: 
//  const VectorFloat& values: (input) vector of constants
//
// return: a bool8 value indicating status
//
// this method sets the constants of the window
//
bool8 Window::setConstants(const VectorFloat& values_a) {
  
  // make sure the correct number of coefficients are specified
  // for the given algorithm
  //
  if (algorithm_d == RECTANGULAR) {
    if (values_a.length() != DEF_RECT_CONSTANTS.length()) {
      return Error::handle(name(), L"setConstants", ERR, __FILE__, __LINE__);
    }
  }
  else if (algorithm_d == BARTLETT) {
    if (values_a.length() != DEF_BART_CONSTANTS.length()) {
      return Error::handle(name(), L"setConstants", ERR, __FILE__, __LINE__);
    }
  }
  else if (algorithm_d == BLACKMAN) {
    if (values_a.length() != DEF_BLCK_CONSTANTS.length()) {
      return Error::handle(name(), L"setConstants", ERR, __FILE__, __LINE__);
    }
  }
  else if (algorithm_d == DOLPH_CHEBYSHEV) {
    if (values_a.length() != DEF_DCHB_CONSTANTS.length()) {
      return Error::handle(name(), L"setConstants", ERR, __FILE__, __LINE__);
    }
  }
  else if (algorithm_d == GAUSSIAN) {
    if (values_a.length() != DEF_GAUS_CONSTANTS.length()) {
      return Error::handle(name(), L"setConstants", ERR, __FILE__, __LINE__);
    }
  }
  else if (algorithm_d == HAMMING) {
    if (values_a.length() != DEF_HAMM_CONSTANTS.length()) {
      return Error::handle(name(), L"setConstants", ERR, __FILE__, __LINE__);
    }
  }
  else if (algorithm_d == HANNING) {
    if (values_a.length() != DEF_HANN_CONSTANTS.length()) {
      return Error::handle(name(), L"setConstants", ERR, __FILE__, __LINE__);
    }
  }
  else if (algorithm_d == KAISER) {
    if (values_a.length() != DEF_KAIS_CONSTANTS.length()) {
      return Error::handle(name(), L"setConstants", ERR, __FILE__, __LINE__);
    }
  }
  else if (algorithm_d == LIFTER) {
    if (values_a.length() != DEF_LIFT_CONSTANTS.length()) {
      return Error::handle(name(), L"setConstants", ERR, __FILE__, __LINE__);
    }
  }
  else if (algorithm_d == CUSTOM) {
    return Error::handle(name(), L"setConstants", ERR, __FILE__, __LINE__);
  }
  else {
    return Error::handle(name(), L"setConstants", Error::ENUM,
			 __FILE__, __LINE__);
  }
  
  // set the constants
  //
  constants_d.assign(values_a);
  
  // make the window invalid
  //
  is_valid_d = false;
  
  // exit gracefully
  //
  return true;
}

// method: generateRectangular
//
// arguments: 
//  VectorFloat& output: (output) the window function
//  const VectorFloat& params: (input) window function parameters
//  int32 num_points: (input) number of points in the window
//
// return: a bool8 value indicating status
//
// this method generates a rectangular window function:
//
//  w(n) = { G   n <= |M|
//         { 0   elsewhere
//
// there is one parameter allowed for this window:
//
//  params(0): G, the gain of the window
//
bool8 Window::generateRectangular(VectorFloat& output_a,
				    const VectorFloat& params_a,
				    int32 num_points_a) const {
  
  // for even or odd number of points, do the same thing
  //
  return output_a.assign(num_points_a, params_a(0));
}

// method: generateBartlett
//
// arguments: 
//  VectorFloat& output: (output) the window function
//  const VectorFloat& params: (input) window function parameters
//  int32 num_points: (input) number of points in the window
//
// return: a bool8 value indicating status
//
// this method generates a Bartlett (triangular) window function:
//
//  w(n) = { G (1 - |n| / (M + 1))  n <= |M|
//         { 0                      elsewhere
//
// where M = numn_points_a / 2. a good reference for this window
// can be found at:
//
//  R.A. Roberts and C.T. Mullis, Digital Signal Processing,
//  Addison Wesley Publishing Company, Reading, Massachusetts, USA,
//  1987, pp. 183-188.
//
// there is one parameter allowed for this window:
//
//  params(0): G, the gain of the window
//
// for an odd number of points, the window is defined as above.
// for an even number of points, an interpolation technique is used.
//
bool8 Window::generateBartlett(VectorFloat& output_a,
				 const VectorFloat& params_a,
				 int32 num_points_a) const {
  
  // compute some constants related to the length of the window
  //
  int32 N = num_points_a;
  int32 Nm1 = N - 1;
  int32 Nm2 = Nm1 - 1;
  float32 scale = (float32)Nm2 / (float32)Nm1;
  
  // compute the half-window length
  //
  int32 M = (N - 1) / 2;
  int32 Mp1 = M + 1;
  
  // create output space
  //
  output_a.setLength(N);
  
  // case 1: an odd number of points
  //
  if ((N % 2) != 0) {
    for (int32 n = -M; n <= M; n++) {
      output_a(n + M) = (float32)params_a(0) * (1 - Integral::abs(n) / (float32)(Mp1));
    }
  }
  
  // case 2: an even number of points
  //  resample an (N-1) point window
  //
  else {
    
    // generate a window for (N-1) points
    // sample it at a new time scale 
    //
    for (int32 n = 0; n < N; n++) {
      float32 scaled_index = -M + scale * n;
      output_a(n) = (float32)params_a(0) *
	(1 - Integral::abs(scaled_index) / (float32)(Mp1));
    }
  }    
  
  // exit gracefully
  //
  return true;
}

// method: generateBlackman
//
// arguments: 
//  VectorFloat& output: (output) the window function
//  const VectorFloat& params: (input) window function parameters
//  int32 num_points: (input) number of points in the window
//
// return: a bool8 value indicating status
//
// this method generates a Blackman window function:
//
//  w(n) = { G (a1 + a2 * cos(2*PI*n / (2*M+1))
//         {       + (0.5 - a1) * cos(4*PI*n / (2*M+1))    n <= |M|
//         { 0                                             elsewhere
//
// a good reference for this window can be found at:
//
//  S.K. Mitra, Digital Signal Processing,
//  McGraw-Hill, Boston, Massachusetts, USA, 2001, pp. 452.
//
// there are three parameters allowed for this window:
//
//  params(0): G, the gain of the window
//  params(1): a1, the DC offset
//  params(2): a2, the first harmonic weight
//
// for an odd number of points, the window is defined as above.
// for an even number of points, an interpolation technique is used.
//
bool8 Window::generateBlackman(VectorFloat& output_a,
				 const VectorFloat& params_a,
				 int32 num_points_a) const {
  
  // compute some constants related to the length of the window
  //
  int32 N = num_points_a;
  int32 Nm1 = N - 1;
  int32 Nm2 = Nm1 - 1;
  float32 scale = (float32)Nm2 / (float32)Nm1;
  
  // compute the half-window length
  //
  int32 M = (N - 1) / 2;
  int32 Mp1 = M + 1;
  
  // create output space
  //
  output_a.setLength(N);
  
  // case 1: an odd number of points
  //
  if ((N % 2) != 0) {
    for (int32 n = -M; n <= M; n++) {
      float32 arg1 = 2.0 * Integral::PI * (float32)(n) / (float32)Mp1;
      float32 arg2 = 2.0 * arg1;
      output_a(n + M) = (float32)params_a(0) *
	((float32)params_a(1) + (float32)params_a(2) * Integral::cos(arg1) +
	 ((float32)params_a(2) - (float32)params_a(1)) * Integral::cos(arg2));
    }
  }
  
  // case 2: an even number of points
  //  resample an (N-1) point window
  //
  else {
    
    // generate a window for (N-1) points
    // sample it at a new time scale 
    //
    for (int32 n = 0; n < N; n++) {
      float32 scaled_index = -M + scale * n;
      float32 arg1 = 2.0 * Integral::PI * scaled_index / (float32)Mp1;
      float32 arg2 = 2.0 * arg1;
      output_a(n) = (float32)params_a(0) *
	((float32)params_a(1) + (float32)params_a(2) * Integral::cos(arg1) +
	 ((float32)params_a(2) - (float32)params_a(1)) * Integral::cos(arg2));
    }
  }
  
  // exit gracefully
  //
  return true;
}

// method: generateDolphChebyshev
//
// arguments: 
//  VectorFloat& output: (output) the window function
//  const VectorFloat& params: (input) window function parameters
//  int32 num_points: (input) number of points in the window
//
// return: a bool8 value indicating status
//
// this method generates a Dolph-Chebyshev window function:
//
//  w(n) = { G/(2M+1) * [1/gamma + 
//         {             2 * sum Cheb(beta*k*PI/(2M+1)) *
//         {             cos(2*n*k*PI/(2M+1))]                  n <= |M|
//         { 0                                                  elsewhere
//
//  there are a couple of important auxiliary equations:
//
//   alpha_s = (2*M * 2.285 * (float32)dw + 16.4) / 2.056
//
//  where dw is the transition bandwidth in normalized frequency.
//  using alpha_s, we compute gamma and beta (needed above):
//
//   gamma = 10 ** (-alpha_s / 20.0)
//   beta = cosh((1/2M) * acosh(1/gamma))
//
// a good reference for this window can be found at:
//
//  S.K. Mitra, Digital Signal Processing,
//  McGraw-Hill, Boston, Massachusetts, USA, 2001, pp. 456.
//
// there are two parameters allowed for this window:
//
//  params(0): G, the gain of the window
//  params(1): dw, the normalized transition bandwidth
//
// for an odd number of points, the window is defined as above.
// for an even number of points, an interpolation technique is used.
//
bool8 Window::generateDolphChebyshev(VectorFloat& output_a,
				       const VectorFloat& params_a,
				       int32 num_points_a) const {
  
  // compute some constants related to the length of the window
  //
  int32 N = num_points_a;
  int32 Nm1 = N - 1;
  int32 Nm2 = Nm1 - 1;
  float32 scale = (float32)Nm2 / (float32)Nm1;
  
  // compute the half-window length
  //
  int32 M = (N - 1) / 2;
  
  // create output space
  //
  output_a.setLength(num_points_a);
  
  // in this approach, the constant represents the transition bandwidth.
  // the window length is supplied by the user. we compute the appropriate
  // constants to satisfy this relationship.
  //
  float32 A = (2*M * 2.285 * (float32)params_a(1) + 16.4) / 2.056;
  float32 gamma = Integral::pow(10.0, -A / Integral::DB_MAG_SCALE_FACTOR);
  float32 gamma_inv = 1.0 / gamma;
  
  // from gamma and M, we compute beta
  //
  float32 beta = Integral::cosh(Integral::acosh(gamma_inv) / (float32)(2*M));
  
  // case 1: an odd number of points
  //
  if ((N % 2) != 0) {
    for (int32 n = -M; n <= M; n++) {
      
      // compute the inner summation
      //
      float32 sum = 0;
      for (int32 k = 1; k < M; k++) {
	float32 arg1 = (float32)k * Integral::PI / (float32)(N);
	float32 arg2 = 2 * n * arg1;
	float32 val;
	Chebyshev::compute(val, beta * Integral::cos(arg1), k);
	sum += val * Integral::cos(arg2);
      }
      
      // save the output value
      //
      output_a(n + M) = params_a(0) / (float32)(N) * (gamma_inv + 2 * sum);
    }
  }
  
  // case 2: an even number of points
  //  use a resampling strategy
  //
  else {
    for (int32 n = 0; n < N; n++) {
      float32 scaled_index = -M + scale * n;
      
      // compute the inner summation
      //
      float32 sum = 0;
      for (int32 k = 1; k < M; k++) {
	float32 arg1 = (float32)k * Integral::PI / (float32)(Nm1);
	float32 arg2 = 2 * scaled_index * arg1;
	float32 val;
	Chebyshev::compute(val, beta * Integral::cos(arg1), k);
	sum += val * Integral::cos(arg2);
      }
      
      // save the output value
      //
      output_a(n) = params_a(0) / (float32)(Nm1) * (gamma_inv + 2 * sum);
    }
  }
  
  // exit gracefully
  //
  return true;
}

// method: generateGaussian
//
// arguments: 
//  VectorFloat& output: (output) the window function
//  const VectorFloat& params: (input) window function parameters
//  int32 num_points: (input) number of points in the window
//
// return: a bool8 value indicating status
//
// this method generates a Gaussian window function:
//
//                            2                  2
//  w(n) = { G * exp (-0.5 * k * (tan(theta_0/2))  ,           n <= |M|
//         { 0                                                 elsewhere
//
// theta_0 is the width of the main lobe in radians, and is related
// to the frequency response of the window.
// a good reference for this window can be found at:
//
//  R.A. Roberts and C.T. Mullis, Digital Signal Processing,
//  Addison Wesley Publishing Company, Reading, Massachusetts, USA,
//  1987, pp. 183-188.
//
// there are two parameters allowed for this window:
//
//  params(0): G, the gain of the window
//  params(1): theta_0, the main lobe width
//
// for an odd number of points, the window is defined as above.
// for an even number of points, an interpolation technique is used.
//
bool8 Window::generateGaussian(VectorFloat& output_a,
				 const VectorFloat& params_a,
				 int32 num_points_a) const {
  
  // compute some constants related to the length of the window
  //
  int32 N = num_points_a;
  int32 Nm1 = N - 1;
  int32 Nm2 = Nm1 - 1;
  float32 scale = (float32)Nm2 / (float32)Nm1;
  
  // compute the half-window length
  //
  int32 M = (N - 1) / 2;
  
  // define some useful constants related to the Gaussian function
  //
  float32 theta_0 = params_a(1);
  float32 arg1 = (float32)(theta_0 / 2.0);
  float32 arg2 = Integral::tan(arg1);
  
  // create output space
  //
  output_a.setLength(num_points_a);
  
  // case 1: an odd number of points
  //
  if ((N % 2) != 0) {
    
    // loop over all elements
    //
    for (int32 n = -M; n <= M; n++) {
      output_a(n + M) = (float32)params_a(0) *
	Integral::exp(-0.5 * (float32)n * (float32)n * arg2 * arg2);
    }
  }
  
  // case 2: an even number of points
  //  use a resampling strategy
  //
  else {
    for (int32 n = 0; n < N; n++) {
      float32 scaled_index = -M + scale * n;
      output_a(n) = (float32)params_a(0) *
	Integral::exp(-0.5 * (float32)scaled_index * (float32)scaled_index *
		      arg2 * arg2);
    }
  }
  
  // exit gracefully
  //
  return true;
}

// method: generateGeneralizedHanning
//
// arguments: 
//  VectorFloat& output: (output) the window function
//  const VectorFloat& params: (input) window function parameters
//  int32 num_points: (input) number of points in the window
//
// return: a bool8 value indicating status
//
// this method generates a generalized Hanning window function:
//
//  w(n) = { G * (alpha - (1 - alpha) cos (2 * PI * n / (N-1)))  0 <= n < N
//         { 0                                                   elsewhere
//
// the constant alpha can be used to implement a variety of windows,
// including a Hamming (alpha = 0.54) and Hanning (alpha = 0.50).
// a good reference for this window can be found at:
//
//  L.R. Rabiner and R.W. Schafer, Digital Processing of Speech Signals,
//  Prentice-Hall, Englewood Cliffs, New Jersey, USA, 1978, pp. 121.
//
// there are two parameters allowed for this window:
//
//  params(0): G, the gain of the window
//  params(1): alpha, the main lobe width
//
// this window works the same way for both odd and even samples.
//
bool8 Window::generateGeneralizedHanning(VectorFloat& output_a,
					   const VectorFloat& params_a,
					   int32 num_points_a) const {
  
  // compute some constants related to the length of the window
  //
  int32 N = num_points_a;
  int32 Nm1 = N - 1;
  
  // create output space
  //
  output_a.setLength(N);
  
  // loop over all elements
  //
  for (int32 n = 0; n < N; n++) {
    float32 arg1 = 2.0 * Integral::PI * (float32)n / (float32)(Nm1);
    output_a(n) = (float32)params_a(0) * ((float32)params_a(1) -
				 (1 - (float32)params_a(1)) * Integral::cos(arg1));
  }
  
  // exit gracefully
  //
  return true;
}

// method: generateKaiser
//
// arguments: 
//  VectorFloat& output: (output) the window function
//  const VectorFloat& params: (input) window function parameters
//  int32 num_points: (input) number of points in the window
//
// return: a bool8 value indicating status
//
// this method generates a Kaiser window function:
//
//                                   2
//  w(n) = { I (beta * sqrt(1 - (n/M) ) / I (beta)             n <= |M|
//            0                            0
//       = { 0                                                 elsewhere
//
//  there are a couple of important auxiliary design equations:
//
//   alpha_s = 2.285 * N * dw + 8
//
//  where dw is the transition bandwidth in normalized radian frequency and
//  N is the window length. from alpha_a, we compute beta using the
//  following heuristically derived equations:
//
//          { 0.1102(alpha_s - 8.7), 		     alpha_s > 50
//   beta = { 0.5842(alpha_s - 21)**0.4 +
//                0.07886(alpha_s - 21)              21 <= alpha_s <= 50
//          { 0                                      alpha_s < 21
//
// a good reference for this window can be found at:
//
//  S.K. Mitra, Digital Signal Processing,
//  McGraw-Hill, Boston, Massachusetts, USA, 2001, pp. 456.
//
// there are two parameters allowed for this window:
//
//  params(0): G, the gain of the window
//  params(1): dw, the normalized transition bandwidth
//
// for an odd number of points, the window is defined as above.
// for an even number of points, an interpolation technique is used.
//
bool8 Window::generateKaiser(VectorFloat& output_a,
			       const VectorFloat& params_a,
			       int32 num_points_a) const {
  
  // compute some constants related to the length of the window
  //
  int32 N = num_points_a;
  int32 Nm1 = N - 1;
  int32 Nm2 = Nm1 - 1;
  float32 scale = (float32)Nm2 / (float32)Nm1;
  
  // compute the half-window length
  //
  int32 M = (N - 1) / 2;
  
  // create output space
  //
  output_a.setLength(N);
  
  // convert the transition bandwidth (delta_w) and window length (M) to
  // an attenuation in dB (alpha_s)
  //
  float32 alpha_s = N * 2.285 * params_a(1) + 8.0;
    
  // convert the attenuation to beta
  //
  float32 beta;

  if (alpha_s >= 50.0) {
    beta = 0.1102 * (alpha_s - 8.7);
  }
  else if ((alpha_s > 21) && (alpha_s < 50)) {
    beta = 0.5842 * Integral::pow(alpha_s - 21.0, 0.4) +
      0.07886 * (alpha_s - 21.0);
  }
  else {
    return Error::handle(name(), L"generateKaiser",
			 ERR_PRM, __FILE__, __LINE__);
  }
  
  // precompute the normalization factor
  //
  float32 num, denom;
  Bessel::compute(denom, beta, 0);
    
  // case 1: an odd number of points
  //
  if ((N % 2) != 0) {
    for (int32 n = -M; n <= M; n++) {

      // compute the argument and then the Bessel function value
      //
      float32 n_M = (float32)n / (float32)M;
      float32 arg1 = beta * Integral::sqrt(1.0 - n_M * n_M);
      Bessel::compute(num, arg1, 0);

      // save the output value
      //
      output_a(n + M) = params_a(0) * num / denom;
    }
  }

  // case 2: an even number of points
  //  use a resampling technique
  //
  else {
    for (int32 n = 0; n < N; n++) {

      // compute the argument and then the Bessel function value
      //
      float32 scaled_index = -M + scale * n;
      float32 n_M = Integral::min(scaled_index / (float32)M, 1.0);
      float32 arg1 = beta * Integral::sqrt(1.0 - n_M * n_M);
      Bessel::compute(num, arg1, 0);

      // save the output value
      //
      output_a(n) = params_a(0) * num / denom;
    }
  }
  // exit gracefully
  //
  return true;
}

// method: generateLifter
//
// arguments: 
//  VectorFloat& output: (output) the window function
//  const VectorFloat& params: (input) window function parameters
//  int32 num_points: (input) number of points in the window
//
// return: a bool8 value indicating status
//
// this method generates a raised sine, often called a lifter,
// window function:
//
//  w(n) = { G * (1 + h * sin(n * PI / L))          1 <= n <= L
//         { 0                                      elsewhere
//
// the constant h is normally set to L/2. note that L != N.
// a good reference for this window can be found at:
//
//  L. Rabiner and B.H. Juang, Fundamentals of Speech Recognition,
//  Prentice-Hall, Englewood Cliffs, New Jersey, USA, 1993, pp. 169.
//
// there are two parameters allowed for this window:
//
//  params(0): G, the gain of the window
//  params(1): L, the raised sine frequency
//  params(2): h, the raised sine scale factor
//
bool8 Window::generateLifter(VectorFloat& output_a,
			       const VectorFloat& params_a,
			       int32 num_points_a) const {
  
  // compute some constants related to the length of the window
  //
  int32 N = num_points_a;
  float32 L = params_a(1);
  float32 h = params_a(2);

  // create output space
  //
  output_a.setLength(N);
  output_a(0) = 0.0;
  
  // loop over all elements
  //
  for (int32 n = 1; n < N; n++) {

    // check the limits
    //
    if (n <= L) {
      float32 arg1 = (float32)n * Integral::PI / L;
      output_a(n) = (float32)params_a(0) * (1 + h * Integral::sin(arg1));
    }
    else {
      output_a(n) = 0;
    }
  }
  
  // exit gracefully
  //
  return true;
}

// method: getLeadingPad
//
// arguments: none 
//
// return: a int32 value indicating leading pad (in units of frames)
//
int32 Window::getLeadingPad() const {

  // in FRAME_INTERNAL mode, we only deal with the current frame
  //
  if (cmode_d == FRAME_INTERNAL) {
    return 0;
  }
  
  // for a right-aligned window, there is no leading pad
  //
  if (alignment_d == RIGHT) {
    return 0;
  }

  // for a left-aligned window, we need to compute the overlap
  // (fixed a previous bug by doing calculations in num samples instead of seconds)
  //
  int32 frame_nsamp = (int32)Integral::round(frame_dur_d * sample_freq_d);
  int32 duration_nsamp = (int32)Integral::round(duration_d * sample_freq_d);
  float32 overlap = (duration_nsamp - frame_nsamp) / (float32)frame_nsamp;
  //float32 overlap = (duration_d - frame_dur_d) / frame_dur_d;

  if (alignment_d == LEFT) {
    return (int32)Integral::ceil(overlap);
  }

  // for a center-aligned window, we need to compute the overlap
  //
  else if (alignment_d == CENTER) {
    return (int32)Integral::ceil(overlap / 2);
  }

  // else: exit ungracefully
  //
  else {
    return Error::handle(name(), L"getLeadingPad",
			 ERR_PRM, __FILE__, __LINE__);
  }
}

// method: getTrailingPad
//
// arguments: none 
//
// return: a int32 value indicating the trailing pad (in units of frames)
//
int32 Window::getTrailingPad() const {

  // in FRAME-INTERNAL mode, we only deal with the current frame
  //
  if (cmode_d == FRAME_INTERNAL) {
    return 0;
  }
  
  // for a right-aligned window, there is no leading pad
  //
  if (alignment_d == LEFT) {
    return 0;
  }

  // for a right-aligned window, we need to compute the overlap
  // (fixed a previous bug by doing calculations in num samples instead of seconds)
  //
  int32 frame_nsamp = (int32)Integral::round(frame_dur_d * sample_freq_d);
  int32 duration_nsamp = (int32)Integral::round(duration_d * sample_freq_d);
  float32 overlap = (duration_nsamp - frame_nsamp) / (float32)frame_nsamp;
  //float32 overlap = (duration_d - frame_dur_d) / frame_dur_d;

  if (alignment_d == RIGHT) {
    return (int32)Integral::ceil(overlap);
  }

  // for a center-aligned window, we need to compute the overlap
  //
  else if (alignment_d == CENTER) {
    return (int32)Integral::ceil(overlap / 2);
  }

  // else: exit ungracefully
  //
  else {
    return Error::handle(name(), L"getTrailingPad", ERR_PRM,
			 __FILE__, __LINE__);
  }
}