Optimizes _celt_autocorr() by using pitch_xcorr()
[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 void celt_fir(const opus_val16 *x,
92          const opus_val16 *num,
93          opus_val16 *y,
94          int N,
95          int ord,
96          opus_val16 *mem)
97 {
98    int i,j;
99
100    for (i=0;i<N;i++)
101    {
102       opus_val32 sum = SHL32(EXTEND32(x[i]), SIG_SHIFT);
103       for (j=0;j<ord;j++)
104       {
105          sum = MAC16_16(sum,num[j],mem[j]);
106       }
107       for (j=ord-1;j>=1;j--)
108       {
109          mem[j]=mem[j-1];
110       }
111       mem[0] = x[i];
112       y[i] = ROUND16(sum, SIG_SHIFT);
113    }
114 }
115
116 void celt_iir(const opus_val32 *x,
117          const opus_val16 *den,
118          opus_val32 *y,
119          int N,
120          int ord,
121          opus_val16 *mem)
122 {
123    int i,j;
124    for (i=0;i<N;i++)
125    {
126       opus_val32 sum = x[i];
127       for (j=0;j<ord;j++)
128       {
129          sum -= MULT16_16(den[j],mem[j]);
130       }
131       for (j=ord-1;j>=1;j--)
132       {
133          mem[j]=mem[j-1];
134       }
135       mem[0] = ROUND16(sum,SIG_SHIFT);
136       y[i] = sum;
137    }
138 }
139
140 void _celt_autocorr(
141                    const opus_val16 *x,   /*  in: [0...n-1] samples x   */
142                    opus_val32       *ac,  /* out: [0...lag-1] ac values */
143                    const opus_val16       *window,
144                    int          overlap,
145                    int          lag,
146                    int          n
147                   )
148 {
149    opus_val32 d;
150    int i;
151    int fastN=n-lag;
152    VARDECL(opus_val16, xx);
153    SAVE_STACK;
154    ALLOC(xx, n, opus_val16);
155    celt_assert(n>0);
156    celt_assert(overlap>=0);
157    for (i=0;i<n;i++)
158       xx[i] = x[i];
159    for (i=0;i<overlap;i++)
160    {
161       xx[i] = MULT16_16_Q15(x[i],window[i]);
162       xx[n-i-1] = MULT16_16_Q15(x[n-i-1],window[i]);
163    }
164 #ifdef FIXED_POINT
165    {
166       opus_val32 ac0;
167       int shift;
168       ac0 = 1+n;
169       if (n&1) ac0 += SHR32(MULT16_16(xx[0],xx[0]),9);
170       for(i=(n&1);i<n;i+=2)
171       {
172          ac0 += SHR32(MULT16_16(xx[i],xx[i]),9);
173          ac0 += SHR32(MULT16_16(xx[i+1],xx[i+1]),9);
174       }
175
176       shift = celt_ilog2(ac0)-30+10;
177       shift = (shift+1)/2;
178       for(i=0;i<n;i++)
179          xx[i] = VSHR32(xx[i], shift);
180    }
181 #endif
182    pitch_xcorr(xx, xx, ac, fastN, lag+1);
183    while (lag>=0)
184    {
185       for (i = lag+fastN, d = 0; i < n; i++)
186          d = MAC16_16(d, xx[i], xx[i-lag]);
187       ac[lag] += d;
188       /*printf ("%f ", ac[lag]);*/
189       lag--;
190    }
191    /*printf ("\n");*/
192    ac[0] += 10;
193
194    RESTORE_STACK;
195 }