summaryrefslogtreecommitdiffstats
path: root/media/libopus/celt/celt_lpc.c
diff options
context:
space:
mode:
Diffstat (limited to 'media/libopus/celt/celt_lpc.c')
-rw-r--r--media/libopus/celt/celt_lpc.c344
1 files changed, 344 insertions, 0 deletions
diff --git a/media/libopus/celt/celt_lpc.c b/media/libopus/celt/celt_lpc.c
new file mode 100644
index 0000000000..f91721bcab
--- /dev/null
+++ b/media/libopus/celt/celt_lpc.c
@@ -0,0 +1,344 @@
+/* Copyright (c) 2009-2010 Xiph.Org Foundation
+ Written by Jean-Marc Valin */
+/*
+ Redistribution and use in source and binary forms, with or without
+ modification, are permitted provided that the following conditions
+ are met:
+
+ - Redistributions of source code must retain the above copyright
+ notice, this list of conditions and the following disclaimer.
+
+ - Redistributions in binary form must reproduce the above copyright
+ notice, this list of conditions and the following disclaimer in the
+ documentation and/or other materials provided with the distribution.
+
+ THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+ ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+ LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+ A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER
+ OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
+ EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
+ PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
+ PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
+ LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
+ NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
+ SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+*/
+
+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#include "celt_lpc.h"
+#include "stack_alloc.h"
+#include "mathops.h"
+#include "pitch.h"
+
+void _celt_lpc(
+ opus_val16 *_lpc, /* out: [0...p-1] LPC coefficients */
+const opus_val32 *ac, /* in: [0...p] autocorrelation values */
+int p
+)
+{
+ int i, j;
+ opus_val32 r;
+ opus_val32 error = ac[0];
+#ifdef FIXED_POINT
+ opus_val32 lpc[LPC_ORDER];
+#else
+ float *lpc = _lpc;
+#endif
+
+ OPUS_CLEAR(lpc, p);
+#ifdef FIXED_POINT
+ if (ac[0] != 0)
+#else
+ if (ac[0] > 1e-10f)
+#endif
+ {
+ for (i = 0; i < p; i++) {
+ /* Sum up this iteration's reflection coefficient */
+ opus_val32 rr = 0;
+ for (j = 0; j < i; j++)
+ rr += MULT32_32_Q31(lpc[j],ac[i - j]);
+ rr += SHR32(ac[i + 1],6);
+ r = -frac_div32(SHL32(rr,6), error);
+ /* Update LPC coefficients and total error */
+ lpc[i] = SHR32(r,6);
+ for (j = 0; j < (i+1)>>1; j++)
+ {
+ opus_val32 tmp1, tmp2;
+ tmp1 = lpc[j];
+ tmp2 = lpc[i-1-j];
+ lpc[j] = tmp1 + MULT32_32_Q31(r,tmp2);
+ lpc[i-1-j] = tmp2 + MULT32_32_Q31(r,tmp1);
+ }
+
+ error = error - MULT32_32_Q31(MULT32_32_Q31(r,r),error);
+ /* Bail out once we get 30 dB gain */
+#ifdef FIXED_POINT
+ if (error<=SHR32(ac[0],10))
+ break;
+#else
+ if (error<=.001f*ac[0])
+ break;
+#endif
+ }
+ }
+#ifdef FIXED_POINT
+ {
+ /* Convert the int32 lpcs to int16 and ensure there are no wrap-arounds.
+ This reuses the logic in silk_LPC_fit() and silk_bwexpander_32(). Any bug
+ fixes should also be applied there. */
+ int iter, idx = 0;
+ opus_val32 maxabs, absval, chirp_Q16, chirp_minus_one_Q16;
+
+ for (iter = 0; iter < 10; iter++) {
+ maxabs = 0;
+ for (i = 0; i < p; i++) {
+ absval = ABS32(lpc[i]);
+ if (absval > maxabs) {
+ maxabs = absval;
+ idx = i;
+ }
+ }
+ maxabs = PSHR32(maxabs, 13); /* Q25->Q12 */
+
+ if (maxabs > 32767) {
+ maxabs = MIN32(maxabs, 163838);
+ chirp_Q16 = QCONST32(0.999, 16) - DIV32(SHL32(maxabs - 32767, 14),
+ SHR32(MULT32_32_32(maxabs, idx + 1), 2));
+ chirp_minus_one_Q16 = chirp_Q16 - 65536;
+
+ /* Apply bandwidth expansion. */
+ for (i = 0; i < p - 1; i++) {
+ lpc[i] = MULT32_32_Q16(chirp_Q16, lpc[i]);
+ chirp_Q16 += PSHR32(MULT32_32_32(chirp_Q16, chirp_minus_one_Q16), 16);
+ }
+ lpc[p - 1] = MULT32_32_Q16(chirp_Q16, lpc[p - 1]);
+ } else {
+ break;
+ }
+ }
+
+ if (iter == 10) {
+ /* If the coeffs still do not fit into the 16 bit range after 10 iterations,
+ fall back to the A(z)=1 filter. */
+ OPUS_CLEAR(lpc, p);
+ _lpc[0] = 4096; /* Q12 */
+ } else {
+ for (i = 0; i < p; i++) {
+ _lpc[i] = EXTRACT16(PSHR32(lpc[i], 13)); /* Q25->Q12 */
+ }
+ }
+ }
+#endif
+}
+
+
+void celt_fir_c(
+ const opus_val16 *x,
+ const opus_val16 *num,
+ opus_val16 *y,
+ int N,
+ int ord,
+ int arch)
+{
+ int i,j;
+ VARDECL(opus_val16, rnum);
+ SAVE_STACK;
+ celt_assert(x != y);
+ ALLOC(rnum, ord, opus_val16);
+ for(i=0;i<ord;i++)
+ rnum[i] = num[ord-i-1];
+ for (i=0;i<N-3;i+=4)
+ {
+ opus_val32 sum[4];
+ sum[0] = SHL32(EXTEND32(x[i ]), SIG_SHIFT);
+ sum[1] = SHL32(EXTEND32(x[i+1]), SIG_SHIFT);
+ sum[2] = SHL32(EXTEND32(x[i+2]), SIG_SHIFT);
+ sum[3] = SHL32(EXTEND32(x[i+3]), SIG_SHIFT);
+ xcorr_kernel(rnum, x+i-ord, sum, ord, arch);
+ y[i ] = SROUND16(sum[0], SIG_SHIFT);
+ y[i+1] = SROUND16(sum[1], SIG_SHIFT);
+ y[i+2] = SROUND16(sum[2], SIG_SHIFT);
+ y[i+3] = SROUND16(sum[3], SIG_SHIFT);
+ }
+ for (;i<N;i++)
+ {
+ opus_val32 sum = SHL32(EXTEND32(x[i]), SIG_SHIFT);
+ for (j=0;j<ord;j++)
+ sum = MAC16_16(sum,rnum[j],x[i+j-ord]);
+ y[i] = SROUND16(sum, SIG_SHIFT);
+ }
+ RESTORE_STACK;
+}
+
+void celt_iir(const opus_val32 *_x,
+ const opus_val16 *den,
+ opus_val32 *_y,
+ int N,
+ int ord,
+ opus_val16 *mem,
+ int arch)
+{
+#ifdef SMALL_FOOTPRINT
+ int i,j;
+ (void)arch;
+ for (i=0;i<N;i++)
+ {
+ opus_val32 sum = _x[i];
+ for (j=0;j<ord;j++)
+ {
+ sum -= MULT16_16(den[j],mem[j]);
+ }
+ for (j=ord-1;j>=1;j--)
+ {
+ mem[j]=mem[j-1];
+ }
+ mem[0] = SROUND16(sum, SIG_SHIFT);
+ _y[i] = sum;
+ }
+#else
+ int i,j;
+ VARDECL(opus_val16, rden);
+ VARDECL(opus_val16, y);
+ SAVE_STACK;
+
+ celt_assert((ord&3)==0);
+ ALLOC(rden, ord, opus_val16);
+ ALLOC(y, N+ord, opus_val16);
+ for(i=0;i<ord;i++)
+ rden[i] = den[ord-i-1];
+ for(i=0;i<ord;i++)
+ y[i] = -mem[ord-i-1];
+ for(;i<N+ord;i++)
+ y[i]=0;
+ for (i=0;i<N-3;i+=4)
+ {
+ /* Unroll by 4 as if it were an FIR filter */
+ opus_val32 sum[4];
+ sum[0]=_x[i];
+ sum[1]=_x[i+1];
+ sum[2]=_x[i+2];
+ sum[3]=_x[i+3];
+ xcorr_kernel(rden, y+i, sum, ord, arch);
+
+ /* Patch up the result to compensate for the fact that this is an IIR */
+ y[i+ord ] = -SROUND16(sum[0],SIG_SHIFT);
+ _y[i ] = sum[0];
+ sum[1] = MAC16_16(sum[1], y[i+ord ], den[0]);
+ y[i+ord+1] = -SROUND16(sum[1],SIG_SHIFT);
+ _y[i+1] = sum[1];
+ sum[2] = MAC16_16(sum[2], y[i+ord+1], den[0]);
+ sum[2] = MAC16_16(sum[2], y[i+ord ], den[1]);
+ y[i+ord+2] = -SROUND16(sum[2],SIG_SHIFT);
+ _y[i+2] = sum[2];
+
+ sum[3] = MAC16_16(sum[3], y[i+ord+2], den[0]);
+ sum[3] = MAC16_16(sum[3], y[i+ord+1], den[1]);
+ sum[3] = MAC16_16(sum[3], y[i+ord ], den[2]);
+ y[i+ord+3] = -SROUND16(sum[3],SIG_SHIFT);
+ _y[i+3] = sum[3];
+ }
+ for (;i<N;i++)
+ {
+ opus_val32 sum = _x[i];
+ for (j=0;j<ord;j++)
+ sum -= MULT16_16(rden[j],y[i+j]);
+ y[i+ord] = SROUND16(sum,SIG_SHIFT);
+ _y[i] = sum;
+ }
+ for(i=0;i<ord;i++)
+ mem[i] = _y[N-i-1];
+ RESTORE_STACK;
+#endif
+}
+
+int _celt_autocorr(
+ const opus_val16 *x, /* in: [0...n-1] samples x */
+ opus_val32 *ac, /* out: [0...lag-1] ac values */
+ const opus_val16 *window,
+ int overlap,
+ int lag,
+ int n,
+ int arch
+ )
+{
+ opus_val32 d;
+ int i, k;
+ int fastN=n-lag;
+ int shift;
+ const opus_val16 *xptr;
+ VARDECL(opus_val16, xx);
+ SAVE_STACK;
+ ALLOC(xx, n, opus_val16);
+ celt_assert(n>0);
+ celt_assert(overlap>=0);
+ if (overlap == 0)
+ {
+ xptr = x;
+ } else {
+ for (i=0;i<n;i++)
+ xx[i] = x[i];
+ for (i=0;i<overlap;i++)
+ {
+ xx[i] = MULT16_16_Q15(x[i],window[i]);
+ xx[n-i-1] = MULT16_16_Q15(x[n-i-1],window[i]);
+ }
+ xptr = xx;
+ }
+ shift=0;
+#ifdef FIXED_POINT
+ {
+ opus_val32 ac0;
+ ac0 = 1+(n<<7);
+ if (n&1) ac0 += SHR32(MULT16_16(xptr[0],xptr[0]),9);
+ for(i=(n&1);i<n;i+=2)
+ {
+ ac0 += SHR32(MULT16_16(xptr[i],xptr[i]),9);
+ ac0 += SHR32(MULT16_16(xptr[i+1],xptr[i+1]),9);
+ }
+
+ shift = celt_ilog2(ac0)-30+10;
+ shift = (shift)/2;
+ if (shift>0)
+ {
+ for(i=0;i<n;i++)
+ xx[i] = PSHR32(xptr[i], shift);
+ xptr = xx;
+ } else
+ shift = 0;
+ }
+#endif
+ celt_pitch_xcorr(xptr, xptr, ac, fastN, lag+1, arch);
+ for (k=0;k<=lag;k++)
+ {
+ for (i = k+fastN, d = 0; i < n; i++)
+ d = MAC16_16(d, xptr[i], xptr[i-k]);
+ ac[k] += d;
+ }
+#ifdef FIXED_POINT
+ shift = 2*shift;
+ if (shift<=0)
+ ac[0] += SHL32((opus_int32)1, -shift);
+ if (ac[0] < 268435456)
+ {
+ int shift2 = 29 - EC_ILOG(ac[0]);
+ for (i=0;i<=lag;i++)
+ ac[i] = SHL32(ac[i], shift2);
+ shift -= shift2;
+ } else if (ac[0] >= 536870912)
+ {
+ int shift2=1;
+ if (ac[0] >= 1073741824)
+ shift2++;
+ for (i=0;i<=lag;i++)
+ ac[i] = SHR32(ac[i], shift2);
+ shift += shift2;
+ }
+#endif
+
+ RESTORE_STACK;
+ return shift;
+}
generated by cgit v1.2.3 (git 2.39.1) at 2024-11-05 05:29:50 +0000