FreeCalypso > hg > gsm-codec-lib

comparison libgsmefr/az_lsp.c @ 53:49dd1ac8e75b

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
libgsmefr: import most *.c files from ETSI source
author Mychaela Falconia <falcon@freecalypso.org>
date Fri, 25 Nov 2022 16:18:21 +0000
parents
children 902bc4b64cc6
comparison
equal deleted inserted replaced
52:988fd7ff514f 53:49dd1ac8e75b
1 /***********************************************************************
2 *
3 * FUNCTION: Az_lsp
4 *
5 * PURPOSE: Compute the LSPs from the LP coefficients (order=10)
6 *
7 * DESCRIPTION:
8 * - The sum and difference filters are computed and divided by
9 * 1+z^{-1} and 1-z^{-1}, respectively.
10 *
11 * f1[i] = a[i] + a[11-i] - f1[i-1] ; i=1,...,5
12 * f2[i] = a[i] - a[11-i] + f2[i-1] ; i=1,...,5
13 *
14 * - The roots of F1(z) and F2(z) are found using Chebyshev polynomial
15 * evaluation. The polynomials are evaluated at 60 points regularly
16 * spaced in the frequency domain. The sign change interval is
17 * subdivided 4 times to better track the root.
18 * The LSPs are found in the cosine domain [1,-1].
19 *
20 * - If less than 10 roots are found, the LSPs from the past frame are
21 * used.
22 *
23 ***********************************************************************/
24
25 #include "typedef.h"
26 #include "basic_op.h"
27 #include "oper_32b.h"
28 #include "count.h"
29 #include "cnst.h"
30
31 #include "grid.tab"
32
33 /* M = LPC order, NC = M/2 */
34
35 #define NC M/2
36
37 /* local function */
38
39 static Word16 Chebps (Word16 x, Word16 f[], Word16 n);
40
41 void Az_lsp (
42 Word16 a[], /* (i) : predictor coefficients */
43 Word16 lsp[], /* (o) : line spectral pairs */
44 Word16 old_lsp[] /* (i) : old lsp[] (in case not found 10 roots) */
45 )
46 {
47 Word16 i, j, nf, ip;
48 Word16 xlow, ylow, xhigh, yhigh, xmid, ymid, xint;
49 Word16 x, y, sign, exp;
50 Word16 *coef;
51 Word16 f1[M / 2 + 1], f2[M / 2 + 1];
52 Word32 t0;
53
54 /*-------------------------------------------------------------*
55 * find the sum and diff. pol. F1(z) and F2(z) *
56 * F1(z) <--- F1(z)/(1+z**-1) & F2(z) <--- F2(z)/(1-z**-1) *
57 * *
58 * f1[0] = 1.0; *
59 * f2[0] = 1.0; *
60 * *
61 * for (i = 0; i< NC; i++) *
62 * { *
63 * f1[i+1] = a[i+1] + a[M-i] - f1[i] ; *
64 * f2[i+1] = a[i+1] - a[M-i] + f2[i] ; *
65 * } *
66 *-------------------------------------------------------------*/
67
68 f1[0] = 1024; move16 (); /* f1[0] = 1.0 */
69 f2[0] = 1024; move16 (); /* f2[0] = 1.0 */
70
71 for (i = 0; i < NC; i++)
72 {
73 t0 = L_mult (a[i + 1], 8192); /* x = (a[i+1] + a[M-i]) >> 2 */
74 t0 = L_mac (t0, a[M - i], 8192);
75 x = extract_h (t0);
76 /* f1[i+1] = a[i+1] + a[M-i] - f1[i] */
77 f1[i + 1] = sub (x, f1[i]);move16 ();
78
79 t0 = L_mult (a[i + 1], 8192); /* x = (a[i+1] - a[M-i]) >> 2 */
80 t0 = L_msu (t0, a[M - i], 8192);
81 x = extract_h (t0);
82 /* f2[i+1] = a[i+1] - a[M-i] + f2[i] */
83 f2[i + 1] = add (x, f2[i]);move16 ();
84 }
85
86 /*-------------------------------------------------------------*
87 * find the LSPs using the Chebychev pol. evaluation *
88 *-------------------------------------------------------------*/
89
90 nf = 0; move16 (); /* number of found frequencies */
91 ip = 0; move16 (); /* indicator for f1 or f2 */
92
93 coef = f1; move16 ();
94
95 xlow = grid[0]; move16 ();
96 ylow = Chebps (xlow, coef, NC);move16 ();
97
98 j = 0;
99 test (); test ();
100 /* while ( (nf < M) && (j < grid_points) ) */
101 while ((sub (nf, M) < 0) && (sub (j, grid_points) < 0))
102 {
103 j++;
104 xhigh = xlow; move16 ();
105 yhigh = ylow; move16 ();
106 xlow = grid[j]; move16 ();
107 ylow = Chebps (xlow, coef, NC);
108 move16 ();
109
110 test ();
111 if (L_mult (ylow, yhigh) <= (Word32) 0L)
112 {
113
114 /* divide 4 times the interval */
115
116 for (i = 0; i < 4; i++)
117 {
118 /* xmid = (xlow + xhigh)/2 */
119 xmid = add (shr (xlow, 1), shr (xhigh, 1));
120 ymid = Chebps (xmid, coef, NC);
121 move16 ();
122
123 test ();
124 if (L_mult (ylow, ymid) <= (Word32) 0L)
125 {
126 yhigh = ymid; move16 ();
127 xhigh = xmid; move16 ();
128 }
129 else
130 {
131 ylow = ymid; move16 ();
132 xlow = xmid; move16 ();
133 }
134 }
135
136 /*-------------------------------------------------------------*
137 * Linear interpolation *
138 * xint = xlow - ylow*(xhigh-xlow)/(yhigh-ylow); *
139 *-------------------------------------------------------------*/
140
141 x = sub (xhigh, xlow);
142 y = sub (yhigh, ylow);
143
144 test ();
145 if (y == 0)
146 {
147 xint = xlow; move16 ();
148 }
149 else
150 {
151 sign = y; move16 ();
152 y = abs_s (y);
153 exp = norm_s (y);
154 y = shl (y, exp);
155 y = div_s ((Word16) 16383, y);
156 t0 = L_mult (x, y);
157 t0 = L_shr (t0, sub (20, exp));
158 y = extract_l (t0); /* y= (xhigh-xlow)/(yhigh-ylow) */
159
160 test ();
161 if (sign < 0)
162 y = negate (y);
163
164 t0 = L_mult (ylow, y);
165 t0 = L_shr (t0, 11);
166 xint = sub (xlow, extract_l (t0)); /* xint = xlow - ylow*y */
167 }
168
169 lsp[nf] = xint; move16 ();
170 xlow = xint; move16 ();
171 nf++;
172
173 test ();
174 if (ip == 0)
175 {
176 ip = 1; move16 ();
177 coef = f2; move16 ();
178 }
179 else
180 {
181 ip = 0; move16 ();
182 coef = f1; move16 ();
183 }
184 ylow = Chebps (xlow, coef, NC);
185 move16 ();
186
187 }
188 test (); test ();
189 }
190
191 /* Check if M roots found */
192
193 test ();
194 if (sub (nf, M) < 0)
195 {
196 for (i = 0; i < M; i++)
197 {
198 lsp[i] = old_lsp[i]; move16 ();
199 }
200
201 }
202 return;
203 }
204
205 /************************************************************************
206 *
207 * FUNCTION: Chebps
208 *
209 * PURPOSE: Evaluates the Chebyshev polynomial series
210 *
211 * DESCRIPTION:
212 * - The polynomial order is n = m/2 = 5
213 * - The polynomial F(z) (F1(z) or F2(z)) is given by
214 * F(w) = 2 exp(-j5w) C(x)
215 * where
216 * C(x) = T_n(x) + f(1)T_n-1(x) + ... +f(n-1)T_1(x) + f(n)/2
217 * and T_m(x) = cos(mw) is the mth order Chebyshev polynomial ( x=cos(w) )
218 * - The function returns the value of C(x) for the input x.
219 *
220 ***********************************************************************/
221
222 static Word16 Chebps (Word16 x, Word16 f[], Word16 n)
223 {
224 Word16 i, cheb;
225 Word16 b0_h, b0_l, b1_h, b1_l, b2_h, b2_l;
226 Word32 t0;
227
228 b2_h = 256; move16 (); /* b2 = 1.0 */
229 b2_l = 0; move16 ();
230
231 t0 = L_mult (x, 512); /* 2*x */
232 t0 = L_mac (t0, f[1], 8192); /* + f[1] */
233 L_Extract (t0, &b1_h, &b1_l); /* b1 = 2*x + f[1] */
234
235 for (i = 2; i < n; i++)
236 {
237 t0 = Mpy_32_16 (b1_h, b1_l, x); /* t0 = 2.0*x*b1 */
238 t0 = L_shl (t0, 1);
239 t0 = L_mac (t0, b2_h, (Word16) 0x8000); /* t0 = 2.0*x*b1 - b2 */
240 t0 = L_msu (t0, b2_l, 1);
241 t0 = L_mac (t0, f[i], 8192); /* t0 = 2.0*x*b1 - b2 + f[i] */
242
243 L_Extract (t0, &b0_h, &b0_l); /* b0 = 2.0*x*b1 - b2 + f[i]*/
244
245 b2_l = b1_l; move16 (); /* b2 = b1; */
246 b2_h = b1_h; move16 ();
247 b1_l = b0_l; move16 (); /* b1 = b0; */
248 b1_h = b0_h; move16 ();
249 }
250
251 t0 = Mpy_32_16 (b1_h, b1_l, x); /* t0 = x*b1; */
252 t0 = L_mac (t0, b2_h, (Word16) 0x8000); /* t0 = x*b1 - b2 */
253 t0 = L_msu (t0, b2_l, 1);
254 t0 = L_mac (t0, f[i], 4096); /* t0 = x*b1 - b2 + f[i]/2 */
255
256 t0 = L_shl (t0, 6);
257
258 cheb = extract_h (t0);
259
260 return (cheb);
261 }