az_lsp.cpp revision b875f69d8e867cb64bd101e66d85a880537c2b72
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/**************************************************************************************** 19Portions 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) 26Permission to distribute, modify and use this file under the standard license 27terms listed above has been obtained from the copyright holder. 28****************************************************************************************/ 29/* 30 31 Pathname: ./audio/gsm-amr/c/src/az_lsp.c 32 Funtions: Chebps 33 Chebps_Wrapper 34 Az_lsp 35 36------------------------------------------------------------------------------ 37 REVISION HISTORY 38 39 Description: Finished first pass of optimization. 40 41 Description: Made changes based on review comments. 42 43 Description: Made input to Chebps_Wrapper consistent with that of Chebps. 44 45 Description: Replaced current Pseudo-code with the UMTS code version 3.2.0. 46 Updated coding template. 47 48 Description: Replaced basic_op.h and oper_32b.h with the header files of the 49 math functions used by the file. 50 51 Description: Made the following changes per comments from Phase 2/3 review: 52 1. Used "-" operator instead of calling sub function in the 53 az_lsp() code. 54 2. Copied detailed function description of az_lsp from the 55 header file. 56 3. Modified local variable definition to one per line. 57 4. Used NC in the definition of f1 and f2 arrays. 58 5. Added curly brackets in the IF statement. 59 60 Description: Changed function interface to pass in a pointer to overflow 61 flag into the function instead of using a global flag. Removed 62 inclusion of unneeded header files. 63 64 Description: For Chebps() and Az_lsp() 65 1. Eliminated unused include files. 66 2. Replaced array addressing by pointers 67 3. Eliminated math operations that unnecessary checked for 68 saturation. 69 4. Eliminated not needed variables 70 5. Eliminated if-else checks for saturation 71 6. Deleted unused function cheps_wraper 72 73 Description: Added casting to eliminate warnings 74 75 76 Who: Date: 77 Description: 78 79------------------------------------------------------------------------------ 80 MODULE DESCRIPTION 81 82 These modules compute the LSPs from the LP coefficients. 83 84------------------------------------------------------------------------------ 85*/ 86 87 88/*---------------------------------------------------------------------------- 89; INCLUDES 90----------------------------------------------------------------------------*/ 91#include "az_lsp.h" 92#include "cnst.h" 93#include "basic_op.h" 94 95/*---------------------------------------------------------------------------- 96; MACROS 97; Define module specific macros here 98----------------------------------------------------------------------------*/ 99 100 101/*---------------------------------------------------------------------------- 102; DEFINES 103; Include all pre-processor statements here. Include conditional 104; compile variables also. 105----------------------------------------------------------------------------*/ 106#define NC M/2 /* M = LPC order, NC = M/2 */ 107 108/*---------------------------------------------------------------------------- 109; LOCAL FUNCTION DEFINITIONS 110; Function Prototype declaration 111----------------------------------------------------------------------------*/ 112 113 114/*---------------------------------------------------------------------------- 115; LOCAL VARIABLE DEFINITIONS 116; Variable declaration - defined here and used outside this module 117----------------------------------------------------------------------------*/ 118 119/* 120------------------------------------------------------------------------------ 121 FUNCTION NAME: Chebps 122------------------------------------------------------------------------------ 123 INPUT AND OUTPUT DEFINITIONS 124 125 Inputs: 126 x = input value (Word16) 127 f = polynomial (Word16) 128 n = polynomial order (Word16) 129 130 pOverflow = pointer to overflow (Flag) 131 132 Outputs: 133 pOverflow -> 1 if the operations in the function resulted in saturation. 134 135 Returns: 136 cheb = Chebyshev polynomial for the input value x.(Word16) 137 138 Global Variables Used: 139 None. 140 141 Local Variables Needed: 142 None. 143 144------------------------------------------------------------------------------ 145 FUNCTION DESCRIPTION 146 147 This module evaluates the Chebyshev polynomial series. 148 - The polynomial order is n = m/2 = 5 149 - The polynomial F(z) (F1(z) or F2(z)) is given by 150 F(w) = 2 exp(-j5w) C(x) 151 where 152 C(x) = T_n(x) + f(1)T_n-1(x) + ... +f(n-1)T_1(x) + f(n)/2 153 and T_m(x) = cos(mw) is the mth order Chebyshev 154 polynomial ( x=cos(w) ) 155 - C(x) for the input x is returned. 156 157------------------------------------------------------------------------------ 158 REQUIREMENTS 159 160 None. 161 162------------------------------------------------------------------------------ 163 REFERENCES 164 165 az_lsp.c, UMTS GSM AMR speech codec, R99 - Version 3.2.0, March 2, 2001 166 167------------------------------------------------------------------------------ 168 PSEUDO-CODE 169 170static Word16 Chebps (Word16 x, 171 Word16 f[], // (n) 172 Word16 n) 173{ 174 Word16 i, cheb; 175 Word16 b0_h, b0_l, b1_h, b1_l, b2_h, b2_l; 176 Word32 t0; 177 178// The reference ETSI code uses a global flag for Overflow. However, in the 179// actual implementation a pointer to Overflow flag is passed in as a 180// parameter to the function. This pointer is passed into all the basic math 181// functions invoked 182 183 b2_h = 256; // b2 = 1.0 184 b2_l = 0; 185 186 t0 = L_mult (x, 512); // 2*x 187 t0 = L_mac (t0, f[1], 8192); // + f[1] 188 L_Extract (t0, &b1_h, &b1_l); // b1 = 2*x + f[1] 189 190 for (i = 2; i < n; i++) 191 { 192 t0 = Mpy_32_16 (b1_h, b1_l, x); // t0 = 2.0*x*b1 193 t0 = L_shl (t0, 1); 194 t0 = L_mac (t0, b2_h, (Word16) 0x8000); // t0 = 2.0*x*b1 - b2 195 t0 = L_msu (t0, b2_l, 1); 196 t0 = L_mac (t0, f[i], 8192); // t0 = 2.0*x*b1 - b2 + f[i] 197 198 L_Extract (t0, &b0_h, &b0_l); // b0 = 2.0*x*b1 - b2 + f[i] 199 200 b2_l = b1_l; // b2 = b1; 201 b2_h = b1_h; 202 b1_l = b0_l; // b1 = b0; 203 b1_h = b0_h; 204 } 205 206 t0 = Mpy_32_16 (b1_h, b1_l, x); // t0 = x*b1; 207 t0 = L_mac (t0, b2_h, (Word16) 0x8000); // t0 = x*b1 - b2 208 t0 = L_msu (t0, b2_l, 1); 209 t0 = L_mac (t0, f[i], 4096); // t0 = x*b1 - b2 + f[i]/2 210 211 t0 = L_shl (t0, 6); 212 213 cheb = extract_h (t0); 214 215 return (cheb); 216} 217 218------------------------------------------------------------------------------ 219 RESOURCES USED [optional] 220 221 When the code is written for a specific target processor the 222 the resources used should be documented below. 223 224 HEAP MEMORY USED: x bytes 225 226 STACK MEMORY USED: x bytes 227 228 CLOCK CYCLES: (cycle count equation for this function) + (variable 229 used to represent cycle count for each subroutine 230 called) 231 where: (cycle count variable) = cycle count for [subroutine 232 name] 233 234------------------------------------------------------------------------------ 235 CAUTION [optional] 236 [State any special notes, constraints or cautions for users of this function] 237 238------------------------------------------------------------------------------ 239*/ 240#ifdef __clang__ 241__attribute__((no_sanitize("integer"))) 242#endif 243static Word16 Chebps(Word16 x, 244 Word16 f[], /* (n) */ 245 Word16 n, 246 Flag *pOverflow) 247{ 248 Word16 i; 249 Word16 cheb; 250 Word16 b1_h; 251 Word16 b1_l; 252 Word32 t0; 253 Word32 L_temp; 254 Word16 *p_f = &f[1]; 255 256 OSCL_UNUSED_ARG(pOverflow); 257 258 /* L_temp = 1.0 */ 259 260 L_temp = 0x01000000L; 261 262 t0 = ((Word32) x << 10) + ((Word32) * (p_f++) << 14); 263 264 /* b1 = t0 = 2*x + f[1] */ 265 266 b1_h = (Word16)(t0 >> 16); 267 b1_l = (Word16)((t0 >> 1) - (b1_h << 15)); 268 269 270 for (i = 2; i < n; i++) 271 { 272 /* t0 = 2.0*x*b1 */ 273 t0 = ((Word32) b1_h * x); 274 t0 += ((Word32) b1_l * x) >> 15; 275 t0 <<= 2; 276 277 /* t0 = 2.0*x*b1 - b2 */ 278 t0 -= L_temp; 279 280 /* t0 = 2.0*x*b1 - b2 + f[i] */ 281 t0 += (Word32) * (p_f++) << 14; 282 283 L_temp = ((Word32) b1_h << 16) + ((Word32) b1_l << 1); 284 285 /* b0 = 2.0*x*b1 - b2 + f[i]*/ 286 b1_h = (Word16)(t0 >> 16); 287 b1_l = (Word16)((t0 >> 1) - (b1_h << 15)); 288 289 } 290 291 /* t0 = x*b1; */ 292 t0 = ((Word32) b1_h * x); 293 t0 += ((Word32) b1_l * x) >> 15; 294 t0 <<= 1; 295 296 297 /* t0 = x*b1 - b2 */ 298 t0 -= L_temp; 299 300 /* t0 = x*b1 - b2 + f[i]/2 */ 301 t0 += (Word32) * (p_f) << 13; 302 303 304 if ((UWord32)(t0 - 0xfe000000L) < (UWord32)0x03ffffffL) 305 { 306 cheb = (Word16)(t0 >> 10); 307 } 308 else 309 { 310 if (t0 > (Word32) 0x01ffffffL) 311 { 312 cheb = MAX_16; 313 314 } 315 else 316 { 317 cheb = MIN_16; 318 } 319 } 320 321 return (cheb); 322} 323 324 325/* 326------------------------------------------------------------------------------ 327 FUNCTION NAME: Az_lsp 328------------------------------------------------------------------------------ 329 INPUT AND OUTPUT DEFINITIONS FOR Az_lsp 330 331 Inputs: 332 a = predictor coefficients (Word16) 333 lsp = line spectral pairs (Word16) 334 old_lsp = old line spectral pairs (Word16) 335 336 pOverflow = pointer to overflow (Flag) 337 338 Outputs: 339 pOverflow -> 1 if the operations in the function resulted in saturation. 340 341 Returns: 342 None. 343 344 Global Variables Used: 345 None. 346 347 Local Variables Needed: 348 None. 349 350------------------------------------------------------------------------------ 351 FUNCTION DESCRIPTION 352 353 This function computes the LSPs from the LP coefficients. 354 355 The sum and difference filters are computed and divided by 1+z^{-1} and 356 1-z^{-1}, respectively. 357 358 f1[i] = a[i] + a[11-i] - f1[i-1] ; i=1,...,5 359 f2[i] = a[i] - a[11-i] + f2[i-1] ; i=1,...,5 360 361 The roots of F1(z) and F2(z) are found using Chebyshev polynomial evaluation. 362 The polynomials are evaluated at 60 points regularly spaced in the 363 frequency domain. The sign change interval is subdivided 4 times to better 364 track the root. The LSPs are found in the cosine domain [1,-1]. 365 366 If less than 10 roots are found, the LSPs from the past frame are used. 367 368------------------------------------------------------------------------------ 369 REQUIREMENTS 370 371 None. 372 373------------------------------------------------------------------------------ 374 REFERENCES 375 376 az_lsp.c, UMTS GSM AMR speech codec, R99 - Version 3.2.0, March 2, 2001 377 378------------------------------------------------------------------------------ 379 PSEUDO-CODE 380 381void Az_lsp ( 382 Word16 a[], // (i) : predictor coefficients (MP1) 383 Word16 lsp[], // (o) : line spectral pairs (M) 384 Word16 old_lsp[] // (i) : old lsp[] (in case not found 10 roots) (M) 385) 386{ 387 Word16 i, j, nf, ip; 388 Word16 xlow, ylow, xhigh, yhigh, xmid, ymid, xint; 389 Word16 x, y, sign, exp; 390 Word16 *coef; 391 Word16 f1[M / 2 + 1], f2[M / 2 + 1]; 392 Word32 t0; 393 394 *-------------------------------------------------------------* 395 * find the sum and diff. pol. F1(z) and F2(z) * 396 * F1(z) <--- F1(z)/(1+z**-1) & F2(z) <--- F2(z)/(1-z**-1) * 397 * * 398 * f1[0] = 1.0; * 399 * f2[0] = 1.0; * 400 * * 401 * for (i = 0; i< NC; i++) * 402 * { * 403 * f1[i+1] = a[i+1] + a[M-i] - f1[i] ; * 404 * f2[i+1] = a[i+1] - a[M-i] + f2[i] ; * 405 * } * 406 *-------------------------------------------------------------* 407 408 f1[0] = 1024; // f1[0] = 1.0 409 f2[0] = 1024; // f2[0] = 1.0 410 411// The reference ETSI code uses a global flag for Overflow. However, in the 412// actual implementation a pointer to Overflow flag is passed in as a 413// parameter to the function. This pointer is passed into all the basic math 414// functions invoked 415 416 for (i = 0; i < NC; i++) 417 { 418 t0 = L_mult (a[i + 1], 8192); // x = (a[i+1] + a[M-i]) >> 2 419 t0 = L_mac (t0, a[M - i], 8192); 420 x = extract_h (t0); 421 // f1[i+1] = a[i+1] + a[M-i] - f1[i] 422 f1[i + 1] = sub (x, f1[i]); 423 424 t0 = L_mult (a[i + 1], 8192); // x = (a[i+1] - a[M-i]) >> 2 425 t0 = L_msu (t0, a[M - i], 8192); 426 x = extract_h (t0); 427 // f2[i+1] = a[i+1] - a[M-i] + f2[i] 428 f2[i + 1] = add (x, f2[i]); 429 } 430 431 *-------------------------------------------------------------* 432 * find the LSPs using the Chebychev pol. evaluation * 433 *-------------------------------------------------------------* 434 435 nf = 0; // number of found frequencies 436 ip = 0; // indicator for f1 or f2 437 438 coef = f1; 439 440 xlow = grid[0]; 441 ylow = Chebps (xlow, coef, NC); 442 443 j = 0; 444 // while ( (nf < M) && (j < grid_points) ) 445 while ((sub (nf, M) < 0) && (sub (j, grid_points) < 0)) 446 { 447 j++; 448 xhigh = xlow; 449 yhigh = ylow; 450 xlow = grid[j]; 451 ylow = Chebps (xlow, coef, NC); 452 453 if (L_mult (ylow, yhigh) <= (Word32) 0L) 454 { 455 456 // divide 4 times the interval 457 458 for (i = 0; i < 4; i++) 459 { 460 // xmid = (xlow + xhigh)/2 461 xmid = add (shr (xlow, 1), shr (xhigh, 1)); 462 ymid = Chebps (xmid, coef, NC); 463 464 if (L_mult (ylow, ymid) <= (Word32) 0L) 465 { 466 yhigh = ymid; 467 xhigh = xmid; 468 } 469 else 470 { 471 ylow = ymid; 472 xlow = xmid; 473 } 474 } 475 476 *-------------------------------------------------------------* 477 * Linear interpolation * 478 * xint = xlow - ylow*(xhigh-xlow)/(yhigh-ylow); * 479 *-------------------------------------------------------------* 480 481 x = sub (xhigh, xlow); 482 y = sub (yhigh, ylow); 483 484 if (y == 0) 485 { 486 xint = xlow; 487 } 488 else 489 { 490 sign = y; 491 y = abs_s (y); 492 exp = norm_s (y); 493 y = shl (y, exp); 494 y = div_s ((Word16) 16383, y); 495 t0 = L_mult (x, y); 496 t0 = L_shr (t0, sub (20, exp)); 497 y = extract_l (t0); // y= (xhigh-xlow)/(yhigh-ylow) 498 499 if (sign < 0) 500 y = negate (y); 501 502 t0 = L_mult (ylow, y); 503 t0 = L_shr (t0, 11); 504 xint = sub (xlow, extract_l (t0)); // xint = xlow - ylow*y 505 } 506 507 lsp[nf] = xint; 508 xlow = xint; 509 nf++; 510 511 if (ip == 0) 512 { 513 ip = 1; 514 coef = f2; 515 } 516 else 517 { 518 ip = 0; 519 coef = f1; 520 } 521 ylow = Chebps (xlow, coef, NC); 522 523 } 524 } 525 526 // Check if M roots found 527 528 if (sub (nf, M) < 0) 529 { 530 for (i = 0; i < M; i++) 531 { 532 lsp[i] = old_lsp[i]; 533 } 534 535 } 536 return; 537} 538 539------------------------------------------------------------------------------ 540 RESOURCES USED [optional] 541 542 When the code is written for a specific target processor the 543 the resources used should be documented below. 544 545 HEAP MEMORY USED: x bytes 546 547 STACK MEMORY USED: x bytes 548 549 CLOCK CYCLES: (cycle count equation for this function) + (variable 550 used to represent cycle count for each subroutine 551 called) 552 where: (cycle count variable) = cycle count for [subroutine 553 name] 554 555------------------------------------------------------------------------------ 556 CAUTION [optional] 557 [State any special notes, constraints or cautions for users of this function] 558 559------------------------------------------------------------------------------ 560*/ 561 562void Az_lsp( 563 Word16 a[], /* (i) : predictor coefficients (MP1) */ 564 Word16 lsp[], /* (o) : line spectral pairs (M) */ 565 Word16 old_lsp[], /* (i) : old lsp[] (in case not found 10 roots) (M) */ 566 Flag *pOverflow /* (i/o): overflow flag */ 567) 568{ 569 Word16 i; 570 Word16 j; 571 Word16 nf; 572 Word16 ip; 573 Word16 xlow; 574 Word16 ylow; 575 Word16 xhigh; 576 Word16 yhigh; 577 Word16 xmid; 578 Word16 ymid; 579 Word16 xint; 580 Word16 x; 581 Word16 y; 582 Word16 sign; 583 Word16 exp; 584 Word16 *coef; 585 Word16 f1[NC + 1]; 586 Word16 f2[NC + 1]; 587 Word32 L_temp1; 588 Word32 L_temp2; 589 Word16 *p_f1 = f1; 590 Word16 *p_f2 = f2; 591 592 /*-------------------------------------------------------------* 593 * find the sum and diff. pol. F1(z) and F2(z) * 594 * F1(z) <--- F1(z)/(1+z**-1) & F2(z) <--- F2(z)/(1-z**-1) * 595 * * 596 * f1[0] = 1.0; * 597 * f2[0] = 1.0; * 598 * * 599 * for (i = 0; i< NC; i++) * 600 * { * 601 * f1[i+1] = a[i+1] + a[M-i] - f1[i] ; * 602 * f2[i+1] = a[i+1] - a[M-i] + f2[i] ; * 603 * } * 604 *-------------------------------------------------------------*/ 605 606 *p_f1 = 1024; /* f1[0] = 1.0 */ 607 *p_f2 = 1024; /* f2[0] = 1.0 */ 608 609 for (i = 0; i < NC; i++) 610 { 611 L_temp1 = (Word32) * (a + i + 1); 612 L_temp2 = (Word32) * (a + M - i); 613 /* x = (a[i+1] + a[M-i]) >> 2 */ 614 x = (Word16)((L_temp1 + L_temp2) >> 2); 615 /* y = (a[i+1] - a[M-i]) >> 2 */ 616 y = (Word16)((L_temp1 - L_temp2) >> 2); 617 /* f1[i+1] = a[i+1] + a[M-i] - f1[i] */ 618 x -= *(p_f1++); 619 *(p_f1) = x; 620 /* f2[i+1] = a[i+1] - a[M-i] + f2[i] */ 621 y += *(p_f2++); 622 *(p_f2) = y; 623 } 624 625 /*-------------------------------------------------------------* 626 * find the LSPs using the Chebychev pol. evaluation * 627 *-------------------------------------------------------------*/ 628 629 nf = 0; /* number of found frequencies */ 630 ip = 0; /* indicator for f1 or f2 */ 631 632 coef = f1; 633 634 xlow = *(grid); 635 ylow = Chebps(xlow, coef, NC, pOverflow); 636 637 j = 0; 638 639 while ((nf < M) && (j < grid_points)) 640 { 641 j++; 642 xhigh = xlow; 643 yhigh = ylow; 644 xlow = *(grid + j); 645 ylow = Chebps(xlow, coef, NC, pOverflow); 646 647 if (((Word32)ylow*yhigh) <= 0) 648 { 649 /* divide 4 times the interval */ 650 for (i = 4; i != 0; i--) 651 { 652 /* xmid = (xlow + xhigh)/2 */ 653 x = xlow >> 1; 654 y = xhigh >> 1; 655 xmid = x + y; 656 657 ymid = Chebps(xmid, coef, NC, pOverflow); 658 659 if (((Word32)ylow*ymid) <= 0) 660 { 661 yhigh = ymid; 662 xhigh = xmid; 663 } 664 else 665 { 666 ylow = ymid; 667 xlow = xmid; 668 } 669 } 670 671 /*-------------------------------------------------------------* 672 * Linear interpolation * 673 * xint = xlow - ylow*(xhigh-xlow)/(yhigh-ylow); * 674 *-------------------------------------------------------------*/ 675 676 x = xhigh - xlow; 677 y = yhigh - ylow; 678 679 if (y == 0) 680 { 681 xint = xlow; 682 } 683 else 684 { 685 sign = y; 686 y = abs_s(y); 687 exp = norm_s(y); 688 y <<= exp; 689 y = div_s((Word16) 16383, y); 690 691 y = ((Word32)x * y) >> (19 - exp); 692 693 if (sign < 0) 694 { 695 y = -y; 696 } 697 698 /* xint = xlow - ylow*y */ 699 xint = xlow - (((Word32) ylow * y) >> 10); 700 } 701 702 *(lsp + nf) = xint; 703 xlow = xint; 704 nf++; 705 706 if (ip == 0) 707 { 708 ip = 1; 709 coef = f2; 710 } 711 else 712 { 713 ip = 0; 714 coef = f1; 715 } 716 717 ylow = Chebps(xlow, coef, NC, pOverflow); 718 719 } 720 } 721 722 /* Check if M roots found */ 723 724 if (nf < M) 725 { 726 for (i = NC; i != 0 ; i--) 727 { 728 *lsp++ = *old_lsp++; 729 *lsp++ = *old_lsp++; 730 } 731 } 732 733} 734 735Word16 Chebps_Wrapper(Word16 x, 736 Word16 f[], /* (n) */ 737 Word16 n, 738 Flag *pOverflow) 739{ 740 return Chebps(x, f, n, pOverflow); 741} 742 743