delaynr.c
上传用户:szhypcb168
上传日期:2007-01-06
资源大小:2187k
文件大小:8k
源码类别:

语音压缩

开发平台:

Unix_Linux

  1. #define NFRAC 6
  2. #define TRUE 1
  3. #define FALSE 0
  4. #define M1 -4
  5. #define M2  3
  6. /* five fractional delays plus zero delay = 6 */
  7. static float frac[NFRAC] = {0.25, 0.33333334, 0.5, 0.66666667, 0.75, 0.0};
  8. static int twelfths[NFRAC] = {3, 4, 6, 8, 9, 0};
  9. /**************************************************************************
  10. *
  11. * NAME
  12. * sinc
  13. *
  14. * FUNCTION
  15. *
  16. *
  17. * SYNOPSIS
  18. * subroutine sinc(arg)
  19. *
  20. *   formal 
  21. *                       data    I/O
  22. *       name            type    type    function
  23. *       -------------------------------------------------------------------
  24. * arg float i
  25. * y(n) float func
  26. *
  27. ***************************************************************************
  28. *
  29. * DESCRIPTION
  30. *
  31. * This interpolating (sinc) function is Hamming windowed to bandlimit
  32. * the signal to reduce aliasing.
  33. *
  34. ***************************************************************************
  35. *
  36. * CALLED BY
  37. *
  38. * pitchvq  psearch
  39. *
  40. * CALLS
  41. *
  42. *
  43. ***************************************************************************/
  44. #include <math.h>
  45. #include <stdio.h>
  46. float
  47. sinc(arg)
  48. float arg;
  49. {
  50.   float pi;
  51.   pi = 4.0 * atan(1.0);
  52.   if (arg == 0.0)
  53.     return(1.0);
  54.   else
  55.     return(sin(pi * arg) / (pi * arg));
  56. /**************************************************************************
  57. *
  58. * NAME
  59. * qd_nr
  60. *
  61. * FUNCTION
  62. *
  63. *
  64. * SYNOPSIS
  65. * subroutine qd_nr(d)
  66. *
  67. *   formal 
  68. *                       data    I/O
  69. *       name            type    type    function
  70. *       -------------------------------------------------------------------
  71. * d float i
  72. * qd_nr(d) int o quantize d function
  73. *
  74. ***************************************************************************
  75. *
  76. * DESCRIPTION
  77. *
  78. * Quantize d function
  79. *
  80. ***************************************************************************
  81. *
  82. * CALLED BY
  83. *
  84. * delay_nr
  85. *
  86. * CALLS
  87. *
  88. *
  89. ***************************************************************************/
  90. int 
  91. qd_nr(d)
  92. float d;
  93. {
  94.   int i, index, ok;
  95.   /* *the second and fourth elements of the data statement below */
  96.   /* when multiplied by 3 must be greater than 1 or 2 respectively. */
  97.   ok = FALSE;
  98.   for (i = 0; i < NFRAC - 1; i++)
  99.   {
  100.     if (fabs(d - frac[i]) < 1.e-2)
  101.     {
  102.       index = i;
  103.       ok = TRUE;
  104.     }
  105.   }
  106.   if (!ok)
  107.   {
  108.     fprintf(stderr, "delay: Invalid pitch delay = %fn", d);
  109.     exit(1);
  110.   }
  111.   return(index);
  112. }
  113. /**************************************************************************
  114. *
  115. * NAME
  116. * delay_nr
  117. *
  118. * FUNCTION
  119. * Time delay a bandlimited signal
  120. * using blockwise recursion (aka "nonrecursive")
  121. *
  122. * SYNOPSIS
  123. * subroutine delay_nr(x, start, n, d, m, y)
  124. *
  125. *   formal 
  126. *                       data    I/O
  127. *       name            type    type    function
  128. *       -------------------------------------------------------------------
  129. * x(n) float i signal input
  130. * start int i beginning of output sequence
  131. * n int i length of input sequence
  132. * d float i fractional pitch
  133. * m int i integer pitch
  134. * y(n) float o delayed input signal
  135. *
  136. ***************************************************************************
  137. *
  138. * DESCRIPTION
  139. *
  140. * Subroutine to time delay a bandlimited signal by resampling the
  141. * reconstructed data (aka sinc interpolation).  The well known
  142. * reconstruction formula is:
  143. *
  144. *              |    M2      sin[pi(t-nT)/T]    M2
  145. *   y(n) = X(t)| = SUM x(n) --------------- = SUM x(n) sinc[(t-nT)/T]
  146. *              |   n=M1         pi(t-nT)/T    n=M1
  147. *              |t=n+d
  148. *
  149. * The interpolating (sinc) function is Hamming windowed to bandlimit
  150. * the signal to reduce aliasing.
  151. *
  152. * Multiple simultaneous time delays may be efficiently calculated
  153. * by polyphase filtering.  Polyphase filters decimate the unused
  154. * filter coefficients.  See Chapter 6 in C&R for details. 
  155. *
  156. ***************************************************************************
  157. *
  158. * CALLED BY
  159. *
  160. * pitchvq   psearch
  161. *
  162. * CALLS
  163. *
  164. * ham
  165. *
  166. ***************************************************************************
  167. *
  168. * REFERENCES
  169. *
  170. * Crochiere & Rabiner, Multirate Digital Signal Processing,
  171. * P-H 1983, Chapters 2 and 6.
  172. *
  173. * Kroon & Atal, "Pitch Predictors with High Temporal Resolution,"
  174. * ICASSP '90, S12.6
  175. *
  176. * CELP Documentation, Version 3.1, Figure 9, 15 Dec '89.
  177. *
  178. **************************************************************************/
  179. #define SIZE (M2 - M1 + 1)
  180. delay_nr(x, start, n, d, m, y)
  181. float x[], d, y[];
  182. int m, n, start;
  183. {
  184.   static float wsinc[SIZE][NFRAC], hwin[12*SIZE+1];
  185.   float pitch, pitch2, pitch3;
  186.   static int first = TRUE;
  187.   int i, j, k, index, k2, k3, mtwo;
  188.   /* Generate Hamming windowed sinc interpolating function for each */
  189.   /* allowable fraction.  The interpolating functions are stored in   */
  190.   /* time reverse order (i.e., delay appears as advance) to align     */
  191.   /* with the adaptive code book v0 array.  The interp filters are:   */
  192.   /*  wsinc[.,0] frac = 1/4 (3/12)       */
  193.   /*  wsinc[.,1] frac = 1/3 (4/12)       */
  194.   /*  . .       */
  195.   /*  wsinc[.,4] frac = 3/4 (9/12)       */
  196.   /*  wsinc[.,5] frac = 0 zero delay for programming ease*/
  197.   if (first)
  198.   {
  199.     ham(hwin, 12*SIZE+1);
  200.     for (i = 0; i < NFRAC; i++)
  201.       for (k = M1, j = 0; j < SIZE; j++)
  202.       {
  203.         wsinc[j][i] = sinc(frac[i] + k) * hwin[12*j+twelfths[i]];
  204.         k++;
  205.       }
  206.     first = FALSE;
  207.   }
  208.   index = qd_nr(d);
  209.   /* *Resample:   */
  210.   /* *If delay is greater than n(=60), then lay down n values at   */
  211.   /* *the appropriate delay.   */
  212.   if (m >= n)
  213.   {
  214.     for (i = 0; i < n + M1; i++)
  215.     {
  216.       y[i] = 0.0;
  217.       for (k = M1, j = 0; j < SIZE; j++)
  218.       {
  219.         y[i] += x[i+k+start-m-1] * wsinc[j][index];
  220.         k++;
  221.       }
  222.     }
  223.     for (i = 0; i < M2; i++)
  224.       x[start+i-1] = y[i];
  225.     for (i = n+M1; i < n; i++)
  226.     {
  227.       y[i] = 0.0;
  228.       for (k = M1, j = 0; j < SIZE; j++)
  229.       {
  230.         y[i] += x[i+k+start-m-1] * wsinc[j][index];
  231.         k++;
  232.       }
  233.     }
  234.   }
  235.   else
  236.   /* *If delay is less than n/2 (=30), then lay down two pitch periods  */
  237.   /* *with the appropriate delay, and part of third pitch period with */
  238.   /* *the appropriate delay. */
  239.   {
  240.     pitch = m + d;
  241.     if (m < n / 2)
  242.     {
  243.       for (i = 0; i < m + M1; i++)
  244.       {
  245.         y[i] = 0.0;
  246.         for (k = M1, j = 0; j < SIZE; j++)
  247.         {
  248.           y[i] += x[i+k+start-m-1] * wsinc[j][index];
  249.           k++;
  250.         }
  251.       }
  252.       for (i = 0; i < M2 + 1; i++)
  253.       {
  254.         x[start+i-1] = y[i];
  255.       }
  256.       for (i = m + M1; i < m; i++)
  257.       {
  258.         y[i] = 0.0;
  259.         for (k = M1, j = 0; j < SIZE; j++)
  260.         {
  261.           y[i] += x[i+k+start-m-1] * wsinc[j][index];
  262.           k++;
  263.         }
  264.        }
  265.       pitch2 = 2 * pitch;
  266.       k2 = (int)(pitch2 - 2 * m);
  267.       mtwo = (int)(pitch2);
  268.       d = pitch2 - (2 * m + k2);
  269.       if (fabs(d) > 0.1 && fabs(d) < 0.9)
  270.         index = qd_nr(d);
  271.       else 
  272.         index = 5;
  273.       for (i = m; i < mtwo; i++)
  274.       {
  275.         y[i] = 0.0;
  276.         for (k = M1, j = 0; j < SIZE; j++)
  277.         {
  278.           y[i] += x[i+k+start-2*m-k2-1] * wsinc[j][index];
  279.           k++;
  280.         }
  281.       }
  282.       pitch3 = 3 * pitch;
  283.       k3 = (int)(pitch3 - (3 * m + k2));
  284.       d = pitch3 - (3 * m + k2 + k3);
  285.       if (fabs(d) > 0.1 && fabs(d) < 0.9)
  286.         index = qd_nr(d);
  287.       else
  288.         index = 5;
  289.       for (i = mtwo; i < n; i++)
  290.       {
  291.         y[i] = 0.0;
  292.          for (k = M1, j = 0; j < SIZE; j++)
  293.          {
  294.            y[i] += x[i+k+start-m-mtwo-k3-k2-1] * wsinc[j][index];
  295.            k++;
  296.          }
  297.       }     
  298.     }
  299.     else
  300.     /* *If delay is greater than n/2 (=30), then lay down one pitch */
  301.     /* *period with the appropriate delay, and part of the second pitch */
  302.     /* *period with appropriate delay.    */
  303.     {
  304.       for (i = 0; i < m + M1; i++)
  305.       {
  306.         y[i] = 0.0;
  307.         for (k = M1, j = 0; j < SIZE; j++)
  308.         {
  309.           y[i] += x[i+k+start-m-1] * wsinc[j][index];
  310.           k++;
  311.         }
  312.       }
  313.       for (i = 0; i < M2 + 1; i++)
  314.         x[start+i-1] = y[i];
  315.       for (i = m + M1; i < m; i++)
  316.       {
  317.         y[i] = 0.0;
  318.         for (k = M1, j = 0; j < SIZE; j++)
  319.         {
  320.           y[i] += x[i+k+start-m-1] * wsinc[j][index];
  321.           k++;
  322.         }
  323.       }
  324.       pitch2 = 2 * pitch;
  325.       k2 = (int)(pitch2 - 2 * m);
  326.       d = pitch2 - (2 * m + k2);
  327.       if (fabs(d) > 0.1 && fabs(d) < 0.9)
  328.         index = qd_nr(d);
  329.       else
  330.         index = 5;
  331.       for (i = m; i < n; i ++)
  332.       {
  333.         y[i] = 0.0;
  334.         for (k = M1, j = 0; j < SIZE; j++)
  335.         {
  336.           y[i] += x[i+k+start-2*m-k2-1] * wsinc[j][index];
  337.           k++;
  338.         }
  339.       }
  340.     }
  341.   }
  342.   for (i = 0; i < M2 + 1; i++)
  343.     x[start+i-1] = 0;
  344. }