1/* 2 ** Copyright 2003-2010, VisualOn, Inc. 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 express or implied. 13 ** See the License for the specific language governing permissions and 14 ** limitations under the License. 15 */ 16 17/*********************************************************************** 18* File: az_isp.c 19* 20* Description: 21*-----------------------------------------------------------------------* 22* Compute the ISPs from the LPC coefficients (order=M) * 23*-----------------------------------------------------------------------* 24* * 25* The ISPs are the roots of the two polynomials F1(z) and F2(z) * 26* defined as * 27* F1(z) = A(z) + z^-m A(z^-1) * 28* and F2(z) = A(z) - z^-m A(z^-1) * 29* * 30* For a even order m=2n, F1(z) has M/2 conjugate roots on the unit * 31* circle and F2(z) has M/2-1 conjugate roots on the unit circle in * 32* addition to two roots at 0 and pi. * 33* * 34* For a 16th order LP analysis, F1(z) and F2(z) can be written as * 35* * 36* F1(z) = (1 + a[M]) PRODUCT (1 - 2 cos(w_i) z^-1 + z^-2 ) * 37* i=0,2,4,6,8,10,12,14 * 38* * 39* F2(z) = (1 - a[M]) (1 - z^-2) PRODUCT (1 - 2 cos(w_i) z^-1 + z^-2 ) * 40* i=1,3,5,7,9,11,13 * 41* * 42* The ISPs are the M-1 frequencies w_i, i=0...M-2 plus the last * 43* predictor coefficient a[M]. * 44*-----------------------------------------------------------------------* 45 46************************************************************************/ 47 48#include "typedef.h" 49#include "basic_op.h" 50#include "oper_32b.h" 51#include "stdio.h" 52#include "grid100.tab" 53 54#define M 16 55#define NC (M/2) 56 57/* local function */ 58static __inline Word16 Chebps2(Word16 x, Word16 f[], Word32 n); 59 60void Az_isp( 61 Word16 a[], /* (i) Q12 : predictor coefficients */ 62 Word16 isp[], /* (o) Q15 : Immittance spectral pairs */ 63 Word16 old_isp[] /* (i) : old isp[] (in case not found M roots) */ 64 ) 65{ 66 Word32 i, j, nf, ip, order; 67 Word16 xlow, ylow, xhigh, yhigh, xmid, ymid, xint; 68 Word16 x, y, sign, exp; 69 Word16 *coef; 70 Word16 f1[NC + 1], f2[NC]; 71 Word32 t0; 72 /*-------------------------------------------------------------* 73 * find the sum and diff polynomials F1(z) and F2(z) * 74 * F1(z) = [A(z) + z^M A(z^-1)] * 75 * F2(z) = [A(z) - z^M A(z^-1)]/(1-z^-2) * 76 * * 77 * for (i=0; i<NC; i++) * 78 * { * 79 * f1[i] = a[i] + a[M-i]; * 80 * f2[i] = a[i] - a[M-i]; * 81 * } * 82 * f1[NC] = 2.0*a[NC]; * 83 * * 84 * for (i=2; i<NC; i++) Divide by (1-z^-2) * 85 * f2[i] += f2[i-2]; * 86 *-------------------------------------------------------------*/ 87 for (i = 0; i < NC; i++) 88 { 89 t0 = a[i] << 15; 90 f1[i] = vo_round(t0 + (a[M - i] << 15)); /* =(a[i]+a[M-i])/2 */ 91 f2[i] = vo_round(t0 - (a[M - i] << 15)); /* =(a[i]-a[M-i])/2 */ 92 } 93 f1[NC] = a[NC]; 94 for (i = 2; i < NC; i++) /* Divide by (1-z^-2) */ 95 f2[i] = add1(f2[i], f2[i - 2]); 96 97 /*---------------------------------------------------------------------* 98 * Find the ISPs (roots of F1(z) and F2(z) ) using the * 99 * Chebyshev polynomial evaluation. * 100 * The roots of F1(z) and F2(z) are alternatively searched. * 101 * We start by finding the first root of F1(z) then we switch * 102 * to F2(z) then back to F1(z) and so on until all roots are found. * 103 * * 104 * - Evaluate Chebyshev pol. at grid points and check for sign change.* 105 * - If sign change track the root by subdividing the interval * 106 * 2 times and ckecking sign change. * 107 *---------------------------------------------------------------------*/ 108 nf = 0; /* number of found frequencies */ 109 ip = 0; /* indicator for f1 or f2 */ 110 coef = f1; 111 order = NC; 112 xlow = vogrid[0]; 113 ylow = Chebps2(xlow, coef, order); 114 j = 0; 115 while ((nf < M - 1) && (j < GRID_POINTS)) 116 { 117 j ++; 118 xhigh = xlow; 119 yhigh = ylow; 120 xlow = vogrid[j]; 121 ylow = Chebps2(xlow, coef, order); 122 if ((ylow * yhigh) <= (Word32) 0) 123 { 124 /* divide 2 times the interval */ 125 for (i = 0; i < 2; i++) 126 { 127 xmid = (xlow >> 1) + (xhigh >> 1); /* xmid = (xlow + xhigh)/2 */ 128 ymid = Chebps2(xmid, coef, order); 129 if ((ylow * ymid) <= (Word32) 0) 130 { 131 yhigh = ymid; 132 xhigh = xmid; 133 } else 134 { 135 ylow = ymid; 136 xlow = xmid; 137 } 138 } 139 /*-------------------------------------------------------------* 140 * Linear interpolation * 141 * xint = xlow - ylow*(xhigh-xlow)/(yhigh-ylow); * 142 *-------------------------------------------------------------*/ 143 x = xhigh - xlow; 144 y = yhigh - ylow; 145 if (y == 0) 146 { 147 xint = xlow; 148 } else 149 { 150 sign = y; 151 y = abs_s(y); 152 exp = norm_s(y); 153 y = y << exp; 154 y = div_s((Word16) 16383, y); 155 t0 = x * y; 156 t0 = (t0 >> (19 - exp)); 157 y = vo_extract_l(t0); /* y= (xhigh-xlow)/(yhigh-ylow) in Q11 */ 158 if (sign < 0) 159 y = -y; 160 t0 = ylow * y; /* result in Q26 */ 161 t0 = (t0 >> 10); /* result in Q15 */ 162 xint = vo_sub(xlow, vo_extract_l(t0)); /* xint = xlow - ylow*y */ 163 } 164 isp[nf] = xint; 165 xlow = xint; 166 nf++; 167 if (ip == 0) 168 { 169 ip = 1; 170 coef = f2; 171 order = NC - 1; 172 } else 173 { 174 ip = 0; 175 coef = f1; 176 order = NC; 177 } 178 ylow = Chebps2(xlow, coef, order); 179 } 180 } 181 /* Check if M-1 roots found */ 182 if(nf < M - 1) 183 { 184 for (i = 0; i < M; i++) 185 { 186 isp[i] = old_isp[i]; 187 } 188 } else 189 { 190 isp[M - 1] = a[M] << 3; /* From Q12 to Q15 with saturation */ 191 } 192 return; 193} 194 195/*--------------------------------------------------------------* 196* function Chebps2: * 197* ~~~~~~~ * 198* Evaluates the Chebishev polynomial series * 199*--------------------------------------------------------------* 200* * 201* The polynomial order is * 202* n = M/2 (M is the prediction order) * 203* The polynomial is given by * 204* C(x) = f(0)T_n(x) + f(1)T_n-1(x) + ... +f(n-1)T_1(x) + f(n)/2 * 205* Arguments: * 206* x: input value of evaluation; x = cos(frequency) in Q15 * 207* f[]: coefficients of the pol. in Q11 * 208* n: order of the pol. * 209* * 210* The value of C(x) is returned. (Satured to +-1.99 in Q14) * 211* * 212*--------------------------------------------------------------*/ 213 214static __inline Word16 Chebps2(Word16 x, Word16 f[], Word32 n) 215{ 216 Word32 i, cheb; 217 Word16 b0_h, b0_l, b1_h, b1_l, b2_h, b2_l; 218 Word32 t0; 219 220 /* Note: All computation are done in Q24. */ 221 222 t0 = f[0] << 13; 223 b2_h = t0 >> 16; 224 b2_l = (t0 & 0xffff)>>1; 225 226 t0 = ((b2_h * x)<<1) + (((b2_l * x)>>15)<<1); 227 t0 <<= 1; 228 t0 += (f[1] << 13); /* + f[1] in Q24 */ 229 230 b1_h = t0 >> 16; 231 b1_l = (t0 & 0xffff) >> 1; 232 233 for (i = 2; i < n; i++) 234 { 235 t0 = ((b1_h * x)<<1) + (((b1_l * x)>>15)<<1); 236 237 t0 += (b2_h * (-16384))<<1; 238 t0 += (f[i] << 12); 239 t0 <<= 1; 240 t0 -= (b2_l << 1); /* t0 = 2.0*x*b1 - b2 + f[i]; */ 241 242 b0_h = t0 >> 16; 243 b0_l = (t0 & 0xffff) >> 1; 244 245 b2_l = b1_l; /* b2 = b1; */ 246 b2_h = b1_h; 247 b1_l = b0_l; /* b1 = b0; */ 248 b1_h = b0_h; 249 } 250 251 t0 = ((b1_h * x)<<1) + (((b1_l * x)>>15)<<1); 252 t0 += (b2_h * (-32768))<<1; /* t0 = x*b1 - b2 */ 253 t0 -= (b2_l << 1); 254 t0 += (f[n] << 12); /* t0 = x*b1 - b2 + f[i]/2 */ 255 256 t0 = L_shl2(t0, 6); /* Q24 to Q30 with saturation */ 257 258 cheb = extract_h(t0); /* Result in Q14 */ 259 260 if (cheb == -32768) 261 { 262 cheb = -32767; /* to avoid saturation in Az_isp */ 263 } 264 return (cheb); 265} 266 267 268 269