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

Routine:
  void FAfiltAP (AFILE *AFpI, AFILE *AFpO, long int NsampO, float h[],
                 int Ncof, long int loffs)

Purpose:
  Filter an audio file with an all-pole filter

Description:
  This routine convolves the data from the input audio file with an all-pole
  filter response.

Parameters:
   -> AFILE *AFpI
      Audio file pointer for the input audio file
   -> AFILE *AFpO
      Audio file pointer for the output audio file
   -> long int NsampO
      Number of output samples to be calculated
   -> float h[]
      Array of Ncof all-pole filter coefficients
   -> long int loffs
      Data offset into the input data for the first output point

Author / revision:
  P. Kabal  Copyright (C) 1996
  $Revision: 1.10 $  $Date: 1996/06/01 02:46:11 $

-------------------------------------------------------------------------*/

static char rcsid[] = "$Id: FAfiltAP.c 1.10 1996/06/01 AFsp-V2R1 $";

#include <stdio.h>
#include <libtsp.h>
#include "FiltAudio.h"

#define MINV(a, b)	(((a) < (b)) ? (a) : (b))
#define MAXV(a, b)	(((a) > (b)) ? (a) : (b))

#define NBUF	5120
#define MAXWUP	1000

void
FAfiltAP (AFpI, AFpO, NsampO, h, Ncof, loffs)

     AFILE *AFpI;
     AFILE *AFpO;
     long int NsampO;
     const float h[];
     int Ncof;
     long int loffs;

{
  float x[NBUF];
  int mem;
  int Nxmax;
  int Nx;
  long int l, k;

/*
   Notes:
   - The input signal d(.) is the data in the file, with d(0) corresponding to
     the first data value in the file.
   - Indexing: l is an offset into d(), referring to sample d(l).
*/

/* Batch processing
   - The data will be processed in batches by reading into a buffer x(.,.).
     The batches of input samples will be of equal size, Nx, except for the
     last batch.  For batch j,
       x(j,l') = d(loffs+j*Nx+l'), for 0 <= l' < Nx,
   - The k'th output point y(k) is calculated at position d(loffs+k), that is
     the start of the impulse response, h(0), is aligned with d(loffs+k).
       y(k) --> h[0] <==> d(l),    where l=loffs+k
                h[0] <==> x(j,l').
   - For batch j=0,
       l = loffs  - pointer to d(loffs),
       l' = 0     - pointer to x(0,0) = d(loffs),
       k = 0      - pointer to y(0).
   - For each batch, k and l' advance by Nx,
       k <- k + Nx,
       l' <- l' + Nx.
   - When the index l' for x(j,l') advances beyond Nx, we bring l' back
     into range by subtracting Nx from it and incrementing the batch number,
*/

/* Buffer allocation
   Let the total buffer size be Nb.  This is allocated to filter memory (mem)
   and the input data (Nx).  The output data will overlay the input data.
*/
  mem = Ncof - 1;
  Nxmax = NBUF - mem;
  if (Nxmax <= 0)
    UThalt ("FAfiltAP: Too many filter sections");
  l = loffs;

/* Warm up points; same as the main loop but with no AFwriteData */
  VRfZero (x, mem);
  l = MINV (loffs, MAXV (0, loffs - MAXWUP));
  while (l < loffs) {
    Nx = (int) MINV (Nxmax, loffs - l);
    AFreadData (AFpI, l, &x[mem], Nx);
    l = l + Nx;
    FIfiltAP (&x[mem], x, Nx, h, Ncof);
    VRfShift (x, mem, Nx);
  }

/* Main processing loop */
  k = 0;
  while (k < NsampO) {

/* Read the input data into the input buffer */
    Nx = (int) MINV (Nxmax, NsampO - k);
    AFreadData (AFpI, l, &x[mem], Nx);
    l = l + Nx;

/* Convolve the input samples with the filter response */
    FIfiltAP (&x[mem], x, Nx, h, Ncof);

/* Write the output data to the output audio file */
    AFwriteData (AFpO, &x[mem], Nx);
    k = k + Nx;

/* Update the filter memory */
    VRfShift (x, mem, Nx);
  }

  return;
}