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