190 lines
		
	
	
		
			4.8 KiB
		
	
	
	
		
			C
		
	
	
	
	
	
		
		
			
		
	
	
			190 lines
		
	
	
		
			4.8 KiB
		
	
	
	
		
			C
		
	
	
	
	
	
| 
								 | 
							
								/* 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 FOUNDATION 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 "plc.h"
							 | 
						||
| 
								 | 
							
								#include "stack_alloc.h"
							 | 
						||
| 
								 | 
							
								#include "mathops.h"
							 | 
						||
| 
								 | 
							
								
							 | 
						||
| 
								 | 
							
								
							 | 
						||
| 
								 | 
							
								
							 | 
						||
| 
								 | 
							
								
							 | 
						||
| 
								 | 
							
								void _celt_lpc(
							 | 
						||
| 
								 | 
							
								      celt_word16       *_lpc, /* out: [0...p-1] LPC coefficients      */
							 | 
						||
| 
								 | 
							
								const celt_word32 *ac,  /* in:  [0...p] autocorrelation values  */
							 | 
						||
| 
								 | 
							
								int          p
							 | 
						||
| 
								 | 
							
								)
							 | 
						||
| 
								 | 
							
								{
							 | 
						||
| 
								 | 
							
								   int i, j;  
							 | 
						||
| 
								 | 
							
								   celt_word32 r;
							 | 
						||
| 
								 | 
							
								   celt_word32 error = ac[0];
							 | 
						||
| 
								 | 
							
								#ifdef FIXED_POINT
							 | 
						||
| 
								 | 
							
								   celt_word32 lpc[LPC_ORDER];
							 | 
						||
| 
								 | 
							
								#else
							 | 
						||
| 
								 | 
							
								   float *lpc = _lpc;
							 | 
						||
| 
								 | 
							
								#endif
							 | 
						||
| 
								 | 
							
								
							 | 
						||
| 
								 | 
							
								   for (i = 0; i < p; i++)
							 | 
						||
| 
								 | 
							
								      lpc[i] = 0;
							 | 
						||
| 
								 | 
							
								   if (ac[0] != 0)
							 | 
						||
| 
								 | 
							
								   {
							 | 
						||
| 
								 | 
							
								      for (i = 0; i < p; i++) {
							 | 
						||
| 
								 | 
							
								         /* Sum up this iteration's reflection coefficient */
							 | 
						||
| 
								 | 
							
								         celt_word32 rr = 0;
							 | 
						||
| 
								 | 
							
								         for (j = 0; j < i; j++)
							 | 
						||
| 
								 | 
							
								            rr += MULT32_32_Q31(lpc[j],ac[i - j]);
							 | 
						||
| 
								 | 
							
								         rr += SHR32(ac[i + 1],3);
							 | 
						||
| 
								 | 
							
								         r = -frac_div32(SHL32(rr,3), error);
							 | 
						||
| 
								 | 
							
								         /*  Update LPC coefficients and total error */
							 | 
						||
| 
								 | 
							
								         lpc[i] = SHR32(r,3);
							 | 
						||
| 
								 | 
							
								         for (j = 0; j < (i+1)>>1; j++)
							 | 
						||
| 
								 | 
							
								         {
							 | 
						||
| 
								 | 
							
								            celt_word32 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
							 | 
						||
| 
								 | 
							
								   for (i=0;i<p;i++)
							 | 
						||
| 
								 | 
							
								      _lpc[i] = ROUND16(lpc[i],16);
							 | 
						||
| 
								 | 
							
								#endif
							 | 
						||
| 
								 | 
							
								}
							 | 
						||
| 
								 | 
							
								
							 | 
						||
| 
								 | 
							
								void fir(const celt_word16 *x,
							 | 
						||
| 
								 | 
							
								         const celt_word16 *num,
							 | 
						||
| 
								 | 
							
								         celt_word16 *y,
							 | 
						||
| 
								 | 
							
								         int N,
							 | 
						||
| 
								 | 
							
								         int ord,
							 | 
						||
| 
								 | 
							
								         celt_word16 *mem)
							 | 
						||
| 
								 | 
							
								{
							 | 
						||
| 
								 | 
							
								   int i,j;
							 | 
						||
| 
								 | 
							
								
							 | 
						||
| 
								 | 
							
								   for (i=0;i<N;i++)
							 | 
						||
| 
								 | 
							
								   {
							 | 
						||
| 
								 | 
							
								      celt_word32 sum = SHL32(EXTEND32(x[i]), SIG_SHIFT);
							 | 
						||
| 
								 | 
							
								      for (j=0;j<ord;j++)
							 | 
						||
| 
								 | 
							
								      {
							 | 
						||
| 
								 | 
							
								         sum += MULT16_16(num[j],mem[j]);
							 | 
						||
| 
								 | 
							
								      }
							 | 
						||
| 
								 | 
							
								      for (j=ord-1;j>=1;j--)
							 | 
						||
| 
								 | 
							
								      {
							 | 
						||
| 
								 | 
							
								         mem[j]=mem[j-1];
							 | 
						||
| 
								 | 
							
								      }
							 | 
						||
| 
								 | 
							
								      mem[0] = x[i];
							 | 
						||
| 
								 | 
							
								      y[i] = ROUND16(sum, SIG_SHIFT);
							 | 
						||
| 
								 | 
							
								   }
							 | 
						||
| 
								 | 
							
								}
							 | 
						||
| 
								 | 
							
								
							 | 
						||
| 
								 | 
							
								void iir(const celt_word32 *x,
							 | 
						||
| 
								 | 
							
								         const celt_word16 *den,
							 | 
						||
| 
								 | 
							
								         celt_word32 *y,
							 | 
						||
| 
								 | 
							
								         int N,
							 | 
						||
| 
								 | 
							
								         int ord,
							 | 
						||
| 
								 | 
							
								         celt_word16 *mem)
							 | 
						||
| 
								 | 
							
								{
							 | 
						||
| 
								 | 
							
								   int i,j;
							 | 
						||
| 
								 | 
							
								   for (i=0;i<N;i++)
							 | 
						||
| 
								 | 
							
								   {
							 | 
						||
| 
								 | 
							
								      celt_word32 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] = ROUND16(sum,SIG_SHIFT);
							 | 
						||
| 
								 | 
							
								      y[i] = sum;
							 | 
						||
| 
								 | 
							
								   }
							 | 
						||
| 
								 | 
							
								}
							 | 
						||
| 
								 | 
							
								
							 | 
						||
| 
								 | 
							
								void _celt_autocorr(
							 | 
						||
| 
								 | 
							
								                   const celt_word16 *x,   /*  in: [0...n-1] samples x   */
							 | 
						||
| 
								 | 
							
								                   celt_word32       *ac,  /* out: [0...lag-1] ac values */
							 | 
						||
| 
								 | 
							
								                   const celt_word16       *window,
							 | 
						||
| 
								 | 
							
								                   int          overlap,
							 | 
						||
| 
								 | 
							
								                   int          lag, 
							 | 
						||
| 
								 | 
							
								                   int          n
							 | 
						||
| 
								 | 
							
								                  )
							 | 
						||
| 
								 | 
							
								{
							 | 
						||
| 
								 | 
							
								   celt_word32 d;
							 | 
						||
| 
								 | 
							
								   int i;
							 | 
						||
| 
								 | 
							
								   VARDECL(celt_word16, xx);
							 | 
						||
| 
								 | 
							
								   SAVE_STACK;
							 | 
						||
| 
								 | 
							
								   ALLOC(xx, n, celt_word16);
							 | 
						||
| 
								 | 
							
								   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]);
							 | 
						||
| 
								 | 
							
								   }
							 | 
						||
