postfil.c
上传用户:touchwatch
上传日期:2007-01-06
资源大小:168k
文件大小:10k
源码类别:

语音压缩

开发平台:

Unix_Linux

  1. /*************************************************************************/
  2. /*                                                                       */
  3. /*                            LD-CELP  G.728                             */
  4. /*                                                                       */
  5. /*    Low-Delay Code Excitation Linear Prediction speech compression.    */
  6. /*                                                                       */
  7. /*    Code edited by Michael Concannon.                                  */
  8. /*    Based on code written by Alex Zatsman, Analog Devices 1993         */
  9. /*                                                                       */
  10. /*************************************************************************/
  11. /* Postfilter of LD-CELP Decoder */
  12. #include "common.h"
  13. #include "prototyp.h"
  14. #include "parm.h"
  15. #include "data.h"
  16. #define ABS(X) ((X)<0.0?(-X):(X))
  17. /*  Parameters from the adapter: */
  18. #define SPORDER 10
  19. #define DECIM 4
  20. #define PMSIZE  (NPWSZ+KPMAX) /* Postfilter Memory SIZE */
  21. #define PDMSIZE (PMSIZE/DECIM)  /* Post. Decim. Memory SIZE */
  22. #define DPERMAX (KPMAX/DECIM) /* Max. Decimated Period */
  23. #define DPERMIN (KPMIN/DECIM) /* Min. Decimated Period */
  24. #define SHIFTSZ (KPMAX+NPWSZ-NFRSZ+IDIM)
  25. static real tap_mem[KPMAX+NPWSZ+IDIM]; /* Post-Filter Memory for syn. sp. */
  26. int  pitch_period = 50;
  27. real pitch_tap=0.0;
  28. real pitch_gain=1;
  29. static real shzscale[SPORDER+1]= /* Precomputed Scales for IIR coefficients */
  30. { 1.0, 0.6500244140625, 0.4224853515625, 0.27459716796875, 
  31.   0.17852783203125, 0.11602783203125, 0.075439453125, 
  32.   0.04901123046875, 0.0318603515625, 0.02069091796875, 0.01348876953125};
  33. static real shpscale[SPORDER+1]= /* Precomputed Scales for FIR Coefficients */
  34. { 1.0, 0.75, 0.5625, 0.421875, 0.31640625, 0.2373046875, 0.177978515625, 
  35.   0.13348388671875, 0.10009765625, 0.0750732421875, 0.05633544921875};
  36. static real
  37.     shpcoef[SPORDER+1], /* Short Term Filter (Poles/IIR) Coefficients */
  38.     shzcoef[SPORDER+1], /* Short Term Filter (Zeros/FIR) Coefficients */
  39.     tiltz,
  40.     fil_mem[PMSIZE]; /* Post-Filter Memory for residual */
  41. static void  longterm(real[], real[]);
  42. static void shortterm(real[], real[]);
  43. /* Compute sum of absolute values of vector V */
  44. static real vec_abs(real v[])
  45. {
  46.     int i;
  47.     real r = ABS(v[0]);
  48.     for(i=1; i<IDIM; i++)
  49. r+=ABS(v[i]);
  50.     return r;
  51. }
  52.     /* Inverse Filter */
  53. void inv_filter(real input[])
  54. {
  55.   int i,j,k;
  56.   static int ip=IDIM;
  57.   static real mem1[SPORDER+NFRSZ];
  58.   
  59.   /** Shift in input into mem1 **/
  60.   for(i=IDIM; i<SPORDER+IDIM; i++)
  61.     mem1[i-IDIM] = mem1[i];
  62.   for(i=0; i<IDIM; i++)
  63.     mem1[SPORDER+i] = input[i];
  64.   for(k=0; k<IDIM; k++) {
  65.     real tmp = mem1[SPORDER+k];
  66.     for(j=1; j<=SPORDER; j++)
  67.       tmp += mem1[SPORDER+k-j]*a10[j];
  68.     fil_mem[PMSIZE-NFRSZ+ip+k] = tmp;
  69.   }
  70.   if(ip == (NFRSZ-IDIM))
  71.     ip = 0;
  72.   else
  73.     ip += IDIM;
  74. }
  75. void
  76. postfilter(real input[], real output[])
  77. {
  78.     int i;
  79.     static
  80.     real
  81. temp[IDIM],  /* Output of long term filter*/
  82. temp2[IDIM],  /* Input of short term filter*/
  83. new_gain, /* Gain of filtered output */
  84. input_gain,  /* Gain of input */
  85. scale; /* Scaling factor for gain preservation */
  86.     static real scalefil=1.0; /* Smoother version of scale */
  87.     static real vec_abs();
  88.     longterm(input, temp);
  89.     shortterm(temp, temp2);
  90.     /* Computed scale for gain preservation */
  91.     new_gain = vec_abs(temp2);
  92.     if (new_gain > 1.0)
  93.     {
  94. input_gain = vec_abs(input);
  95. scale = input_gain/new_gain;
  96.     }
  97.     else 
  98. scale = 1.0;
  99.     
  100.     /* Smooth out scale, then scale the output */
  101.     for(i=0; i<IDIM; i++) {
  102. scalefil = AGCFAC * scalefil + (1.0 - AGCFAC)*scale;
  103. output[i] = scalefil * temp2[i];
  104.     }
  105. }
  106. static void
  107. longterm(real input[], real output[])
  108. {
  109.     int i;
  110.     real out;
  111.     static real lmemory[KPMAX];
  112.     /* Add weighted pitch_period-delayed signal */
  113.     for(i=0; i<IDIM; i++) {
  114. out = pitch_tap * lmemory[KPMAX+i-pitch_period];
  115. out += input[i];
  116. output[i] = pitch_gain*out;
  117.     }
  118.     
  119.     /* Shift-in input to lmemory */
  120.     for (i=0; i<KPMAX-IDIM; i++)
  121. lmemory[i] = lmemory[i+IDIM];
  122.     for(i=0; i<IDIM; i++)
  123. lmemory[KPMAX-IDIM+i] = input[i];
  124. }
  125. /*
  126.   Again, memories (shpmem, shzmem) are in reverse order, 
  127.   i.e. [0] is the oldest.
  128.  */
  129. static void
  130. shortterm(real input[], real output[])
  131. {
  132.     int k,j;
  133.     
  134.     static real shpmem[SPORDER], shzmem[SPORDER];
  135.     for(k=0; k<IDIM; k++) {
  136.       
  137.       /* FIR Part */
  138.       
  139.       real in = input[k], out;
  140.       out = in;
  141.       for(j=SPORDER-1; j>=1; j--) {
  142. out += shzmem[j]*shzcoef[j+1];
  143. shzmem[j] = shzmem[j-1];
  144.       }
  145.       out += shzmem[0] * shzcoef[1];
  146.       shzmem[0] = in;
  147.       
  148.       /* IIR Part */
  149.       
  150. for(j=SPORDER-1; j>=1; j--) {
  151.     out -= shpmem[j]*shpcoef[j+1];
  152.     shpmem[j] = shpmem[j-1];
  153. }
  154. out -= shpmem[0] * shpcoef[1];
  155. shpmem[0] = out;
  156. output[k] = out+tiltz*shpmem[1];
  157.     }
  158. }
  159. /*********************************/
  160. /***** Postfilter Adapter ********/
  161. /*********************************/
  162. static int extract_pitch();
  163. void
  164. psf_adapter (real frame[])
  165. {
  166.     pitch_period = extract_pitch();
  167.     /** Compute Pitch Tap **/
  168.     {
  169.       int i;
  170.       real corr=0.0, corr_per=0.0;
  171.       /** Shift old memory **/
  172.       for(i=0;i<SHIFTSZ;i++)
  173. tap_mem[i] = tap_mem[i+NFRSZ];
  174.       /** Shift new frame into memory **/
  175.       for(i=0;i<NFRSZ;i++)
  176. tap_mem[SHIFTSZ+i] = frame[i];
  177.       for(i=KPMAX-pitch_period; i<(KPMAX-pitch_period+NPWSZ); i++) {
  178.   corr     += tap_mem[i] * tap_mem[i];
  179.   corr_per += tap_mem[i] * tap_mem[i+pitch_period];
  180.       }
  181.       if REALZEROP(corr)
  182.   pitch_tap = 0.0;
  183.       else
  184.   pitch_tap = corr_per/corr;
  185.     }
  186.     
  187.     /** Compute Long Term Coefficients **/
  188.     
  189.     {
  190.       if (pitch_tap > 1)
  191.   pitch_tap = 1.0;
  192.       if (pitch_tap < PPFTH)
  193.   pitch_tap = 0.0;
  194.       pitch_tap = PPFZCF * pitch_tap;
  195.       pitch_gain = 1.0/(1.0+pitch_tap);
  196.     }
  197. };
  198. void compute_sh_coeff()
  199. {
  200.   /** Compute Short Term Coefficients **/
  201.   {
  202.     int i;
  203.     for(i=1; i<=SPORDER; i++) {
  204. shzcoef[i] = shzscale[i]*a10[i];
  205.   shpcoef[i] = shpscale[i]*a10[i];
  206.       }
  207.       tiltz = TILTF * k10;
  208.   }
  209. }
  210. static int best_period (real buffer[], int buflen,
  211. int pmin, int pmax)
  212. {
  213.     int i, per, best_per = -1;
  214.     real best_corr = -(BIG);
  215.     for(per = pmin; per<pmax; per++) {
  216. real corr = 0.0;
  217. for(i=pmax; i<buflen; i++)
  218.     corr += buffer[i] * buffer[i-per];
  219. if (corr > best_corr) {
  220.     best_corr = corr;
  221.     best_per  = per;
  222. }     
  223.     }
  224.     return best_per;
  225. }
  226. #define DCFRSZ NFRSZ/DECIM /* size of decimated frame */
  227. static int
  228. extract_pitch()
  229. {
  230.     int
  231. i, j, k,
  232. best_per=KPMAX,  /* Best Period (undecimated) */
  233. best_dper = KPMAX/DECIM, /* Best Decimated Period */
  234. best_old_per=KPMAX,  /* Best Old Period */
  235. permin,  /* Limits for search of best period */
  236. permax;
  237.     real
  238. best_corr=-(BIG), best_old_corr=-(BIG), tap0=0.0, tap1=0.0;
  239.     static int old_per = (KPMIN+KPMAX)>>1;
  240.     static real
  241. fil_decim_mem[PDMSIZE],
  242. fil_out_mem[NFRSZ+DECIM];
  243. #define FIL_DECIM_MEM(I) fil_decim_mem[I]
  244. #define FIL_OUT_MEM(I)   fil_out_mem[I]
  245.     
  246.     /** Shift decimated filtered output */
  247.     for(i=DCFRSZ; i<PDMSIZE; i++)
  248. FIL_DECIM_MEM(i-DCFRSZ) = FIL_DECIM_MEM(i);
  249.   
  250.     /* Filter and  decimate  input */
  251.   {
  252.       int decim_phase = 0, dk;
  253.       for (k = 0, dk=0; k<NFRSZ; k++)
  254.       {
  255.   real tmp;
  256.   tmp = fil_mem[PMSIZE-NFRSZ+k] - A1 * FIL_OUT_MEM(2)
  257.                  - A2 * FIL_OUT_MEM(1)
  258.                  - A3 * FIL_OUT_MEM(0);
  259.   decim_phase++;
  260.   if (decim_phase == 4) {
  261.       FIL_DECIM_MEM(PDMSIZE-DCFRSZ+dk) = 
  262.   B0 * tmp
  263. + B1 * FIL_OUT_MEM(2)
  264. + B2 * FIL_OUT_MEM(1)
  265. + B3 * FIL_OUT_MEM(0);
  266.       decim_phase = 0;
  267.       dk++;
  268.   }
  269.   FIL_OUT_MEM(0) = FIL_OUT_MEM(1);
  270.   FIL_OUT_MEM(1) = FIL_OUT_MEM(2);
  271.   FIL_OUT_MEM(2) = tmp;
  272.       }
  273.   }
  274.     /* Find best Correlation in decimated domain: */
  275.     best_dper = best_period(fil_decim_mem, PDMSIZE, DPERMIN, DPERMAX);
  276.     /* Now fine-tune best correlation on undecimated  domain */
  277.     
  278.     permin = best_dper * DECIM - DECIM + 1;
  279.     permax = best_dper * DECIM + DECIM - 1;
  280.     if (permax > KPMAX)
  281. permax = KPMAX;
  282.     if (permin < KPMIN)
  283. permin = KPMIN;
  284.     
  285.   {
  286.       int per;
  287.       best_corr = -(BIG);
  288.       for(per = permin; per<=permax; per++) {
  289.   real corr = 0.0;
  290.   for(i=1,j=(per+1);   i<= NPWSZ;   i++,j++)
  291.       corr += fil_mem[PMSIZE-i]*fil_mem[PMSIZE-j];
  292.   if (corr > best_corr) {
  293.       best_corr = corr;
  294.       best_per = per;
  295.   }
  296.       }
  297.   }
  298.     
  299.     /** If we are not exceeding old period by too much, we have a real
  300.        period and not a multiple */
  301.     
  302.     permax = old_per + KPDELTA;
  303.     if (best_per <= permax)
  304. goto done;
  305.     /** Now compute best period around the old period **/
  306.     
  307.     permin = old_per - KPDELTA;
  308.     if (permin<KPMIN) permin = KPMIN;
  309.   {
  310.       int per;
  311.       best_old_corr = -(BIG);
  312.       for(per = permin; per<=permax; per++) {
  313.   real corr = 0.0;
  314.   for(i=1,j=(per+1);
  315.       i<=NPWSZ;
  316.       i++,j++)
  317.       corr += fil_mem[PMSIZE-i]*fil_mem[PMSIZE-j];
  318.   if (corr > best_old_corr) {
  319.       best_old_corr = corr;
  320.       best_old_per = per;
  321.   }
  322.       }
  323.   }
  324.     
  325. /***** Compute the tap ****/
  326.   {
  327.       real s0=0.0, s1=0.0;
  328.       for(i=1; i<=NPWSZ; i++) {
  329. s0 += fil_mem[PMSIZE-i-best_per] * fil_mem[PMSIZE-i-best_per];
  330. s1 += fil_mem[PMSIZE-i-best_old_per] * fil_mem[PMSIZE-i-best_old_per];
  331.       }
  332.       if (! REALZEROP(s0))
  333.   tap0 = best_corr/s0;
  334.       if (! REALZEROP(s1))
  335.   tap1 = best_old_corr/s1;
  336.       tap0 = CLIPP(tap0, 0.0, 1.0);
  337.       tap1 = CLIPP(tap1, 0.0, 1.0);
  338.       if (tap1 > TAPTH * tap0)
  339.   best_per = best_old_per;
  340.   }
  341.   done:
  342.     /** Shift fil_mem **/
  343.     for(i=NFRSZ; i<PMSIZE; i++)
  344. fil_mem[i-NFRSZ] = fil_mem[i];
  345.     old_per = best_per;
  346.     return best_per;
  347. }
  348. void
  349. init_postfilter()
  350. {
  351.     int i;
  352.     shzscale[0] = shpscale[0] = 1.0;
  353.     for (i=1; i<=SPORDER; i++)
  354.     {
  355. shzscale[i] = SPFZCF * shzscale[i-1];
  356. shpscale[i] = SPFPCF * shpscale[i-1];
  357.     }
  358. }