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