opencore/codecs_v2/audio/gsm_amr/amr_nb/common/src/az_lsp.cpp

   1 /* ------------------------------------------------------------------
   2  * Copyright (C) 1998-2009 PacketVideo
   3  *
   4  * Licensed under the Apache License, Version 2.0 (the "License");
   5  * you may not use this file except in compliance with the License.
   6  * You may obtain a copy of the License at
   7  *
   8  *      http://www.apache.org/licenses/LICENSE-2.0
   9  *
  10  * Unless required by applicable law or agreed to in writing, software
  11  * distributed under the License is distributed on an "AS IS" BASIS,
  12  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either
  13  * express or implied.
  14  * See the License for the specific language governing permissions
  15  * and limitations under the License.
  16  * -------------------------------------------------------------------
  17  */
  18 /****************************************************************************************
  19 Portions of this file are derived from the following 3GPP standard:
  20
  21     3GPP TS 26.073
  22     ANSI-C code for the Adaptive Multi-Rate (AMR) speech codec
  23     Available from http://www.3gpp.org
  24
  25 (C) 2004, 3GPP Organizational Partners (ARIB, ATIS, CCSA, ETSI, TTA, TTC)
  26 Permission to distribute, modify and use this file under the standard license
  27 terms listed above has been obtained from the copyright holder.
  28 ****************************************************************************************/
  29 /*
  30
  31  Filename: az_lsp.cpp
  32  Funtions: Chebps
  33            Chebps_Wrapper
  34            Az_lsp
  35
  36 ------------------------------------------------------------------------------
  37  MODULE DESCRIPTION
  38
  39  These modules compute the LSPs from the LP coefficients.
  40
  41 ------------------------------------------------------------------------------
  42 */
  43
  44
  45 /*----------------------------------------------------------------------------
  46 ; INCLUDES
  47 ----------------------------------------------------------------------------*/
  48 #include "az_lsp.h"
  49 #include "cnst.h"
  50 #include "basic_op.h"
  51
  52 /*----------------------------------------------------------------------------
  53 ; MACROS
  54 ; Define module specific macros here
  55 ----------------------------------------------------------------------------*/
  56
  57
  58 /*----------------------------------------------------------------------------
  59 ; DEFINES
  60 ; Include all pre-processor statements here. Include conditional
  61 ; compile variables also.
  62 ----------------------------------------------------------------------------*/
  63 #define NC   M/2                  /* M = LPC order, NC = M/2 */
  64
  65 /*----------------------------------------------------------------------------
  66 ; LOCAL FUNCTION DEFINITIONS
  67 ; Function Prototype declaration
  68 ----------------------------------------------------------------------------*/
  69
  70
  71 /*----------------------------------------------------------------------------
  72 ; LOCAL VARIABLE DEFINITIONS
  73 ; Variable declaration - defined here and used outside this module
  74 ----------------------------------------------------------------------------*/
  75
  76 /*
  77 ------------------------------------------------------------------------------
  78  FUNCTION NAME: Chebps
  79 ------------------------------------------------------------------------------
  80  INPUT AND OUTPUT DEFINITIONS
  81
  82  Inputs:
  83     x = input value (Word16)
  84     f = polynomial (Word16)
  85     n = polynomial order (Word16)
  86
  87     pOverflow = pointer to overflow (Flag)
  88
  89  Outputs:
  90     pOverflow -> 1 if the operations in the function resulted in saturation.
  91
  92  Returns:
  93     cheb = Chebyshev polynomial for the input value x.(Word16)
  94
  95  Global Variables Used:
  96     None.
  97
  98  Local Variables Needed:
  99     None.
 100
 101 ------------------------------------------------------------------------------
 102  FUNCTION DESCRIPTION
 103
 104  This module evaluates the Chebyshev polynomial series.
 105     - The polynomial order is   n = m/2 = 5
 106     - The polynomial F(z) (F1(z) or F2(z)) is given by
 107         F(w) = 2 exp(-j5w) C(x)
 108         where
 109         C(x) = T_n(x) + f(1)T_n-1(x) + ... +f(n-1)T_1(x) + f(n)/2
 110         and T_m(x) = cos(mw) is the mth order Chebyshev
 111         polynomial ( x=cos(w) )
 112     - C(x) for the input x is returned.
 113
 114 ------------------------------------------------------------------------------
 115  REQUIREMENTS
 116
 117  None.
 118
 119 ------------------------------------------------------------------------------
 120  REFERENCES
 121
 122  az_lsp.c, UMTS GSM AMR speech codec, R99 - Version 3.2.0, March 2, 2001
 123
 124 ------------------------------------------------------------------------------
 125  PSEUDO-CODE
 126
 127 static Word16 Chebps (Word16 x,
 128                       Word16 f[], // (n)
 129                       Word16 n)
 130 {
 131     Word16 i, cheb;
 132     Word16 b0_h, b0_l, b1_h, b1_l, b2_h, b2_l;
 133     Word32 t0;
 134
 135 // The reference ETSI code uses a global flag for Overflow. However, in the
 136 // actual implementation a pointer to Overflow flag is passed in as a
 137 // parameter to the function. This pointer is passed into all the basic math
 138 // functions invoked
 139
 140     b2_h = 256; // b2 = 1.0
 141     b2_l = 0;
 142
 143     t0 = L_mult (x, 512);          // 2*x
 144     t0 = L_mac (t0, f[1], 8192);   // + f[1]
 145     L_Extract (t0, &b1_h, &b1_l);  // b1 = 2*x + f[1]
 146
 147     for (i = 2; i < n; i++)
 148     {
 149         t0 = Mpy_32_16 (b1_h, b1_l, x);         // t0 = 2.0*x*b1
 150         t0 = L_shl (t0, 1);
 151         t0 = L_mac (t0, b2_h, (Word16) 0x8000); // t0 = 2.0*x*b1 - b2
 152         t0 = L_msu (t0, b2_l, 1);
 153         t0 = L_mac (t0, f[i], 8192);            // t0 = 2.0*x*b1 - b2 + f[i]
 154
 155         L_Extract (t0, &b0_h, &b0_l);           // b0 = 2.0*x*b1 - b2 + f[i]
 156
 157         b2_l = b1_l; // b2 = b1;
 158         b2_h = b1_h;
 159         b1_l = b0_l; // b1 = b0;
 160         b1_h = b0_h;
 161     }
 162
 163     t0 = Mpy_32_16 (b1_h, b1_l, x);             // t0 = x*b1;
 164     t0 = L_mac (t0, b2_h, (Word16) 0x8000);     // t0 = x*b1 - b2
 165     t0 = L_msu (t0, b2_l, 1);
 166     t0 = L_mac (t0, f[i], 4096);                // t0 = x*b1 - b2 + f[i]/2
 167
 168     t0 = L_shl (t0, 6);
 169
 170     cheb = extract_h (t0);
 171
 172     return (cheb);
 173 }
 174
 175 ------------------------------------------------------------------------------
 176  CAUTION [optional]
 177  [State any special notes, constraints or cautions for users of this function]
 178
 179 ------------------------------------------------------------------------------
 180 */
 181
 182 static Word16 Chebps(Word16 x,
 183                      Word16 f[], /* (n) */
 184                      Word16 n,
 185                      Flag *pOverflow)
 186 {
 187     Word16 i;
 188     Word16 cheb;
 189     Word16 b1_h;
 190     Word16 b1_l;
 191     Word32 t0;
 192     Word32 L_temp;
 193     Word16 *p_f = &f[1];
 194
 195     OSCL_UNUSED_ARG(pOverflow);
 196
 197     /* L_temp = 1.0 */
 198
 199     L_temp = 0x01000000L;
 200
 201     t0 = ((Word32) x << 10) + ((Word32) * (p_f++) << 14);
 202
 203     /* b1 = t0 = 2*x + f[1]  */
 204
 205     b1_h = (Word16)(t0 >> 16);
 206     b1_l = (Word16)((t0 >> 1) - (b1_h << 15));
 207
 208
 209     for (i = 2; i < n; i++)
 210     {
 211         /* t0 = 2.0*x*b1    */
 212         t0  = ((Word32) b1_h * x);
 213         t0 += ((Word32) b1_l * x) >> 15;
 214         t0 <<= 2;
 215
 216         /* t0 = 2.0*x*b1 - b2   */
 217         t0 -= L_temp;
 218
 219         /* t0 = 2.0*x*b1 - b2 + f[i] */
 220         t0 += (Word32) * (p_f++) << 14;
 221
 222         L_temp = ((Word32) b1_h << 16) + ((Word32) b1_l << 1);
 223
 224         /* b0 = 2.0*x*b1 - b2 + f[i]*/
 225         b1_h = (Word16)(t0 >> 16);
 226         b1_l = (Word16)((t0 >> 1) - (b1_h << 15));
 227
 228     }
 229
 230     /* t0 = x*b1; */
 231     t0  = ((Word32) b1_h * x);
 232     t0 += ((Word32) b1_l * x) >> 15;
 233     t0 <<= 1;
 234
 235
 236     /* t0 = x*b1 - b2   */
 237     t0 -= L_temp;
 238
 239     /* t0 = x*b1 - b2 + f[i]/2 */
 240     t0 += (Word32) * (p_f) << 13;
 241
 242
 243     if ((UWord32)(t0 + 33554432) < 67108863)
 244     {
 245         cheb = (Word16)(t0 >> 10);
 246     }
 247     else
 248     {
 249         if (t0 > (Word32) 0x01ffffffL)
 250         {
 251             cheb = MAX_16;
 252
 253         }
 254         else
 255         {
 256             cheb = MIN_16;
 257         }
 258     }
 259
 260     return (cheb);
 261 }
 262
 263
 264 /*
 265 ------------------------------------------------------------------------------
 266  FUNCTION NAME: Az_lsp
 267 ------------------------------------------------------------------------------
 268  INPUT AND OUTPUT DEFINITIONS FOR Az_lsp
 269
 270  Inputs:
 271     a = predictor coefficients (Word16)
 272     lsp = line spectral pairs (Word16)
 273     old_lsp = old line spectral pairs (Word16)
 274
 275     pOverflow = pointer to overflow (Flag)
 276
 277  Outputs:
 278     pOverflow -> 1 if the operations in the function resulted in saturation.
 279
 280  Returns:
 281     None.
 282
 283  Global Variables Used:
 284     None.
 285
 286  Local Variables Needed:
 287     None.
 288
 289 ------------------------------------------------------------------------------
 290  FUNCTION DESCRIPTION
 291
 292  This function computes the LSPs from the LP coefficients.
 293
 294  The sum and difference filters are computed and divided by 1+z^{-1} and
 295  1-z^{-1}, respectively.
 296
 297      f1[i] = a[i] + a[11-i] - f1[i-1] ;   i=1,...,5
 298      f2[i] = a[i] - a[11-i] + f2[i-1] ;   i=1,...,5
 299
 300  The roots of F1(z) and F2(z) are found using Chebyshev polynomial evaluation.
 301  The polynomials are evaluated at 60 points regularly spaced in the
 302  frequency domain. The sign change interval is subdivided 4 times to better
 303  track the root. The LSPs are found in the cosine domain [1,-1].
 304
 305  If less than 10 roots are found, the LSPs from the past frame are used.
 306
 307 ------------------------------------------------------------------------------
 308  REQUIREMENTS
 309
 310  None.
 311
 312 ------------------------------------------------------------------------------
 313  REFERENCES
 314
 315  az_lsp.c, UMTS GSM AMR speech codec, R99 - Version 3.2.0, March 2, 2001
 316
 317 ------------------------------------------------------------------------------
 318  PSEUDO-CODE
 319
 320 void Az_lsp (
 321     Word16 a[],         // (i)  : predictor coefficients (MP1)
 322     Word16 lsp[],       // (o)  : line spectral pairs (M)
 323     Word16 old_lsp[]    // (i)  : old lsp[] (in case not found 10 roots) (M)
 324 )
 325 {
 326     Word16 i, j, nf, ip;
 327     Word16 xlow, ylow, xhigh, yhigh, xmid, ymid, xint;
 328     Word16 x, y, sign, exp;
 329     Word16 *coef;
 330     Word16 f1[M / 2 + 1], f2[M / 2 + 1];
 331     Word32 t0;
 332
 333      *-------------------------------------------------------------*
 334      *  find the sum and diff. pol. F1(z) and F2(z)                *
 335      *    F1(z) <--- F1(z)/(1+z**-1) & F2(z) <--- F2(z)/(1-z**-1)  *
 336      *                                                             *
 337      * f1[0] = 1.0;                                                *
 338      * f2[0] = 1.0;                                                *
 339      *                                                             *
 340      * for (i = 0; i< NC; i++)                                     *
 341      * {                                                           *
 342      *   f1[i+1] = a[i+1] + a[M-i] - f1[i] ;                       *
 343      *   f2[i+1] = a[i+1] - a[M-i] + f2[i] ;                       *
 344      * }                                                           *
 345      *-------------------------------------------------------------*
 346
 347     f1[0] = 1024; // f1[0] = 1.0
 348     f2[0] = 1024; // f2[0] = 1.0
 349
 350 // The reference ETSI code uses a global flag for Overflow. However, in the
 351 // actual implementation a pointer to Overflow flag is passed in as a
 352 // parameter to the function. This pointer is passed into all the basic math
 353 // functions invoked
 354
 355     for (i = 0; i < NC; i++)
 356     {
 357         t0 = L_mult (a[i + 1], 8192);   // x = (a[i+1] + a[M-i]) >> 2
 358         t0 = L_mac (t0, a[M - i], 8192);
 359         x = extract_h (t0);
 360         // f1[i+1] = a[i+1] + a[M-i] - f1[i]
 361         f1[i + 1] = sub (x, f1[i]);
 362
 363         t0 = L_mult (a[i + 1], 8192);   // x = (a[i+1] - a[M-i]) >> 2
 364         t0 = L_msu (t0, a[M - i], 8192);
 365         x = extract_h (t0);
 366         // f2[i+1] = a[i+1] - a[M-i] + f2[i]
 367         f2[i + 1] = add (x, f2[i]);
 368     }
 369
 370      *-------------------------------------------------------------*
 371      * find the LSPs using the Chebychev pol. evaluation           *
 372      *-------------------------------------------------------------*
 373
 374     nf = 0; // number of found frequencies
 375     ip = 0; // indicator for f1 or f2
 376
 377     coef = f1;
 378
 379     xlow = grid[0];
 380     ylow = Chebps (xlow, coef, NC);
 381
 382     j = 0;
 383     // while ( (nf < M) && (j < grid_points) )
 384     while ((sub (nf, M) < 0) && (sub (j, grid_points) < 0))
 385     {
 386         j++;
 387         xhigh = xlow;
 388         yhigh = ylow;
 389         xlow = grid[j];
 390         ylow = Chebps (xlow, coef, NC);
 391
 392         if (L_mult (ylow, yhigh) <= (Word32) 0L)
 393         {
 394
 395             // divide 4 times the interval
 396
 397             for (i = 0; i < 4; i++)
 398             {
 399                 // xmid = (xlow + xhigh)/2
 400                 xmid = add (shr (xlow, 1), shr (xhigh, 1));
 401                 ymid = Chebps (xmid, coef, NC);
 402
 403                 if (L_mult (ylow, ymid) <= (Word32) 0L)
 404                 {
 405                     yhigh = ymid;
 406                     xhigh = xmid;
 407                 }
 408                 else
 409                 {
 410                     ylow = ymid;
 411                     xlow = xmid;
 412                 }
 413             }
 414
 415              *-------------------------------------------------------------*
 416              * Linear interpolation                                        *
 417              *    xint = xlow - ylow*(xhigh-xlow)/(yhigh-ylow);            *
 418              *-------------------------------------------------------------*
 419
 420             x = sub (xhigh, xlow);
 421             y = sub (yhigh, ylow);
 422
 423             if (y == 0)
 424             {
 425                 xint = xlow;
 426             }
 427             else
 428             {
 429                 sign = y;
 430                 y = abs_s (y);
 431                 exp = norm_s (y);
 432                 y = shl (y, exp);
 433                 y = div_s ((Word16) 16383, y);
 434                 t0 = L_mult (x, y);
 435                 t0 = L_shr (t0, sub (20, exp));
 436                 y = extract_l (t0);     // y= (xhigh-xlow)/(yhigh-ylow)
 437
 438                 if (sign < 0)
 439                     y = negate (y);
 440
 441                 t0 = L_mult (ylow, y);
 442                 t0 = L_shr (t0, 11);
 443                 xint = sub (xlow, extract_l (t0)); // xint = xlow - ylow*y
 444             }
 445
 446             lsp[nf] = xint;
 447             xlow = xint;
 448             nf++;
 449
 450             if (ip == 0)
 451             {
 452                 ip = 1;
 453                 coef = f2;
 454             }
 455             else
 456             {
 457                 ip = 0;
 458                 coef = f1;
 459             }
 460             ylow = Chebps (xlow, coef, NC);
 461
 462         }
 463     }
 464
 465     // Check if M roots found
 466
 467     if (sub (nf, M) < 0)
 468     {
 469         for (i = 0; i < M; i++)
 470         {
 471             lsp[i] = old_lsp[i];
 472         }
 473
 474     }
 475     return;
 476 }
 477
 478 ------------------------------------------------------------------------------
 479  CAUTION [optional]
 480  [State any special notes, constraints or cautions for users of this function]
 481
 482 ------------------------------------------------------------------------------
 483 */
 484
 485 void Az_lsp(
 486     Word16 a[],         /* (i)  : predictor coefficients (MP1)               */
 487     Word16 lsp[],       /* (o)  : line spectral pairs (M)                    */
 488     Word16 old_lsp[],   /* (i)  : old lsp[] (in case not found 10 roots) (M) */
 489     Flag   *pOverflow   /* (i/o): overflow flag                              */
 490 )
 491 {
 492     register Word16 i;
 493     register Word16 j;
 494     register Word16 nf;
 495     register Word16 ip;
 496     Word16 xlow;
 497     Word16 ylow;
 498     Word16 xhigh;
 499     Word16 yhigh;
 500     Word16 xmid;
 501     Word16 ymid;
 502     Word16 xint;
 503     Word16 x;
 504     Word16 y;
 505     Word16 sign;
 506     Word16 exp;
 507     Word16 *coef;
 508     Word16 f1[NC + 1];
 509     Word16 f2[NC + 1];
 510     Word32 L_temp1;
 511     Word32 L_temp2;
 512     Word16 *p_f1 = f1;
 513     Word16 *p_f2 = f2;
 514
 515     /*-------------------------------------------------------------*
 516      *  find the sum and diff. pol. F1(z) and F2(z)                *
 517      *    F1(z) <--- F1(z)/(1+z**-1) & F2(z) <--- F2(z)/(1-z**-1)  *
 518      *                                                             *
 519      * f1[0] = 1.0;                                                *
 520      * f2[0] = 1.0;                                                *
 521      *                                                             *
 522      * for (i = 0; i< NC; i++)                                     *
 523      * {                                                           *
 524      *   f1[i+1] = a[i+1] + a[M-i] - f1[i] ;                       *
 525      *   f2[i+1] = a[i+1] - a[M-i] + f2[i] ;                       *
 526      * }                                                           *
 527      *-------------------------------------------------------------*/
 528
 529     *p_f1 = 1024;                       /* f1[0] = 1.0 */
 530     *p_f2 = 1024;                       /* f2[0] = 1.0 */
 531
 532     for (i = 0; i < NC; i++)
 533     {
 534         L_temp1 = (Word32) * (a + i + 1);
 535         L_temp2 = (Word32) * (a + M - i);
 536         /* x = (a[i+1] + a[M-i]) >> 2  */
 537         x = (Word16)((L_temp1 + L_temp2) >> 2);
 538         /* y = (a[i+1] - a[M-i]) >> 2 */
 539         y = (Word16)((L_temp1 - L_temp2) >> 2);
 540         /* f1[i+1] = a[i+1] + a[M-i] - f1[i] */
 541         x -= *(p_f1++);
 542         *(p_f1) = x;
 543         /* f2[i+1] = a[i+1] - a[M-i] + f2[i] */
 544         y += *(p_f2++);
 545         *(p_f2) = y;
 546     }
 547
 548     /*-------------------------------------------------------------*
 549      * find the LSPs using the Chebychev pol. evaluation           *
 550      *-------------------------------------------------------------*/
 551
 552     nf = 0;                         /* number of found frequencies */
 553     ip = 0;                         /* indicator for f1 or f2      */
 554
 555     coef = f1;
 556
 557     xlow = *(grid);
 558     ylow = Chebps(xlow, coef, NC, pOverflow);
 559
 560     j = 0;
 561
 562     while ((nf < M) && (j < grid_points))
 563     {
 564         j++;
 565         xhigh = xlow;
 566         yhigh = ylow;
 567         xlow = *(grid + j);
 568         ylow = Chebps(xlow, coef, NC, pOverflow);
 569
 570         if (((Word32)ylow*yhigh) <= 0)
 571         {
 572             /* divide 4 times the interval */
 573             for (i = 4; i != 0; i--)
 574             {
 575                 /* xmid = (xlow + xhigh)/2 */
 576                 x = xlow >> 1;
 577                 y = xhigh >> 1;
 578                 xmid = x + y;
 579
 580                 ymid = Chebps(xmid, coef, NC, pOverflow);
 581
 582                 if (((Word32)ylow*ymid) <= 0)
 583                 {
 584                     yhigh = ymid;
 585                     xhigh = xmid;
 586                 }
 587                 else
 588                 {
 589                     ylow = ymid;
 590                     xlow = xmid;
 591                 }
 592             }
 593
 594             /*-------------------------------------------------------------*
 595              * Linear interpolation                                        *
 596              *    xint = xlow - ylow*(xhigh-xlow)/(yhigh-ylow);            *
 597              *-------------------------------------------------------------*/
 598
 599             x = xhigh - xlow;
 600             y = yhigh - ylow;
 601
 602             if (y == 0)
 603             {
 604                 xint = xlow;
 605             }
 606             else
 607             {
 608                 sign = y;
 609                 y = abs_s(y);
 610                 exp = norm_s(y);
 611                 y <<= exp;
 612                 y = div_s((Word16) 16383, y);
 613
 614                 y = ((Word32)x * y) >> (19 - exp);
 615
 616                 if (sign < 0)
 617                 {
 618                     y = -y;
 619                 }
 620
 621                 /* xint = xlow - ylow*y */
 622                 xint = xlow - (((Word32) ylow * y) >> 10);
 623             }
 624
 625             *(lsp + nf) = xint;
 626             xlow = xint;
 627             nf++;
 628
 629             if (ip == 0)
 630             {
 631                 ip = 1;
 632                 coef = f2;
 633             }
 634             else
 635             {
 636                 ip = 0;
 637                 coef = f1;
 638             }
 639
 640             ylow = Chebps(xlow, coef, NC, pOverflow);
 641
 642         }
 643     }
 644
 645     /* Check if M roots found */
 646
 647     if (nf < M)
 648     {
 649         for (i = NC; i != 0 ; i--)
 650         {
 651             *lsp++ = *old_lsp++;
 652             *lsp++ = *old_lsp++;
 653         }
 654     }
 655
 656 }
 657
 658 Word16 Chebps_Wrapper(Word16 x,
 659                       Word16 f[], /* (n) */
 660                       Word16 n,
 661                       Flag *pOverflow)
 662 {
 663     return Chebps(x, f, n, pOverflow);
 664 }
 665