FreeCalypso > hg > gsm-codec-lib

comparison libgsmefr/levinson.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 0e41ca9ebf45
comparison
equal deleted inserted replaced
52:988fd7ff514f 53:49dd1ac8e75b
1 /*************************************************************************
2 *
3 * FUNCTION: Levinson()
4 *
5 * PURPOSE: Levinson-Durbin algorithm in double precision. To compute the
6 * LP filter parameters from the speech autocorrelations.
7 *
8 * DESCRIPTION:
9 * R[i] autocorrelations.
10 * A[i] filter coefficients.
11 * K reflection coefficients.
12 * Alpha prediction gain.
13 *
14 * Initialisation:
15 * A[0] = 1
16 * K = -R[1]/R[0]
17 * A[1] = K
18 * Alpha = R[0] * (1-K**2]
19 *
20 * Do for i = 2 to M
21 *
22 * S = SUM ( R[j]*A[i-j] ,j=1,i-1 ) + R[i]
23 *
24 * K = -S / Alpha
25 *
26 * An[j] = A[j] + K*A[i-j] for j=1 to i-1
27 * where An[i] = new A[i]
28 * An[i]=K
29 *
30 * Alpha=Alpha * (1-K**2)
31 *
32 * END
33 *
34 *************************************************************************/
35
36 #include "typedef.h"
37 #include "basic_op.h"
38 #include "oper_32b.h"
39 #include "count.h"
40
41 /* Lpc order == 10 */
42
43 #define M 10
44
45 /* Last A(z) for case of unstable filter */
46
47 Word16 old_A[M + 1];
48
49 void Levinson (
50 Word16 Rh[], /* (i) : Rh[m+1] Vector of autocorrelations (msb) */
51 Word16 Rl[], /* (i) : Rl[m+1] Vector of autocorrelations (lsb) */
52 Word16 A[], /* (o) : A[m] LPC coefficients (m = 10) */
53 Word16 rc[] /* (o) : rc[4] First 4 reflection coefficients */
54
55 )
56 {
57 Word16 i, j;
58 Word16 hi, lo;
59 Word16 Kh, Kl; /* reflexion coefficient; hi and lo */
60 Word16 alp_h, alp_l, alp_exp; /* Prediction gain; hi lo and exponent */
61 Word16 Ah[M + 1], Al[M + 1]; /* LPC coef. in double prec. */
62 Word16 Anh[M + 1], Anl[M + 1];/* LPC coef.for next iteration in double
63 prec. */
64 Word32 t0, t1, t2; /* temporary variable */
65
66 /* K = A[1] = -R[1] / R[0] */
67
68 t1 = L_Comp (Rh[1], Rl[1]);
69 t2 = L_abs (t1); /* abs R[1] */
70 t0 = Div_32 (t2, Rh[0], Rl[0]); /* R[1]/R[0] */
71 test ();
72 if (t1 > 0)
73 t0 = L_negate (t0); /* -R[1]/R[0] */
74 L_Extract (t0, &Kh, &Kl); /* K in DPF */
75
76 rc[0] = round (t0); move16 ();
77
78 t0 = L_shr (t0, 4); /* A[1] in */
79 L_Extract (t0, &Ah[1], &Al[1]); /* A[1] in DPF */
80
81 /* Alpha = R[0] * (1-K**2) */
82
83 t0 = Mpy_32 (Kh, Kl, Kh, Kl); /* K*K */
84 t0 = L_abs (t0); /* Some case <0 !! */
85 t0 = L_sub ((Word32) 0x7fffffffL, t0); /* 1 - K*K */
86 L_Extract (t0, &hi, &lo); /* DPF format */
87 t0 = Mpy_32 (Rh[0], Rl[0], hi, lo); /* Alpha in */
88
89 /* Normalize Alpha */
90
91 alp_exp = norm_l (t0);
92 t0 = L_shl (t0, alp_exp);
93 L_Extract (t0, &alp_h, &alp_l); /* DPF format */
94
95 /*--------------------------------------*
96 * ITERATIONS I=2 to M *
97 *--------------------------------------*/
98
99 for (i = 2; i <= M; i++)
100 {
101 /* t0 = SUM ( R[j]*A[i-j] ,j=1,i-1 ) + R[i] */
102
103 t0 = 0; move32 ();
104 for (j = 1; j < i; j++)
105 {
106 t0 = L_add (t0, Mpy_32 (Rh[j], Rl[j], Ah[i - j], Al[i - j]));
107 }
108 t0 = L_shl (t0, 4);
109
110 t1 = L_Comp (Rh[i], Rl[i]);
111 t0 = L_add (t0, t1); /* add R[i] */
112
113 /* K = -t0 / Alpha */
114
115 t1 = L_abs (t0);
116 t2 = Div_32 (t1, alp_h, alp_l); /* abs(t0)/Alpha */
117 test ();
118 if (t0 > 0)
119 t2 = L_negate (t2); /* K =-t0/Alpha */
120 t2 = L_shl (t2, alp_exp); /* denormalize; compare to Alpha */
121 L_Extract (t2, &Kh, &Kl); /* K in DPF */
122
123 test ();
124 if (sub (i, 5) < 0)
125 {
126 rc[i - 1] = round (t2); move16 ();
127 }
128 /* Test for unstable filter. If unstable keep old A(z) */
129
130 test ();
131 if (sub (abs_s (Kh), 32750) > 0)
132 {
133 for (j = 0; j <= M; j++)
134 {
135 A[j] = old_A[j]; move16 ();
136 }
137
138 for (j = 0; j < 4; j++)
139 {
140 rc[j] = 0; move16 ();
141 }
142
143 return;
144 }
145 /*------------------------------------------*
146 * Compute new LPC coeff. -> An[i] *
147 * An[j]= A[j] + K*A[i-j] , j=1 to i-1 *
148 * An[i]= K *
149 *------------------------------------------*/
150
151 for (j = 1; j < i; j++)
152 {
153 t0 = Mpy_32 (Kh, Kl, Ah[i - j], Al[i - j]);
154 t0 = L_mac (t0, Ah[j], 32767);
155 t0 = L_mac (t0, Al[j], 1);
156 L_Extract (t0, &Anh[j], &Anl[j]);
157 }
158 t2 = L_shr (t2, 4);
159 L_Extract (t2, &Anh[i], &Anl[i]);
160
161 /* Alpha = Alpha * (1-K**2) */
162
163 t0 = Mpy_32 (Kh, Kl, Kh, Kl); /* K*K */
164 t0 = L_abs (t0); /* Some case <0 !! */
165 t0 = L_sub ((Word32) 0x7fffffffL, t0); /* 1 - K*K */
166 L_Extract (t0, &hi, &lo); /* DPF format */
167 t0 = Mpy_32 (alp_h, alp_l, hi, lo);
168
169 /* Normalize Alpha */
170
171 j = norm_l (t0);
172 t0 = L_shl (t0, j);
173 L_Extract (t0, &alp_h, &alp_l); /* DPF format */
174 alp_exp = add (alp_exp, j); /* Add normalization to
175 alp_exp */
176
177 /* A[j] = An[j] */
178
179 for (j = 1; j <= i; j++)
180 {
181 Ah[j] = Anh[j]; move16 ();
182 Al[j] = Anl[j]; move16 ();
183 }
184 }
185
186 A[0] = 4096; move16 ();
187 for (i = 1; i <= M; i++)
188 {
189 t0 = L_Comp (Ah[i], Al[i]);
190 old_A[i] = A[i] = round (L_shl (t0, 1));move16 (); move16 ();
191 }
192
193 return;
194 }