1 /*
2 * Copyright (c) 2011 The WebRTC project authors. All Rights Reserved.
3 *
4 * Use of this source code is governed by a BSD-style license
5 * that can be found in the LICENSE file in the root of the source
6 * tree. An additional intellectual property rights grant can be found
7 * in the file PATENTS. All contributing project authors may
8 * be found in the AUTHORS file in the root of the source tree.
9 */
10
11
12 /*
13 * This file contains the function WebRtcSpl_LevinsonDurbin().
14 * The description header can be found in signal_processing_library.h
15 *
16 */
17
18 #include "signal_processing_library.h"
19
20 #define SPL_LEVINSON_MAXORDER 20
21
WebRtcSpl_LevinsonDurbin(WebRtc_Word32 * R,WebRtc_Word16 * A,WebRtc_Word16 * K,WebRtc_Word16 order)22 WebRtc_Word16 WebRtcSpl_LevinsonDurbin(WebRtc_Word32 *R, WebRtc_Word16 *A, WebRtc_Word16 *K,
23 WebRtc_Word16 order)
24 {
25 WebRtc_Word16 i, j;
26 // Auto-correlation coefficients in high precision
27 WebRtc_Word16 R_hi[SPL_LEVINSON_MAXORDER + 1], R_low[SPL_LEVINSON_MAXORDER + 1];
28 // LPC coefficients in high precision
29 WebRtc_Word16 A_hi[SPL_LEVINSON_MAXORDER + 1], A_low[SPL_LEVINSON_MAXORDER + 1];
30 // LPC coefficients for next iteration
31 WebRtc_Word16 A_upd_hi[SPL_LEVINSON_MAXORDER + 1], A_upd_low[SPL_LEVINSON_MAXORDER + 1];
32 // Reflection coefficient in high precision
33 WebRtc_Word16 K_hi, K_low;
34 // Prediction gain Alpha in high precision and with scale factor
35 WebRtc_Word16 Alpha_hi, Alpha_low, Alpha_exp;
36 WebRtc_Word16 tmp_hi, tmp_low;
37 WebRtc_Word32 temp1W32, temp2W32, temp3W32;
38 WebRtc_Word16 norm;
39
40 // Normalize the autocorrelation R[0]...R[order+1]
41
42 norm = WebRtcSpl_NormW32(R[0]);
43
44 for (i = order; i >= 0; i--)
45 {
46 temp1W32 = WEBRTC_SPL_LSHIFT_W32(R[i], norm);
47 // Put R in hi and low format
48 R_hi[i] = (WebRtc_Word16)WEBRTC_SPL_RSHIFT_W32(temp1W32, 16);
49 R_low[i] = (WebRtc_Word16)WEBRTC_SPL_RSHIFT_W32((temp1W32
50 - WEBRTC_SPL_LSHIFT_W32((WebRtc_Word32)R_hi[i], 16)), 1);
51 }
52
53 // K = A[1] = -R[1] / R[0]
54
55 temp2W32 = WEBRTC_SPL_LSHIFT_W32((WebRtc_Word32)R_hi[1],16)
56 + WEBRTC_SPL_LSHIFT_W32((WebRtc_Word32)R_low[1],1); // R[1] in Q31
57 temp3W32 = WEBRTC_SPL_ABS_W32(temp2W32); // abs R[1]
58 temp1W32 = WebRtcSpl_DivW32HiLow(temp3W32, R_hi[0], R_low[0]); // abs(R[1])/R[0] in Q31
59 // Put back the sign on R[1]
60 if (temp2W32 > 0)
61 {
62 temp1W32 = -temp1W32;
63 }
64
65 // Put K in hi and low format
66 K_hi = (WebRtc_Word16)WEBRTC_SPL_RSHIFT_W32(temp1W32, 16);
67 K_low = (WebRtc_Word16)WEBRTC_SPL_RSHIFT_W32((temp1W32
68 - WEBRTC_SPL_LSHIFT_W32((WebRtc_Word32)K_hi, 16)), 1);
69
70 // Store first reflection coefficient
71 K[0] = K_hi;
72
73 temp1W32 = WEBRTC_SPL_RSHIFT_W32(temp1W32, 4); // A[1] in Q27
74
75 // Put A[1] in hi and low format
76 A_hi[1] = (WebRtc_Word16)WEBRTC_SPL_RSHIFT_W32(temp1W32, 16);
77 A_low[1] = (WebRtc_Word16)WEBRTC_SPL_RSHIFT_W32((temp1W32
78 - WEBRTC_SPL_LSHIFT_W32((WebRtc_Word32)A_hi[1], 16)), 1);
79
80 // Alpha = R[0] * (1-K^2)
81
82 temp1W32 = (((WEBRTC_SPL_MUL_16_16(K_hi, K_low) >> 14) + WEBRTC_SPL_MUL_16_16(K_hi, K_hi))
83 << 1); // temp1W32 = k^2 in Q31
84
85 temp1W32 = WEBRTC_SPL_ABS_W32(temp1W32); // Guard against <0
86 temp1W32 = (WebRtc_Word32)0x7fffffffL - temp1W32; // temp1W32 = (1 - K[0]*K[0]) in Q31
87
88 // Store temp1W32 = 1 - K[0]*K[0] on hi and low format
89 tmp_hi = (WebRtc_Word16)WEBRTC_SPL_RSHIFT_W32(temp1W32, 16);
90 tmp_low = (WebRtc_Word16)WEBRTC_SPL_RSHIFT_W32((temp1W32
91 - WEBRTC_SPL_LSHIFT_W32((WebRtc_Word32)tmp_hi, 16)), 1);
92
93 // Calculate Alpha in Q31
94 temp1W32 = ((WEBRTC_SPL_MUL_16_16(R_hi[0], tmp_hi)
95 + (WEBRTC_SPL_MUL_16_16(R_hi[0], tmp_low) >> 15)
96 + (WEBRTC_SPL_MUL_16_16(R_low[0], tmp_hi) >> 15)) << 1);
97
98 // Normalize Alpha and put it in hi and low format
99
100 Alpha_exp = WebRtcSpl_NormW32(temp1W32);
101 temp1W32 = WEBRTC_SPL_LSHIFT_W32(temp1W32, Alpha_exp);
102 Alpha_hi = (WebRtc_Word16)WEBRTC_SPL_RSHIFT_W32(temp1W32, 16);
103 Alpha_low = (WebRtc_Word16)WEBRTC_SPL_RSHIFT_W32((temp1W32
104 - WEBRTC_SPL_LSHIFT_W32((WebRtc_Word32)Alpha_hi, 16)), 1);
105
106 // Perform the iterative calculations in the Levinson-Durbin algorithm
107
108 for (i = 2; i <= order; i++)
109 {
110 /* ----
111 temp1W32 = R[i] + > R[j]*A[i-j]
112 /
113 ----
114 j=1..i-1
115 */
116
117 temp1W32 = 0;
118
119 for (j = 1; j < i; j++)
120 {
121 // temp1W32 is in Q31
122 temp1W32 += ((WEBRTC_SPL_MUL_16_16(R_hi[j], A_hi[i-j]) << 1)
123 + (((WEBRTC_SPL_MUL_16_16(R_hi[j], A_low[i-j]) >> 15)
124 + (WEBRTC_SPL_MUL_16_16(R_low[j], A_hi[i-j]) >> 15)) << 1));
125 }
126
127 temp1W32 = WEBRTC_SPL_LSHIFT_W32(temp1W32, 4);
128 temp1W32 += (WEBRTC_SPL_LSHIFT_W32((WebRtc_Word32)R_hi[i], 16)
129 + WEBRTC_SPL_LSHIFT_W32((WebRtc_Word32)R_low[i], 1));
130
131 // K = -temp1W32 / Alpha
132 temp2W32 = WEBRTC_SPL_ABS_W32(temp1W32); // abs(temp1W32)
133 temp3W32 = WebRtcSpl_DivW32HiLow(temp2W32, Alpha_hi, Alpha_low); // abs(temp1W32)/Alpha
134
135 // Put the sign of temp1W32 back again
136 if (temp1W32 > 0)
137 {
138 temp3W32 = -temp3W32;
139 }
140
141 // Use the Alpha shifts from earlier to de-normalize
142 norm = WebRtcSpl_NormW32(temp3W32);
143 if ((Alpha_exp <= norm) || (temp3W32 == 0))
144 {
145 temp3W32 = WEBRTC_SPL_LSHIFT_W32(temp3W32, Alpha_exp);
146 } else
147 {
148 if (temp3W32 > 0)
149 {
150 temp3W32 = (WebRtc_Word32)0x7fffffffL;
151 } else
152 {
153 temp3W32 = (WebRtc_Word32)0x80000000L;
154 }
155 }
156
157 // Put K on hi and low format
158 K_hi = (WebRtc_Word16)WEBRTC_SPL_RSHIFT_W32(temp3W32, 16);
159 K_low = (WebRtc_Word16)WEBRTC_SPL_RSHIFT_W32((temp3W32
160 - WEBRTC_SPL_LSHIFT_W32((WebRtc_Word32)K_hi, 16)), 1);
161
162 // Store Reflection coefficient in Q15
163 K[i - 1] = K_hi;
164
165 // Test for unstable filter.
166 // If unstable return 0 and let the user decide what to do in that case
167
168 if ((WebRtc_Word32)WEBRTC_SPL_ABS_W16(K_hi) > (WebRtc_Word32)32750)
169 {
170 return 0; // Unstable filter
171 }
172
173 /*
174 Compute updated LPC coefficient: Anew[i]
175 Anew[j]= A[j] + K*A[i-j] for j=1..i-1
176 Anew[i]= K
177 */
178
179 for (j = 1; j < i; j++)
180 {
181 // temp1W32 = A[j] in Q27
182 temp1W32 = WEBRTC_SPL_LSHIFT_W32((WebRtc_Word32)A_hi[j],16)
183 + WEBRTC_SPL_LSHIFT_W32((WebRtc_Word32)A_low[j],1);
184
185 // temp1W32 += K*A[i-j] in Q27
186 temp1W32 += ((WEBRTC_SPL_MUL_16_16(K_hi, A_hi[i-j])
187 + (WEBRTC_SPL_MUL_16_16(K_hi, A_low[i-j]) >> 15)
188 + (WEBRTC_SPL_MUL_16_16(K_low, A_hi[i-j]) >> 15)) << 1);
189
190 // Put Anew in hi and low format
191 A_upd_hi[j] = (WebRtc_Word16)WEBRTC_SPL_RSHIFT_W32(temp1W32, 16);
192 A_upd_low[j] = (WebRtc_Word16)WEBRTC_SPL_RSHIFT_W32((temp1W32
193 - WEBRTC_SPL_LSHIFT_W32((WebRtc_Word32)A_upd_hi[j], 16)), 1);
194 }
195
196 // temp3W32 = K in Q27 (Convert from Q31 to Q27)
197 temp3W32 = WEBRTC_SPL_RSHIFT_W32(temp3W32, 4);
198
199 // Store Anew in hi and low format
200 A_upd_hi[i] = (WebRtc_Word16)WEBRTC_SPL_RSHIFT_W32(temp3W32, 16);
201 A_upd_low[i] = (WebRtc_Word16)WEBRTC_SPL_RSHIFT_W32((temp3W32
202 - WEBRTC_SPL_LSHIFT_W32((WebRtc_Word32)A_upd_hi[i], 16)), 1);
203
204 // Alpha = Alpha * (1-K^2)
205
206 temp1W32 = (((WEBRTC_SPL_MUL_16_16(K_hi, K_low) >> 14)
207 + WEBRTC_SPL_MUL_16_16(K_hi, K_hi)) << 1); // K*K in Q31
208
209 temp1W32 = WEBRTC_SPL_ABS_W32(temp1W32); // Guard against <0
210 temp1W32 = (WebRtc_Word32)0x7fffffffL - temp1W32; // 1 - K*K in Q31
211
212 // Convert 1- K^2 in hi and low format
213 tmp_hi = (WebRtc_Word16)WEBRTC_SPL_RSHIFT_W32(temp1W32, 16);
214 tmp_low = (WebRtc_Word16)WEBRTC_SPL_RSHIFT_W32((temp1W32
215 - WEBRTC_SPL_LSHIFT_W32((WebRtc_Word32)tmp_hi, 16)), 1);
216
217 // Calculate Alpha = Alpha * (1-K^2) in Q31
218 temp1W32 = ((WEBRTC_SPL_MUL_16_16(Alpha_hi, tmp_hi)
219 + (WEBRTC_SPL_MUL_16_16(Alpha_hi, tmp_low) >> 15)
220 + (WEBRTC_SPL_MUL_16_16(Alpha_low, tmp_hi) >> 15)) << 1);
221
222 // Normalize Alpha and store it on hi and low format
223
224 norm = WebRtcSpl_NormW32(temp1W32);
225 temp1W32 = WEBRTC_SPL_LSHIFT_W32(temp1W32, norm);
226
227 Alpha_hi = (WebRtc_Word16)WEBRTC_SPL_RSHIFT_W32(temp1W32, 16);
228 Alpha_low = (WebRtc_Word16)WEBRTC_SPL_RSHIFT_W32((temp1W32
229 - WEBRTC_SPL_LSHIFT_W32((WebRtc_Word32)Alpha_hi, 16)), 1);
230
231 // Update the total normalization of Alpha
232 Alpha_exp = Alpha_exp + norm;
233
234 // Update A[]
235
236 for (j = 1; j <= i; j++)
237 {
238 A_hi[j] = A_upd_hi[j];
239 A_low[j] = A_upd_low[j];
240 }
241 }
242
243 /*
244 Set A[0] to 1.0 and store the A[i] i=1...order in Q12
245 (Convert from Q27 and use rounding)
246 */
247
248 A[0] = 4096;
249
250 for (i = 1; i <= order; i++)
251 {
252 // temp1W32 in Q27
253 temp1W32 = WEBRTC_SPL_LSHIFT_W32((WebRtc_Word32)A_hi[i], 16)
254 + WEBRTC_SPL_LSHIFT_W32((WebRtc_Word32)A_low[i], 1);
255 // Round and store upper word
256 A[i] = (WebRtc_Word16)WEBRTC_SPL_RSHIFT_W32((temp1W32<<1)+(WebRtc_Word32)32768, 16);
257 }
258 return 1; // Stable filters
259 }
260