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