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
_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 OPUS_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 = -frac_div32(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_c(const opus_val16 * x,const opus_val16 * num,opus_val16 * y,int N,int ord,int arch)91 void celt_fir_c(
92 const opus_val16 *x,
93 const opus_val16 *num,
94 opus_val16 *y,
95 int N,
96 int ord,
97 int arch)
98 {
99 int i,j;
100 VARDECL(opus_val16, rnum);
101 SAVE_STACK;
102 celt_assert(x != y);
103 ALLOC(rnum, ord, opus_val16);
104 for(i=0;i<ord;i++)
105 rnum[i] = num[ord-i-1];
106 for (i=0;i<N-3;i+=4)
107 {
108 opus_val32 sum[4];
109 sum[0] = SHL32(EXTEND32(x[i ]), SIG_SHIFT);
110 sum[1] = SHL32(EXTEND32(x[i+1]), SIG_SHIFT);
111 sum[2] = SHL32(EXTEND32(x[i+2]), SIG_SHIFT);
112 sum[3] = SHL32(EXTEND32(x[i+3]), SIG_SHIFT);
113 xcorr_kernel(rnum, x+i-ord, sum, ord, arch);
114 y[i ] = ROUND16(sum[0], SIG_SHIFT);
115 y[i+1] = ROUND16(sum[1], SIG_SHIFT);
116 y[i+2] = ROUND16(sum[2], SIG_SHIFT);
117 y[i+3] = ROUND16(sum[3], SIG_SHIFT);
118 }
119 for (;i<N;i++)
120 {
121 opus_val32 sum = SHL32(EXTEND32(x[i]), SIG_SHIFT);
122 for (j=0;j<ord;j++)
123 sum = MAC16_16(sum,rnum[j],x[i+j-ord]);
124 y[i] = ROUND16(sum, SIG_SHIFT);
125 }
126 RESTORE_STACK;
127 }
128
celt_iir(const opus_val32 * _x,const opus_val16 * den,opus_val32 * _y,int N,int ord,opus_val16 * mem,int arch)129 void celt_iir(const opus_val32 *_x,
130 const opus_val16 *den,
131 opus_val32 *_y,
132 int N,
133 int ord,
134 opus_val16 *mem,
135 int arch)
136 {
137 #ifdef SMALL_FOOTPRINT
138 int i,j;
139 (void)arch;
140 for (i=0;i<N;i++)
141 {
142 opus_val32 sum = _x[i];
143 for (j=0;j<ord;j++)
144 {
145 sum -= MULT16_16(den[j],mem[j]);
146 }
147 for (j=ord-1;j>=1;j--)
148 {
149 mem[j]=mem[j-1];
150 }
151 mem[0] = SROUND16(sum, SIG_SHIFT);
152 _y[i] = sum;
153 }
154 #else
155 int i,j;
156 VARDECL(opus_val16, rden);
157 VARDECL(opus_val16, y);
158 SAVE_STACK;
159
160 celt_assert((ord&3)==0);
161 ALLOC(rden, ord, opus_val16);
162 ALLOC(y, N+ord, opus_val16);
163 for(i=0;i<ord;i++)
164 rden[i] = den[ord-i-1];
165 for(i=0;i<ord;i++)
166 y[i] = -mem[ord-i-1];
167 for(;i<N+ord;i++)
168 y[i]=0;
169 for (i=0;i<N-3;i+=4)
170 {
171 /* Unroll by 4 as if it were an FIR filter */
172 opus_val32 sum[4];
173 sum[0]=_x[i];
174 sum[1]=_x[i+1];
175 sum[2]=_x[i+2];
176 sum[3]=_x[i+3];
177 xcorr_kernel(rden, y+i, sum, ord, arch);
178
179 /* Patch up the result to compensate for the fact that this is an IIR */
180 y[i+ord ] = -SROUND16(sum[0],SIG_SHIFT);
181 _y[i ] = sum[0];
182 sum[1] = MAC16_16(sum[1], y[i+ord ], den[0]);
183 y[i+ord+1] = -SROUND16(sum[1],SIG_SHIFT);
184 _y[i+1] = sum[1];
185 sum[2] = MAC16_16(sum[2], y[i+ord+1], den[0]);
186 sum[2] = MAC16_16(sum[2], y[i+ord ], den[1]);
187 y[i+ord+2] = -SROUND16(sum[2],SIG_SHIFT);
188 _y[i+2] = sum[2];
189
190 sum[3] = MAC16_16(sum[3], y[i+ord+2], den[0]);
191 sum[3] = MAC16_16(sum[3], y[i+ord+1], den[1]);
192 sum[3] = MAC16_16(sum[3], y[i+ord ], den[2]);
193 y[i+ord+3] = -SROUND16(sum[3],SIG_SHIFT);
194 _y[i+3] = sum[3];
195 }
196 for (;i<N;i++)
197 {
198 opus_val32 sum = _x[i];
199 for (j=0;j<ord;j++)
200 sum -= MULT16_16(rden[j],y[i+j]);
201 y[i+ord] = SROUND16(sum,SIG_SHIFT);
202 _y[i] = sum;
203 }
204 for(i=0;i<ord;i++)
205 mem[i] = _y[N-i-1];
206 RESTORE_STACK;
207 #endif
208 }
209
_celt_autocorr(const opus_val16 * x,opus_val32 * ac,const opus_val16 * window,int overlap,int lag,int n,int arch)210 int _celt_autocorr(
211 const opus_val16 *x, /* in: [0...n-1] samples x */
212 opus_val32 *ac, /* out: [0...lag-1] ac values */
213 const opus_val16 *window,
214 int overlap,
215 int lag,
216 int n,
217 int arch
218 )
219 {
220 opus_val32 d;
221 int i, k;
222 int fastN=n-lag;
223 int shift;
224 const opus_val16 *xptr;
225 VARDECL(opus_val16, xx);
226 SAVE_STACK;
227 ALLOC(xx, n, opus_val16);
228 celt_assert(n>0);
229 celt_assert(overlap>=0);
230 if (overlap == 0)
231 {
232 xptr = x;
233 } else {
234 for (i=0;i<n;i++)
235 xx[i] = x[i];
236 for (i=0;i<overlap;i++)
237 {
238 xx[i] = MULT16_16_Q15(x[i],window[i]);
239 xx[n-i-1] = MULT16_16_Q15(x[n-i-1],window[i]);
240 }
241 xptr = xx;
242 }
243 shift=0;
244 #ifdef FIXED_POINT
245 {
246 opus_val32 ac0;
247 ac0 = 1+(n<<7);
248 if (n&1) ac0 += SHR32(MULT16_16(xptr[0],xptr[0]),9);
249 for(i=(n&1);i<n;i+=2)
250 {
251 ac0 += SHR32(MULT16_16(xptr[i],xptr[i]),9);
252 ac0 += SHR32(MULT16_16(xptr[i+1],xptr[i+1]),9);
253 }
254
255 shift = celt_ilog2(ac0)-30+10;
256 shift = (shift)/2;
257 if (shift>0)
258 {
259 for(i=0;i<n;i++)
260 xx[i] = PSHR32(xptr[i], shift);
261 xptr = xx;
262 } else
263 shift = 0;
264 }
265 #endif
266 celt_pitch_xcorr(xptr, xptr, ac, fastN, lag+1, arch);
267 for (k=0;k<=lag;k++)
268 {
269 for (i = k+fastN, d = 0; i < n; i++)
270 d = MAC16_16(d, xptr[i], xptr[i-k]);
271 ac[k] += d;
272 }
273 #ifdef FIXED_POINT
274 shift = 2*shift;
275 if (shift<=0)
276 ac[0] += SHL32((opus_int32)1, -shift);
277 if (ac[0] < 268435456)
278 {
279 int shift2 = 29 - EC_ILOG(ac[0]);
280 for (i=0;i<=lag;i++)
281 ac[i] = SHL32(ac[i], shift2);
282 shift -= shift2;
283 } else if (ac[0] >= 536870912)
284 {
285 int shift2=1;
286 if (ac[0] >= 1073741824)
287 shift2++;
288 for (i=0;i<=lag;i++)
289 ac[i] = SHR32(ac[i], shift2);
290 shift += shift2;
291 }
292 #endif
293
294 RESTORE_STACK;
295 return shift;
296 }
297