#include "nrutil.h"
#include <stdio.h>


/* ================Bilinear transformed cepstrum coefficients ======= */
/* Input :
          l  - number of cepstrum coefficients
	  c  - cepstrum vector, c[1],c[2],...,c[l]
	  m  - Order of bilinear transform
	  a  - Mel-warping constant -1<a<1 
  Output:
          cc - Bilinear transformed cepstrum coefficients, cc[0],..,cc[m-1]
----------------------------------------------------------------------- */
void BLT_CEPS (l, c, m, cc, a)
int l,m;
float c[],cc[],a;
{
    int i,k,n;
    int mCoefs;                   /* Number of mel-warped cepstrum coeffs */
    float       *d0, *d1;         /* Temporary arrays */
    float *cep,                   /* Copy of cepstrum coefficients */
          prev,prev1,             /* Temporary delayed variables */
	  *inp,*out;              /* Input/output vectors of k-stage */
    float inp1;                   /* Delayed input */
    float out1;                   /* Delayed output */
    float   small = 0.1e-20,
           gain;                  /* gain factor, 1-a^2 */

    mCoefs=m-1;         /* Number of mel-warped cepstrum coeffs */
    n=l-1;              /* Number of cepstrum coefficients, 0..l-1 */
    cep = vector(0,n);  /* Cepstrum vector */
    inp = vector(0,n);  /* Input vector used in the 2,3,..m stage */
    out = vector(0,n);  /* Output vector used in the 2,3,..m stage */
    d0  = vector(0,n);  /* Temporary array used in stages 0 and 1  */
    d1  = vector(0,n);  /* Temporary array used in stages 0 and 1 */

    for (i=0;i<=n;++i)
      cep[i]=c[i+1]; 

   /* -------------- 0th stage ------------ */
   /*                      1 
              cep(i)-->  ----------  -->  d0(i)  
                         1 - a z^-1		       
    ----------------------------------------- */

    prev=cep[n];
    for (i=n-1; i>0; --i) {
      if (i != (n-1)) prev  = d0[i+1];
      d0[i] = cep[i] + a*prev;
      if (fabs(d0[i]) < small)
	d0[i]=0.0;
    }
    cc[0] = cep[0] + a*d0[1];  /* -- 0th mel-warped cepstrum coef */

   /* -------------- 1st stage ------------ */
   /*                    z^-1 
             d0(i)-->  ----------  ---> d1(i)
                       1 - a z^-1		       
    ---------------------------------------- */

    prev=cep[n];
    prev1=0.0;
    for (i=n-1; i>0; --i) {
      if (i!=(n-1)) prev1  = d1[i+1];
      d1[i] = prev + a*prev1;
      if (fabs(d1[i])<small) d1[i]=0.0;
      prev  = d0[i];
    }
    
    cc[1] = d0[1] + (a * d1[1]);   /* -- 1st mel-warped cepstrum coef */


    if (mCoefs == 1){
      gain = 1.0 - a * a;
      for (i = 1; i <= mCoefs; i++)
	cc[i] *= gain;
      return;
    }

    /* - Initialize input to 2nd stage to be the output from 1st stage */
    /**/
    for (i=1;i<n;++i)
      inp[i]=d1[i];

   /* -------------- * k-th stage ( 2 through m )  ------------ */
   /*                   z^-1 - a
             inp(i)-->  ----------  ---> out(i)
                       1 - a z^-1		       
    ----------------------------------------------------------- */
     

    for (k = 2; k <= mCoefs; k++)
    {
        inp1=0.0;
	out1=0.0;
        for (i=n-1; i>0; i--) {
	  if (i!= (n-1)) { out1=out[i+1];
			   inp1=inp[i+1];
			 }
	  out[i] = inp1 + a*(out1-inp[i]);
	  if (fabs(out[i]) < small) out[i]=0.0;
	}

	cc[k] = inp[1] + a * (out[1]-cc[k-1]);  /* kth mel-warped cepstrum */

	/* -- Feed the output as an input to the next stage -- */
	/**/
	for (i=1;i<n; ++i)
	  inp[i]=out[i];   
    }

    /* --- adjustment gain factor (1.0 -a**2) from the first stage---- */
    /**/
    gain = 1.0 - a * a;
    for (i = 1; i <= mCoefs; i++)
	cc[i] *= gain;

  free_vector(cep,0,n);
  free_vector(inp,0,n);
  free_vector(out,0,n);
  free_vector(d0,0,n);
  free_vector(d1,0,n);


 
}