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