celp.c

this the source code of addio compression standard CELP. Also, it is optimizied for the execution sp

/**************************************************************************
*                                                                         *
*       CELP Voice Decoder --- A fast version                             *
*       Base on the original Version 3.2c.                                *
*       Copyright CMU   1994.5.
*                                                                         *
*                                                                         *
**************************************************************************/
#define TRUE		1
#define FALSE		0
#define STREAMBITS	144
#define CODELENGTH1	15
#define CODELENGTH2	11
#define PARITYLENGTH	(CODELENGTH1 - CODELENGTH2)
#define SYNDRUN		100
#define OMEGA		0.994127  /* Bandwidth expansion for LPC analysis (15 Hz) */
#define ALPHA		0.8  /* Bandwidth expansion for postfilter */
#define BETA		0.5  /* Bandwidth expansion for postfilter */
#define mmax(A,B)        ((A)>(B)?(A):(B))
#define mmin(A,B)        ((A)<(B)?(A):(B))
/*  #define nint(F)          (((F)>0)?(int)(F+0.5):(int)(F-0.5))  */
#include <stdio.h>
#include <math.h>
#include <strings.h>
#include "ccsub.h"

int cbgbits = 5, cbindex = 0, gindex = 0, idb = 0, ncsize = 512, no = 10;
int nseg = 0, pindex = 0, frame = 0, tauptr = 0, minptr = 0, plevel1 = 0;
int plevel2 = 0, pbits[MAXNP + 2] = {8, 6, 5, 0, 0};
int mxsw = 1;

float bb[MAXNP + 1], e0[MAXLP];
float fc[MAXNO + 1], fcn[MAXNO + 1], fci[MAXNO + 1];
float h[MAXLP], gamma2 = 0.8, prewt = 0.0;

  /* *read adaptive code book index (pitch delay) file		*/
float pdelay[MAXPD] = 
{
#include "pdelay.h"
};

  /* *load stochastic code book vector	file			*/
float x[MAXCODE] = 
{
#include "codebook.h"
};

char ptype[10] = "max2", cbgtype[10] = "log";
char pstype[10] = "hier";

main(argc, argv)
int argc;
char *argv[];

