/* ------------------------------------------------------------------ * Copyright (C) 1998-2009 PacketVideo * * Licensed under the Apache License, Version 2.0 (the "License"); * you may not use this file except in compliance with the License. * You may obtain a copy of the License at * * http://www.apache.org/licenses/LICENSE-2.0 * * Unless required by applicable law or agreed to in writing, software * distributed under the License is distributed on an "AS IS" BASIS, * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either * express or implied. * See the License for the specific language governing permissions * and limitations under the License. * ------------------------------------------------------------------- */ /**************************************************************************************** Portions of this file are derived from the following 3GPP standard: 3GPP TS 26.073 ANSI-C code for the Adaptive Multi-Rate (AMR) speech codec Available from http://www.3gpp.org (C) 2004, 3GPP Organizational Partners (ARIB, ATIS, CCSA, ETSI, TTA, TTC) Permission to distribute, modify and use this file under the standard license terms listed above has been obtained from the copyright holder. ****************************************************************************************/ /* Pathname: ./audio/gsm-amr/c/src/az_lsp.c Funtions: Chebps Chebps_Wrapper Az_lsp ------------------------------------------------------------------------------ REVISION HISTORY Description: Finished first pass of optimization. Description: Made changes based on review comments. Description: Made input to Chebps_Wrapper consistent with that of Chebps. Description: Replaced current Pseudo-code with the UMTS code version 3.2.0. Updated coding template. Description: Replaced basic_op.h and oper_32b.h with the header files of the math functions used by the file. Description: Made the following changes per comments from Phase 2/3 review: 1. Used "-" operator instead of calling sub function in the az_lsp() code. 2. Copied detailed function description of az_lsp from the header file. 3. Modified local variable definition to one per line. 4. Used NC in the definition of f1 and f2 arrays. 5. Added curly brackets in the IF statement. Description: Changed function interface to pass in a pointer to overflow flag into the function instead of using a global flag. Removed inclusion of unneeded header files. Description: For Chebps() and Az_lsp() 1. Eliminated unused include files. 2. Replaced array addressing by pointers 3. Eliminated math operations that unnecessary checked for saturation. 4. Eliminated not needed variables 5. Eliminated if-else checks for saturation 6. Deleted unused function cheps_wraper Description: Added casting to eliminate warnings Who: Date: Description: ------------------------------------------------------------------------------ MODULE DESCRIPTION These modules compute the LSPs from the LP coefficients. ------------------------------------------------------------------------------ */ /*---------------------------------------------------------------------------- ; INCLUDES ----------------------------------------------------------------------------*/ #include "az_lsp.h" #include "cnst.h" #include "basic_op.h" /*---------------------------------------------------------------------------- ; MACROS ; Define module specific macros here ----------------------------------------------------------------------------*/ /*---------------------------------------------------------------------------- ; DEFINES ; Include all pre-processor statements here. Include conditional ; compile variables also. ----------------------------------------------------------------------------*/ #define NC M/2 /* M = LPC order, NC = M/2 */ /*---------------------------------------------------------------------------- ; LOCAL FUNCTION DEFINITIONS ; Function Prototype declaration ----------------------------------------------------------------------------*/ /*---------------------------------------------------------------------------- ; LOCAL VARIABLE DEFINITIONS ; Variable declaration - defined here and used outside this module ----------------------------------------------------------------------------*/ /* ------------------------------------------------------------------------------ FUNCTION NAME: Chebps ------------------------------------------------------------------------------ INPUT AND OUTPUT DEFINITIONS Inputs: x = input value (Word16) f = polynomial (Word16) n = polynomial order (Word16) pOverflow = pointer to overflow (Flag) Outputs: pOverflow -> 1 if the operations in the function resulted in saturation. Returns: cheb = Chebyshev polynomial for the input value x.(Word16) Global Variables Used: None. Local Variables Needed: None. ------------------------------------------------------------------------------ FUNCTION DESCRIPTION This module evaluates the Chebyshev polynomial series. - The polynomial order is n = m/2 = 5 - The polynomial F(z) (F1(z) or F2(z)) is given by F(w) = 2 exp(-j5w) C(x) where C(x) = T_n(x) + f(1)T_n-1(x) + ... +f(n-1)T_1(x) + f(n)/2 and T_m(x) = cos(mw) is the mth order Chebyshev polynomial ( x=cos(w) ) - C(x) for the input x is returned. ------------------------------------------------------------------------------ REQUIREMENTS None. ------------------------------------------------------------------------------ REFERENCES az_lsp.c, UMTS GSM AMR speech codec, R99 - Version 3.2.0, March 2, 2001 ------------------------------------------------------------------------------ PSEUDO-CODE static Word16 Chebps (Word16 x, Word16 f[], // (n) Word16 n) { Word16 i, cheb; Word16 b0_h, b0_l, b1_h, b1_l, b2_h, b2_l; Word32 t0; // The reference ETSI code uses a global flag for Overflow. However, in the // actual implementation a pointer to Overflow flag is passed in as a // parameter to the function. This pointer is passed into all the basic math // functions invoked b2_h = 256; // b2 = 1.0 b2_l = 0; t0 = L_mult (x, 512); // 2*x t0 = L_mac (t0, f[1], 8192); // + f[1] L_Extract (t0, &b1_h, &b1_l); // b1 = 2*x + f[1] for (i = 2; i < n; i++) { t0 = Mpy_32_16 (b1_h, b1_l, x); // t0 = 2.0*x*b1 t0 = L_shl (t0, 1); t0 = L_mac (t0, b2_h, (Word16) 0x8000); // t0 = 2.0*x*b1 - b2 t0 = L_msu (t0, b2_l, 1); t0 = L_mac (t0, f[i], 8192); // t0 = 2.0*x*b1 - b2 + f[i] L_Extract (t0, &b0_h, &b0_l); // b0 = 2.0*x*b1 - b2 + f[i] b2_l = b1_l; // b2 = b1; b2_h = b1_h; b1_l = b0_l; // b1 = b0; b1_h = b0_h; } t0 = Mpy_32_16 (b1_h, b1_l, x); // t0 = x*b1; t0 = L_mac (t0, b2_h, (Word16) 0x8000); // t0 = x*b1 - b2 t0 = L_msu (t0, b2_l, 1); t0 = L_mac (t0, f[i], 4096); // t0 = x*b1 - b2 + f[i]/2 t0 = L_shl (t0, 6); cheb = extract_h (t0); return (cheb); } ------------------------------------------------------------------------------ RESOURCES USED [optional] When the code is written for a specific target processor the the resources used should be documented below. HEAP MEMORY USED: x bytes STACK MEMORY USED: x bytes CLOCK CYCLES: (cycle count equation for this function) + (variable used to represent cycle count for each subroutine called) where: (cycle count variable) = cycle count for [subroutine name] ------------------------------------------------------------------------------ CAUTION [optional] [State any special notes, constraints or cautions for users of this function] ------------------------------------------------------------------------------ */ static Word16 Chebps(Word16 x, Word16 f[], /* (n) */ Word16 n, Flag *pOverflow) { Word16 i; Word16 cheb; Word16 b1_h; Word16 b1_l; Word32 t0; Word32 L_temp; Word16 *p_f = &f[1]; OSCL_UNUSED_ARG(pOverflow); /* L_temp = 1.0 */ L_temp = 0x01000000L; t0 = ((Word32) x << 10) + ((Word32) * (p_f++) << 14); /* b1 = t0 = 2*x + f[1] */ b1_h = (Word16)(t0 >> 16); b1_l = (Word16)((t0 >> 1) - (b1_h << 15)); for (i = 2; i < n; i++) { /* t0 = 2.0*x*b1 */ t0 = ((Word32) b1_h * x); t0 += ((Word32) b1_l * x) >> 15; t0 <<= 2; /* t0 = 2.0*x*b1 - b2 */ t0 -= L_temp; /* t0 = 2.0*x*b1 - b2 + f[i] */ t0 += (Word32) * (p_f++) << 14; L_temp = ((Word32) b1_h << 16) + ((Word32) b1_l << 1); /* b0 = 2.0*x*b1 - b2 + f[i]*/ b1_h = (Word16)(t0 >> 16); b1_l = (Word16)((t0 >> 1) - (b1_h << 15)); } /* t0 = x*b1; */ t0 = ((Word32) b1_h * x); t0 += ((Word32) b1_l * x) >> 15; t0 <<= 1; /* t0 = x*b1 - b2 */ t0 -= L_temp; /* t0 = x*b1 - b2 + f[i]/2 */ t0 += (Word32) * (p_f) << 13; if ((UWord32)(t0 - 0xfe000000L) < (UWord32)0x03ffffffL) { cheb = (Word16)(t0 >> 10); } else { if (t0 > (Word32) 0x01ffffffL) { cheb = MAX_16; } else { cheb = MIN_16; } } return (cheb); } /* ------------------------------------------------------------------------------ FUNCTION NAME: Az_lsp ------------------------------------------------------------------------------ INPUT AND OUTPUT DEFINITIONS FOR Az_lsp Inputs: a = predictor coefficients (Word16) lsp = line spectral pairs (Word16) old_lsp = old line spectral pairs (Word16) pOverflow = pointer to overflow (Flag) Outputs: pOverflow -> 1 if the operations in the function resulted in saturation. Returns: None. Global Variables Used: None. Local Variables Needed: None. ------------------------------------------------------------------------------ FUNCTION DESCRIPTION This function computes the LSPs from the LP coefficients. The sum and difference filters are computed and divided by 1+z^{-1} and 1-z^{-1}, respectively. f1[i] = a[i] + a[11-i] - f1[i-1] ; i=1,...,5 f2[i] = a[i] - a[11-i] + f2[i-1] ; i=1,...,5 The roots of F1(z) and F2(z) are found using Chebyshev polynomial evaluation. The polynomials are evaluated at 60 points regularly spaced in the frequency domain. The sign change interval is subdivided 4 times to better track the root. The LSPs are found in the cosine domain [1,-1]. If less than 10 roots are found, the LSPs from the past frame are used. ------------------------------------------------------------------------------ REQUIREMENTS None. ------------------------------------------------------------------------------ REFERENCES az_lsp.c, UMTS GSM AMR speech codec, R99 - Version 3.2.0, March 2, 2001 ------------------------------------------------------------------------------ PSEUDO-CODE void Az_lsp ( Word16 a[], // (i) : predictor coefficients (MP1) Word16 lsp[], // (o) : line spectral pairs (M) Word16 old_lsp[] // (i) : old lsp[] (in case not found 10 roots) (M) ) { Word16 i, j, nf, ip; Word16 xlow, ylow, xhigh, yhigh, xmid, ymid, xint; Word16 x, y, sign, exp; Word16 *coef; Word16 f1[M / 2 + 1], f2[M / 2 + 1]; Word32 t0; *-------------------------------------------------------------* * find the sum and diff. pol. F1(z) and F2(z) * * F1(z) <--- F1(z)/(1+z**-1) & F2(z) <--- F2(z)/(1-z**-1) * * * * f1[0] = 1.0; * * f2[0] = 1.0; * * * * for (i = 0; i< NC; i++) * * { * * f1[i+1] = a[i+1] + a[M-i] - f1[i] ; * * f2[i+1] = a[i+1] - a[M-i] + f2[i] ; * * } * *-------------------------------------------------------------* f1[0] = 1024; // f1[0] = 1.0 f2[0] = 1024; // f2[0] = 1.0 // The reference ETSI code uses a global flag for Overflow. However, in the // actual implementation a pointer to Overflow flag is passed in as a // parameter to the function. This pointer is passed into all the basic math // functions invoked for (i = 0; i < NC; i++) { t0 = L_mult (a[i + 1], 8192); // x = (a[i+1] + a[M-i]) >> 2 t0 = L_mac (t0, a[M - i], 8192); x = extract_h (t0); // f1[i+1] = a[i+1] + a[M-i] - f1[i] f1[i + 1] = sub (x, f1[i]); t0 = L_mult (a[i + 1], 8192); // x = (a[i+1] - a[M-i]) >> 2 t0 = L_msu (t0, a[M - i], 8192); x = extract_h (t0); // f2[i+1] = a[i+1] - a[M-i] + f2[i] f2[i + 1] = add (x, f2[i]); } *-------------------------------------------------------------* * find the LSPs using the Chebychev pol. evaluation * *-------------------------------------------------------------* nf = 0; // number of found frequencies ip = 0; // indicator for f1 or f2 coef = f1; xlow = grid[0]; ylow = Chebps (xlow, coef, NC); j = 0; // while ( (nf < M) && (j < grid_points) ) while ((sub (nf, M) < 0) && (sub (j, grid_points) < 0)) { j++; xhigh = xlow; yhigh = ylow; xlow = grid[j]; ylow = Chebps (xlow, coef, NC); if (L_mult (ylow, yhigh) <= (Word32) 0L) { // divide 4 times the interval for (i = 0; i < 4; i++) { // xmid = (xlow + xhigh)/2 xmid = add (shr (xlow, 1), shr (xhigh, 1)); ymid = Chebps (xmid, coef, NC); if (L_mult (ylow, ymid) <= (Word32) 0L) { yhigh = ymid; xhigh = xmid; } else { ylow = ymid; xlow = xmid; } } *-------------------------------------------------------------* * Linear interpolation * * xint = xlow - ylow*(xhigh-xlow)/(yhigh-ylow); * *-------------------------------------------------------------* x = sub (xhigh, xlow); y = sub (yhigh, ylow); if (y == 0) { xint = xlow; } else { sign = y; y = abs_s (y); exp = norm_s (y); y = shl (y, exp); y = div_s ((Word16) 16383, y); t0 = L_mult (x, y); t0 = L_shr (t0, sub (20, exp)); y = extract_l (t0); // y= (xhigh-xlow)/(yhigh-ylow) if (sign < 0) y = negate (y); t0 = L_mult (ylow, y); t0 = L_shr (t0, 11); xint = sub (xlow, extract_l (t0)); // xint = xlow - ylow*y } lsp[nf] = xint; xlow = xint; nf++; if (ip == 0) { ip = 1; coef = f2; } else { ip = 0; coef = f1; } ylow = Chebps (xlow, coef, NC); } } // Check if M roots found if (sub (nf, M) < 0) { for (i = 0; i < M; i++) { lsp[i] = old_lsp[i]; } } return; } ------------------------------------------------------------------------------ RESOURCES USED [optional] When the code is written for a specific target processor the the resources used should be documented below. HEAP MEMORY USED: x bytes STACK MEMORY USED: x bytes CLOCK CYCLES: (cycle count equation for this function) + (variable used to represent cycle count for each subroutine called) where: (cycle count variable) = cycle count for [subroutine name] ------------------------------------------------------------------------------ CAUTION [optional] [State any special notes, constraints or cautions for users of this function] ------------------------------------------------------------------------------ */ void Az_lsp( Word16 a[], /* (i) : predictor coefficients (MP1) */ Word16 lsp[], /* (o) : line spectral pairs (M) */ Word16 old_lsp[], /* (i) : old lsp[] (in case not found 10 roots) (M) */ Flag *pOverflow /* (i/o): overflow flag */ ) { Word16 i; Word16 j; Word16 nf; Word16 ip; Word16 xlow; Word16 ylow; Word16 xhigh; Word16 yhigh; Word16 xmid; Word16 ymid; Word16 xint; Word16 x; Word16 y; Word16 sign; Word16 exp; Word16 *coef; Word16 f1[NC + 1]; Word16 f2[NC + 1]; Word32 L_temp1; Word32 L_temp2; Word16 *p_f1 = f1; Word16 *p_f2 = f2; /*-------------------------------------------------------------* * find the sum and diff. pol. F1(z) and F2(z) * * F1(z) <--- F1(z)/(1+z**-1) & F2(z) <--- F2(z)/(1-z**-1) * * * * f1[0] = 1.0; * * f2[0] = 1.0; * * * * for (i = 0; i< NC; i++) * * { * * f1[i+1] = a[i+1] + a[M-i] - f1[i] ; * * f2[i+1] = a[i+1] - a[M-i] + f2[i] ; * * } * *-------------------------------------------------------------*/ *p_f1 = 1024; /* f1[0] = 1.0 */ *p_f2 = 1024; /* f2[0] = 1.0 */ for (i = 0; i < NC; i++) { L_temp1 = (Word32) * (a + i + 1); L_temp2 = (Word32) * (a + M - i); /* x = (a[i+1] + a[M-i]) >> 2 */ x = (Word16)((L_temp1 + L_temp2) >> 2); /* y = (a[i+1] - a[M-i]) >> 2 */ y = (Word16)((L_temp1 - L_temp2) >> 2); /* f1[i+1] = a[i+1] + a[M-i] - f1[i] */ x -= *(p_f1++); *(p_f1) = x; /* f2[i+1] = a[i+1] - a[M-i] + f2[i] */ y += *(p_f2++); *(p_f2) = y; } /*-------------------------------------------------------------* * find the LSPs using the Chebychev pol. evaluation * *-------------------------------------------------------------*/ nf = 0; /* number of found frequencies */ ip = 0; /* indicator for f1 or f2 */ coef = f1; xlow = *(grid); ylow = Chebps(xlow, coef, NC, pOverflow); j = 0; while ((nf < M) && (j < grid_points)) { j++; xhigh = xlow; yhigh = ylow; xlow = *(grid + j); ylow = Chebps(xlow, coef, NC, pOverflow); if (((Word32)ylow*yhigh) <= 0) { /* divide 4 times the interval */ for (i = 4; i != 0; i--) { /* xmid = (xlow + xhigh)/2 */ x = xlow >> 1; y = xhigh >> 1; xmid = x + y; ymid = Chebps(xmid, coef, NC, pOverflow); if (((Word32)ylow*ymid) <= 0) { yhigh = ymid; xhigh = xmid; } else { ylow = ymid; xlow = xmid; } } /*-------------------------------------------------------------* * Linear interpolation * * xint = xlow - ylow*(xhigh-xlow)/(yhigh-ylow); * *-------------------------------------------------------------*/ x = xhigh - xlow; y = yhigh - ylow; if (y == 0) { xint = xlow; } else { sign = y; y = abs_s(y); exp = norm_s(y); y <<= exp; y = div_s((Word16) 16383, y); y = ((Word32)x * y) >> (19 - exp); if (sign < 0) { y = -y; } /* xint = xlow - ylow*y */ xint = xlow - (((Word32) ylow * y) >> 10); } *(lsp + nf) = xint; xlow = xint; nf++; if (ip == 0) { ip = 1; coef = f2; } else { ip = 0; coef = f1; } ylow = Chebps(xlow, coef, NC, pOverflow); } } /* Check if M roots found */ if (nf < M) { for (i = NC; i != 0 ; i--) { *lsp++ = *old_lsp++; *lsp++ = *old_lsp++; } } } Word16 Chebps_Wrapper(Word16 x, Word16 f[], /* (n) */ Word16 n, Flag *pOverflow) { return Chebps(x, f, n, pOverflow); }