/*
 *  These functions calculate the values for quadrature mirror filter
 *  coefficient arrays.
 *
 * Copyright (C) 1991--94 Wickerhauser Consulting.  All Rights Reserved.
 * May be freely copied for noncommercial use.  See
 * ``Adapted Wavelet Analysis from Theory to Software'' ISBN 1-56881-041-5
 * by Mladen Victor Wickerhauser [AK Peters, Ltd., Wellesley, Mass., 1994]
 * 
 */

#include <assert.h>
#include <stdlib.h>
#include "real.h"
#include "common.h"
#include "qf.h"		/* Data structs and function prototypes */
#include "oqfs.h"	/* Orthogonal filter coefficients.   */

/********************************************************************
 * mkpqf()
 *
 * Function `mkpqf()' returns a pointer to a `pqf' data structure
 *  containing the coefficients, length and the pre-periodized filter
 *  named by the input parameters.
 *
 * Calling sequence:
 *       mkpqf( coefs, alpha, omega, flags )
 * 
 * Inputs:
 *	(const real *)coefs     These are the filter coefficients,
 *				  assumed to be given in the array
 *				  `coefs[0],...,coefs[omega-alpha]'
 *				  are the only nonzero values.
 *
 *      (int)alpha      These are to be the first and last valid
 *      (int)omega        indices of the pqf->f struct member.
 *
 *	(int)flags	This is reserved for later uses, such as 
 *			  indicating when to generate a full QF
 *			  sequence from just one symmetric half.
 *
 * Return value:
 *	(pqf *)mkpqf  The return value is a pointer to a newly 
 *			  allocated pqf struct containing `coefs[]',
 *			  `alpha', `omega', and the preperiodized
 *			  version of  `coefs[]'.
 * Assumptions:
 *     (1) Conventional indexing: `alpha <= 0 <= omega'.
 */
extern pqf *
  mkpqf(
	const real *coefs,	/* Original filter coefficients. */
	int alpha,		/* Least valid index of `?->f[]'.    */
	int omega,		/* Greatest valid index of `?->f[]'. */
	int flags)		/* Reserved for future use. */
{
  pqf *qm;
  int M;

  assert(alpha <= 0);		/* Conventional indexing. */
  assert(0 <= omega);		/* Conventional indexing. */

  qm = (pqf *)calloc(1,sizeof(pqf)); assert(qm);
  qm->alpha = alpha;
  qm->omega = omega;
  qm->f = coefs-alpha;
  M = IFH( omega-alpha );
  qm->fp = (real *)calloc( PQFL(M), sizeof(real)); assert(qm->fp);
  qfcirc( qm->fp, qm->f, alpha, omega );

  qm->center = coe( qm->f, alpha, omega );
  qm->deviation = lphdev( qm->f, alpha, omega );
  return(qm);
}

/********************************************************************
 * qf()
 *
 * Function `qf()' returns a pointer to a `pqf' data structure
 *  containing the coefficients, name, length and kind of the filter
 *  named by the input parameters.
 *
 * Calling sequence:
 *       qf( name, range, kind )
 * 
 * Inputs:
 *	(const char *)name     This is the name of the filter, specified as
 *				 at least the first letter of "Beylkin",
 *				 "Coifman", "Vaidyanathan", or "Daubechies",
 *  				 written as a string (enclosed in " marks).
 *  				 Case is ignored, only the first letter 
 *				 matters, and "Standard" is an alias for "D".
 *
 *      (int)range      This is the length of the requested filter.  Allowed
 *				values depend on what `name' is.
 *
 *      (int)kind	If kind==LOW_PASS_QF, qf() returns the summing or
 *			  low-pass filter `pqf' structure.
 *			If kind==HIGH_PASS_QF, qf() returns the differencing
 *			  or high-pass filter `pqf' structure.
 *
 * Return value:
 *	(pqf *)qf    If a `name'd filter of the requested `range' and
 *			  `kind' is listed, then the return value is
 *			  a pointer to a newly-allocated pqf struct
 *			  specifying a  conventionally-indexed filter
 *			  to be used for convolution-decimation.
 * 			  Otherwise, the returned pointer is NULL.
 * 
 */