{

  int cbi[MAXLL / MAXL], framesnr = 0;
  int framedm = 0, framedm2 = 0, i, j, k, l = 60, lun[10];
  int ll = 240, lp = 60, nn, np = 1, nrec1 = 1, nrec2 = 1, nrec3 = 1;
  int findex[MAXNO];
  int error = 0, total = 0, flag;
  int pdtabi[MAXPD], sync = 1;

  short iarf[MAXLL], npf[MAXLL], pf[MAXLL];

  float cbg[MAXLL / MAXL], pgs[MAXLL / MAXL];
  float sold[MAXLL], snew[MAXLL], ssub[MAXLL], v[MAXLL];
  float vdecoded[MAXLL], rcn[MAXNO], hamw[MAXLL], hamws[MAXL], dpa[MAXPA];
  float dps[MAXPA], newfreq[MAXNO], unqfreq[MAXNO], lsp[MAXLL / MAXL][MAXNO];
  float dppa[MAXPA], dpps[MAXPA]; 
  float scale = 1.0, decodedgain, dm[9], dm2[9], taus[4], descale = 1.0;

  /* *load pitch delay coding tables for bit assignment		*/
  /* *pdencode.h for encoding, pddecode.h for decoding 		*/
  static int pdencode[MAXPD] = 
  {
#include "pdencode.h"
  };

  static float pddecode[MAXPD] = 
  {
#include "pddecode.h"
  };

  /* *filter coefficients for 2nd order 100 Hz HPF with 60 Hz notch:	*/
  /* static float ahpf[3] = {1.0, -1.99778, 1.0};  */
  /* static float bhpf[3] = {1.0, -1.88 0.89};  */

  /* *filter coefficients for 2nd order Butterworth 100 Hz HPF:		*/
  static float ahpf[3] = {0.946, -1.892, 0.946};
  static float bhpf[3] = {1.0, -1.889033, 0.8948743};

  /* *filter coefficients for 2nd order Butterworth 275 Hz HPF:		*/
  /* static float ahpfo[3] = {0.858, -1.716, 0.858}; */
  /* static float bhpfo[3] = {1.0, -1.696452, 0.7368054}; */

  static float ahpfo[3] = {0.946, -1.892, 0.946};
  static float bhpfo[3] = {1.0, -1.889033, 0.8948743};


/*
		*bit stream 
*/

  int cbbits = 9, pointer, bitpointer, bitsum1, bitsum2;
  int mask[STREAMBITS], ssum, pstream[STREAMBITS]; 
  static int sbits[MAXNO] = {3, 4, 4, 4, 4, 3, 3, 3, 3, 3};
  short stream[STREAMBITS], savestream[STREAMBITS];
  char line[38];

/*
		*filter memories (should be maxno+1)
*/

  static float dhpf1[3], dhpf2[3], dsa[MAXNO+1], dss[MAXNO+1];
  static float dhpf1o[3], dhpf2o[3];
  static float dp1[MAXNO+1], dp2[MAXNO+1], dp3[2];
  static float ip, op, sumsnr, sumdm[10], sumdm2[10];

/*
		*error control coding parameters:
*/
  float ber = 0.0, realerror, syndavg = 0.0;
  int codeword[CODELENGTH1], hmatrix[CODELENGTH1];
  int syndrometable[CODELENGTH1], paritybit, twoerror, protect;
  int syndrome, eccbits = 5;	

  /* *load bit protection vector 				*/
  static int bitprotect[CODELENGTH2] = 
  {
#include "bitprot.h"
  };

  /* *load bit permutation vector				 */
  static int bitpermute[STREAMBITS] = 
  {
#include "bitperm.h"
  };

  char tempstr[82];

  static char ifile[82];
  static char ofile[82];
  static char stype[12]="kang";

  FILE *fopen(), *fp25;
  FILE  *fp125, *fp126, *fp127;

  int nrec9 = 1;

  /* *start CELP                                         */

//strcpy(ifile, *++argv);
//strcpy(ofile, *++argv);
strcpy(ifile, argv[1]);
strcpy(ofile, argv[2]);




  /* *intialize mask						*/

  for (i = 0; i < STREAMBITS; i++)
  {
    mask[i] = 0;
  }

  /* ********************* initialize********************	 */

  /* *number of codewords/LPC frame				 */

  nn = ll / l;

  /* *dimension of d1a and d1b???				 */

  idb = MMAX + MAXNP - 1 + l;
  plevel1 = 1 << pbits[0];

  /* *levels of delta tau					 */

  plevel2 = 1 << pbits[1];

  /* *number of bits per subframe				 */

  bitsum1 = cbbits + cbgbits + pbits[0] + pbits[2];
  bitsum2 = cbbits + cbgbits + pbits[1] + pbits[2];

  /* *enable/disable error control coding			 */

  protect = TRUE;

  /* *for double error detecting FEC codes (NOT USED)		*/

  twoerror = FALSE;

  /* *intialize arrays						*/

  for (i = 0; i < MAXLP; i++) h[i] = e0[i] = 0.0;
  for (i = 0; i < MAXLL; i++) sold[i] = 0.0;
  for (i = 0; i < STREAMBITS; i++) stream[i] = savestream[i] = 0;

  /* *start nseg at 0 to do pitch on odd segments 		*/
  /*   (nseg is incremented before csub) 			*/

  nseg = 0;

  /* *generate matrix for error control coding			 */

  matrixgen(CODELENGTH1, CODELENGTH2, hmatrix, syndrometable);


  /* *** open and define files					 */

  fp125 = fopen(strcat(strcpy(tempstr, ofile), ".spd"), "w");
  if (fp125 == NULL)
  {
    perror("celp: Error opening the output.spd file");
    exit(0);
  }

  fp126 = fopen(strcat(strcpy(tempstr, ofile), "npf.spd"), "w");
  if (fp126 == NULL)
  {
    perror("celp: Error opening the output.spd file");
    exit(0);
  }

  fp127 = fopen(strcat(strcpy(tempstr, ofile), "hpf.spd"), "w");
  if (fp127 == NULL)
  {
    perror("celp: Error opening the output.spd file");
    exit(0);
  }


  /* *bit stream channel file					*/
  fp25 = fopen(ifile, "r");
  if (fp25 == NULL)
  {
    perror("celp: Error opening the channel file");
    exit(0);
  }


  /* *generate pdtabi for delta delay coding			*/
  for (i = 0; i < MAXPD; i++)
  {
    pdtabi[pdencode[i]] = i;
  }

  /* ......................... m a i n  l o o p ........................ */

  /* *** CHANNEL ..................................................	 */

  /* *read in channel file (if synthesis only)				 */
  while (fscanf(fp25, "%s", line) != EOF)
  {
    pointer = -1; 
    gethex(STREAMBITS, pstream, line);
    frame++;

    /* *display a propeller (rotating bar) once per frame		 */

//    mark(0);

    /* *unpermute bitstream 						 */

    for (i = 0; i < STREAMBITS; i++)
      stream[bitpermute[i] - 1] = pstream[i];

    /* *** corrupt the channel with bit errors 				 */

    biterror(ber, mask, stream, STREAMBITS, &error, &total);



    /* *** SYNTHESIS ..........................................		 */

    /* *unpack parity bits						 */

    if (protect)
    {
      pointer = pointer - PARITYLENGTH - 2;

#ifndef ANALY
      pointer = 138;
#endif

      for (i = 0; i < PARITYLENGTH; i++)
	unpack(stream, 1, &codeword[CODELENGTH2 + i], &pointer);

      /* *extract code word from stream array				 */

      for (i = 0; i < CODELENGTH2; i++)
	codeword[i] = stream[bitprotect[i] - 1];

      /* repack Bisnu bit (remains constant for now)			 */

      codeword[10] = 0;

      /* *Hamming decode						 */

      decodeham(CODELENGTH1, hmatrix, syndrometable, paritybit, codeword,
		 &twoerror, &syndrome);

      /* *disable parity check (if parity not used)			 */
      twoerror = FALSE;

      /* *bit error rate estimator (running avg of bad syndromes)	 */

      if (syndrome != 0)
        syndrome = 1;
      syndavg = (1.0 - (1.0 / SYNDRUN)) * syndavg + (1.0 / SYNDRUN) * (float) syndrome;

      /* *repack protected bits						 */

      for (i = 0; i < CODELENGTH2; i++)
	stream[bitprotect[i] - 1] = codeword[i];

      /* *frame repeat if two errors detected in code word		 */

      if (twoerror)
	printf("celp: two errors have occured in frame %d\n", frame);

    }

    pointer = -1;

    /* *unpack data stream						 */

    for (i = 0; i < no; i++)
      unpack(stream, sbits[i], &findex[i], &pointer);

    /* *decode lsp's							 */

    lspdecode34(findex, no, newfreq);

    /* *interpolate spectrum lsp's for nn subframes			 */

    intsynth(newfreq, nn, lsp, twoerror, syndavg);


    /* *decode all code book and pitch parameters			 */

    bitpointer = pointer;
    dcodtau(pbits[0], pbits[1], bitsum1, bitsum2, &bitpointer, nn, stream, pddecode, pdtabi, taus);
    dcodpg(pbits[2], bitsum1, bitsum2, &bitpointer, nn, stream, pgs);
    dcodcbi(cbbits, bitsum1, bitsum2, &bitpointer, nn, stream, cbi);
    dcodcbg(cbgbits, bitsum1, bitsum2, &bitpointer, nn, stream, cbg);

    /* *** synthesize each subframe					 */

    k = 0;
    for (i = 0; i < nn; i++)
    {
      nseg++;


      /* *decode values for subframe					 */

      cbindex = cbi[i];
      decodedgain = cbg[i];
      if (protect)
	smoothcbgain(&decodedgain, twoerror, syndavg, cbg, i + 1);

      /* *code book synthesis						 */

      vdecode(decodedgain, l, &vdecoded[k]);

      if (protect)
	smoothtau(&taus[i], twoerror, syndavg, taus[2], i + 1);
      bb[0] = taus[i];
      bb[2] = pgs[i];
      if (protect)
	smoothpgain(&bb[2], twoerror, syndavg, pgs, i + 1);


      /* *pitch synthesis						 */

      pitchvq(&vdecoded[k], l, dps, idb, bb, "long");

      /* *pitch pre filter (synthesis)					 */

      if (prewt != 0.0) 
        prefilt(&vdecoded[k], l, dpps);

      /* convert lsp's to pc's 						 */

      lsptopc(&lsp[i][0], fci);	

      /* lpc synthesis	 						 */ 

      polefilt(fci, no, dss, &vdecoded[k], l);	



      /* *** write nonpostfiltered output speech disk files		 */

      for (j = 0; j < l; j++)
      {
        vdecoded[k + j] = descale * vdecoded[k + j];
	npf[k + j] = nint(mmax(-32768., mmin(32767., vdecoded[k + j])));
      }

      /* *write npf output speech file "ofilenpf"			 */

      if ( fwrite(&npf[k], sizeof(short), l, fp126)  != l)
      {
        fprintf(stderr, "celp: *** Error writing ofilenpf.spd\n");
        exit(1);
      }


      /* *** post filtering						 */

      postfilt(&vdecoded[k], l, ALPHA, BETA, &ip, &op, dp1, dp2, dp3);

      /* *** test for output speech clipping				 */

      while (clip(&vdecoded[k], l))
      {

      /* frame repeat & recall synthesizer? 				 */
      /* or scale vdecoded?						 */

	printf("celp: Clipping detected at frame %d\n", frame);
        for (j = 0; j < l; j++)
          vdecoded[k + j] = 0.05 * vdecoded[k + j]; 
      }

      /* *** write postfiltered output speech disk files		 */

      for (j = 0; j < l; j++)
	pf[k + j] = nint(mmax(-32768., mmin(32767., vdecoded[k + j])));

      /* *write output speech file "ofile"				 */

      if ( fwrite(&pf[k], sizeof(short), l, fp125)  != l)
      {
        fprintf(stderr, "celp: *** Error writing ofile.spd\n");
        exit(1);
      }


      /* *high pass filter output speech	 */

/*
      zerofilt_S(ahpfo, 2, dhpf1o, &vdecoded[k], l);
      polefilt_S(bhpfo, 2, dhpf2o, &vdecoded[k], l);

      for (j = 0; j < l; j++)
        pf[k + j] = nint(mmax(-32768., mmin(32767., vdecoded[k + j])));


      if ( fwrite(&pf[k], sizeof(short), l, fp127)  != l)
      {
        fprintf(stderr, "celp: *** Error writing ofilehpf.spd\n");
        exit(1);
      }
*/
      k += l;
    }

  } /* end main while loop; either the analyzer loop, analyzer and       */
    /* sungraph loop or the synthesizer only loop	 		 */


  /* ...........e n d  m a i n  l o o p ...............................  */


  fclose(fp25);
  fclose(fp125);
  fclose(fp126);
  fclose(fp127);
  exit(0);
}