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