levinson.c revision e2e838afcf03e603a41a0455846eaf9614537c16
1/* 2 ** Copyright 2003-2010, VisualOn, Inc. 3 ** 4 ** Licensed under the Apache License, Version 2.0 (the "License"); 5 ** you may not use this file except in compliance with the License. 6 ** You may obtain a copy of the License at 7 ** 8 ** http://www.apache.org/licenses/LICENSE-2.0 9 ** 10 ** Unless required by applicable law or agreed to in writing, software 11 ** distributed under the License is distributed on an "AS IS" BASIS, 12 ** WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 13 ** See the License for the specific language governing permissions and 14 ** limitations under the License. 15 */ 16 17/*********************************************************************** 18* File: levinson.c * 19* * 20* Description:LEVINSON-DURBIN algorithm in double precision * 21* * 22************************************************************************/ 23/*---------------------------------------------------------------------------* 24 * LEVINSON.C * 25 *---------------------------------------------------------------------------* 26 * * 27 * LEVINSON-DURBIN algorithm in double precision * 28 * * 29 * * 30 * Algorithm * 31 * * 32 * R[i] autocorrelations. * 33 * A[i] filter coefficients. * 34 * K reflection coefficients. * 35 * Alpha prediction gain. * 36 * * 37 * Initialization: * 38 * A[0] = 1 * 39 * K = -R[1]/R[0] * 40 * A[1] = K * 41 * Alpha = R[0] * (1-K**2] * 42 * * 43 * Do for i = 2 to M * 44 * * 45 * S = SUM ( R[j]*A[i-j] ,j=1,i-1 ) + R[i] * 46 * * 47 * K = -S / Alpha * 48 * * 49 * An[j] = A[j] + K*A[i-j] for j=1 to i-1 * 50 * where An[i] = new A[i] * 51 * An[i]=K * 52 * * 53 * Alpha=Alpha * (1-K**2) * 54 * * 55 * END * 56 * * 57 * Remarks on the dynamics of the calculations. * 58 * * 59 * The numbers used are in double precision in the following format : * 60 * A = AH <<16 + AL<<1. AH and AL are 16 bit signed integers. * 61 * Since the LSB's also contain a sign bit, this format does not * 62 * correspond to standard 32 bit integers. We use this format since * 63 * it allows fast execution of multiplications and divisions. * 64 * * 65 * "DPF" will refer to this special format in the following text. * 66 * See oper_32b.c * 67 * * 68 * The R[i] were normalized in routine AUTO (hence, R[i] < 1.0). * 69 * The K[i] and Alpha are theoretically < 1.0. * 70 * The A[i], for a sampling frequency of 8 kHz, are in practice * 71 * always inferior to 16.0. * 72 * * 73 * These characteristics allow straigthforward fixed-point * 74 * implementation. We choose to represent the parameters as * 75 * follows : * 76 * * 77 * R[i] Q31 +- .99.. * 78 * K[i] Q31 +- .99.. * 79 * Alpha Normalized -> mantissa in Q31 plus exponent * 80 * A[i] Q27 +- 15.999.. * 81 * * 82 * The additions are performed in 32 bit. For the summation used * 83 * to calculate the K[i], we multiply numbers in Q31 by numbers * 84 * in Q27, with the result of the multiplications in Q27, * 85 * resulting in a dynamic of +- 16. This is sufficient to avoid * 86 * overflow, since the final result of the summation is * 87 * necessarily < 1.0 as both the K[i] and Alpha are * 88 * theoretically < 1.0. * 89 *___________________________________________________________________________*/ 90#include "typedef.h" 91#include "basic_op.h" 92#include "oper_32b.h" 93#include "acelp.h" 94 95#define M 16 96#define NC (M/2) 97 98void Init_Levinson( 99 Word16 * mem /* output :static memory (18 words) */ 100 ) 101{ 102 Set_zero(mem, 18); /* old_A[0..M-1] = 0, old_rc[0..1] = 0 */ 103 return; 104} 105 106 107void Levinson( 108 Word16 Rh[], /* (i) : Rh[M+1] Vector of autocorrelations (msb) */ 109 Word16 Rl[], /* (i) : Rl[M+1] Vector of autocorrelations (lsb) */ 110 Word16 A[], /* (o) Q12 : A[M] LPC coefficients (m = 16) */ 111 Word16 rc[], /* (o) Q15 : rc[M] Reflection coefficients. */ 112 Word16 * mem /* (i/o) :static memory (18 words) */ 113 ) 114{ 115 Word32 i, j; 116 Word16 hi, lo; 117 Word16 Kh, Kl; /* reflection coefficient; hi and lo */ 118 Word16 alp_h, alp_l, alp_exp; /* Prediction gain; hi lo and exponent */ 119 Word16 Ah[M + 1], Al[M + 1]; /* LPC coef. in double prec. */ 120 Word16 Anh[M + 1], Anl[M + 1]; /* LPC coef.for next iteration in double prec. */ 121 Word32 t0, t1, t2; /* temporary variable */ 122 Word16 *old_A, *old_rc; 123 124 /* Last A(z) for case of unstable filter */ 125 old_A = mem; 126 old_rc = mem + M; 127 128 /* K = A[1] = -R[1] / R[0] */ 129 130 t1 = ((Rh[1] << 16) + (Rl[1] << 1)); /* R[1] in Q31 */ 131 t2 = L_abs(t1); /* abs R[1] */ 132 t0 = Div_32(t2, Rh[0], Rl[0]); /* R[1]/R[0] in Q31 */ 133 if (t1 > 0) 134 t0 = -t0; /* -R[1]/R[0] */ 135 136 Kh = t0 >> 16; 137 Kl = (t0 & 0xffff)>>1; 138 rc[0] = Kh; 139 t0 = (t0 >> 4); /* A[1] in Q27 */ 140 141 Ah[1] = t0 >> 16; 142 Al[1] = (t0 & 0xffff)>>1; 143 144 /* Alpha = R[0] * (1-K**2) */ 145 t0 = Mpy_32(Kh, Kl, Kh, Kl); /* K*K in Q31 */ 146 t0 = L_abs(t0); /* Some case <0 !! */ 147 t0 = vo_L_sub((Word32) 0x7fffffffL, t0); /* 1 - K*K in Q31 */ 148 149 hi = t0 >> 16; 150 lo = (t0 & 0xffff)>>1; 151 152 t0 = Mpy_32(Rh[0], Rl[0], hi, lo); /* Alpha in Q31 */ 153 154 /* Normalize Alpha */ 155 alp_exp = norm_l(t0); 156 t0 = (t0 << alp_exp); 157 158 alp_h = t0 >> 16; 159 alp_l = (t0 & 0xffff)>>1; 160 /*--------------------------------------* 161 * ITERATIONS I=2 to M * 162 *--------------------------------------*/ 163 for (i = 2; i <= M; i++) 164 { 165 /* t0 = SUM ( R[j]*A[i-j] ,j=1,i-1 ) + R[i] */ 166 t0 = 0; 167 for (j = 1; j < i; j++) 168 t0 = vo_L_add(t0, Mpy_32(Rh[j], Rl[j], Ah[i - j], Al[i - j])); 169 170 t0 = t0 << 4; /* result in Q27 -> convert to Q31 */ 171 /* No overflow possible */ 172 t1 = ((Rh[i] << 16) + (Rl[i] << 1)); 173 t0 = vo_L_add(t0, t1); /* add R[i] in Q31 */ 174 175 /* K = -t0 / Alpha */ 176 t1 = L_abs(t0); 177 t2 = Div_32(t1, alp_h, alp_l); /* abs(t0)/Alpha */ 178 if (t0 > 0) 179 t2 = -t2; /* K =-t0/Alpha */ 180 t2 = (t2 << alp_exp); /* denormalize; compare to Alpha */ 181 182 Kh = t2 >> 16; 183 Kl = (t2 & 0xffff)>>1; 184 185 rc[i - 1] = Kh; 186 /* Test for unstable filter. If unstable keep old A(z) */ 187 if (abs_s(Kh) > 32750) 188 { 189 A[0] = 4096; /* Ai[0] not stored (always 1.0) */ 190 for (j = 0; j < M; j++) 191 { 192 A[j + 1] = old_A[j]; 193 } 194 rc[0] = old_rc[0]; /* only two rc coefficients are needed */ 195 rc[1] = old_rc[1]; 196 return; 197 } 198 /*------------------------------------------* 199 * Compute new LPC coeff. -> An[i] * 200 * An[j]= A[j] + K*A[i-j] , j=1 to i-1 * 201 * An[i]= K * 202 *------------------------------------------*/ 203 for (j = 1; j < i; j++) 204 { 205 t0 = Mpy_32(Kh, Kl, Ah[i - j], Al[i - j]); 206 t0 = vo_L_add(t0, ((Ah[j] << 16) + (Al[j] << 1))); 207 Anh[j] = t0 >> 16; 208 Anl[j] = (t0 & 0xffff)>>1; 209 } 210 t2 = (t2 >> 4); /* t2 = K in Q31 ->convert to Q27 */ 211 212 VO_L_Extract(t2, &Anh[i], &Anl[i]); /* An[i] in Q27 */ 213 214 /* Alpha = Alpha * (1-K**2) */ 215 t0 = Mpy_32(Kh, Kl, Kh, Kl); /* K*K in Q31 */ 216 t0 = L_abs(t0); /* Some case <0 !! */ 217 t0 = vo_L_sub((Word32) 0x7fffffffL, t0); /* 1 - K*K in Q31 */ 218 hi = t0 >> 16; 219 lo = (t0 & 0xffff)>>1; 220 t0 = Mpy_32(alp_h, alp_l, hi, lo); /* Alpha in Q31 */ 221 222 /* Normalize Alpha */ 223 j = norm_l(t0); 224 t0 = (t0 << j); 225 alp_h = t0 >> 16; 226 alp_l = (t0 & 0xffff)>>1; 227 alp_exp += j; /* Add normalization to alp_exp */ 228 229 /* A[j] = An[j] */ 230 for (j = 1; j <= i; j++) 231 { 232 Ah[j] = Anh[j]; 233 Al[j] = Anl[j]; 234 } 235 } 236 /* Truncate A[i] in Q27 to Q12 with rounding */ 237 A[0] = 4096; 238 for (i = 1; i <= M; i++) 239 { 240 t0 = (Ah[i] << 16) + (Al[i] << 1); 241 old_A[i - 1] = A[i] = vo_round((t0 << 1)); 242 } 243 old_rc[0] = rc[0]; 244 old_rc[1] = rc[1]; 245 246 return; 247} 248 249 250 251