1 /*
2 ** Copyright 2003-2010, VisualOn, Inc.
3 **
4 ** Licensed under the Apache License, Version 2.0 (the "License");
5 ** you may not use this file except in compliance with the License.
6 ** You may obtain a copy of the License at
7 **
8 ** http://www.apache.org/licenses/LICENSE-2.0
9 **
10 ** Unless required by applicable law or agreed to in writing, software
11 ** distributed under the License is distributed on an "AS IS" BASIS,
12 ** WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
13 ** See the License for the specific language governing permissions and
14 ** limitations under the License.
15 */
16
17 /***********************************************************************
18 * File: levinson.c *
19 * *
20 * Description:LEVINSON-DURBIN algorithm in double precision *
21 * *
22 ************************************************************************/
23 /*---------------------------------------------------------------------------*
24 * LEVINSON.C *
25 *---------------------------------------------------------------------------*
26 * *
27 * LEVINSON-DURBIN algorithm in double precision *
28 * *
29 * *
30 * Algorithm *
31 * *
32 * R[i] autocorrelations. *
33 * A[i] filter coefficients. *
34 * K reflection coefficients. *
35 * Alpha prediction gain. *
36 * *
37 * Initialization: *
38 * A[0] = 1 *
39 * K = -R[1]/R[0] *
40 * A[1] = K *
41 * Alpha = R[0] * (1-K**2] *
42 * *
43 * Do for i = 2 to M *
44 * *
45 * S = SUM ( R[j]*A[i-j] ,j=1,i-1 ) + R[i] *
46 * *
47 * K = -S / Alpha *
48 * *
49 * An[j] = A[j] + K*A[i-j] for j=1 to i-1 *
50 * where An[i] = new A[i] *
51 * An[i]=K *
52 * *
53 * Alpha=Alpha * (1-K**2) *
54 * *
55 * END *
56 * *
57 * Remarks on the dynamics of the calculations. *
58 * *
59 * The numbers used are in double precision in the following format : *
60 * A = AH <<16 + AL<<1. AH and AL are 16 bit signed integers. *
61 * Since the LSB's also contain a sign bit, this format does not *
62 * correspond to standard 32 bit integers. We use this format since *
63 * it allows fast execution of multiplications and divisions. *
64 * *
65 * "DPF" will refer to this special format in the following text. *
66 * See oper_32b.c *
67 * *
68 * The R[i] were normalized in routine AUTO (hence, R[i] < 1.0). *
69 * The K[i] and Alpha are theoretically < 1.0. *
70 * The A[i], for a sampling frequency of 8 kHz, are in practice *
71 * always inferior to 16.0. *
72 * *
73 * These characteristics allow straigthforward fixed-point *
74 * implementation. We choose to represent the parameters as *
75 * follows : *
76 * *
77 * R[i] Q31 +- .99.. *
78 * K[i] Q31 +- .99.. *
79 * Alpha Normalized -> mantissa in Q31 plus exponent *
80 * A[i] Q27 +- 15.999.. *
81 * *
82 * The additions are performed in 32 bit. For the summation used *
83 * to calculate the K[i], we multiply numbers in Q31 by numbers *
84 * in Q27, with the result of the multiplications in Q27, *
85 * resulting in a dynamic of +- 16. This is sufficient to avoid *
86 * overflow, since the final result of the summation is *
87 * necessarily < 1.0 as both the K[i] and Alpha are *
88 * theoretically < 1.0. *
89 *___________________________________________________________________________*/
90 #include "typedef.h"
91 #include "basic_op.h"
92 #include "oper_32b.h"
93 #include "acelp.h"
94
95 #define M 16
96 #define NC (M/2)
97
Init_Levinson(Word16 * mem)98 void Init_Levinson(
99 Word16 * mem /* output :static memory (18 words) */
100 )
101 {
102 Set_zero(mem, 18); /* old_A[0..M-1] = 0, old_rc[0..1] = 0 */
103 return;
104 }
105
106
Levinson(Word16 Rh[],Word16 Rl[],Word16 A[],Word16 rc[],Word16 * mem)107 void Levinson(
108 Word16 Rh[], /* (i) : Rh[M+1] Vector of autocorrelations (msb) */
109 Word16 Rl[], /* (i) : Rl[M+1] Vector of autocorrelations (lsb) */
110 Word16 A[], /* (o) Q12 : A[M] LPC coefficients (m = 16) */
111 Word16 rc[], /* (o) Q15 : rc[M] Reflection coefficients. */
112 Word16 * mem /* (i/o) :static memory (18 words) */
113 )
114 {
115 Word32 i, j;
116 Word16 hi, lo;
117 Word16 Kh, Kl; /* reflection coefficient; hi and lo */
118 Word16 alp_h, alp_l, alp_exp; /* Prediction gain; hi lo and exponent */
119 Word16 Ah[M + 1], Al[M + 1]; /* LPC coef. in double prec. */
120 Word16 Anh[M + 1], Anl[M + 1]; /* LPC coef.for next iteration in double prec. */
121 Word32 t0, t1, t2; /* temporary variable */
122 Word16 *old_A, *old_rc;
123
124 /* Last A(z) for case of unstable filter */
125 old_A = mem;
126 old_rc = mem + M;
127
128 /* K = A[1] = -R[1] / R[0] */
129
130 t1 = ((Rh[1] << 16) + (Rl[1] << 1)); /* R[1] in Q31 */
131 t2 = L_abs(t1); /* abs R[1] */
132 t0 = Div_32(t2, Rh[0], Rl[0]); /* R[1]/R[0] in Q31 */
133 if (t1 > 0)
134 t0 = -t0; /* -R[1]/R[0] */
135
136 Kh = t0 >> 16;
137 Kl = (t0 & 0xffff)>>1;
138 rc[0] = Kh;
139 t0 = (t0 >> 4); /* A[1] in Q27 */
140
141 Ah[1] = t0 >> 16;
142 Al[1] = (t0 & 0xffff)>>1;
143
144 /* Alpha = R[0] * (1-K**2) */
145 t0 = Mpy_32(Kh, Kl, Kh, Kl); /* K*K in Q31 */
146 t0 = L_abs(t0); /* Some case <0 !! */
147 t0 = vo_L_sub((Word32) 0x7fffffffL, t0); /* 1 - K*K in Q31 */
148
149 hi = t0 >> 16;
150 lo = (t0 & 0xffff)>>1;
151
152 t0 = Mpy_32(Rh[0], Rl[0], hi, lo); /* Alpha in Q31 */
153
154 /* Normalize Alpha */
155 alp_exp = norm_l(t0);
156 t0 = (t0 << alp_exp);
157
158 alp_h = t0 >> 16;
159 alp_l = (t0 & 0xffff)>>1;
160 /*--------------------------------------*
161 * ITERATIONS I=2 to M *
162 *--------------------------------------*/
163 for (i = 2; i <= M; i++)
164 {
165 /* t0 = SUM ( R[j]*A[i-j] ,j=1,i-1 ) + R[i] */
166 t0 = 0;
167 for (j = 1; j < i; j++)
168 t0 = vo_L_add(t0, Mpy_32(Rh[j], Rl[j], Ah[i - j], Al[i - j]));
169
170 t0 = t0 << 4; /* result in Q27 -> convert to Q31 */
171 /* No overflow possible */
172 t1 = ((Rh[i] << 16) + (Rl[i] << 1));
173 t0 = vo_L_add(t0, t1); /* add R[i] in Q31 */
174
175 /* K = -t0 / Alpha */
176 t1 = L_abs(t0);
177 t2 = Div_32(t1, alp_h, alp_l); /* abs(t0)/Alpha */
178 if (t0 > 0)
179 t2 = -t2; /* K =-t0/Alpha */
180 t2 = (t2 << alp_exp); /* denormalize; compare to Alpha */
181
182 Kh = t2 >> 16;
183 Kl = (t2 & 0xffff)>>1;
184
185 rc[i - 1] = Kh;
186 /* Test for unstable filter. If unstable keep old A(z) */
187 if (abs_s(Kh) > 32750)
188 {
189 A[0] = 4096; /* Ai[0] not stored (always 1.0) */
190 for (j = 0; j < M; j++)
191 {
192 A[j + 1] = old_A[j];
193 }
194 rc[0] = old_rc[0]; /* only two rc coefficients are needed */
195 rc[1] = old_rc[1];
196 return;
197 }
198 /*------------------------------------------*
199 * Compute new LPC coeff. -> An[i] *
200 * An[j]= A[j] + K*A[i-j] , j=1 to i-1 *
201 * An[i]= K *
202 *------------------------------------------*/
203 for (j = 1; j < i; j++)
204 {
205 t0 = Mpy_32(Kh, Kl, Ah[i - j], Al[i - j]);
206 t0 = vo_L_add(t0, ((Ah[j] << 16) + (Al[j] << 1)));
207 Anh[j] = t0 >> 16;
208 Anl[j] = (t0 & 0xffff)>>1;
209 }
210 t2 = (t2 >> 4); /* t2 = K in Q31 ->convert to Q27 */
211
212 VO_L_Extract(t2, &Anh[i], &Anl[i]); /* An[i] in Q27 */
213
214 /* Alpha = Alpha * (1-K**2) */
215 t0 = Mpy_32(Kh, Kl, Kh, Kl); /* K*K in Q31 */
216 t0 = L_abs(t0); /* Some case <0 !! */
217 t0 = vo_L_sub((Word32) 0x7fffffffL, t0); /* 1 - K*K in Q31 */
218 hi = t0 >> 16;
219 lo = (t0 & 0xffff)>>1;
220 t0 = Mpy_32(alp_h, alp_l, hi, lo); /* Alpha in Q31 */
221
222 /* Normalize Alpha */
223 j = norm_l(t0);
224 t0 = (t0 << j);
225 alp_h = t0 >> 16;
226 alp_l = (t0 & 0xffff)>>1;
227 alp_exp += j; /* Add normalization to alp_exp */
228
229 /* A[j] = An[j] */
230 for (j = 1; j <= i; j++)
231 {
232 Ah[j] = Anh[j];
233 Al[j] = Anl[j];
234 }
235 }
236 /* Truncate A[i] in Q27 to Q12 with rounding */
237 A[0] = 4096;
238 for (i = 1; i <= M; i++)
239 {
240 t0 = (Ah[i] << 16) + (Al[i] << 1);
241 old_A[i - 1] = A[i] = vo_round((t0 << 1));
242 }
243 old_rc[0] = rc[0];
244 old_rc[1] = rc[1];
245
246 return;
247 }
248
249
250
251