cannam@154: /* Copyright (c) 2009-2010 Xiph.Org Foundation
cannam@154:    Written by Jean-Marc Valin */
cannam@154: /*
cannam@154:    Redistribution and use in source and binary forms, with or without
cannam@154:    modification, are permitted provided that the following conditions
cannam@154:    are met:
cannam@154: 
cannam@154:    - Redistributions of source code must retain the above copyright
cannam@154:    notice, this list of conditions and the following disclaimer.
cannam@154: 
cannam@154:    - Redistributions in binary form must reproduce the above copyright
cannam@154:    notice, this list of conditions and the following disclaimer in the
cannam@154:    documentation and/or other materials provided with the distribution.
cannam@154: 
cannam@154:    THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
cannam@154:    ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
cannam@154:    LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
cannam@154:    A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER
cannam@154:    OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
cannam@154:    EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
cannam@154:    PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
cannam@154:    PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
cannam@154:    LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
cannam@154:    NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
cannam@154:    SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
cannam@154: */
cannam@154: 
cannam@154: #ifdef HAVE_CONFIG_H
cannam@154: #include "config.h"
cannam@154: #endif
cannam@154: 
cannam@154: #include "celt_lpc.h"
cannam@154: #include "stack_alloc.h"
cannam@154: #include "mathops.h"
cannam@154: #include "pitch.h"
cannam@154: 
cannam@154: void _celt_lpc(
cannam@154:       opus_val16       *_lpc, /* out: [0...p-1] LPC coefficients      */
cannam@154: const opus_val32 *ac,  /* in:  [0...p] autocorrelation values  */
cannam@154: int          p
cannam@154: )
cannam@154: {
cannam@154:    int i, j;
cannam@154:    opus_val32 r;
cannam@154:    opus_val32 error = ac[0];
cannam@154: #ifdef FIXED_POINT
cannam@154:    opus_val32 lpc[LPC_ORDER];
cannam@154: #else
cannam@154:    float *lpc = _lpc;
cannam@154: #endif
cannam@154: 
cannam@154:    OPUS_CLEAR(lpc, p);
cannam@154:    if (ac[0] != 0)
cannam@154:    {
cannam@154:       for (i = 0; i < p; i++) {
cannam@154:          /* Sum up this iteration's reflection coefficient */
cannam@154:          opus_val32 rr = 0;
cannam@154:          for (j = 0; j < i; j++)
cannam@154:             rr += MULT32_32_Q31(lpc[j],ac[i - j]);
cannam@154:          rr += SHR32(ac[i + 1],3);
cannam@154:          r = -frac_div32(SHL32(rr,3), error);
cannam@154:          /*  Update LPC coefficients and total error */
cannam@154:          lpc[i] = SHR32(r,3);
cannam@154:          for (j = 0; j < (i+1)>>1; j++)
cannam@154:          {
cannam@154:             opus_val32 tmp1, tmp2;
cannam@154:             tmp1 = lpc[j];
cannam@154:             tmp2 = lpc[i-1-j];
cannam@154:             lpc[j]     = tmp1 + MULT32_32_Q31(r,tmp2);
cannam@154:             lpc[i-1-j] = tmp2 + MULT32_32_Q31(r,tmp1);
cannam@154:          }
cannam@154: 
cannam@154:          error = error - MULT32_32_Q31(MULT32_32_Q31(r,r),error);
cannam@154:          /* Bail out once we get 30 dB gain */
cannam@154: #ifdef FIXED_POINT
cannam@154:          if (error<SHR32(ac[0],10))
cannam@154:             break;
cannam@154: #else
cannam@154:          if (error<.001f*ac[0])
cannam@154:             break;
cannam@154: #endif
cannam@154:       }
cannam@154:    }
cannam@154: #ifdef FIXED_POINT
cannam@154:    for (i=0;i<p;i++)
cannam@154:       _lpc[i] = ROUND16(lpc[i],16);
cannam@154: #endif
cannam@154: }
cannam@154: 
cannam@154: 
cannam@154: void celt_fir_c(
cannam@154:          const opus_val16 *x,
cannam@154:          const opus_val16 *num,
cannam@154:          opus_val16 *y,
cannam@154:          int N,
cannam@154:          int ord,
cannam@154:          int arch)
cannam@154: {
cannam@154:    int i,j;
cannam@154:    VARDECL(opus_val16, rnum);
cannam@154:    SAVE_STACK;
cannam@154:    celt_assert(x != y);
cannam@154:    ALLOC(rnum, ord, opus_val16);
cannam@154:    for(i=0;i<ord;i++)
cannam@154:       rnum[i] = num[ord-i-1];
cannam@154:    for (i=0;i<N-3;i+=4)
cannam@154:    {
cannam@154:       opus_val32 sum[4];
cannam@154:       sum[0] = SHL32(EXTEND32(x[i  ]), SIG_SHIFT);
cannam@154:       sum[1] = SHL32(EXTEND32(x[i+1]), SIG_SHIFT);
cannam@154:       sum[2] = SHL32(EXTEND32(x[i+2]), SIG_SHIFT);
cannam@154:       sum[3] = SHL32(EXTEND32(x[i+3]), SIG_SHIFT);
cannam@154:       xcorr_kernel(rnum, x+i-ord, sum, ord, arch);
cannam@154:       y[i  ] = ROUND16(sum[0], SIG_SHIFT);
cannam@154:       y[i+1] = ROUND16(sum[1], SIG_SHIFT);
cannam@154:       y[i+2] = ROUND16(sum[2], SIG_SHIFT);
cannam@154:       y[i+3] = ROUND16(sum[3], SIG_SHIFT);
cannam@154:    }
cannam@154:    for (;i<N;i++)
cannam@154:    {
cannam@154:       opus_val32 sum = SHL32(EXTEND32(x[i]), SIG_SHIFT);
cannam@154:       for (j=0;j<ord;j++)
cannam@154:          sum = MAC16_16(sum,rnum[j],x[i+j-ord]);
cannam@154:       y[i] = ROUND16(sum, SIG_SHIFT);
cannam@154:    }
cannam@154:    RESTORE_STACK;
cannam@154: }
cannam@154: 
cannam@154: void celt_iir(const opus_val32 *_x,
cannam@154:          const opus_val16 *den,
cannam@154:          opus_val32 *_y,
cannam@154:          int N,
cannam@154:          int ord,
cannam@154:          opus_val16 *mem,
cannam@154:          int arch)
cannam@154: {
cannam@154: #ifdef SMALL_FOOTPRINT
cannam@154:    int i,j;
cannam@154:    (void)arch;
cannam@154:    for (i=0;i<N;i++)
cannam@154:    {
cannam@154:       opus_val32 sum = _x[i];
cannam@154:       for (j=0;j<ord;j++)
cannam@154:       {
cannam@154:          sum -= MULT16_16(den[j],mem[j]);
cannam@154:       }
cannam@154:       for (j=ord-1;j>=1;j--)
cannam@154:       {
cannam@154:          mem[j]=mem[j-1];
cannam@154:       }
cannam@154:       mem[0] = SROUND16(sum, SIG_SHIFT);
cannam@154:       _y[i] = sum;
cannam@154:    }
cannam@154: #else
cannam@154:    int i,j;
cannam@154:    VARDECL(opus_val16, rden);
cannam@154:    VARDECL(opus_val16, y);
cannam@154:    SAVE_STACK;
cannam@154: 
cannam@154:    celt_assert((ord&3)==0);
cannam@154:    ALLOC(rden, ord, opus_val16);
cannam@154:    ALLOC(y, N+ord, opus_val16);
cannam@154:    for(i=0;i<ord;i++)
cannam@154:       rden[i] = den[ord-i-1];
cannam@154:    for(i=0;i<ord;i++)
cannam@154:       y[i] = -mem[ord-i-1];
cannam@154:    for(;i<N+ord;i++)
cannam@154:       y[i]=0;
cannam@154:    for (i=0;i<N-3;i+=4)
cannam@154:    {
cannam@154:       /* Unroll by 4 as if it were an FIR filter */
cannam@154:       opus_val32 sum[4];
cannam@154:       sum[0]=_x[i];
cannam@154:       sum[1]=_x[i+1];
cannam@154:       sum[2]=_x[i+2];
cannam@154:       sum[3]=_x[i+3];
cannam@154:       xcorr_kernel(rden, y+i, sum, ord, arch);
cannam@154: 
cannam@154:       /* Patch up the result to compensate for the fact that this is an IIR */
cannam@154:       y[i+ord  ] = -SROUND16(sum[0],SIG_SHIFT);
cannam@154:       _y[i  ] = sum[0];
cannam@154:       sum[1] = MAC16_16(sum[1], y[i+ord  ], den[0]);
cannam@154:       y[i+ord+1] = -SROUND16(sum[1],SIG_SHIFT);
cannam@154:       _y[i+1] = sum[1];
cannam@154:       sum[2] = MAC16_16(sum[2], y[i+ord+1], den[0]);
cannam@154:       sum[2] = MAC16_16(sum[2], y[i+ord  ], den[1]);
cannam@154:       y[i+ord+2] = -SROUND16(sum[2],SIG_SHIFT);
cannam@154:       _y[i+2] = sum[2];
cannam@154: 
cannam@154:       sum[3] = MAC16_16(sum[3], y[i+ord+2], den[0]);
cannam@154:       sum[3] = MAC16_16(sum[3], y[i+ord+1], den[1]);
cannam@154:       sum[3] = MAC16_16(sum[3], y[i+ord  ], den[2]);
cannam@154:       y[i+ord+3] = -SROUND16(sum[3],SIG_SHIFT);
cannam@154:       _y[i+3] = sum[3];
cannam@154:    }
cannam@154:    for (;i<N;i++)
cannam@154:    {
cannam@154:       opus_val32 sum = _x[i];
cannam@154:       for (j=0;j<ord;j++)
cannam@154:          sum -= MULT16_16(rden[j],y[i+j]);
cannam@154:       y[i+ord] = SROUND16(sum,SIG_SHIFT);
cannam@154:       _y[i] = sum;
cannam@154:    }
cannam@154:    for(i=0;i<ord;i++)
cannam@154:       mem[i] = _y[N-i-1];
cannam@154:    RESTORE_STACK;
cannam@154: #endif
cannam@154: }
cannam@154: 
cannam@154: int _celt_autocorr(
cannam@154:                    const opus_val16 *x,   /*  in: [0...n-1] samples x   */
cannam@154:                    opus_val32       *ac,  /* out: [0...lag-1] ac values */
cannam@154:                    const opus_val16       *window,
cannam@154:                    int          overlap,
cannam@154:                    int          lag,
cannam@154:                    int          n,
cannam@154:                    int          arch
cannam@154:                   )
cannam@154: {
cannam@154:    opus_val32 d;
cannam@154:    int i, k;
cannam@154:    int fastN=n-lag;
cannam@154:    int shift;
cannam@154:    const opus_val16 *xptr;
cannam@154:    VARDECL(opus_val16, xx);
cannam@154:    SAVE_STACK;
cannam@154:    ALLOC(xx, n, opus_val16);
cannam@154:    celt_assert(n>0);
cannam@154:    celt_assert(overlap>=0);
cannam@154:    if (overlap == 0)
cannam@154:    {
cannam@154:       xptr = x;
cannam@154:    } else {
cannam@154:       for (i=0;i<n;i++)
cannam@154:          xx[i] = x[i];
cannam@154:       for (i=0;i<overlap;i++)
cannam@154:       {
cannam@154:          xx[i] = MULT16_16_Q15(x[i],window[i]);
cannam@154:          xx[n-i-1] = MULT16_16_Q15(x[n-i-1],window[i]);
cannam@154:       }
cannam@154:       xptr = xx;
cannam@154:    }
cannam@154:    shift=0;
cannam@154: #ifdef FIXED_POINT
cannam@154:    {
cannam@154:       opus_val32 ac0;
cannam@154:       ac0 = 1+(n<<7);
cannam@154:       if (n&1) ac0 += SHR32(MULT16_16(xptr[0],xptr[0]),9);
cannam@154:       for(i=(n&1);i<n;i+=2)
cannam@154:       {
cannam@154:          ac0 += SHR32(MULT16_16(xptr[i],xptr[i]),9);
cannam@154:          ac0 += SHR32(MULT16_16(xptr[i+1],xptr[i+1]),9);
cannam@154:       }
cannam@154: 
cannam@154:       shift = celt_ilog2(ac0)-30+10;
cannam@154:       shift = (shift)/2;
cannam@154:       if (shift>0)
cannam@154:       {
cannam@154:          for(i=0;i<n;i++)
cannam@154:             xx[i] = PSHR32(xptr[i], shift);
cannam@154:          xptr = xx;
cannam@154:       } else
cannam@154:          shift = 0;
cannam@154:    }
cannam@154: #endif
cannam@154:    celt_pitch_xcorr(xptr, xptr, ac, fastN, lag+1, arch);
cannam@154:    for (k=0;k<=lag;k++)
cannam@154:    {
cannam@154:       for (i = k+fastN, d = 0; i < n; i++)
cannam@154:          d = MAC16_16(d, xptr[i], xptr[i-k]);
cannam@154:       ac[k] += d;
cannam@154:    }
cannam@154: #ifdef FIXED_POINT
cannam@154:    shift = 2*shift;
cannam@154:    if (shift<=0)
cannam@154:       ac[0] += SHL32((opus_int32)1, -shift);
cannam@154:    if (ac[0] < 268435456)
cannam@154:    {
cannam@154:       int shift2 = 29 - EC_ILOG(ac[0]);
cannam@154:       for (i=0;i<=lag;i++)
cannam@154:          ac[i] = SHL32(ac[i], shift2);
cannam@154:       shift -= shift2;
cannam@154:    } else if (ac[0] >= 536870912)
cannam@154:    {
cannam@154:       int shift2=1;
cannam@154:       if (ac[0] >= 1073741824)
cannam@154:          shift2++;
cannam@154:       for (i=0;i<=lag;i++)
cannam@154:          ac[i] = SHR32(ac[i], shift2);
cannam@154:       shift += shift2;
cannam@154:    }
cannam@154: #endif
cannam@154: 
cannam@154:    RESTORE_STACK;
cannam@154:    return shift;
cannam@154: }