/*-------------- Telecommunications & Signal Processing Lab ---------------
                             McGill University

Routine:
  void RSKaiserLPF (float h[], int N, double Fc, double alpha, double Gain,
                    double Woffs, double Wspan)

Purpose:
  Generate a Kaiser windowed lowpass filter

Description:
  This routine calculates the impulse response of a lowpass filter which is
  obtained by applying a Kaiser window to a sin(x)/x function.  Consider the
  continuous-time version of this response
                           sin(2 pi Fc t)
    h(t) = G * Kw(t/T, a)  -------------- ,
                                pi t
  where Kw(t/T,a) is a Kaiser window with parameter a, spanning -T <= t <= T.
  The cutoff of the lowpass filter is Fc.  The gain of the filter is set to
  approximate G (Gain) in the passband.

  References:
  J. F. Kaiser, "Nonrecursive digital filter design using the I0-sinh window
  function", Proc. 1974 IEEE Int. Symp. on Circuits and Syst., pp. 20-23,
  April 1974.

  A. Antoniou, Digital Filters: Analysis, Design and Applications, 2nd ed.,
  McGraw-Hill, 1993.

Parameters:
  <-  float h[]
      Array containing the impulse response samples.  The reference point in
      the array is at offset (N-1)/2 - soffs.
   -> int N
      Number of impulse response samples
   -> double Fc
      Lowpass filter cutoff (normalized frequency, 0 to 1/2)
   -> double alpha
      Kaiser window parameter
   -> double Gain
      Low frequency gain
   -> double Woffs
      Offset of the first sample from the beginning of the non-zero part of
      the window
   -> double Wspan
      Span of the non-zero part of the Kaiser window in samples (2T in the
      formulation above)

Author / revision:
  P. Kabal  Copyright (C) 1996
  $Revision: 1.5 $  $Date: 1996/07/16 14:32:12 $

-------------------------------------------------------------------------*/

static char rcsid[] = "$Id: RSKaiserLPF.c 1.5 1996/07/16 AFsp-V2R1 $";

#include <libtsp.h>
#include <libtsp/nucleus.h>
#include "ResampAudio.h"

void
RSKaiserLPF (h, N, Fc, alpha, Gain, Woffs, Wspan)

     float h[];
     int N;
     double Fc;
     double alpha;
     double Gain;
     double Woffs;
     double Wspan;

{
  int i, k;
  double t, u, T, Ga;

  /* Multiply the Kaiser window by the sin(x)/x response */
  /* Use symmetry to reduce computations */
  Ga = 2.0 * Fc * Gain;
  T = 0.5 * Wspan;
  for (i = 0, k = N-1; i <= k; ++i, --k) {
    t = i - T + Woffs;
    u = k - T + Woffs;
    h[i] = Ga * FNsinc (2.0 * Fc * t) * FIxKaiser (t / T, alpha);
    if (-t == u)
      h[k] = h[i];
    else
      h[k] = Ga * FNsinc (2.0 * Fc * u) * FIxKaiser (u / T, alpha);
  }

  return;
}