Correct many whitespace errors under libcelt/ and remove
[opus.git] / libcelt / plc.c
1 /* Copyright (c) 2009-2010 Xiph.Org Foundation
2    Written by Jean-Marc Valin */
3 /*
4    Redistribution and use in source and binary forms, with or without
5    modification, are permitted provided that the following conditions
6    are met:
7
8    - Redistributions of source code must retain the above copyright
9    notice, this list of conditions and the following disclaimer.
10
11    - Redistributions in binary form must reproduce the above copyright
12    notice, this list of conditions and the following disclaimer in the
13    documentation and/or other materials provided with the distribution.
14
15    THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
16    ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
17    LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
18    A PARTICULAR PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE FOUNDATION OR
19    CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
20    EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
21    PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
22    PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
23    LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
24    NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
25    SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
26 */
27
28 #ifdef HAVE_CONFIG_H
29 #include "config.h"
30 #endif
31
32 #include "plc.h"
33 #include "stack_alloc.h"
34 #include "mathops.h"
35
36
37
38
39 void _celt_lpc(
40       opus_val16       *_lpc, /* out: [0...p-1] LPC coefficients      */
41 const opus_val32 *ac,  /* in:  [0...p] autocorrelation values  */
42 int          p
43 )
44 {
45    int i, j;
46    opus_val32 r;
47    opus_val32 error = ac[0];
48 #ifdef FIXED_POINT
49    opus_val32 lpc[LPC_ORDER];
50 #else
51    float *lpc = _lpc;
52 #endif
53
54    for (i = 0; i < p; i++)
55       lpc[i] = 0;
56    if (ac[0] != 0)
57    {
58       for (i = 0; i < p; i++) {
59          /* Sum up this iteration's reflection coefficient */
60          opus_val32 rr = 0;
61          for (j = 0; j < i; j++)
62             rr += MULT32_32_Q31(lpc[j],ac[i - j]);
63          rr += SHR32(ac[i + 1],3);
64          r = -frac_div32(SHL32(rr,3), error);
65          /*  Update LPC coefficients and total error */
66          lpc[i] = SHR32(r,3);
67          for (j = 0; j < (i+1)>>1; j++)
68          {
69             opus_val32 tmp1, tmp2;
70             tmp1 = lpc[j];
71             tmp2 = lpc[i-1-j];
72             lpc[j]     = tmp1 + MULT32_32_Q31(r,tmp2);
73             lpc[i-1-j] = tmp2 + MULT32_32_Q31(r,tmp1);
74          }
75
76          error = error - MULT32_32_Q31(MULT32_32_Q31(r,r),error);
77          /* Bail out once we get 30 dB gain */
78 #ifdef FIXED_POINT
79          if (error<SHR32(ac[0],10))
80             break;
81 #else
82          if (error<.001f*ac[0])
83             break;
84 #endif
85       }
86    }
87 #ifdef FIXED_POINT
88    for (i=0;i<p;i++)
89       _lpc[i] = ROUND16(lpc[i],16);
90 #endif
91 }
92
93 void fir(const opus_val16 *x,
94          const opus_val16 *num,
95          opus_val16 *y,
96          int N,
97          int ord,
98          opus_val16 *mem)
99 {
100    int i,j;
101
102    for (i=0;i<N;i++)
103    {
104       opus_val32 sum = SHL32(EXTEND32(x[i]), SIG_SHIFT);
105       for (j=0;j<ord;j++)
106       {
107          sum += MULT16_16(num[j],mem[j]);
108       }
109       for (j=ord-1;j>=1;j--)
110       {
111          mem[j]=mem[j-1];
112       }
113       mem[0] = x[i];
114       y[i] = ROUND16(sum, SIG_SHIFT);
115    }
116 }
117
118 void iir(const opus_val32 *x,
119          const opus_val16 *den,
120          opus_val32 *y,
121          int N,
122          int ord,
123          opus_val16 *mem)
124 {
125    int i,j;
126    for (i=0;i<N;i++)
127    {
128       opus_val32 sum = x[i];
129       for (j=0;j<ord;j++)
130       {
131          sum -= MULT16_16(den[j],mem[j]);
132       }
133       for (j=ord-1;j>=1;j--)
134       {
135          mem[j]=mem[j-1];
136       }
137       mem[0] = ROUND16(sum,SIG_SHIFT);
138       y[i] = sum;
139    }
140 }
141
142 void _celt_autocorr(
143                    const opus_val16 *x,   /*  in: [0...n-1] samples x   */
144                    opus_val32       *ac,  /* out: [0...lag-1] ac values */
145                    const opus_val16       *window,
146                    int          overlap,
147                    int          lag,
148                    int          n
149                   )
150 {
151    opus_val32 d;
152    int i;
153    VARDECL(opus_val16, xx);
154    SAVE_STACK;
155    ALLOC(xx, n, opus_val16);
156    for (i=0;i<n;i++)
157       xx[i] = x[i];
158    for (i=0;i<overlap;i++)
159    {
160       xx[i] = MULT16_16_Q15(x[i],window[i]);
161       xx[n-i-1] = MULT16_16_Q15(x[n-i-1],window[i]);
162    }
163 #ifdef FIXED_POINT
164    {
165       opus_val32 ac0=0;
166       int shift;
167       for(i=0;i<n;i++)
168          ac0 += SHR32(MULT16_16(xx[i],xx[i]),8);
169       ac0 += 1+n;
170
171       shift = celt_ilog2(ac0)-30+9;
172       shift = (shift+1)/2;
173       for(i=0;i<n;i++)
174          xx[i] = VSHR32(xx[i], shift);
175    }
176 #endif
177    while (lag>=0)
178    {
179       for (i = lag, d = 0; i < n; i++)
180          d += xx[i] * xx[i-lag];
181       ac[lag] = d;
182       /*printf ("%f ", ac[lag]);*/
183       lag--;
184    }
185    /*printf ("\n");*/
186    ac[0] += 10;
187
188    RESTORE_STACK;
189 }