Improve OPUS_GET_LSB_DEPTH documentation.
[opus.git] / celt / celt_lpc.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 COPYRIGHT OWNER
19    OR 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 "celt_lpc.h"
33 #include "stack_alloc.h"
34 #include "mathops.h"
35 #include "pitch.h"
36
37 void _celt_lpc(
38       opus_val16       *_lpc, /* out: [0...p-1] LPC coefficients      */
39 const opus_val32 *ac,  /* in:  [0...p] autocorrelation values  */
40 int          p
41 )
42 {
43    int i, j;
44    opus_val32 r;
45    opus_val32 error = ac[0];
46 #ifdef FIXED_POINT
47    opus_val32 lpc[LPC_ORDER];
48 #else
49    float *lpc = _lpc;
50 #endif
51
52    for (i = 0; i < p; i++)
53       lpc[i] = 0;
54    if (ac[0] != 0)
55    {
56       for (i = 0; i < p; i++) {
57          /* Sum up this iteration's reflection coefficient */
58          opus_val32 rr = 0;
59          for (j = 0; j < i; j++)
60             rr += MULT32_32_Q31(lpc[j],ac[i - j]);
61          rr += SHR32(ac[i + 1],3);
62          r = -frac_div32(SHL32(rr,3), error);
63          /*  Update LPC coefficients and total error */
64          lpc[i] = SHR32(r,3);
65          for (j = 0; j < (i+1)>>1; j++)
66          {
67             opus_val32 tmp1, tmp2;
68             tmp1 = lpc[j];
69             tmp2 = lpc[i-1-j];
70             lpc[j]     = tmp1 + MULT32_32_Q31(r,tmp2);
71             lpc[i-1-j] = tmp2 + MULT32_32_Q31(r,tmp1);
72          }
73
74          error = error - MULT32_32_Q31(MULT32_32_Q31(r,r),error);
75          /* Bail out once we get 30 dB gain */
76 #ifdef FIXED_POINT
77          if (error<SHR32(ac[0],10))
78             break;
79 #else
80          if (error<.001f*ac[0])
81             break;
82 #endif
83       }
84    }
85 #ifdef FIXED_POINT
86    for (i=0;i<p;i++)
87       _lpc[i] = ROUND16(lpc[i],16);
88 #endif
89 }
90
91
92 void celt_fir_c(
93          const opus_val16 *_x,
94          const opus_val16 *num,
95          opus_val16 *_y,
96          int N,
97          int ord,
98          opus_val16 *mem,
99          int arch)
100 {
101    int i,j;
102    VARDECL(opus_val16, rnum);
103    VARDECL(opus_val16, x);
104    SAVE_STACK;
105
106    ALLOC(rnum, ord, opus_val16);
107    ALLOC(x, N+ord, opus_val16);
108    for(i=0;i<ord;i++)
109       rnum[i] = num[ord-i-1];
110    for(i=0;i<ord;i++)
111       x[i] = mem[ord-i-1];
112    for (i=0;i<N;i++)
113       x[i+ord]=_x[i];
114    for(i=0;i<ord;i++)
115       mem[i] = _x[N-i-1];
116 #ifdef SMALL_FOOTPRINT
117    (void)arch;
118    for (i=0;i<N;i++)
119    {
120       opus_val32 sum = SHL32(EXTEND32(_x[i]), SIG_SHIFT);
121       for (j=0;j<ord;j++)
122       {
123          sum = MAC16_16(sum,rnum[j],x[i+j]);
124       }
125       _y[i] = SATURATE16(PSHR32(sum, SIG_SHIFT));
126    }
127 #else
128    for (i=0;i<N-3;i+=4)
129    {
130       opus_val32 sum[4]={0,0,0,0};
131       xcorr_kernel(rnum, x+i, sum, ord, arch);
132       _y[i  ] = SATURATE16(ADD32(EXTEND32(_x[i  ]), PSHR32(sum[0], SIG_SHIFT)));
133       _y[i+1] = SATURATE16(ADD32(EXTEND32(_x[i+1]), PSHR32(sum[1], SIG_SHIFT)));
134       _y[i+2] = SATURATE16(ADD32(EXTEND32(_x[i+2]), PSHR32(sum[2], SIG_SHIFT)));
135       _y[i+3] = SATURATE16(ADD32(EXTEND32(_x[i+3]), PSHR32(sum[3], SIG_SHIFT)));
136    }
137    for (;i<N;i++)
138    {
139       opus_val32 sum = 0;
140       for (j=0;j<ord;j++)
141          sum = MAC16_16(sum,rnum[j],x[i+j]);
142       _y[i] = SATURATE16(ADD32(EXTEND32(_x[i]), PSHR32(sum, SIG_SHIFT)));
143    }
144 #endif
145    RESTORE_STACK;
146 }
147
148 void celt_iir(const opus_val32 *_x,
149          const opus_val16 *den,
150          opus_val32 *_y,
151          int N,
152          int ord,
153          opus_val16 *mem,
154          int arch)
155 {
156 #ifdef SMALL_FOOTPRINT
157    int i,j;
158    (void)arch;
159    for (i=0;i<N;i++)
160    {
161       opus_val32 sum = _x[i];
162       for (j=0;j<ord;j++)
163       {
164          sum -= MULT16_16(den[j],mem[j]);
165       }
166       for (j=ord-1;j>=1;j--)
167       {
168          mem[j]=mem[j-1];
169       }
170       mem[0] = ROUND16(sum,SIG_SHIFT);
171       _y[i] = sum;
172    }
173 #else
174    int i,j;
175    VARDECL(opus_val16, rden);
176    VARDECL(opus_val16, y);
177    SAVE_STACK;
178
179    celt_assert((ord&3)==0);
180    ALLOC(rden, ord, opus_val16);
181    ALLOC(y, N+ord, opus_val16);
182    for(i=0;i<ord;i++)
183       rden[i] = den[ord-i-1];
184    for(i=0;i<ord;i++)
185       y[i] = -mem[ord-i-1];
186    for(;i<N+ord;i++)
187       y[i]=0;
188    for (i=0;i<N-3;i+=4)
189    {
190       /* Unroll by 4 as if it were an FIR filter */
191       opus_val32 sum[4];
192       sum[0]=_x[i];
193       sum[1]=_x[i+1];
194       sum[2]=_x[i+2];
195       sum[3]=_x[i+3];
196       xcorr_kernel(rden, y+i, sum, ord, arch);
197
198       /* Patch up the result to compensate for the fact that this is an IIR */
199       y[i+ord  ] = -ROUND16(sum[0],SIG_SHIFT);
200       _y[i  ] = sum[0];
201       sum[1] = MAC16_16(sum[1], y[i+ord  ], den[0]);
202       y[i+ord+1] = -ROUND16(sum[1],SIG_SHIFT);
203       _y[i+1] = sum[1];
204       sum[2] = MAC16_16(sum[2], y[i+ord+1], den[0]);
205       sum[2] = MAC16_16(sum[2], y[i+ord  ], den[1]);
206       y[i+ord+2] = -ROUND16(sum[2],SIG_SHIFT);
207       _y[i+2] = sum[2];
208
209       sum[3] = MAC16_16(sum[3], y[i+ord+2], den[0]);
210       sum[3] = MAC16_16(sum[3], y[i+ord+1], den[1]);
211       sum[3] = MAC16_16(sum[3], y[i+ord  ], den[2]);
212       y[i+ord+3] = -ROUND16(sum[3],SIG_SHIFT);
213       _y[i+3] = sum[3];
214    }
215    for (;i<N;i++)
216    {
217       opus_val32 sum = _x[i];
218       for (j=0;j<ord;j++)
219          sum -= MULT16_16(rden[j],y[i+j]);
220       y[i+ord] = ROUND16(sum,SIG_SHIFT);
221       _y[i] = sum;
222    }
223    for(i=0;i<ord;i++)
224       mem[i] = _y[N-i-1];
225    RESTORE_STACK;
226 #endif
227 }
228
229 int _celt_autocorr(
230                    const opus_val16 *x,   /*  in: [0...n-1] samples x   */
231                    opus_val32       *ac,  /* out: [0...lag-1] ac values */
232                    const opus_val16       *window,
233                    int          overlap,
234                    int          lag,
235                    int          n,
236                    int          arch
237                   )
238 {
239    opus_val32 d;
240    int i, k;
241    int fastN=n-lag;
242    int shift;
243    const opus_val16 *xptr;
244    VARDECL(opus_val16, xx);
245    SAVE_STACK;
246    ALLOC(xx, n, opus_val16);
247    celt_assert(n>0);
248    celt_assert(overlap>=0);
249    if (overlap == 0)
250    {
251       xptr = x;
252    } else {
253       for (i=0;i<n;i++)
254          xx[i] = x[i];
255       for (i=0;i<overlap;i++)
256       {
257          xx[i] = MULT16_16_Q15(x[i],window[i]);
258          xx[n-i-1] = MULT16_16_Q15(x[n-i-1],window[i]);
259       }
260       xptr = xx;
261    }
262    shift=0;
263 #ifdef FIXED_POINT
264    {
265       opus_val32 ac0;
266       ac0 = 1+(n<<7);
267       if (n&1) ac0 += SHR32(MULT16_16(xptr[0],xptr[0]),9);
268       for(i=(n&1);i<n;i+=2)
269       {
270          ac0 += SHR32(MULT16_16(xptr[i],xptr[i]),9);
271          ac0 += SHR32(MULT16_16(xptr[i+1],xptr[i+1]),9);
272       }
273
274       shift = celt_ilog2(ac0)-30+10;
275       shift = (shift)/2;
276       if (shift>0)
277       {
278          for(i=0;i<n;i++)
279             xx[i] = PSHR32(xptr[i], shift);
280          xptr = xx;
281       } else
282          shift = 0;
283    }
284 #endif
285    celt_pitch_xcorr(xptr, xptr, ac, fastN, lag+1, arch);
286    for (k=0;k<=lag;k++)
287    {
288       for (i = k+fastN, d = 0; i < n; i++)
289          d = MAC16_16(d, xptr[i], xptr[i-k]);
290       ac[k] += d;
291    }
292 #ifdef FIXED_POINT
293    shift = 2*shift;
294    if (shift<=0)
295       ac[0] += SHL32((opus_int32)1, -shift);
296    if (ac[0] < 268435456)
297    {
298       int shift2 = 29 - EC_ILOG(ac[0]);
299       for (i=0;i<=lag;i++)
300          ac[i] = SHL32(ac[i], shift2);
301       shift -= shift2;
302    } else if (ac[0] >= 536870912)
303    {
304       int shift2=1;
305       if (ac[0] >= 1073741824)
306          shift2++;
307       for (i=0;i<=lag;i++)
308          ac[i] = SHR32(ac[i], shift2);
309       shift += shift2;
310    }
311 #endif
312
313    RESTORE_STACK;
314    return shift;
315 }