media/libopus/celt/celt_lpc.c

Wed, 31 Dec 2014 06:09:35 +0100

author
Michael Schloh von Bennewitz <michael@schloh.com>
date
Wed, 31 Dec 2014 06:09:35 +0100
changeset 0
6474c204b198
permissions
-rw-r--r--

Cloned upstream origin tor-browser at tor-browser-31.3.0esr-4.5-1-build1
revision ID fc1c9ff7c1b2defdbc039f12214767608f46423f for hacking purpose.

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

mercurial