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 */
58 static __inline Word16 Chebps2(Word16 x, Word16 f[], Word32 n);
59
Az_isp(Word16 a[],Word16 isp[],Word16 old_isp[])60 void 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
Chebps2(Word16 x,Word16 f[],Word32 n)214 static __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