PITCH_A.C
上传用户:meifeng08
上传日期:2013-06-18
资源大小:5304k
文件大小:18k
源码类别:

语音压缩

开发平台:

C/C++

  1. /*
  2.    ITU-T G.729A Speech Coder    ANSI-C Source Code
  3.    Version 1.1    Last modified: September 1996
  4.    Copyright (c) 1996,
  5.    AT&T, France Telecom, NTT, Universite de Sherbrooke
  6.    All rights reserved.
  7. */
  8. /*---------------------------------------------------------------------------*
  9.  * Pitch related functions                                                   *
  10.  * ~~~~~~~~~~~~~~~~~~~~~~~                                                   *
  11.  *---------------------------------------------------------------------------*/
  12. #include "typedef.h"
  13. #include "basic_op.h"
  14. #include "oper_32b.h"
  15. #include "ld8a.h"
  16. #include "tab_ld8a.h"
  17. /*---------------------------------------------------------------------------*
  18.  * Function  Pitch_ol_fast                                                   *
  19.  * ~~~~~~~~~~~~~~~~~~~~~~~                                                   *
  20.  * Compute the open loop pitch lag. (fast version)                           *
  21.  *                                                                           *
  22.  *---------------------------------------------------------------------------*/
  23. Word16 Pitch_ol_fast(  /* output: open loop pitch lag                        */
  24.    Word16 signal[],    /* input : signal used to compute the open loop pitch */
  25.                        /*     signal[-pit_max] to signal[-1] should be known */
  26.    Word16   pit_max,   /* input : maximum pitch lag                          */
  27.    Word16   L_frame    /* input : length of frame to compute pitch           */
  28. )
  29. {
  30.   Word16  i, j;
  31.   Word16  max1, max2, max3;
  32.   Word16  max_h, max_l, ener_h, ener_l;
  33.   Word16  T1, T2, T3;
  34.   Word16  *p, *p1;
  35.   Word32  max, sum, L_temp;
  36.   /* Scaled signal */
  37.   Word16 scaled_signal[L_FRAME+PIT_MAX];
  38.   Word16 *scal_sig;
  39.   scal_sig = &scaled_signal[pit_max];
  40.   /*--------------------------------------------------------*
  41.    *  Verification for risk of overflow.                    *
  42.    *--------------------------------------------------------*/
  43.    Overflow = 0;
  44.    sum = 0;
  45.    for(i= -pit_max; i< L_frame; i+=2)
  46.      sum = L_mac(sum, signal[i], signal[i]);
  47.   /*--------------------------------------------------------*
  48.    * Scaling of input signal.                               *
  49.    *                                                        *
  50.    *   if Overflow        -> scal_sig[i] = signal[i]>>3     *
  51.    *   else if sum < 1^20 -> scal_sig[i] = signal[i]<<3     *
  52.    *   else               -> scal_sig[i] = signal[i]        *
  53.    *--------------------------------------------------------*/
  54.    if(Overflow == 1)
  55.    {
  56.      for(i=-pit_max; i<L_frame; i++)
  57.      {
  58.        scal_sig[i] = shr(signal[i], 3);
  59.      }
  60.    }
  61.    else {
  62.      L_temp = L_sub(sum, (Word32)1048576L);
  63.      if ( L_temp < (Word32)0 )  /* if (sum < 2^20) */
  64.      {
  65.         for(i=-pit_max; i<L_frame; i++)
  66.         {
  67.           scal_sig[i] = shl(signal[i], 3);
  68.         }
  69.      }
  70.      else
  71.      {
  72.        for(i=-pit_max; i<L_frame; i++)
  73.        {
  74.          scal_sig[i] = signal[i];
  75.        }
  76.      }
  77.    }
  78.    /*--------------------------------------------------------------------*
  79.     *  The pitch lag search is divided in three sections.                *
  80.     *  Each section cannot have a pitch multiple.                        *
  81.     *  We find a maximum for each section.                               *
  82.     *  We compare the maxima of each section by favoring small lag.      *
  83.     *                                                                    *
  84.     *  First section:  lag delay = 20 to 39                              *
  85.     *  Second section: lag delay = 40 to 79                              *
  86.     *  Third section:  lag delay = 80 to 143                             *
  87.     *--------------------------------------------------------------------*/
  88.     /* First section */
  89.     max = MIN_32;
  90.     T1  = 20;    /* Only to remove warning from some compilers */
  91.     for (i = 20; i < 40; i++) {
  92.         p  = scal_sig;
  93.         p1 = &scal_sig[-i];
  94.         sum = 0;
  95.         for (j=0; j<L_frame; j+=2, p+=2, p1+=2)
  96.             sum = L_mac(sum, *p, *p1);
  97.         L_temp = L_sub(sum, max);
  98.         if (L_temp > 0) { max = sum; T1 = i;   }
  99.     }
  100.     /* compute energy of maximum */
  101.     sum = 1;                   /* to avoid division by zero */
  102.     p = &scal_sig[-T1];
  103.     for(i=0; i<L_frame; i+=2, p+=2)
  104.         sum = L_mac(sum, *p, *p);
  105.     /* max1 = max/sqrt(energy)                  */
  106.     /* This result will always be on 16 bits !! */
  107.     sum = Inv_sqrt(sum);            /* 1/sqrt(energy),    result in Q30 */
  108.     L_Extract(max, &max_h, &max_l);
  109.     L_Extract(sum, &ener_h, &ener_l);
  110.     sum  = Mpy_32(max_h, max_l, ener_h, ener_l);
  111.     max1 = extract_l(sum);
  112.     /* Second section */
  113.     max = MIN_32;
  114.     T2  = 40;    /* Only to remove warning from some compilers */
  115.     for (i = 40; i < 80; i++) {
  116.         p  = scal_sig;
  117.         p1 = &scal_sig[-i];
  118.         sum = 0;
  119.         for (j=0; j<L_frame; j+=2, p+=2, p1+=2)
  120.             sum = L_mac(sum, *p, *p1);
  121.         L_temp = L_sub(sum, max);
  122.         if (L_temp > 0) { max = sum; T2 = i;   }
  123.     }
  124.     /* compute energy of maximum */
  125.     sum = 1;                   /* to avoid division by zero */
  126.     p = &scal_sig[-T2];
  127.     for(i=0; i<L_frame; i+=2, p+=2)
  128.         sum = L_mac(sum, *p, *p);
  129.     /* max2 = max/sqrt(energy)                  */
  130.     /* This result will always be on 16 bits !! */
  131.     sum = Inv_sqrt(sum);            /* 1/sqrt(energy),    result in Q30 */
  132.     L_Extract(max, &max_h, &max_l);
  133.     L_Extract(sum, &ener_h, &ener_l);
  134.     sum  = Mpy_32(max_h, max_l, ener_h, ener_l);
  135.     max2 = extract_l(sum);
  136.     /* Third section */
  137.     max = MIN_32;
  138.     T3  = 80;    /* Only to remove warning from some compilers */
  139.     for (i = 80; i < 143; i+=2) {
  140.         p  = scal_sig;
  141.         p1 = &scal_sig[-i];
  142.         sum = 0;
  143.         for (j=0; j<L_frame; j+=2, p+=2, p1+=2)
  144.             sum = L_mac(sum, *p, *p1);
  145.         L_temp = L_sub(sum, max);
  146.         if (L_temp > 0) { max = sum; T3 = i;   }
  147.     }
  148.      /* Test around max3 */
  149.      i = T3;
  150.      p  = scal_sig;
  151.      p1 = &scal_sig[-(i+1)];
  152.      sum = 0;
  153.      for (j=0; j<L_frame; j+=2, p+=2, p1+=2)
  154.          sum = L_mac(sum, *p, *p1);
  155.      L_temp = L_sub(sum, max);
  156.      if (L_temp > 0) { max = sum; T3 = i+(Word16)1;   }
  157.      p  = scal_sig;
  158.      p1 = &scal_sig[-(i-1)];
  159.      sum = 0;
  160.      for (j=0; j<L_frame; j+=2, p+=2, p1+=2)
  161.          sum = L_mac(sum, *p, *p1);
  162.      L_temp = L_sub(sum, max);
  163.      if (L_temp > 0) { max = sum; T3 = i-(Word16)1;   }
  164.     /* compute energy of maximum */
  165.     sum = 1;                   /* to avoid division by zero */
  166.     p = &scal_sig[-T3];
  167.     for(i=0; i<L_frame; i+=2, p+=2)
  168.         sum = L_mac(sum, *p, *p);
  169.     /* max1 = max/sqrt(energy)                  */
  170.     /* This result will always be on 16 bits !! */
  171.     sum = Inv_sqrt(sum);            /* 1/sqrt(energy),    result in Q30 */
  172.     L_Extract(max, &max_h, &max_l);
  173.     L_Extract(sum, &ener_h, &ener_l);
  174.     sum  = Mpy_32(max_h, max_l, ener_h, ener_l);
  175.     max3 = extract_l(sum);
  176.    /*-----------------------*
  177.     * Test for multiple.    *
  178.     *-----------------------*/
  179.     /* if( abs(T2*2 - T3) < 5)  */
  180.     /*    max2 += max3 * 0.25;  */
  181.     i = sub(shl(T2,1), T3);
  182.     j = sub(abs_s(i), 5);
  183.     if(j < 0)
  184.       max2 = add(max2, shr(max3, 2));
  185.     /* if( abs(T2*3 - T3) < 7)  */
  186.     /*    max2 += max3 * 0.25;  */
  187.     i = add(i, T2);
  188.     j = sub(abs_s(i), 7);
  189.     if(j < 0)
  190.       max2 = add(max2, shr(max3, 2));
  191.     /* if( abs(T1*2 - T2) < 5)  */
  192.     /*    max1 += max2 * 0.20;  */
  193.     i = sub(shl(T1,1), T2);
  194.     j = sub(abs_s(i), 5);
  195.     if(j < 0)
  196.       max1 = add(max1, mult(max2, 6554));
  197.     /* if( abs(T1*3 - T2) < 7)  */
  198.     /*    max1 += max2 * 0.20;  */
  199.     i = add(i, T1);
  200.     j = sub(abs_s(i), 7);
  201.     if(j < 0)
  202.       max1 = add(max1, mult(max2, 6554));
  203.    /*--------------------------------------------------------------------*
  204.     * Compare the 3 sections maxima.                                     *
  205.     *--------------------------------------------------------------------*/
  206.     if( sub(max1, max2) < 0 ) {max1 = max2; T1 = T2;  }
  207.     if( sub(max1, max3) <0 )  {T1 = T3; }
  208.     return T1;
  209. }
  210. /*--------------------------------------------------------------------------*
  211.  *  Function  Dot_Product()                                                 *
  212.  *  ~~~~~~~~~~~~~~~~~~~~~~                                                  *
  213.  *--------------------------------------------------------------------------*/
  214. Word32 Dot_Product(      /* (o)   :Result of scalar product. */
  215.        Word16   x[],     /* (i)   :First vector.             */
  216.        Word16   y[],     /* (i)   :Second vector.            */
  217.        Word16   lg       /* (i)   :Number of point.          */
  218. )
  219. {
  220.   Word16 i;
  221.   Word32 sum;
  222.   sum = 0;
  223.   for(i=0; i<lg; i++)
  224.     sum = L_mac(sum, x[i], y[i]);
  225.   return sum;
  226. }
  227. /*--------------------------------------------------------------------------*
  228.  *  Function  Pitch_fr3_fast()                                              *
  229.  *  ~~~~~~~~~~~~~~~~~~~~~~~~~~                                              *
  230.  * Fast version of the pitch close loop.                                    *
  231.  *--------------------------------------------------------------------------*/
  232. Word16 Pitch_fr3_fast(/* (o)     : pitch period.                          */
  233.   Word16 exc[],       /* (i)     : excitation buffer                      */
  234.   Word16 xn[],        /* (i)     : target vector                          */
  235.   Word16 h[],         /* (i) Q12 : impulse response of filters.           */
  236.   Word16 L_subfr,     /* (i)     : Length of subframe                     */
  237.   Word16 t0_min,      /* (i)     : minimum value in the searched range.   */
  238.   Word16 t0_max,      /* (i)     : maximum value in the searched range.   */
  239.   Word16 i_subfr,     /* (i)     : indicator for first subframe.          */
  240.   Word16 *pit_frac    /* (o)     : chosen fraction.                       */
  241. )
  242. {
  243.   Word16 t, t0;
  244.   Word16 Dn[L_SUBFR];
  245.   Word16 exc_tmp[L_SUBFR];
  246.   Word32 max, corr, L_temp;
  247.  /*-----------------------------------------------------------------*
  248.   * Compute correlation of target vector with impulse response.     *
  249.   *-----------------------------------------------------------------*/
  250.   Cor_h_X(h, xn, Dn);
  251.  /*-----------------------------------------------------------------*
  252.   * Find maximum integer delay.                                     *
  253.   *-----------------------------------------------------------------*/
  254.   max = MIN_32;
  255.   t0 = t0_min; /* Only to remove warning from some compilers */
  256.   for(t=t0_min; t<=t0_max; t++)
  257.   {
  258.     corr = Dot_Product(Dn, &exc[-t], L_subfr);
  259.     L_temp = L_sub(corr, max);
  260.     if(L_temp > 0) {max = corr; t0 = t;  }
  261.   }
  262.  /*-----------------------------------------------------------------*
  263.   * Test fractions.                                                 *
  264.   *-----------------------------------------------------------------*/
  265.   /* Fraction 0 */
  266.   Pred_lt_3(exc, t0, 0, L_subfr);
  267.   max = Dot_Product(Dn, exc, L_subfr);
  268.   *pit_frac = 0;
  269.   /* If first subframe and lag > 84 do not search fractional pitch */
  270.   if( (i_subfr == 0) && (sub(t0, 84) > 0) )
  271.     return t0;
  272.   Copy(exc, exc_tmp, L_subfr);
  273.   /* Fraction -1/3 */
  274.   Pred_lt_3(exc, t0, -1, L_subfr);
  275.   corr = Dot_Product(Dn, exc, L_subfr);
  276.   L_temp = L_sub(corr, max);
  277.   if(L_temp > 0) {
  278.      max = corr;
  279.      *pit_frac = -1;
  280.      Copy(exc, exc_tmp, L_subfr);
  281.   }
  282.   /* Fraction +1/3 */
  283.   Pred_lt_3(exc, t0, 1, L_subfr);
  284.   corr = Dot_Product(Dn, exc, L_subfr);
  285.   L_temp = L_sub(corr, max);
  286.   if(L_temp > 0) {
  287.      max = corr;
  288.      *pit_frac =  1;
  289.   }
  290.   else
  291.     Copy(exc_tmp, exc, L_subfr);
  292.   return t0;
  293. }
  294. /*---------------------------------------------------------------------*
  295.  * Function  G_pitch:                                                  *
  296.  *           ~~~~~~~~                                                  *
  297.  *---------------------------------------------------------------------*
  298.  * Compute correlations <xn,y1> and <y1,y1> to use in gains quantizer. *
  299.  * Also compute the gain of pitch. Result in Q14                       *
  300.  *  if (gain < 0)  gain =0                                             *
  301.  *  if (gain >1.2) gain =1.2                                           *
  302.  *---------------------------------------------------------------------*/
  303. Word16 G_pitch(      /* (o) Q14 : Gain of pitch lag saturated to 1.2       */
  304.   Word16 xn[],       /* (i)     : Pitch target.                            */
  305.   Word16 y1[],       /* (i)     : Filtered adaptive codebook.              */
  306.   Word16 g_coeff[],  /* (i)     : Correlations need for gain quantization. */
  307.   Word16 L_subfr     /* (i)     : Length of subframe.                      */
  308. )
  309. {
  310.    Word16 i;
  311.    Word16 xy, yy, exp_xy, exp_yy, gain;
  312.    Word32 s;
  313.    Word16 scaled_y1[L_SUBFR];
  314.    /* divide "y1[]" by 4 to avoid overflow */
  315.    for(i=0; i<L_subfr; i++)
  316.      scaled_y1[i] = shr(y1[i], 2);
  317.    /* Compute scalar product <y1[],y1[]> */
  318.    Overflow = 0;
  319.    s = 1;                    /* Avoid case of all zeros */
  320.    for(i=0; i<L_subfr; i++)
  321.      s = L_mac(s, y1[i], y1[i]);
  322.    if (Overflow == 0) {
  323.      exp_yy = norm_l(s);
  324.      yy     = round( L_shl(s, exp_yy) );
  325.    }
  326.    else {
  327.      s = 1;                  /* Avoid case of all zeros */
  328.      for(i=0; i<L_subfr; i++)
  329.        s = L_mac(s, scaled_y1[i], scaled_y1[i]);
  330.      exp_yy = norm_l(s);
  331.      yy     = round( L_shl(s, exp_yy) );
  332.      exp_yy = sub(exp_yy, 4);
  333.    }
  334.    /* Compute scalar product <xn[],y1[]> */
  335.    Overflow = 0;
  336.    s = 0;
  337.    for(i=0; i<L_subfr; i++)
  338.      s = L_mac(s, xn[i], y1[i]);
  339.    if (Overflow == 0) {
  340.      exp_xy = norm_l(s);
  341.      xy     = round( L_shl(s, exp_xy) );
  342.    }
  343.    else {
  344.      s = 0;
  345.      for(i=0; i<L_subfr; i++)
  346.        s = L_mac(s, xn[i], scaled_y1[i]);
  347.      exp_xy = norm_l(s);
  348.      xy     = round( L_shl(s, exp_xy) );
  349.      exp_xy = sub(exp_xy, 2);
  350.    }
  351.    g_coeff[0] = yy;
  352.    g_coeff[1] = sub(15, exp_yy);
  353.    g_coeff[2] = xy;
  354.    g_coeff[3] = sub(15, exp_xy);
  355.    /* If (xy <= 0) gain = 0 */
  356.    if (xy <= 0)
  357.    {
  358.       g_coeff[3] = -15;   /* Force exp_xy to -15 = (15-30) */
  359.       return( (Word16) 0);
  360.    }
  361.    /* compute gain = xy/yy */
  362.    xy = shr(xy, 1);             /* Be sure xy < yy */
  363.    gain = div_s( xy, yy);
  364.    i = sub(exp_xy, exp_yy);
  365.    gain = shr(gain, i);         /* saturation if > 1.99 in Q14 */
  366.    /* if(gain >1.2) gain = 1.2  in Q14 */
  367.    if( sub(gain, 19661) > 0)
  368.    {
  369.      gain = 19661;
  370.    }
  371.    return(gain);
  372. }
  373. /*----------------------------------------------------------------------*
  374.  *    Function Enc_lag3                                                 *
  375.  *             ~~~~~~~~                                                 *
  376.  *   Encoding of fractional pitch lag with 1/3 resolution.              *
  377.  *----------------------------------------------------------------------*
  378.  * The pitch range for the first subframe is divided as follows:        *
  379.  *   19 1/3  to   84 2/3   resolution 1/3                               *
  380.  *   85      to   143      resolution 1                                 *
  381.  *                                                                      *
  382.  * The period in the first subframe is encoded with 8 bits.             *
  383.  * For the range with fractions:                                        *
  384.  *   index = (T-19)*3 + frac - 1;   where T=[19..85] and frac=[-1,0,1]  *
  385.  * and for the integer only range                                       *
  386.  *   index = (T - 85) + 197;        where T=[86..143]                   *
  387.  *----------------------------------------------------------------------*
  388.  * For the second subframe a resolution of 1/3 is always used, and the  *
  389.  * search range is relative to the lag in the first subframe.           *
  390.  * If t0 is the lag in the first subframe then                          *
  391.  *  t_min=t0-5   and  t_max=t0+4   and  the range is given by           *
  392.  *       t_min - 2/3   to  t_max + 2/3                                  *
  393.  *                                                                      *
  394.  * The period in the 2nd subframe is encoded with 5 bits:               *
  395.  *   index = (T-(t_min-1))*3 + frac - 1;    where T[t_min-1 .. t_max+1] *
  396.  *----------------------------------------------------------------------*/
  397. Word16 Enc_lag3(     /* output: Return index of encoding */
  398.   Word16 T0,         /* input : Pitch delay              */
  399.   Word16 T0_frac,    /* input : Fractional pitch delay   */
  400.   Word16 *T0_min,    /* in/out: Minimum search delay     */
  401.   Word16 *T0_max,    /* in/out: Maximum search delay     */
  402.   Word16 pit_min,    /* input : Minimum pitch delay      */
  403.   Word16 pit_max,    /* input : Maximum pitch delay      */
  404.   Word16 pit_flag    /* input : Flag for 1st subframe    */
  405. )
  406. {
  407.   Word16 index, i;
  408.   if (pit_flag == 0)   /* if 1st subframe */
  409.   {
  410.     /* encode pitch delay (with fraction) */
  411.     if (sub(T0, 85) <= 0)
  412.     {
  413.       /* index = t0*3 - 58 + t0_frac   */
  414.       i = add(add(T0, T0), T0);
  415.       index = add(sub(i, 58), T0_frac);
  416.     }
  417.     else {
  418.       index = add(T0, 112);
  419.     }
  420.     /* find T0_min and T0_max for second subframe */
  421.     *T0_min = sub(T0, 5);
  422.     if (sub(*T0_min, pit_min) < 0)
  423.     {
  424.       *T0_min = pit_min;
  425.     }
  426.     *T0_max = add(*T0_min, 9);
  427.     if (sub(*T0_max, pit_max) > 0)
  428.     {
  429.       *T0_max = pit_max;
  430.       *T0_min = sub(*T0_max, 9);
  431.     }
  432.   }
  433.   else      /* if second subframe */
  434.   {
  435.     /* i = t0 - t0_min;               */
  436.     /* index = i*3 + 2 + t0_frac;     */
  437.     i = sub(T0, *T0_min);
  438.     i = add(add(i, i), i);
  439.     index = add(add(i, 2), T0_frac);
  440.   }
  441.   return index;
  442. }