extern pqf *
  qf(				/* Coefficient struct or NULL pointer. */
     const char *name,		/* String "B", "C", "D", or "V". */
     int         range,		/* Length of coefficient array. */
     int         kind)		/* LOW_PASS_QF or HIGH_PASS_QF. */
{
  pqf *qm;

  /* These macros call `mkpqf()' with the right arguments: */
#define MKMKPQF(n, sd)  mkpqf( n##sd##oqf, n##sd##alpha, n##sd##omega, 0 )
#define SETSDPQF(n, k)  ( k==LOW_PASS_QF ) ? MKMKPQF(n, s) : MKMKPQF(n, d)

  qm = 0;			/* Return NULL pointer if unsuccessful. */

  /* Choose an orthogonal filter based on the first letter of
   * the filter name: `name' argument starts with  B, C, D, V
   */
  switch( name[0] ) {

  case 'b':
  case 'B':
    /* Beylkin filters.
     *
     * Here are the coefficients for Gregory BEYLKIN'S wave packets.
     */
    switch( range )
      {
      case 18: qm = SETSDPQF( b18, kind ); break;
      default: break;		/* Fall out: the length is unavailable. */
      }
    break;			/* Beylkin type. */
    
  case 'c':
  case 'C':
    /* Coiflet filters.
     *
     * Here are the coefficients for COIFLETS with respectively
     * 2, 4, 6, 8 and 10 moments vanishing for phi.  Filter Q has
     * range/3 moments vanishing.  Filter P has range/3 moments
     * vanishing with the appropriate shift.
     */
    switch( range ) 
      {
      case 6:  qm = SETSDPQF( c06, kind ); break;
      case 12: qm = SETSDPQF( c12, kind ); break;
      case 18: qm = SETSDPQF( c18, kind ); break;
      case 24: qm = SETSDPQF( c24, kind ); break;
      case 30: qm = SETSDPQF( c30, kind ); break;
      default: break;		/* Fall out: the length is unavailable. */
      }
    break;			/* Coifman type. */
    
  case 's':
  case 'S':			/* STANDARD filters, in old terminology. */
  case 'd':
  case 'D':
    /* Daubechies filters.
     *
     * Initialize quadrature mirror filters P,Q of length 'range' with 
     * smooth limits and minimal phase, as in DAUBECHIES:
     */
    switch( range )
      {
      case 2:  qm = SETSDPQF( d02, kind ); break;
      case 4:  qm = SETSDPQF( d04, kind ); break;
      case 6:  qm = SETSDPQF( d06, kind ); break;
      case 8:  qm = SETSDPQF( d08, kind ); break;
      case 10: qm = SETSDPQF( d10, kind ); break;
      case 12: qm = SETSDPQF( d12, kind ); break;
      case 14: qm = SETSDPQF( d14, kind ); break;
      case 16: qm = SETSDPQF( d16, kind ); break;
      case 18: qm = SETSDPQF( d18, kind ); break;
      case 20: qm = SETSDPQF( d20, kind ); break;
      default: break;		/* Fall out: the length is unavailable. */
      }
    break;			/* Standard or Daubechies type. */
    
  case 'v':
  case 'V':
    /* Vaidyanathan filters 
     * 
     * The following coefficients correspond to one of the filters
     * constructed by Vaidyanathan (Filter #24B from his paper
     * in IEEE Trans. ASSP Vol.36, pp 81-94 (1988). These coefficients
     * give an exact reconstruction scheme, but they don't satisfy ANY
     * moment condition (not even the normalization : the sum of the c_n
     * is close to 1, but not equal to 1). The function one obtains after
     * iterating the reconstruction procedure 9 times looks continuous,
     * but is of course not differentiable. It would be interesting to
     * see how such a filter performs. It has been optimized, for its
     * filter length, for the standard requirements that speech people
     * impose.
     */
    switch( range )
      {
      case 24:  qm = SETSDPQF( v24, kind ); break;
      default: break;		/* Fall out: the length is unavailable. */
      }
    break;			/* Vaidyanathan type. */
    
  default: break;		/* Fall out: the filter is unavailable. */
  }
  return(qm);
}

/*******************************************************************
 * coe()
 *
 * Compute the center-of-energy for a sequence.
 *
 * Basic algorithm:
 *   center =  ( Sum_k k*in[x]^2 ) / (Sum_k in[k]^2 )
 *
 *  Calling sequence:
 *	coe( in, alpha, omega )
 *
 *  Inputs:
 *	(const real *)in	The sequence is `in[alpha],...,in[omega]'
 *
 *      (int)alpha      	These are to be the first and last
 *      (int)omega		  valid indices of the array `in[]'.
 *
 *  Output:
 *	(real)coe		This is in the interval `[alpha, omega]'.
 *
 *  Assumptions:
 *     (1) Conventional indexing: `alpha <= 0 <= omega'.
 */