| 
								 | 
							
								#ifdef FIXED_POINT
							 | 
						||
| 
								 | 
							
								   {
							 | 
						||
| 
								 | 
							
								      celt_word32 ac0=0;
							 | 
						||
| 
								 | 
							
								      int shift;
							 | 
						||
| 
								 | 
							
								      for(i=0;i<n;i++)
							 | 
						||
| 
								 | 
							
								         ac0 += SHR32(MULT16_16(xx[i],xx[i]),8);
							 | 
						||
| 
								 | 
							
								      ac0 += 1+n;
							 | 
						||
| 
								 | 
							
								
							 | 
						||
| 
								 | 
							
								      shift = celt_ilog2(ac0)-30+9;
							 | 
						||
| 
								 | 
							
								      shift = (shift+1)/2;
							 | 
						||
| 
								 | 
							
								      for(i=0;i<n;i++)
							 | 
						||
| 
								 | 
							
								         xx[i] = VSHR32(xx[i], shift);
							 | 
						||
| 
								 | 
							
								   }
							 | 
						||
| 
								 | 
							
								#endif
							 | 
						||
| 
								 | 
							
								   while (lag>=0)
							 | 
						||
| 
								 | 
							
								   {
							 | 
						||
| 
								 | 
							
								      for (i = lag, d = 0; i < n; i++) 
							 | 
						||
| 
								 | 
							
								         d += xx[i] * xx[i-lag];
							 | 
						||
| 
								 | 
							
								      ac[lag] = d;
							 | 
						||
| 
								 | 
							
								      /*printf ("%f ", ac[lag]);*/
							 | 
						||
| 
								 | 
							
								      lag--;
							 | 
						||
| 
								 | 
							
								   }
							 | 
						||
| 
								 | 
							
								   /*printf ("\n");*/
							 | 
						||
| 
								 | 
							
								   ac[0] += 10;
							 | 
						||
| 
								 | 
							
								
							 | 
						||
| 
								 | 
							
								   RESTORE_STACK;
							 | 
						||
| 
								 | 
							
								}
							 |