1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/media/libopus/celt/celt_lpc.c Wed Dec 31 06:09:35 2014 +0100
1.3 @@ -0,0 +1,309 @@
1.4 +/* Copyright (c) 2009-2010 Xiph.Org Foundation
1.5 + Written by Jean-Marc Valin */
1.6 +/*
1.7 + Redistribution and use in source and binary forms, with or without
1.8 + modification, are permitted provided that the following conditions
1.9 + are met:
1.10 +
1.11 + - Redistributions of source code must retain the above copyright
1.12 + notice, this list of conditions and the following disclaimer.
1.13 +
1.14 + - Redistributions in binary form must reproduce the above copyright
1.15 + notice, this list of conditions and the following disclaimer in the
1.16 + documentation and/or other materials provided with the distribution.
1.17 +
1.18 + THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
1.19 + ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
1.20 + LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
1.21 + A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER
1.22 + OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
1.23 + EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
1.24 + PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
1.25 + PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
1.26 + LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
1.27 + NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
1.28 + SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
1.29 +*/
1.30 +
1.31 +#ifdef HAVE_CONFIG_H
1.32 +#include "config.h"
1.33 +#endif
1.34 +
1.35 +#include "celt_lpc.h"
1.36 +#include "stack_alloc.h"
1.37 +#include "mathops.h"
1.38 +#include "pitch.h"
1.39 +
1.40 +void _celt_lpc(
1.41 + opus_val16 *_lpc, /* out: [0...p-1] LPC coefficients */
1.42 +const opus_val32 *ac, /* in: [0...p] autocorrelation values */
1.43 +int p
1.44 +)
1.45 +{
1.46 + int i, j;
1.47 + opus_val32 r;
1.48 + opus_val32 error = ac[0];
1.49 +#ifdef FIXED_POINT
1.50 + opus_val32 lpc[LPC_ORDER];
1.51 +#else
1.52 + float *lpc = _lpc;
1.53 +#endif
1.54 +
1.55 + for (i = 0; i < p; i++)
1.56 + lpc[i] = 0;
1.57 + if (ac[0] != 0)
1.58 + {
1.59 + for (i = 0; i < p; i++) {
1.60 + /* Sum up this iteration's reflection coefficient */
1.61 + opus_val32 rr = 0;
1.62 + for (j = 0; j < i; j++)
1.63 + rr += MULT32_32_Q31(lpc[j],ac[i - j]);
1.64 + rr += SHR32(ac[i + 1],3);
1.65 + r = -frac_div32(SHL32(rr,3), error);
1.66 + /* Update LPC coefficients and total error */
1.67 + lpc[i] = SHR32(r,3);
1.68 + for (j = 0; j < (i+1)>>1; j++)
1.69 + {
1.70 + opus_val32 tmp1, tmp2;
1.71 + tmp1 = lpc[j];
1.72 + tmp2 = lpc[i-1-j];
1.73 + lpc[j] = tmp1 + MULT32_32_Q31(r,tmp2);
1.74 + lpc[i-1-j] = tmp2 + MULT32_32_Q31(r,tmp1);
1.75 + }
1.76 +
1.77 + error = error - MULT32_32_Q31(MULT32_32_Q31(r,r),error);
1.78 + /* Bail out once we get 30 dB gain */
1.79 +#ifdef FIXED_POINT
1.80 + if (error<SHR32(ac[0],10))
1.81 + break;
1.82 +#else
1.83 + if (error<.001f*ac[0])
1.84 + break;
1.85 +#endif
1.86 + }
1.87 + }
1.88 +#ifdef FIXED_POINT
1.89 + for (i=0;i<p;i++)
1.90 + _lpc[i] = ROUND16(lpc[i],16);
1.91 +#endif
1.92 +}
1.93 +
1.94 +void celt_fir(const opus_val16 *_x,
1.95 + const opus_val16 *num,
1.96 + opus_val16 *_y,
1.97 + int N,
1.98 + int ord,
1.99 + opus_val16 *mem)
1.100 +{
1.101 + int i,j;
1.102 + VARDECL(opus_val16, rnum);
1.103 + VARDECL(opus_val16, x);
1.104 + SAVE_STACK;
1.105 +
1.106 + ALLOC(rnum, ord, opus_val16);
1.107 + ALLOC(x, N+ord, opus_val16);
1.108 + for(i=0;i<ord;i++)
1.109 + rnum[i] = num[ord-i-1];
1.110 + for(i=0;i<ord;i++)
1.111 + x[i] = mem[ord-i-1];
1.112 + for (i=0;i<N;i++)
1.113 + x[i+ord]=_x[i];
1.114 + for(i=0;i<ord;i++)
1.115 + mem[i] = _x[N-i-1];
1.116 +#ifdef SMALL_FOOTPRINT
1.117 + for (i=0;i<N;i++)
1.118 + {
1.119 + opus_val32 sum = SHL32(EXTEND32(_x[i]), SIG_SHIFT);
1.120 + for (j=0;j<ord;j++)
1.121 + {
1.122 + sum = MAC16_16(sum,rnum[j],x[i+j]);
1.123 + }
1.124 + _y[i] = SATURATE16(PSHR32(sum, SIG_SHIFT));
1.125 + }
1.126 +#else
1.127 + for (i=0;i<N-3;i+=4)
1.128 + {
1.129 + opus_val32 sum[4]={0,0,0,0};
1.130 + xcorr_kernel(rnum, x+i, sum, ord);
1.131 + _y[i ] = SATURATE16(ADD32(EXTEND32(_x[i ]), PSHR32(sum[0], SIG_SHIFT)));
1.132 + _y[i+1] = SATURATE16(ADD32(EXTEND32(_x[i+1]), PSHR32(sum[1], SIG_SHIFT)));
1.133 + _y[i+2] = SATURATE16(ADD32(EXTEND32(_x[i+2]), PSHR32(sum[2], SIG_SHIFT)));
1.134 + _y[i+3] = SATURATE16(ADD32(EXTEND32(_x[i+3]), PSHR32(sum[3], SIG_SHIFT)));
1.135 + }
1.136 + for (;i<N;i++)
1.137 + {
1.138 + opus_val32 sum = 0;
1.139 + for (j=0;j<ord;j++)
1.140 + sum = MAC16_16(sum,rnum[j],x[i+j]);
1.141 + _y[i] = SATURATE16(ADD32(EXTEND32(_x[i]), PSHR32(sum, SIG_SHIFT)));
1.142 + }
1.143 +#endif
1.144 + RESTORE_STACK;
1.145 +}
1.146 +
1.147 +void celt_iir(const opus_val32 *_x,
1.148 + const opus_val16 *den,
1.149 + opus_val32 *_y,
1.150 + int N,
1.151 + int ord,
1.152 + opus_val16 *mem)
1.153 +{
1.154 +#ifdef SMALL_FOOTPRINT
1.155 + int i,j;
1.156 + for (i=0;i<N;i++)
1.157 + {
1.158 + opus_val32 sum = _x[i];
1.159 + for (j=0;j<ord;j++)
1.160 + {
1.161 + sum -= MULT16_16(den[j],mem[j]);
1.162 + }
1.163 + for (j=ord-1;j>=1;j--)
1.164 + {
1.165 + mem[j]=mem[j-1];
1.166 + }
1.167 + mem[0] = ROUND16(sum,SIG_SHIFT);
1.168 + _y[i] = sum;
1.169 + }
1.170 +#else
1.171 + int i,j;
1.172 + VARDECL(opus_val16, rden);
1.173 + VARDECL(opus_val16, y);
1.174 + SAVE_STACK;
1.175 +
1.176 + celt_assert((ord&3)==0);
1.177 + ALLOC(rden, ord, opus_val16);
1.178 + ALLOC(y, N+ord, opus_val16);
1.179 + for(i=0;i<ord;i++)
1.180 + rden[i] = den[ord-i-1];
1.181 + for(i=0;i<ord;i++)
1.182 + y[i] = -mem[ord-i-1];
1.183 + for(;i<N+ord;i++)
1.184 + y[i]=0;
1.185 + for (i=0;i<N-3;i+=4)
1.186 + {
1.187 + /* Unroll by 4 as if it were an FIR filter */
1.188 + opus_val32 sum[4];
1.189 + sum[0]=_x[i];
1.190 + sum[1]=_x[i+1];
1.191 + sum[2]=_x[i+2];
1.192 + sum[3]=_x[i+3];
1.193 + xcorr_kernel(rden, y+i, sum, ord);
1.194 +
1.195 + /* Patch up the result to compensate for the fact that this is an IIR */
1.196 + y[i+ord ] = -ROUND16(sum[0],SIG_SHIFT);
1.197 + _y[i ] = sum[0];
1.198 + sum[1] = MAC16_16(sum[1], y[i+ord ], den[0]);
1.199 + y[i+ord+1] = -ROUND16(sum[1],SIG_SHIFT);
1.200 + _y[i+1] = sum[1];
1.201 + sum[2] = MAC16_16(sum[2], y[i+ord+1], den[0]);
1.202 + sum[2] = MAC16_16(sum[2], y[i+ord ], den[1]);
1.203 + y[i+ord+2] = -ROUND16(sum[2],SIG_SHIFT);
1.204 + _y[i+2] = sum[2];
1.205 +
1.206 + sum[3] = MAC16_16(sum[3], y[i+ord+2], den[0]);
1.207 + sum[3] = MAC16_16(sum[3], y[i+ord+1], den[1]);
1.208 + sum[3] = MAC16_16(sum[3], y[i+ord ], den[2]);
1.209 + y[i+ord+3] = -ROUND16(sum[3],SIG_SHIFT);
1.210 + _y[i+3] = sum[3];
1.211 + }
1.212 + for (;i<N;i++)
1.213 + {
1.214 + opus_val32 sum = _x[i];
1.215 + for (j=0;j<ord;j++)
1.216 + sum -= MULT16_16(rden[j],y[i+j]);
1.217 + y[i+ord] = ROUND16(sum,SIG_SHIFT);
1.218 + _y[i] = sum;
1.219 + }
1.220 + for(i=0;i<ord;i++)
1.221 + mem[i] = _y[N-i-1];
1.222 + RESTORE_STACK;
1.223 +#endif
1.224 +}
1.225 +
1.226 +int _celt_autocorr(
1.227 + const opus_val16 *x, /* in: [0...n-1] samples x */
1.228 + opus_val32 *ac, /* out: [0...lag-1] ac values */
1.229 + const opus_val16 *window,
1.230 + int overlap,
1.231 + int lag,
1.232 + int n,
1.233 + int arch
1.234 + )
1.235 +{
1.236 + opus_val32 d;
1.237 + int i, k;
1.238 + int fastN=n-lag;
1.239 + int shift;
1.240 + const opus_val16 *xptr;
1.241 + VARDECL(opus_val16, xx);
1.242 + SAVE_STACK;
1.243 + ALLOC(xx, n, opus_val16);
1.244 + celt_assert(n>0);
1.245 + celt_assert(overlap>=0);
1.246 + if (overlap == 0)
1.247 + {
1.248 + xptr = x;
1.249 + } else {
1.250 + for (i=0;i<n;i++)
1.251 + xx[i] = x[i];
1.252 + for (i=0;i<overlap;i++)
1.253 + {
1.254 + xx[i] = MULT16_16_Q15(x[i],window[i]);
1.255 + xx[n-i-1] = MULT16_16_Q15(x[n-i-1],window[i]);
1.256 + }
1.257 + xptr = xx;
1.258 + }
1.259 + shift=0;
1.260 +#ifdef FIXED_POINT
1.261 + {
1.262 + opus_val32 ac0;
1.263 + ac0 = 1+(n<<7);
1.264 + if (n&1) ac0 += SHR32(MULT16_16(xptr[0],xptr[0]),9);
1.265 + for(i=(n&1);i<n;i+=2)
1.266 + {
1.267 + ac0 += SHR32(MULT16_16(xptr[i],xptr[i]),9);
1.268 + ac0 += SHR32(MULT16_16(xptr[i+1],xptr[i+1]),9);
1.269 + }
1.270 +
1.271 + shift = celt_ilog2(ac0)-30+10;
1.272 + shift = (shift)/2;
1.273 + if (shift>0)
1.274 + {
1.275 + for(i=0;i<n;i++)
1.276 + xx[i] = PSHR32(xptr[i], shift);
1.277 + xptr = xx;
1.278 + } else
1.279 + shift = 0;
1.280 + }
1.281 +#endif
1.282 + celt_pitch_xcorr(xptr, xptr, ac, fastN, lag+1, arch);
1.283 + for (k=0;k<=lag;k++)
1.284 + {
1.285 + for (i = k+fastN, d = 0; i < n; i++)
1.286 + d = MAC16_16(d, xptr[i], xptr[i-k]);
1.287 + ac[k] += d;
1.288 + }
1.289 +#ifdef FIXED_POINT
1.290 + shift = 2*shift;
1.291 + if (shift<=0)
1.292 + ac[0] += SHL32((opus_int32)1, -shift);
1.293 + if (ac[0] < 268435456)
1.294 + {
1.295 + int shift2 = 29 - EC_ILOG(ac[0]);
1.296 + for (i=0;i<=lag;i++)
1.297 + ac[i] = SHL32(ac[i], shift2);
1.298 + shift -= shift2;
1.299 + } else if (ac[0] >= 536870912)
1.300 + {
1.301 + int shift2=1;
1.302 + if (ac[0] >= 1073741824)
1.303 + shift2++;
1.304 + for (i=0;i<=lag;i++)
1.305 + ac[i] = SHR32(ac[i], shift2);
1.306 + shift += shift2;
1.307 + }
1.308 +#endif
1.309 +
1.310 + RESTORE_STACK;
1.311 + return shift;
1.312 +}