extern real
  coe(				/* Center of energy. */
      const real *in,		/* Sequence of coefficients. */
      int alpha,		/* First valid `in[]' index.   */
      int omega)		/* Last valid `in[]' index.    */
{
  real center, energy;
  int k;

  assert(alpha <= 0);		/* Conventional indexing. */
  assert(0 <= omega);		/* Conventional indexing. */

  center = 0;  energy = 0;
  for( k=alpha; k<=omega; k++)
    {
      center += k*in[k]*in[k];
      energy += in[k]*in[k];
    }
  if( energy>0 ) center /= energy;
  return(center);
}

/*******************************************************************
 * lphdev()
 *
 * Compute the maximum deviation from linear phase of the
 * convolution-decimation operator associated to a sequence.
 * Basic algorithm:
 *                    Sum_{x>0} (-1)^x Sum_k k*in[k-x]*in[k+x]
 *   deviation =  2 * ----------------------------------------
 *                               Sum_k in[k]^2
 *
 *  Calling sequence:
 *	lphdev( in, alpha, omega )
 *
 *  Inputs:
 *	(const real *)in  The sequence is `in[alpha],...,in[omega]'
 *
 *      (int)alpha      These are to be the first and last valid
 *      (int)omega        indices of the `pqf->f[]' struct member.
 *
 *  Output:
 *	(real)lphdev		This is the absolute value of the maximum.
 *
 * Assumptions:
 *     (1) Conventional indexing: `alpha <= 0 <= omega'.
 */
extern real
  lphdev(			/* Center of energy. */
	 const real *in,	/* Sequence of coefficients. */
	 int alpha,		/* First valid `in[]' index.   */
	 int omega)		/* Last valid `in[]' index.    */
{
  real energy, fx, deviation;
  int k, x, sgn;

  assert(alpha <= 0);		/* Conventional indexing. */
  assert(0 <= omega);		/* Conventional indexing. */

  /* First compute the sum of the squares of the sequence elements: */
  energy = 0;
  for( k=alpha; k<omega; k++)    energy += in[k]*in[k];

  /* If the sum of the squares in nonzero, compute the deviation: */
  deviation = 0;
  if( energy>0 )
    {
      sgn= -1;
      for( x=1;  x <= (omega-alpha)/2;  x++ )
	{
	  fx = 0;
	  for( k=x+alpha; k<=omega-x; k++ )
	    {
	      fx += k*in[k-x]*in[k+x];
	    }
	  deviation += sgn*fx;
	  sgn = -sgn;
	}
      deviation = absval(deviation);
      deviation /= energy;	/* Normalize. */
      deviation *= 2;		/* Final factor from the formula. */
    }
  /* If `energy==0' then `deviation' is trivially 0 */

  return(deviation);
}

/***************************************************************
 * periodize()
 *
 * Periodize an array into a shorter-length array.
 *
 * Calling sequence and basic algorithm:
 *
 * periodize(FQ, Q, F, ALPHA, OMEGA)
 *   For K = 0 to Q-1
 *      Let FQ[K] = 0
 *      Let J = (ALPHA-K)/Q
 *      While K+J*Q <= OMEGA
 *         FQ[K] += F[K+J*Q]
 *         J += 1
 *
 * Inputs:
 *   (real *)fq       Preallocated array of length `q'.
 *   (int)q           This is the period of `fq[]'.
 *   (const real *)f  This array holds the original function. 
 *   (int)alpha       These are to be the first and last valid
 *   (int)omega          indices of the array `f[]'.
 *
 * Output:
 *   (real *)fq       The array `fq[]' is assigned as follows:
 *                      fq[k] = f[k+j0*q]+...+f[k+j1*q],
 *				k = 0, 1, ..., q-1;
 *				j0 = ceil[(alpha-k)/q];
 *				j1 = floor[(omega-k)/q].
 *
 * Assumptions:
 *	(1) `omega-alpha >= q > 0',
 *	(2) `alpha <= 0 <= omega'
 */
