FreeCalypso > hg > gsm-codec-lib

comparison libtwamr/az_lsp.c @ 253:54f6bc41ed10

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
libtwamr: integrate a* modules
author Mychaela Falconia <falcon@freecalypso.org>
date Fri, 05 Apr 2024 06:08:15 +0000
parents
children 07f936338de1
comparison
equal deleted inserted replaced
252:57b4053559ff 253:54f6bc41ed10
1 /*
2 ********************************************************************************
3 *
4 * GSM AMR-NB speech codec R98 Version 7.6.0 December 12, 2001
5 * R99 Version 3.3.0
6 * REL-4 Version 4.1.0
7 *
8 ********************************************************************************
9 *
10 * File : az_lsp.c
11 * Purpose : Compute the LSPs from the LP coefficients
12 *
13 ********************************************************************************
14 */
15 /*
16 ********************************************************************************
17 * MODULE INCLUDE FILE AND VERSION ID
18 ********************************************************************************
19 */
20 #include "namespace.h"
21 #include "az_lsp.h"
22 const char az_lsp_id[] = "@(#)$Id $" az_lsp_h;
23 /*
24 ********************************************************************************
25 * INCLUDE FILES
26 ********************************************************************************
27 */
28 #include "typedef.h"
29 #include "basic_op.h"
30 #include "oper_32b.h"
31 #include "no_count.h"
32 #include "cnst.h"
33
34 /*
35 ********************************************************************************
36 * LOCAL VARIABLES AND TABLES
37 ********************************************************************************
38 */
39 #include "grid.tab"
40 #define NC M/2 /* M = LPC order, NC = M/2 */
41
42 /*
43 ********************************************************************************
44 * LOCAL PROGRAM CODE
45 ********************************************************************************
46 */
47 /*
48 **************************************************************************
49 *
50 * Function : Chebps
51 * Purpose : Evaluates the Chebyshev polynomial series
52 * Description : - The polynomial order is n = m/2 = 5
53 * - The polynomial F(z) (F1(z) or F2(z)) is given by
54 * F(w) = 2 exp(-j5w) C(x)
55 * where
56 * C(x) = T_n(x) + f(1)T_n-1(x) + ... +f(n-1)T_1(x) + f(n)/2
57 * and T_m(x) = cos(mw) is the mth order Chebyshev
58 * polynomial ( x=cos(w) )
59 * Returns : C(x) for the input x.
60 *
61 **************************************************************************
62 */
63 static Word16 Chebps (Word16 x,
64 Word16 f[], /* (n) */
65 Word16 n)
66 {
67 Word16 i, cheb;
68 Word16 b0_h, b0_l, b1_h, b1_l, b2_h, b2_l;
69 Word32 t0;
70
71 b2_h = 256; move16 (); /* b2 = 1.0 */
72 b2_l = 0; move16 ();
73
74 t0 = L_mult (x, 512); /* 2*x */
75 t0 = L_mac (t0, f[1], 8192); /* + f[1] */
76 L_Extract (t0, &b1_h, &b1_l); /* b1 = 2*x + f[1] */
77
78 for (i = 2; i < n; i++)
79 {
80 t0 = Mpy_32_16 (b1_h, b1_l, x); /* t0 = 2.0*x*b1 */
81 t0 = L_shl (t0, 1);
82 t0 = L_mac (t0, b2_h, (Word16) 0x8000); /* t0 = 2.0*x*b1 - b2 */
83 t0 = L_msu (t0, b2_l, 1);
84 t0 = L_mac (t0, f[i], 8192); /* t0 = 2.0*x*b1 - b2 + f[i] */
85
86 L_Extract (t0, &b0_h, &b0_l); /* b0 = 2.0*x*b1 - b2 + f[i]*/
87
88 b2_l = b1_l; move16 (); /* b2 = b1; */
89 b2_h = b1_h; move16 ();
90 b1_l = b0_l; move16 (); /* b1 = b0; */
91 b1_h = b0_h; move16 ();
92 }
93
94 t0 = Mpy_32_16 (b1_h, b1_l, x); /* t0 = x*b1; */
95 t0 = L_mac (t0, b2_h, (Word16) 0x8000); /* t0 = x*b1 - b2 */
96 t0 = L_msu (t0, b2_l, 1);
97 t0 = L_mac (t0, f[i], 4096); /* t0 = x*b1 - b2 + f[i]/2 */
98
99 t0 = L_shl (t0, 6);
100
101 cheb = extract_h (t0);
102
103 return (cheb);
104 }
105
106 /*
107 ********************************************************************************
108 * PUBLIC PROGRAM CODE
109 ********************************************************************************
110 */
111 /*
112 **************************************************************************
113 *
114 * Function : Az_lsp
115 * Purpose : Compute the LSPs from the LP coefficients
116 *
117 **************************************************************************
118 */
119 void Az_lsp (
120 Word16 a[], /* (i) : predictor coefficients (MP1) */
121 Word16 lsp[], /* (o) : line spectral pairs (M) */
122 Word16 old_lsp[] /* (i) : old lsp[] (in case not found 10 roots) (M) */
123 )
124 {
125 Word16 i, j, nf, ip;
126 Word16 xlow, ylow, xhigh, yhigh, xmid, ymid, xint;
127 Word16 x, y, sign, exp;
128 Word16 *coef;
129 Word16 f1[M / 2 + 1], f2[M / 2 + 1];
130 Word32 t0;
131
132 /*-------------------------------------------------------------*
133 * find the sum and diff. pol. F1(z) and F2(z) *
134 * F1(z) <--- F1(z)/(1+z**-1) & F2(z) <--- F2(z)/(1-z**-1) *
135 * *
136 * f1[0] = 1.0; *
137 * f2[0] = 1.0; *
138 * *
139 * for (i = 0; i< NC; i++) *
140 * { *
141 * f1[i+1] = a[i+1] + a[M-i] - f1[i] ; *
142 * f2[i+1] = a[i+1] - a[M-i] + f2[i] ; *
143 * } *
144 *-------------------------------------------------------------*/
145
146 f1[0] = 1024; move16 (); /* f1[0] = 1.0 */
147 f2[0] = 1024; move16 (); /* f2[0] = 1.0 */
148
149 for (i = 0; i < NC; i++)
150 {
151 t0 = L_mult (a[i + 1], 8192); /* x = (a[i+1] + a[M-i]) >> 2 */
152 t0 = L_mac (t0, a[M - i], 8192);
153 x = extract_h (t0);
154 /* f1[i+1] = a[i+1] + a[M-i] - f1[i] */
155 f1[i + 1] = sub (x, f1[i]);move16 ();
156
157 t0 = L_mult (a[i + 1], 8192); /* x = (a[i+1] - a[M-i]) >> 2 */
158 t0 = L_msu (t0, a[M - i], 8192);
159 x = extract_h (t0);
160 /* f2[i+1] = a[i+1] - a[M-i] + f2[i] */
161 f2[i + 1] = add (x, f2[i]);move16 ();
162 }
163
164 /*-------------------------------------------------------------*
165 * find the LSPs using the Chebychev pol. evaluation *
166 *-------------------------------------------------------------*/
167
168 nf = 0; move16 (); /* number of found frequencies */
169 ip = 0; move16 (); /* indicator for f1 or f2 */
170
171 coef = f1; move16 ();
172
173 xlow = grid[0]; move16 ();
174 ylow = Chebps (xlow, coef, NC);move16 ();
175
176 j = 0;
177 test (); test ();
178 /* while ( (nf < M) && (j < grid_points) ) */
179 while ((sub (nf, M) < 0) && (sub (j, grid_points) < 0))
180 {
181 j++;
182 xhigh = xlow; move16 ();
183 yhigh = ylow; move16 ();
184 xlow = grid[j]; move16 ();
185 ylow = Chebps (xlow, coef, NC);
186 move16 ();
187
188 test ();
189 if (L_mult (ylow, yhigh) <= (Word32) 0L)
190 {
191
192 /* divide 4 times the interval */
193
194 for (i = 0; i < 4; i++)
195 {
196 /* xmid = (xlow + xhigh)/2 */
197 xmid = add (shr (xlow, 1), shr (xhigh, 1));
198 ymid = Chebps (xmid, coef, NC);
199 move16 ();
200
201 test ();
202 if (L_mult (ylow, ymid) <= (Word32) 0L)
203 {
204 yhigh = ymid; move16 ();
205 xhigh = xmid; move16 ();
206 }
207 else
208 {
209 ylow = ymid; move16 ();
210 xlow = xmid; move16 ();
211 }
212 }
213
214 /*-------------------------------------------------------------*
215 * Linear interpolation *
216 * xint = xlow - ylow*(xhigh-xlow)/(yhigh-ylow); *
217 *-------------------------------------------------------------*/
218
219 x = sub (xhigh, xlow);
220 y = sub (yhigh, ylow);
221
222 test ();
223 if (y == 0)
224 {
225 xint = xlow; move16 ();
226 }
227 else
228 {
229 sign = y; move16 ();
230 y = abs_s (y);
231 exp = norm_s (y);
232 y = shl (y, exp);
233 y = div_s ((Word16) 16383, y);
234 t0 = L_mult (x, y);
235 t0 = L_shr (t0, sub (20, exp));
236 y = extract_l (t0); /* y= (xhigh-xlow)/(yhigh-ylow) */
237
238 test ();
239 if (sign < 0)
240 y = negate (y);
241
242 t0 = L_mult (ylow, y);
243 t0 = L_shr (t0, 11);
244 xint = sub (xlow, extract_l (t0)); /* xint = xlow - ylow*y */
245 }
246
247 lsp[nf] = xint; move16 ();
248 xlow = xint; move16 ();
249 nf++;
250
251 test ();
252 if (ip == 0)
253 {
254 ip = 1; move16 ();
255 coef = f2; move16 ();
256 }
257 else
258 {
259 ip = 0; move16 ();
260 coef = f1; move16 ();
261 }
262 ylow = Chebps (xlow, coef, NC);
263 move16 ();
264
265 }
266 test (); test ();
267 }
268
269 /* Check if M roots found */
270
271 test ();
272 if (sub (nf, M) < 0)
273 {
274 for (i = 0; i < M; i++)
275 {
276 lsp[i] = old_lsp[i]; move16 ();
277 }
278
279 }
280 return;
281 }
282