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 for (i = 0; i < p; i++) 53 lpc[i] = 0; 54 if (ac[0] != 0) 55 { 56 for (i = 0; i < p; i++) { 57 /* Sum up this iteration's reflection coefficient */ 58 opus_val32 rr = 0; 59 for (j = 0; j < i; j++) 60 rr += MULT32_32_Q31(lpc[j],ac[i - j]); 61 rr += SHR32(ac[i + 1],3); 62 r = -frac_div32(SHL32(rr,3), error); 63 /* Update LPC coefficients and total error */ 64 lpc[i] = SHR32(r,3); 65 for (j = 0; j < (i+1)>>1; j++) 66 { 67 opus_val32 tmp1, tmp2; 68 tmp1 = lpc[j]; 69 tmp2 = lpc[i-1-j]; 70 lpc[j] = tmp1 + MULT32_32_Q31(r,tmp2); 71 lpc[i-1-j] = tmp2 + MULT32_32_Q31(r,tmp1); 72 } 73 74 error = error - MULT32_32_Q31(MULT32_32_Q31(r,r),error); 75 /* Bail out once we get 30 dB gain */ 76 #ifdef FIXED_POINT 77 if (error<SHR32(ac[0],10)) 78 break; 79 #else 80 if (error<.001f*ac[0]) 81 break; 82 #endif 83 } 84 } 85 #ifdef FIXED_POINT 86 for (i=0;i<p;i++) 87 _lpc[i] = ROUND16(lpc[i],16); 88 #endif 89 } 90 91 void celt_fir(const opus_val16 *_x, 92 const opus_val16 *num, 93 opus_val16 *_y, 94 int N, 95 int ord, 96 opus_val16 *mem) 97 { 98 int i,j; 99 VARDECL(opus_val16, rnum); 100 VARDECL(opus_val16, x); 101 SAVE_STACK; 102 103 ALLOC(rnum, ord, opus_val16); 104 ALLOC(x, N+ord, opus_val16); 105 for(i=0;i<ord;i++) 106 rnum[i] = num[ord-i-1]; 107 for(i=0;i<ord;i++) 108 x[i] = mem[ord-i-1]; 109 for (i=0;i<N;i++) 110 x[i+ord]=_x[i]; 111 for(i=0;i<ord;i++) 112 mem[i] = _x[N-i-1]; 113 #ifdef SMALL_FOOTPRINT 114 for (i=0;i<N;i++) 115 { 116 opus_val32 sum = SHL32(EXTEND32(_x[i]), SIG_SHIFT); 117 for (j=0;j<ord;j++) 118 { 119 sum = MAC16_16(sum,rnum[j],x[i+j]); 120 } 121 _y[i] = SATURATE16(PSHR32(sum, SIG_SHIFT)); 122 } 123 #else 124 for (i=0;i<N-3;i+=4) 125 { 126 opus_val32 sum[4]={0,0,0,0}; 127 xcorr_kernel(rnum, x+i, sum, ord); 128 _y[i ] = SATURATE16(ADD32(EXTEND32(_x[i ]), PSHR32(sum[0], SIG_SHIFT))); 129 _y[i+1] = SATURATE16(ADD32(EXTEND32(_x[i+1]), PSHR32(sum[1], SIG_SHIFT))); 130 _y[i+2] = SATURATE16(ADD32(EXTEND32(_x[i+2]), PSHR32(sum[2], SIG_SHIFT))); 131 _y[i+3] = SATURATE16(ADD32(EXTEND32(_x[i+3]), PSHR32(sum[3], SIG_SHIFT))); 132 } 133 for (;i<N;i++) 134 { 135 opus_val32 sum = 0; 136 for (j=0;j<ord;j++) 137 sum = MAC16_16(sum,rnum[j],x[i+j]); 138 _y[i] = SATURATE16(ADD32(EXTEND32(_x[i]), PSHR32(sum, SIG_SHIFT))); 139 } 140 #endif 141 RESTORE_STACK; 142 } 143 144 void celt_iir(const opus_val32 *_x, 145 const opus_val16 *den, 146 opus_val32 *_y, 147 int N, 148 int ord, 149 opus_val16 *mem) 150 { 151 #ifdef SMALL_FOOTPRINT 152 int i,j; 153 for (i=0;i<N;i++) 154 { 155 opus_val32 sum = _x[i]; 156 for (j=0;j<ord;j++) 157 { 158 sum -= MULT16_16(den[j],mem[j]); 159 } 160 for (j=ord-1;j>=1;j--) 161 { 162 mem[j]=mem[j-1]; 163 } 164 mem[0] = ROUND16(sum,SIG_SHIFT); 165 _y[i] = sum; 166 } 167 #else 168 int i,j; 169 VARDECL(opus_val16, rden); 170 VARDECL(opus_val16, y); 171 SAVE_STACK; 172 173 celt_assert((ord&3)==0); 174 ALLOC(rden, ord, opus_val16); 175 ALLOC(y, N+ord, opus_val16); 176 for(i=0;i<ord;i++) 177 rden[i] = den[ord-i-1]; 178 for(i=0;i<ord;i++) 179 y[i] = -mem[ord-i-1]; 180 for(;i<N+ord;i++) 181 y[i]=0; 182 for (i=0;i<N-3;i+=4) 183 { 184 /* Unroll by 4 as if it were an FIR filter */ 185 opus_val32 sum[4]; 186 sum[0]=_x[i]; 187 sum[1]=_x[i+1]; 188 sum[2]=_x[i+2]; 189 sum[3]=_x[i+3]; 190 xcorr_kernel(rden, y+i, sum, ord); 191 192 /* Patch up the result to compensate for the fact that this is an IIR */ 193 y[i+ord ] = -ROUND16(sum[0],SIG_SHIFT); 194 _y[i ] = sum[0]; 195 sum[1] = MAC16_16(sum[1], y[i+ord ], den[0]); 196 y[i+ord+1] = -ROUND16(sum[1],SIG_SHIFT); 197 _y[i+1] = sum[1]; 198 sum[2] = MAC16_16(sum[2], y[i+ord+1], den[0]); 199 sum[2] = MAC16_16(sum[2], y[i+ord ], den[1]); 200 y[i+ord+2] = -ROUND16(sum[2],SIG_SHIFT); 201 _y[i+2] = sum[2]; 202 203 sum[3] = MAC16_16(sum[3], y[i+ord+2], den[0]); 204 sum[3] = MAC16_16(sum[3], y[i+ord+1], den[1]); 205 sum[3] = MAC16_16(sum[3], y[i+ord ], den[2]); 206 y[i+ord+3] = -ROUND16(sum[3],SIG_SHIFT); 207 _y[i+3] = sum[3]; 208 } 209 for (;i<N;i++) 210 { 211 opus_val32 sum = _x[i]; 212 for (j=0;j<ord;j++) 213 sum -= MULT16_16(rden[j],y[i+j]); 214 y[i+ord] = ROUND16(sum,SIG_SHIFT); 215 _y[i] = sum; 216 } 217 for(i=0;i<ord;i++) 218 mem[i] = _y[N-i-1]; 219 RESTORE_STACK; 220 #endif 221 } 222 223 int _celt_autocorr( 224 const opus_val16 *x, /* in: [0...n-1] samples x */ 225 opus_val32 *ac, /* out: [0...lag-1] ac values */ 226 const opus_val16 *window, 227 int overlap, 228 int lag, 229 int n, 230 int arch 231 ) 232 { 233 opus_val32 d; 234 int i, k; 235 int fastN=n-lag; 236 int shift; 237 const opus_val16 *xptr; 238 VARDECL(opus_val16, xx); 239 SAVE_STACK; 240 ALLOC(xx, n, opus_val16); 241 celt_assert(n>0); 242 celt_assert(overlap>=0); 243 if (overlap == 0) 244 { 245 xptr = x; 246 } else { 247 for (i=0;i<n;i++) 248 xx[i] = x[i]; 249 for (i=0;i<overlap;i++) 250 { 251 xx[i] = MULT16_16_Q15(x[i],window[i]); 252 xx[n-i-1] = MULT16_16_Q15(x[n-i-1],window[i]); 253 } 254 xptr = xx; 255 } 256 shift=0; 257 #ifdef FIXED_POINT 258 { 259 opus_val32 ac0; 260 ac0 = 1+(n<<7); 261 if (n&1) ac0 += SHR32(MULT16_16(xptr[0],xptr[0]),9); 262 for(i=(n&1);i<n;i+=2) 263 { 264 ac0 += SHR32(MULT16_16(xptr[i],xptr[i]),9); 265 ac0 += SHR32(MULT16_16(xptr[i+1],xptr[i+1]),9); 266 } 267 268 shift = celt_ilog2(ac0)-30+10; 269 shift = (shift)/2; 270 if (shift>0) 271 { 272 for(i=0;i<n;i++) 273 xx[i] = PSHR32(xptr[i], shift); 274 xptr = xx; 275 } else 276 shift = 0; 277 } 278 #endif 279 celt_pitch_xcorr(xptr, xptr, ac, fastN, lag+1, arch); 280 for (k=0;k<=lag;k++) 281 { 282 for (i = k+fastN, d = 0; i < n; i++) 283 d = MAC16_16(d, xptr[i], xptr[i-k]); 284 ac[k] += d; 285 } 286 #ifdef FIXED_POINT 287 shift = 2*shift; 288 if (shift<=0) 289 ac[0] += SHL32((opus_int32)1, -shift); 290 if (ac[0] < 268435456) 291 { 292 int shift2 = 29 - EC_ILOG(ac[0]); 293 for (i=0;i<=lag;i++) 294 ac[i] = SHL32(ac[i], shift2); 295 shift -= shift2; 296 } else if (ac[0] >= 536870912) 297 { 298 int shift2=1; 299 if (ac[0] >= 1073741824) 300 shift2++; 301 for (i=0;i<=lag;i++) 302 ac[i] = SHR32(ac[i], shift2); 303 shift += shift2; 304 } 305 #endif 306 307 RESTORE_STACK; 308 return shift; 309 } 310
