FreeCalypso > hg > gsm-codec-lib
comparison libtwamr/levinson.c @ 385:c713061b6edf
libtwamr: integrate levinson.c
author | Mychaela Falconia <falcon@freecalypso.org> |
---|---|
date | Mon, 06 May 2024 06:01:56 +0000 |
parents | |
children |
comparison
equal
deleted
inserted
replaced
384:a8dab7028e4d | 385:c713061b6edf |
---|---|
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 : levinson.c | |
11 * Purpose : Levinson-Durbin algorithm in double precision. | |
12 * : To compute the LP filter parameters from the | |
13 * : speech autocorrelations. | |
14 * | |
15 ***************************************************************************** | |
16 */ | |
17 | |
18 /* | |
19 ***************************************************************************** | |
20 * MODULE INCLUDE FILE AND VERSION ID | |
21 ***************************************************************************** | |
22 */ | |
23 #include "namespace.h" | |
24 #include "levinson.h" | |
25 | |
26 /* | |
27 ***************************************************************************** | |
28 * INCLUDE FILES | |
29 ***************************************************************************** | |
30 */ | |
31 #include "typedef.h" | |
32 #include "basic_op.h" | |
33 #include "oper_32b.h" | |
34 #include "no_count.h" | |
35 #include "cnst.h" | |
36 | |
37 /* | |
38 ***************************************************************************** | |
39 * LOCAL VARIABLES AND TABLES | |
40 ***************************************************************************** | |
41 */ | |
42 /*---------------------------------------------------------------* | |
43 * Constants (defined in "cnst.h") * | |
44 *---------------------------------------------------------------* | |
45 * M : LPC order | |
46 *---------------------------------------------------------------*/ | |
47 | |
48 /* | |
49 ***************************************************************************** | |
50 * PUBLIC PROGRAM CODE | |
51 ***************************************************************************** | |
52 */ | |
53 | |
54 /************************************************************************* | |
55 * | |
56 * Function: Levinson_reset | |
57 * Purpose: Initializes state memory to zero | |
58 * | |
59 ************************************************************************** | |
60 */ | |
61 void Levinson_reset (LevinsonState *state) | |
62 { | |
63 Word16 i; | |
64 | |
65 state->old_A[0] = 4096; | |
66 for(i = 1; i < M + 1; i++) | |
67 state->old_A[i] = 0; | |
68 } | |
69 | |
70 /************************************************************************* | |
71 * | |
72 * FUNCTION: Levinson() | |
73 * | |
74 * PURPOSE: Levinson-Durbin algorithm in double precision. To compute the | |
75 * LP filter parameters from the speech autocorrelations. | |
76 * | |
77 * DESCRIPTION: | |
78 * R[i] autocorrelations. | |
79 * A[i] filter coefficients. | |
80 * K reflection coefficients. | |
81 * Alpha prediction gain. | |
82 * | |
83 * Initialisation: | |
84 * A[0] = 1 | |
85 * K = -R[1]/R[0] | |
86 * A[1] = K | |
87 * Alpha = R[0] * (1-K**2] | |
88 * | |
89 * Do for i = 2 to M | |
90 * | |
91 * S = SUM ( R[j]*A[i-j] ,j=1,i-1 ) + R[i] | |
92 * | |
93 * K = -S / Alpha | |
94 * | |
95 * An[j] = A[j] + K*A[i-j] for j=1 to i-1 | |
96 * where An[i] = new A[i] | |
97 * An[i]=K | |
98 * | |
99 * Alpha=Alpha * (1-K**2) | |
100 * | |
101 * END | |
102 * | |
103 *************************************************************************/ | |
104 int Levinson ( | |
105 LevinsonState *st, | |
106 Word16 Rh[], /* i : Rh[m+1] Vector of autocorrelations (msb) */ | |
107 Word16 Rl[], /* i : Rl[m+1] Vector of autocorrelations (lsb) */ | |
108 Word16 A[], /* o : A[m] LPC coefficients (m = 10) */ | |
109 Word16 rc[] /* o : rc[4] First 4 reflection coefficients */ | |
110 ) | |
111 { | |
112 Word16 i, j; | |
113 Word16 hi, lo; | |
114 Word16 Kh, Kl; /* reflexion coefficient; hi and lo */ | |
115 Word16 alp_h, alp_l, alp_exp; /* Prediction gain; hi lo and exponent */ | |
116 Word16 Ah[M + 1], Al[M + 1]; /* LPC coef. in double prec. */ | |
117 Word16 Anh[M + 1], Anl[M + 1];/* LPC coef.for next iteration in double | |
118 prec. */ | |
119 Word32 t0, t1, t2; /* temporary variable */ | |
120 | |
121 /* K = A[1] = -R[1] / R[0] */ | |
122 | |
123 t1 = L_Comp (Rh[1], Rl[1]); | |
124 t2 = L_abs (t1); /* abs R[1] */ | |
125 t0 = Div_32 (t2, Rh[0], Rl[0]); /* R[1]/R[0] */ | |
126 test (); | |
127 if (t1 > 0) | |
128 t0 = L_negate (t0); /* -R[1]/R[0] */ | |
129 L_Extract (t0, &Kh, &Kl); /* K in DPF */ | |
130 | |
131 rc[0] = round (t0); move16 (); | |
132 | |
133 t0 = L_shr (t0, 4); /* A[1] in */ | |
134 L_Extract (t0, &Ah[1], &Al[1]); /* A[1] in DPF */ | |
135 | |
136 /* Alpha = R[0] * (1-K**2) */ | |
137 | |
138 t0 = Mpy_32 (Kh, Kl, Kh, Kl); /* K*K */ | |
139 t0 = L_abs (t0); /* Some case <0 !! */ | |
140 t0 = L_sub ((Word32) 0x7fffffffL, t0); /* 1 - K*K */ | |
141 L_Extract (t0, &hi, &lo); /* DPF format */ | |
142 t0 = Mpy_32 (Rh[0], Rl[0], hi, lo); /* Alpha in */ | |
143 | |
144 /* Normalize Alpha */ | |
145 | |
146 alp_exp = norm_l (t0); | |
147 t0 = L_shl (t0, alp_exp); | |
148 L_Extract (t0, &alp_h, &alp_l); /* DPF format */ | |
149 | |
150 /*--------------------------------------* | |
151 * ITERATIONS I=2 to M * | |
152 *--------------------------------------*/ | |
153 | |
154 for (i = 2; i <= M; i++) | |
155 { | |
156 /* t0 = SUM ( R[j]*A[i-j] ,j=1,i-1 ) + R[i] */ | |
157 | |
158 t0 = 0; move32 (); | |
159 for (j = 1; j < i; j++) | |
160 { | |
161 t0 = L_add (t0, Mpy_32 (Rh[j], Rl[j], Ah[i - j], Al[i - j])); | |
162 } | |
163 t0 = L_shl (t0, 4); | |
164 | |
165 t1 = L_Comp (Rh[i], Rl[i]); | |
166 t0 = L_add (t0, t1); /* add R[i] */ | |
167 | |
168 /* K = -t0 / Alpha */ | |
169 | |
170 t1 = L_abs (t0); | |
171 t2 = Div_32 (t1, alp_h, alp_l); /* abs(t0)/Alpha */ | |
172 test (); | |
173 if (t0 > 0) | |
174 t2 = L_negate (t2); /* K =-t0/Alpha */ | |
175 t2 = L_shl (t2, alp_exp); /* denormalize; compare to Alpha */ | |
176 L_Extract (t2, &Kh, &Kl); /* K in DPF */ | |
177 | |
178 test (); | |
179 if (sub (i, 5) < 0) | |
180 { | |
181 rc[i - 1] = round (t2); move16 (); | |
182 } | |
183 /* Test for unstable filter. If unstable keep old A(z) */ | |
184 | |
185 test (); | |
186 if (sub (abs_s (Kh), 32750) > 0) | |
187 { | |
188 for (j = 0; j <= M; j++) | |
189 { | |
190 A[j] = st->old_A[j]; move16 (); | |
191 } | |
192 | |
193 for (j = 0; j < 4; j++) | |
194 { | |
195 rc[j] = 0; move16 (); | |
196 } | |
197 | |
198 return 0; | |
199 } | |
200 /*------------------------------------------* | |
201 * Compute new LPC coeff. -> An[i] * | |
202 * An[j]= A[j] + K*A[i-j] , j=1 to i-1 * | |
203 * An[i]= K * | |
204 *------------------------------------------*/ | |
205 | |
206 for (j = 1; j < i; j++) | |
207 { | |
208 t0 = Mpy_32 (Kh, Kl, Ah[i - j], Al[i - j]); | |
209 t0 = L_add(t0, L_Comp(Ah[j], Al[j])); | |
210 L_Extract (t0, &Anh[j], &Anl[j]); | |
211 } | |
212 t2 = L_shr (t2, 4); | |
213 L_Extract (t2, &Anh[i], &Anl[i]); | |
214 | |
215 /* Alpha = Alpha * (1-K**2) */ | |
216 | |
217 t0 = Mpy_32 (Kh, Kl, Kh, Kl); /* K*K */ | |
218 t0 = L_abs (t0); /* Some case <0 !! */ | |
219 t0 = L_sub ((Word32) 0x7fffffffL, t0); /* 1 - K*K */ | |
220 L_Extract (t0, &hi, &lo); /* DPF format */ | |
221 t0 = Mpy_32 (alp_h, alp_l, hi, lo); | |
222 | |
223 /* Normalize Alpha */ | |
224 | |
225 j = norm_l (t0); | |
226 t0 = L_shl (t0, j); | |
227 L_Extract (t0, &alp_h, &alp_l); /* DPF format */ | |
228 alp_exp = add (alp_exp, j); /* Add normalization to | |
229 alp_exp */ | |
230 | |
231 /* A[j] = An[j] */ | |
232 | |
233 for (j = 1; j <= i; j++) | |
234 { | |
235 Ah[j] = Anh[j]; move16 (); | |
236 Al[j] = Anl[j]; move16 (); | |
237 } | |
238 } | |
239 | |
240 A[0] = 4096; move16 (); | |
241 for (i = 1; i <= M; i++) | |
242 { | |
243 t0 = L_Comp (Ah[i], Al[i]); | |
244 st->old_A[i] = A[i] = round (L_shl (t0, 1));move16 (); move16 (); | |
245 } | |
246 | |
247 return 0; | |
248 } |