Line data Source code
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 0 : 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 0 : opus_val32 error = ac[0];
46 : #ifdef FIXED_POINT
47 : opus_val32 lpc[LPC_ORDER];
48 : #else
49 0 : float *lpc = _lpc;
50 : #endif
51 :
52 0 : OPUS_CLEAR(lpc, p);
53 0 : if (ac[0] != 0)
54 : {
55 0 : for (i = 0; i < p; i++) {
56 : /* Sum up this iteration's reflection coefficient */
57 0 : opus_val32 rr = 0;
58 0 : for (j = 0; j < i; j++)
59 0 : rr += MULT32_32_Q31(lpc[j],ac[i - j]);
60 0 : rr += SHR32(ac[i + 1],3);
61 0 : r = -frac_div32(SHL32(rr,3), error);
62 : /* Update LPC coefficients and total error */
63 0 : lpc[i] = SHR32(r,3);
64 0 : for (j = 0; j < (i+1)>>1; j++)
65 : {
66 : opus_val32 tmp1, tmp2;
67 0 : tmp1 = lpc[j];
68 0 : tmp2 = lpc[i-1-j];
69 0 : lpc[j] = tmp1 + MULT32_32_Q31(r,tmp2);
70 0 : lpc[i-1-j] = tmp2 + MULT32_32_Q31(r,tmp1);
71 : }
72 :
73 0 : 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 0 : if (error<.001f*ac[0])
80 0 : 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 0 : }
89 :
90 :
91 0 : 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 :
103 0 : ALLOC(rnum, ord, opus_val16);
104 0 : for(i=0;i<ord;i++)
105 0 : rnum[i] = num[ord-i-1];
106 0 : for (i=0;i<N-3;i+=4)
107 : {
108 : opus_val32 sum[4];
109 0 : sum[0] = SHL32(EXTEND32(x[i ]), SIG_SHIFT);
110 0 : sum[1] = SHL32(EXTEND32(x[i+1]), SIG_SHIFT),
111 0 : sum[2] = SHL32(EXTEND32(x[i+2]), SIG_SHIFT);
112 0 : sum[3] = SHL32(EXTEND32(x[i+3]), SIG_SHIFT);
113 0 : xcorr_kernel(rnum, x+i-ord, sum, ord, arch);
114 0 : y[i ] = ROUND16(sum[0], SIG_SHIFT);
115 0 : y[i+1] = ROUND16(sum[1], SIG_SHIFT);
116 0 : y[i+2] = ROUND16(sum[2], SIG_SHIFT);
117 0 : y[i+3] = ROUND16(sum[3], SIG_SHIFT);
118 : }
119 0 : for (;i<N;i++)
120 : {
121 0 : opus_val32 sum = SHL32(EXTEND32(x[i]), SIG_SHIFT);
122 0 : for (j=0;j<ord;j++)
123 0 : sum = MAC16_16(sum,rnum[j],x[i+j-ord]);
124 0 : y[i] = ROUND16(sum, SIG_SHIFT);
125 : }
126 : RESTORE_STACK;
127 0 : }
128 :
129 0 : 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 0 : celt_assert((ord&3)==0);
161 0 : ALLOC(rden, ord, opus_val16);
162 0 : ALLOC(y, N+ord, opus_val16);
163 0 : for(i=0;i<ord;i++)
164 0 : rden[i] = den[ord-i-1];
165 0 : for(i=0;i<ord;i++)
166 0 : y[i] = -mem[ord-i-1];
167 0 : for(;i<N+ord;i++)
168 0 : y[i]=0;
169 0 : 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 0 : sum[0]=_x[i];
174 0 : sum[1]=_x[i+1];
175 0 : sum[2]=_x[i+2];
176 0 : sum[3]=_x[i+3];
177 0 : 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 0 : y[i+ord ] = -SROUND16(sum[0],SIG_SHIFT);
181 0 : _y[i ] = sum[0];
182 0 : sum[1] = MAC16_16(sum[1], y[i+ord ], den[0]);
183 0 : y[i+ord+1] = -SROUND16(sum[1],SIG_SHIFT);
184 0 : _y[i+1] = sum[1];
185 0 : sum[2] = MAC16_16(sum[2], y[i+ord+1], den[0]);
186 0 : sum[2] = MAC16_16(sum[2], y[i+ord ], den[1]);
187 0 : y[i+ord+2] = -SROUND16(sum[2],SIG_SHIFT);
188 0 : _y[i+2] = sum[2];
189 :
190 0 : sum[3] = MAC16_16(sum[3], y[i+ord+2], den[0]);
191 0 : sum[3] = MAC16_16(sum[3], y[i+ord+1], den[1]);
192 0 : sum[3] = MAC16_16(sum[3], y[i+ord ], den[2]);
193 0 : y[i+ord+3] = -SROUND16(sum[3],SIG_SHIFT);
194 0 : _y[i+3] = sum[3];
195 : }
196 0 : for (;i<N;i++)
197 : {
198 0 : opus_val32 sum = _x[i];
199 0 : for (j=0;j<ord;j++)
200 0 : sum -= MULT16_16(rden[j],y[i+j]);
201 0 : y[i+ord] = SROUND16(sum,SIG_SHIFT);
202 0 : _y[i] = sum;
203 : }
204 0 : for(i=0;i<ord;i++)
205 0 : mem[i] = _y[N-i-1];
206 : RESTORE_STACK;
207 : #endif
208 0 : }
209 :
210 0 : 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 0 : int fastN=n-lag;
223 : int shift;
224 : const opus_val16 *xptr;
225 : VARDECL(opus_val16, xx);
226 : SAVE_STACK;
227 0 : ALLOC(xx, n, opus_val16);
228 0 : celt_assert(n>0);
229 0 : celt_assert(overlap>=0);
230 0 : if (overlap == 0)
231 : {
232 0 : xptr = x;
233 : } else {
234 0 : for (i=0;i<n;i++)
235 0 : xx[i] = x[i];
236 0 : for (i=0;i<overlap;i++)
237 : {
238 0 : xx[i] = MULT16_16_Q15(x[i],window[i]);
239 0 : xx[n-i-1] = MULT16_16_Q15(x[n-i-1],window[i]);
240 : }
241 0 : xptr = xx;
242 : }
243 0 : 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 0 : celt_pitch_xcorr(xptr, xptr, ac, fastN, lag+1, arch);
267 0 : for (k=0;k<=lag;k++)
268 : {
269 0 : for (i = k+fastN, d = 0; i < n; i++)
270 0 : d = MAC16_16(d, xptr[i], xptr[i-k]);
271 0 : 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 0 : return shift;
296 : }
|
