FreeCalypso > hg > gsm-codec-lib

comparison libtwamr/levinson.c @ 385:c713061b6edf

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
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 }