static void
  periodize(
	    real *fq,        /* Preallocated output array. */
	    int q,           /* Length of `fq[]'.          */
	    const real *f,   /* Input array.               */
	    int alpha,       /* First valid `f[]' index.   */
	    int omega)       /* Last valid `f[]' index.    */
{
  int j, k;

  assert(q>0);			/* Nontrivial period.     */
  assert(q<=omega-alpha);	/* Periodization needed.  */
  assert(alpha <= 0);		/* Conventional indexing. */
  assert(0 <= omega);		/* Conventional indexing. */

  for(k=0; k<q; k++)
    {
      fq[k] = 0;
      j = (alpha-k)/q;		/* `j==ceil((alpha-k)/q)' */
      while( (k+j*q) <= omega )
	{
	  fq[k] += f[k+j*q];
	  j++;
	}
    }
  return;
}

/***************************************************************
 *  qfcirc()
 *
 * This function computes the periodized filter coefficients
 * associated to the input array of quadrature mirror filter
 * coefficients, for lengths 2, 4, 6,...  It then fills an array
 * with these periodized coefficients, duplicating some of them
 * to facilitate circular convolution.  If the input array is
 *          f[alpha], f[alpha+1], ..., f[omega],
 * and M is the largest integer satisfying `2M<=omega-alpha',
 * then the output array is the concatenation of:
 *
 * Start     Contents                                         End
 * -----  -------------------------------------------------- -----
 *   0            f2[0],f2[1],f2[0],f2[1];                     3
 *   4    f4[0],f4[1],f4[2],f4[3],f4[0],f4[1],f4[2],f4[3];     11
 *  12    f6[0],f6[1],f6[2],...,f6[5],f6[0],f6[1],...,f6[5];   23
 *  24    f8[0],f8[1],f8[2],...,f8[7],f8[0],f8[1],...,f8[7];   39
 *   .          ...         ...                   ...      ;    .
 * S(m-1) f{2m}[0],...,f{2m}[2m-1],f{2m}[0],...,f{2m}[2m-1]; S(m)-1
 *   .          ...         ...             ...            ;    .
 * S(M-1) f{2M}[0],...,f{2M}[2M-1],f{2M}[0],...,f{2M}[2M-1]; S(M)-1
 *
 * where S(m)= 4+8+16+...+ 4(m-1) = 2m(m-1) gives the
 * starting index of the m-th segment.  
 *
 * The "midpoint" or "origin" of the m-th periodic subarray is: 
 *		PQFO(m) = S(m-1)+2m = 2m*m.
 *
 * The total length of concatenated periodic subarrays `1,...,M' is
 *		PQFL(M) = S(M-1)+4M = 2M*(M+1).
 *
 *  Calling sequence and basic algorithm:
 *
 * qfcirc(FP, F, ALPHA, OMEGA)
 *   Let M = IFH(omega-alpha)
 *   For N = 1 to M
 *      Let Q = 2*N
 *      Let FQ = FP + PQFO(N)
 *      periodize( FQ, Q, F, ALPHA, OMEGA )
 *      For K = 0 to Q-1
 *         Let FQ[K-Q] = FQ[K]
 *
 * Inputs:
 *      (real *)fp	This must point to a preallocated array
 *			  with at least `PQFL(M)' memory 
 *			  locations.  This function call fills
 *			  those locations, beginning with 0.
 *
 *	(const real *)f	This is an array of quadrature mirror
 *			  filter coefficients.
 *
 *      (int)alpha      These are to be the first and last valid
 *      (int)omega        indices of the array `f[]'.
 *
 * Output:
 *      (real *)fp	This array is filled with periodized,
 *			  partially duplicated arrays with 
 *			  lengths 4, 8, ..., 4M.
 *
 * Assumptions:
 *     (1) Conventional indexing: `alpha <= 0 <= omega'.
 */
extern void
  qfcirc(
	 real *fp,	   /* Preallocated output array  */
	 const real *f,	   /* Filter coefficients        */
	 int alpha,	   /* First valid `f[]' index    */
	 int omega)	   /* Last valid `f[]' index     */
{
  int k, m, M, q;
  real *fq;

  assert(alpha <= 0);		/* Conventional indexing */
  assert(0 <= omega);		/* Conventional indexing */

  M = IFH(omega-alpha);		/* Max with `2*M <= omega-alpha' */

  for( m=1; m<=M; m++ )
    {
      q  = 2*m;
      fq = fp+PQFO(m);
      periodize( fq, q, f, alpha, omega);
      for( k=0; k<q; k++ )
	fq[k-q] = fq[k];
    }
  return;
}