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 22int16_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