az_lsp.cpp revision 4f1efc098cb5791c3e9f483f2af84aef70d2d0a0
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 241static Word16 Chebps(Word16 x, 242 Word16 f[], /* (n) */ 243 Word16 n, 244 Flag *pOverflow) 245{ 246 Word16 i; 247 Word16 cheb; 248 Word16 b1_h; 249 Word16 b1_l; 250 Word32 t0; 251 Word32 L_temp; 252 Word16 *p_f = &f[1]; 253 254 OSCL_UNUSED_ARG(pOverflow); 255 256 /* L_temp = 1.0 */ 257 258 L_temp = 0x01000000L; 259 260 t0 = ((Word32) x << 10) + ((Word32) * (p_f++) << 14); 261 262 /* b1 = t0 = 2*x + f[1] */ 263 264 b1_h = (Word16)(t0 >> 16); 265 b1_l = (Word16)((t0 >> 1) - (b1_h << 15)); 266 267 268 for (i = 2; i < n; i++) 269 { 270 /* t0 = 2.0*x*b1 */ 271 t0 = ((Word32) b1_h * x); 272 t0 += ((Word32) b1_l * x) >> 15; 273 t0 <<= 2; 274 275 /* t0 = 2.0*x*b1 - b2 */ 276 t0 -= L_temp; 277 278 /* t0 = 2.0*x*b1 - b2 + f[i] */ 279 t0 += (Word32) * (p_f++) << 14; 280 281 L_temp = ((Word32) b1_h << 16) + ((Word32) b1_l << 1); 282 283 /* b0 = 2.0*x*b1 - b2 + f[i]*/ 284 b1_h = (Word16)(t0 >> 16); 285 b1_l = (Word16)((t0 >> 1) - (b1_h << 15)); 286 287 } 288 289 /* t0 = x*b1; */ 290 t0 = ((Word32) b1_h * x); 291 t0 += ((Word32) b1_l * x) >> 15; 292 t0 <<= 1; 293 294 295 /* t0 = x*b1 - b2 */ 296 t0 -= L_temp; 297 298 /* t0 = x*b1 - b2 + f[i]/2 */ 299 t0 += (Word32) * (p_f) << 13; 300 301 302 if ((UWord32)(t0 - 0xfe000000L) < 0x01ffffffL - 0xfe000000L) 303 { 304 cheb = (Word16)(t0 >> 10); 305 } 306 else 307 { 308 if (t0 > (Word32) 0x01ffffffL) 309 { 310 cheb = MAX_16; 311 312 } 313 else 314 { 315 cheb = MIN_16; 316 } 317 } 318 319 return (cheb); 320} 321 322 323/* 324------------------------------------------------------------------------------ 325 FUNCTION NAME: Az_lsp 326------------------------------------------------------------------------------ 327 INPUT AND OUTPUT DEFINITIONS FOR Az_lsp 328 329 Inputs: 330 a = predictor coefficients (Word16) 331 lsp = line spectral pairs (Word16) 332 old_lsp = old line spectral pairs (Word16) 333 334 pOverflow = pointer to overflow (Flag) 335 336 Outputs: 337 pOverflow -> 1 if the operations in the function resulted in saturation. 338 339 Returns: 340 None. 341 342 Global Variables Used: 343 None. 344 345 Local Variables Needed: 346 None. 347 348------------------------------------------------------------------------------ 349 FUNCTION DESCRIPTION 350 351 This function computes the LSPs from the LP coefficients. 352 353 The sum and difference filters are computed and divided by 1+z^{-1} and 354 1-z^{-1}, respectively. 355 356 f1[i] = a[i] + a[11-i] - f1[i-1] ; i=1,...,5 357 f2[i] = a[i] - a[11-i] + f2[i-1] ; i=1,...,5 358 359 The roots of F1(z) and F2(z) are found using Chebyshev polynomial evaluation. 360 The polynomials are evaluated at 60 points regularly spaced in the 361 frequency domain. The sign change interval is subdivided 4 times to better 362 track the root. The LSPs are found in the cosine domain [1,-1]. 363 364 If less than 10 roots are found, the LSPs from the past frame are used. 365 366------------------------------------------------------------------------------ 367 REQUIREMENTS 368 369 None. 370 371------------------------------------------------------------------------------ 372 REFERENCES 373 374 az_lsp.c, UMTS GSM AMR speech codec, R99 - Version 3.2.0, March 2, 2001 375 376------------------------------------------------------------------------------ 377 PSEUDO-CODE 378 379void Az_lsp ( 380 Word16 a[], // (i) : predictor coefficients (MP1) 381 Word16 lsp[], // (o) : line spectral pairs (M) 382 Word16 old_lsp[] // (i) : old lsp[] (in case not found 10 roots) (M) 383) 384{ 385 Word16 i, j, nf, ip; 386 Word16 xlow, ylow, xhigh, yhigh, xmid, ymid, xint; 387 Word16 x, y, sign, exp; 388 Word16 *coef; 389 Word16 f1[M / 2 + 1], f2[M / 2 + 1]; 390 Word32 t0; 391 392 *-------------------------------------------------------------* 393 * find the sum and diff. pol. F1(z) and F2(z) * 394 * F1(z) <--- F1(z)/(1+z**-1) & F2(z) <--- F2(z)/(1-z**-1) * 395 * * 396 * f1[0] = 1.0; * 397 * f2[0] = 1.0; * 398 * * 399 * for (i = 0; i< NC; i++) * 400 * { * 401 * f1[i+1] = a[i+1] + a[M-i] - f1[i] ; * 402 * f2[i+1] = a[i+1] - a[M-i] + f2[i] ; * 403 * } * 404 *-------------------------------------------------------------* 405 406 f1[0] = 1024; // f1[0] = 1.0 407 f2[0] = 1024; // f2[0] = 1.0 408 409// The reference ETSI code uses a global flag for Overflow. However, in the 410// actual implementation a pointer to Overflow flag is passed in as a 411// parameter to the function. This pointer is passed into all the basic math 412// functions invoked 413 414 for (i = 0; i < NC; i++) 415 { 416 t0 = L_mult (a[i + 1], 8192); // x = (a[i+1] + a[M-i]) >> 2 417 t0 = L_mac (t0, a[M - i], 8192); 418 x = extract_h (t0); 419 // f1[i+1] = a[i+1] + a[M-i] - f1[i] 420 f1[i + 1] = sub (x, f1[i]); 421 422 t0 = L_mult (a[i + 1], 8192); // x = (a[i+1] - a[M-i]) >> 2 423 t0 = L_msu (t0, a[M - i], 8192); 424 x = extract_h (t0); 425 // f2[i+1] = a[i+1] - a[M-i] + f2[i] 426 f2[i + 1] = add (x, f2[i]); 427 } 428 429 *-------------------------------------------------------------* 430 * find the LSPs using the Chebychev pol. evaluation * 431 *-------------------------------------------------------------* 432 433 nf = 0; // number of found frequencies 434 ip = 0; // indicator for f1 or f2 435 436 coef = f1; 437 438 xlow = grid[0]; 439 ylow = Chebps (xlow, coef, NC); 440 441 j = 0; 442 // while ( (nf < M) && (j < grid_points) ) 443 while ((sub (nf, M) < 0) && (sub (j, grid_points) < 0)) 444 { 445 j++; 446 xhigh = xlow; 447 yhigh = ylow; 448 xlow = grid[j]; 449 ylow = Chebps (xlow, coef, NC); 450 451 if (L_mult (ylow, yhigh) <= (Word32) 0L) 452 { 453 454 // divide 4 times the interval 455 456 for (i = 0; i < 4; i++) 457 { 458 // xmid = (xlow + xhigh)/2 459 xmid = add (shr (xlow, 1), shr (xhigh, 1)); 460 ymid = Chebps (xmid, coef, NC); 461 462 if (L_mult (ylow, ymid) <= (Word32) 0L) 463 { 464 yhigh = ymid; 465 xhigh = xmid; 466 } 467 else 468 { 469 ylow = ymid; 470 xlow = xmid; 471 } 472 } 473 474 *-------------------------------------------------------------* 475 * Linear interpolation * 476 * xint = xlow - ylow*(xhigh-xlow)/(yhigh-ylow); * 477 *-------------------------------------------------------------* 478 479 x = sub (xhigh, xlow); 480 y = sub (yhigh, ylow); 481 482 if (y == 0) 483 { 484 xint = xlow; 485 } 486 else 487 { 488 sign = y; 489 y = abs_s (y); 490 exp = norm_s (y); 491 y = shl (y, exp); 492 y = div_s ((Word16) 16383, y); 493 t0 = L_mult (x, y); 494 t0 = L_shr (t0, sub (20, exp)); 495 y = extract_l (t0); // y= (xhigh-xlow)/(yhigh-ylow) 496 497 if (sign < 0) 498 y = negate (y); 499 500 t0 = L_mult (ylow, y); 501 t0 = L_shr (t0, 11); 502 xint = sub (xlow, extract_l (t0)); // xint = xlow - ylow*y 503 } 504 505 lsp[nf] = xint; 506 xlow = xint; 507 nf++; 508 509 if (ip == 0) 510 { 511 ip = 1; 512 coef = f2; 513 } 514 else 515 { 516 ip = 0; 517 coef = f1; 518 } 519 ylow = Chebps (xlow, coef, NC); 520 521 } 522 } 523 524 // Check if M roots found 525 526 if (sub (nf, M) < 0) 527 { 528 for (i = 0; i < M; i++) 529 { 530 lsp[i] = old_lsp[i]; 531 } 532 533 } 534 return; 535} 536 537------------------------------------------------------------------------------ 538 RESOURCES USED [optional] 539 540 When the code is written for a specific target processor the 541 the resources used should be documented below. 542 543 HEAP MEMORY USED: x bytes 544 545 STACK MEMORY USED: x bytes 546 547 CLOCK CYCLES: (cycle count equation for this function) + (variable 548 used to represent cycle count for each subroutine 549 called) 550 where: (cycle count variable) = cycle count for [subroutine 551 name] 552 553------------------------------------------------------------------------------ 554 CAUTION [optional] 555 [State any special notes, constraints or cautions for users of this function] 556 557------------------------------------------------------------------------------ 558*/ 559 560void Az_lsp( 561 Word16 a[], /* (i) : predictor coefficients (MP1) */ 562 Word16 lsp[], /* (o) : line spectral pairs (M) */ 563 Word16 old_lsp[], /* (i) : old lsp[] (in case not found 10 roots) (M) */ 564 Flag *pOverflow /* (i/o): overflow flag */ 565) 566{ 567 register Word16 i; 568 register Word16 j; 569 register Word16 nf; 570 register Word16 ip; 571 Word16 xlow; 572 Word16 ylow; 573 Word16 xhigh; 574 Word16 yhigh; 575 Word16 xmid; 576 Word16 ymid; 577 Word16 xint; 578 Word16 x; 579 Word16 y; 580 Word16 sign; 581 Word16 exp; 582 Word16 *coef; 583 Word16 f1[NC + 1]; 584 Word16 f2[NC + 1]; 585 Word32 L_temp1; 586 Word32 L_temp2; 587 Word16 *p_f1 = f1; 588 Word16 *p_f2 = f2; 589 590 /*-------------------------------------------------------------* 591 * find the sum and diff. pol. F1(z) and F2(z) * 592 * F1(z) <--- F1(z)/(1+z**-1) & F2(z) <--- F2(z)/(1-z**-1) * 593 * * 594 * f1[0] = 1.0; * 595 * f2[0] = 1.0; * 596 * * 597 * for (i = 0; i< NC; i++) * 598 * { * 599 * f1[i+1] = a[i+1] + a[M-i] - f1[i] ; * 600 * f2[i+1] = a[i+1] - a[M-i] + f2[i] ; * 601 * } * 602 *-------------------------------------------------------------*/ 603 604 *p_f1 = 1024; /* f1[0] = 1.0 */ 605 *p_f2 = 1024; /* f2[0] = 1.0 */ 606 607 for (i = 0; i < NC; i++) 608 { 609 L_temp1 = (Word32) * (a + i + 1); 610 L_temp2 = (Word32) * (a + M - i); 611 /* x = (a[i+1] + a[M-i]) >> 2 */ 612 x = (Word16)((L_temp1 + L_temp2) >> 2); 613 /* y = (a[i+1] - a[M-i]) >> 2 */ 614 y = (Word16)((L_temp1 - L_temp2) >> 2); 615 /* f1[i+1] = a[i+1] + a[M-i] - f1[i] */ 616 x -= *(p_f1++); 617 *(p_f1) = x; 618 /* f2[i+1] = a[i+1] - a[M-i] + f2[i] */ 619 y += *(p_f2++); 620 *(p_f2) = y; 621 } 622 623 /*-------------------------------------------------------------* 624 * find the LSPs using the Chebychev pol. evaluation * 625 *-------------------------------------------------------------*/ 626 627 nf = 0; /* number of found frequencies */ 628 ip = 0; /* indicator for f1 or f2 */ 629 630 coef = f1; 631 632 xlow = *(grid); 633 ylow = Chebps(xlow, coef, NC, pOverflow); 634 635 j = 0; 636 637 while ((nf < M) && (j < grid_points)) 638 { 639 j++; 640 xhigh = xlow; 641 yhigh = ylow; 642 xlow = *(grid + j); 643 ylow = Chebps(xlow, coef, NC, pOverflow); 644 645 if (((Word32)ylow*yhigh) <= 0) 646 { 647 /* divide 4 times the interval */ 648 for (i = 4; i != 0; i--) 649 { 650 /* xmid = (xlow + xhigh)/2 */ 651 x = xlow >> 1; 652 y = xhigh >> 1; 653 xmid = x + y; 654 655 ymid = Chebps(xmid, coef, NC, pOverflow); 656 657 if (((Word32)ylow*ymid) <= 0) 658 { 659 yhigh = ymid; 660 xhigh = xmid; 661 } 662 else 663 { 664 ylow = ymid; 665 xlow = xmid; 666 } 667 } 668 669 /*-------------------------------------------------------------* 670 * Linear interpolation * 671 * xint = xlow - ylow*(xhigh-xlow)/(yhigh-ylow); * 672 *-------------------------------------------------------------*/ 673 674 x = xhigh - xlow; 675 y = yhigh - ylow; 676 677 if (y == 0) 678 { 679 xint = xlow; 680 } 681 else 682 { 683 sign = y; 684 y = abs_s(y); 685 exp = norm_s(y); 686 y <<= exp; 687 y = div_s((Word16) 16383, y); 688 689 y = ((Word32)x * y) >> (19 - exp); 690 691 if (sign < 0) 692 { 693 y = -y; 694 } 695 696 /* xint = xlow - ylow*y */ 697 xint = xlow - (((Word32) ylow * y) >> 10); 698 } 699 700 *(lsp + nf) = xint; 701 xlow = xint; 702 nf++; 703 704 if (ip == 0) 705 { 706 ip = 1; 707 coef = f2; 708 } 709 else 710 { 711 ip = 0; 712 coef = f1; 713 } 714 715 ylow = Chebps(xlow, coef, NC, pOverflow); 716 717 } 718 } 719 720 /* Check if M roots found */ 721 722 if (nf < M) 723 { 724 for (i = NC; i != 0 ; i--) 725 { 726 *lsp++ = *old_lsp++; 727 *lsp++ = *old_lsp++; 728 } 729 } 730 731} 732 733Word16 Chebps_Wrapper(Word16 x, 734 Word16 f[], /* (n) */ 735 Word16 n, 736 Flag *pOverflow) 737{ 738 return Chebps(x, f, n, pOverflow); 739} 740 741