1 /* Copyright (c) 2009-2010 Xiph.Org Foundation 2 Written by Jean-Marc Valin */ 3 /* 4 Redistribution and use in source and binary forms, with or without 5 modification, are permitted provided that the following conditions 6 are met: 7 8 - Redistributions of source code must retain the above copyright 9 notice, this list of conditions and the following disclaimer. 10 11 - Redistributions in binary form must reproduce the above copyright 12 notice, this list of conditions and the following disclaimer in the 13 documentation and/or other materials provided with the distribution. 14 15 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS 16 ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT 17 LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR 18 A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER 19 OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, 20 EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, 21 PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR 22 PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF 23 LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING 24 NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS 25 SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. 26 */ 27 28 #ifdef HAVE_CONFIG_H 29 #include "config.h" 30 #endif 31 32 #include "celt_lpc.h" 33 #include "stack_alloc.h" 34 #include "mathops.h" 35 #include "pitch.h" 36 37 void _celt_lpc( 38 opus_val16 *_lpc, /* out: [0...p-1] LPC coefficients */ 39 const opus_val32 *ac, /* in: [0...p] autocorrelation values */ 40 int p 41 ) 42 { 43 int i, j; 44 opus_val32 r; 45 opus_val32 error = ac[0]; 46 #ifdef FIXED_POINT 47 opus_val32 lpc[LPC_ORDER]; 48 #else 49 float *lpc = _lpc; 50 #endif 51 52 OPUS_CLEAR(lpc, p); 53 if (ac[0] != 0) 54 { 55 for (i = 0; i < p; i++) { 56 /* Sum up this iteration's reflection coefficient */ 57 opus_val32 rr = 0; 58 for (j = 0; j < i; j++) 59 rr += MULT32_32_Q31(lpc[j],ac[i - j]); 60 rr += SHR32(ac[i + 1],3); 61 r = -frac_div32(SHL32(rr,3), error); 62 /* Update LPC coefficients and total error */ 63 lpc[i] = SHR32(r,3); 64 for (j = 0; j < (i+1)>>1; j++) 65 { 66 opus_val32 tmp1, tmp2; 67 tmp1 = lpc[j]; 68 tmp2 = lpc[i-1-j]; 69 lpc[j] = tmp1 + MULT32_32_Q31(r,tmp2); 70 lpc[i-1-j] = tmp2 + MULT32_32_Q31(r,tmp1); 71 } 72 73 error = error - MULT32_32_Q31(MULT32_32_Q31(r,r),error); 74 /* Bail out once we get 30 dB gain */ 75 #ifdef FIXED_POINT 76 if (error<SHR32(ac[0],10)) 77 break; 78 #else 79 if (error<.001f*ac[0]) 80 break; 81 #endif 82 } 83 } 84 #ifdef FIXED_POINT 85 for (i=0;i<p;i++) 86 _lpc[i] = ROUND16(lpc[i],16); 87 #endif 88 } 89 90 91 void celt_fir_c( 92 const opus_val16 *x, 93 const opus_val16 *num, 94 opus_val16 *y, 95 int N, 96 int ord, 97 int arch) 98 { 99 int i,j; 100 VARDECL(opus_val16, rnum); 101 SAVE_STACK; 102 103 ALLOC(rnum, ord, opus_val16); 104 for(i=0;i<ord;i++) 105 rnum[i] = num[ord-i-1]; 106 for (i=0;i<N-3;i+=4) 107 { 108 opus_val32 sum[4]; 109 sum[0] = SHL32(EXTEND32(x[i ]), SIG_SHIFT); 110 sum[1] = SHL32(EXTEND32(x[i+1]), SIG_SHIFT), 111 sum[2] = SHL32(EXTEND32(x[i+2]), SIG_SHIFT); 112 sum[3] = SHL32(EXTEND32(x[i+3]), SIG_SHIFT); 113 xcorr_kernel(rnum, x+i-ord, sum, ord, arch); 114 y[i ] = ROUND16(sum[0], SIG_SHIFT); 115 y[i+1] = ROUND16(sum[1], SIG_SHIFT); 116 y[i+2] = ROUND16(sum[2], SIG_SHIFT); 117 y[i+3] = ROUND16(sum[3], SIG_SHIFT); 118 } 119 for (;i<N;i++) 120 { 121 opus_val32 sum = SHL32(EXTEND32(x[i]), SIG_SHIFT); 122 for (j=0;j<ord;j++) 123 sum = MAC16_16(sum,rnum[j],x[i+j-ord]); 124 y[i] = ROUND16(sum, SIG_SHIFT); 125 } 126 RESTORE_STACK; 127 } 128 129 void celt_iir(const opus_val32 *_x, 130 const opus_val16 *den, 131 opus_val32 *_y, 132 int N, 133 int ord, 134 opus_val16 *mem, 135 int arch) 136 { 137 #ifdef SMALL_FOOTPRINT 138 int i,j; 139 (void)arch; 140 for (i=0;i<N;i++) 141 { 142 opus_val32 sum = _x[i]; 143 for (j=0;j<ord;j++) 144 { 145 sum -= MULT16_16(den[j],mem[j]); 146 } 147 for (j=ord-1;j>=1;j--) 148 { 149 mem[j]=mem[j-1]; 150 } 151 mem[0] = SROUND16(sum, SIG_SHIFT); 152 _y[i] = sum; 153 } 154 #else 155 int i,j; 156 VARDECL(opus_val16, rden); 157 VARDECL(opus_val16, y); 158 SAVE_STACK; 159 160 celt_assert((ord&3)==0); 161 ALLOC(rden, ord, opus_val16); 162 ALLOC(y, N+ord, opus_val16); 163 for(i=0;i<ord;i++) 164 rden[i] = den[ord-i-1]; 165 for(i=0;i<ord;i++) 166 y[i] = -mem[ord-i-1]; 167 for(;i<N+ord;i++) 168 y[i]=0; 169 for (i=0;i<N-3;i+=4) 170 { 171 /* Unroll by 4 as if it were an FIR filter */ 172 opus_val32 sum[4]; 173 sum[0]=_x[i]; 174 sum[1]=_x[i+1]; 175 sum[2]=_x[i+2]; 176 sum[3]=_x[i+3]; 177 xcorr_kernel(rden, y+i, sum, ord, arch); 178 179 /* Patch up the result to compensate for the fact that this is an IIR */ 180 y[i+ord ] = -SROUND16(sum[0],SIG_SHIFT); 181 _y[i ] = sum[0]; 182 sum[1] = MAC16_16(sum[1], y[i+ord ], den[0]); 183 y[i+ord+1] = -SROUND16(sum[1],SIG_SHIFT); 184 _y[i+1] = sum[1]; 185 sum[2] = MAC16_16(sum[2], y[i+ord+1], den[0]); 186 sum[2] = MAC16_16(sum[2], y[i+ord ], den[1]); 187 y[i+ord+2] = -SROUND16(sum[2],SIG_SHIFT); 188 _y[i+2] = sum[2]; 189 190 sum[3] = MAC16_16(sum[3], y[i+ord+2], den[0]); 191 sum[3] = MAC16_16(sum[3], y[i+ord+1], den[1]); 192 sum[3] = MAC16_16(sum[3], y[i+ord ], den[2]); 193 y[i+ord+3] = -SROUND16(sum[3],SIG_SHIFT); 194 _y[i+3] = sum[3]; 195 } 196 for (;i<N;i++) 197 { 198 opus_val32 sum = _x[i]; 199 for (j=0;j<ord;j++) 200 sum -= MULT16_16(rden[j],y[i+j]); 201 y[i+ord] = SROUND16(sum,SIG_SHIFT); 202 _y[i] = sum; 203 } 204 for(i=0;i<ord;i++) 205 mem[i] = _y[N-i-1]; 206 RESTORE_STACK; 207 #endif 208 } 209 210 int _celt_autocorr( 211 const opus_val16 *x, /* in: [0...n-1] samples x */ 212 opus_val32 *ac, /* out: [0...lag-1] ac values */ 213 const opus_val16 *window, 214 int overlap, 215 int lag, 216 int n, 217 int arch 218 ) 219 { 220 opus_val32 d; 221 int i, k; 222 int fastN=n-lag; 223 int shift; 224 const opus_val16 *xptr; 225 VARDECL(opus_val16, xx); 226 SAVE_STACK; 227 ALLOC(xx, n, opus_val16); 228 celt_assert(n>0); 229 celt_assert(overlap>=0); 230 if (overlap == 0) 231 { 232 xptr = x; 233 } else { 234 for (i=0;i<n;i++) 235 xx[i] = x[i]; 236 for (i=0;i<overlap;i++) 237 { 238 xx[i] = MULT16_16_Q15(x[i],window[i]); 239 xx[n-i-1] = MULT16_16_Q15(x[n-i-1],window[i]); 240 } 241 xptr = xx; 242 } 243 shift=0; 244 #ifdef FIXED_POINT 245 { 246 opus_val32 ac0; 247 ac0 = 1+(n<<7); 248 if (n&1) ac0 += SHR32(MULT16_16(xptr[0],xptr[0]),9); 249 for(i=(n&1);i<n;i+=2) 250 { 251 ac0 += SHR32(MULT16_16(xptr[i],xptr[i]),9); 252 ac0 += SHR32(MULT16_16(xptr[i+1],xptr[i+1]),9); 253 } 254 255 shift = celt_ilog2(ac0)-30+10; 256 shift = (shift)/2; 257 if (shift>0) 258 { 259 for(i=0;i<n;i++) 260 xx[i] = PSHR32(xptr[i], shift); 261 xptr = xx; 262 } else 263 shift = 0; 264 } 265 #endif 266 celt_pitch_xcorr(xptr, xptr, ac, fastN, lag+1, arch); 267 for (k=0;k<=lag;k++) 268 { 269 for (i = k+fastN, d = 0; i < n; i++) 270 d = MAC16_16(d, xptr[i], xptr[i-k]); 271 ac[k] += d; 272 } 273 #ifdef FIXED_POINT 274 shift = 2*shift; 275 if (shift<=0) 276 ac[0] += SHL32((opus_int32)1, -shift); 277 if (ac[0] < 268435456) 278 { 279 int shift2 = 29 - EC_ILOG(ac[0]); 280 for (i=0;i<=lag;i++) 281 ac[i] = SHL32(ac[i], shift2); 282 shift -= shift2; 283 } else if (ac[0] >= 536870912) 284 { 285 int shift2=1; 286 if (ac[0] >= 1073741824) 287 shift2++; 288 for (i=0;i<=lag;i++) 289 ac[i] = SHR32(ac[i], shift2); 290 shift += shift2; 291 } 292 #endif 293 294 RESTORE_STACK; 295 return shift; 296 } 297
