1 /* Copyright (c) 2007-2008 CSIRO
2 Copyright (c) 2007-2009 Xiph.Org Foundation
3 Copyright (c) 2008-2009 Gregory Maxwell
4 Written by Jean-Marc Valin and Gregory Maxwell */
5 /*
6 Redistribution and use in source and binary forms, with or without
7 modification, are permitted provided that the following conditions
8 are met:
9
10 - Redistributions of source code must retain the above copyright
11 notice, this list of conditions and the following disclaimer.
12
13 - Redistributions in binary form must reproduce the above copyright
14 notice, this list of conditions and the following disclaimer in the
15 documentation and/or other materials provided with the distribution.
16
17 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
18 ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
19 LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
20 A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER
21 OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
22 EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
23 PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
24 PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
25 LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
26 NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
27 SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
28 */
29
30 #ifdef HAVE_CONFIG_H
31 #include "config.h"
32 #endif
33
34 #include <math.h>
35 #include "bands.h"
36 #include "modes.h"
37 #include "vq.h"
38 #include "cwrs.h"
39 #include "stack_alloc.h"
40 #include "os_support.h"
41 #include "mathops.h"
42 #include "rate.h"
43
celt_lcg_rand(opus_uint32 seed)44 opus_uint32 celt_lcg_rand(opus_uint32 seed)
45 {
46 return 1664525 * seed + 1013904223;
47 }
48
49 /* This is a cos() approximation designed to be bit-exact on any platform. Bit exactness
50 with this approximation is important because it has an impact on the bit allocation */
bitexact_cos(opus_int16 x)51 static opus_int16 bitexact_cos(opus_int16 x)
52 {
53 opus_int32 tmp;
54 opus_int16 x2;
55 tmp = (4096+((opus_int32)(x)*(x)))>>13;
56 celt_assert(tmp<=32767);
57 x2 = tmp;
58 x2 = (32767-x2) + FRAC_MUL16(x2, (-7651 + FRAC_MUL16(x2, (8277 + FRAC_MUL16(-626, x2)))));
59 celt_assert(x2<=32766);
60 return 1+x2;
61 }
62
bitexact_log2tan(int isin,int icos)63 static int bitexact_log2tan(int isin,int icos)
64 {
65 int lc;
66 int ls;
67 lc=EC_ILOG(icos);
68 ls=EC_ILOG(isin);
69 icos<<=15-lc;
70 isin<<=15-ls;
71 return (ls-lc)*(1<<11)
72 +FRAC_MUL16(isin, FRAC_MUL16(isin, -2597) + 7932)
73 -FRAC_MUL16(icos, FRAC_MUL16(icos, -2597) + 7932);
74 }
75
76 #ifdef FIXED_POINT
77 /* Compute the amplitude (sqrt energy) in each of the bands */
compute_band_energies(const CELTMode * m,const celt_sig * X,celt_ener * bandE,int end,int C,int M)78 void compute_band_energies(const CELTMode *m, const celt_sig *X, celt_ener *bandE, int end, int C, int M)
79 {
80 int i, c, N;
81 const opus_int16 *eBands = m->eBands;
82 N = M*m->shortMdctSize;
83 c=0; do {
84 for (i=0;i<end;i++)
85 {
86 int j;
87 opus_val32 maxval=0;
88 opus_val32 sum = 0;
89
90 j=M*eBands[i]; do {
91 maxval = MAX32(maxval, X[j+c*N]);
92 maxval = MAX32(maxval, -X[j+c*N]);
93 } while (++j<M*eBands[i+1]);
94
95 if (maxval > 0)
96 {
97 int shift = celt_ilog2(maxval)-10;
98 j=M*eBands[i]; do {
99 sum = MAC16_16(sum, EXTRACT16(VSHR32(X[j+c*N],shift)),
100 EXTRACT16(VSHR32(X[j+c*N],shift)));
101 } while (++j<M*eBands[i+1]);
102 /* We're adding one here to ensure the normalized band isn't larger than unity norm */
103 bandE[i+c*m->nbEBands] = EPSILON+VSHR32(EXTEND32(celt_sqrt(sum)),-shift);
104 } else {
105 bandE[i+c*m->nbEBands] = EPSILON;
106 }
107 /*printf ("%f ", bandE[i+c*m->nbEBands]);*/
108 }
109 } while (++c<C);
110 /*printf ("\n");*/
111 }
112
113 /* Normalise each band such that the energy is one. */
normalise_bands(const CELTMode * m,const celt_sig * OPUS_RESTRICT freq,celt_norm * OPUS_RESTRICT X,const celt_ener * bandE,int end,int C,int M)114 void normalise_bands(const CELTMode *m, const celt_sig * OPUS_RESTRICT freq, celt_norm * OPUS_RESTRICT X, const celt_ener *bandE, int end, int C, int M)
115 {
116 int i, c, N;
117 const opus_int16 *eBands = m->eBands;
118 N = M*m->shortMdctSize;
119 c=0; do {
120 i=0; do {
121 opus_val16 g;
122 int j,shift;
123 opus_val16 E;
124 shift = celt_zlog2(bandE[i+c*m->nbEBands])-13;
125 E = VSHR32(bandE[i+c*m->nbEBands], shift);
126 g = EXTRACT16(celt_rcp(SHL32(E,3)));
127 j=M*eBands[i]; do {
128 X[j+c*N] = MULT16_16_Q15(VSHR32(freq[j+c*N],shift-1),g);
129 } while (++j<M*eBands[i+1]);
130 } while (++i<end);
131 } while (++c<C);
132 }
133
134 #else /* FIXED_POINT */
135 /* Compute the amplitude (sqrt energy) in each of the bands */
compute_band_energies(const CELTMode * m,const celt_sig * X,celt_ener * bandE,int end,int C,int M)136 void compute_band_energies(const CELTMode *m, const celt_sig *X, celt_ener *bandE, int end, int C, int M)
137 {
138 int i, c, N;
139 const opus_int16 *eBands = m->eBands;
140 N = M*m->shortMdctSize;
141 c=0; do {
142 for (i=0;i<end;i++)
143 {
144 int j;
145 opus_val32 sum = 1e-27f;
146 for (j=M*eBands[i];j<M*eBands[i+1];j++)
147 sum += X[j+c*N]*X[j+c*N];
148 bandE[i+c*m->nbEBands] = celt_sqrt(sum);
149 /*printf ("%f ", bandE[i+c*m->nbEBands]);*/
150 }
151 } while (++c<C);
152 /*printf ("\n");*/
153 }
154
155 /* Normalise each band such that the energy is one. */
normalise_bands(const CELTMode * m,const celt_sig * OPUS_RESTRICT freq,celt_norm * OPUS_RESTRICT X,const celt_ener * bandE,int end,int C,int M)156 void normalise_bands(const CELTMode *m, const celt_sig * OPUS_RESTRICT freq, celt_norm * OPUS_RESTRICT X, const celt_ener *bandE, int end, int C, int M)
157 {
158 int i, c, N;
159 const opus_int16 *eBands = m->eBands;
160 N = M*m->shortMdctSize;
161 c=0; do {
162 for (i=0;i<end;i++)
163 {
164 int j;
165 opus_val16 g = 1.f/(1e-27f+bandE[i+c*m->nbEBands]);
166 for (j=M*eBands[i];j<M*eBands[i+1];j++)
167 X[j+c*N] = freq[j+c*N]*g;
168 }
169 } while (++c<C);
170 }
171
172 #endif /* FIXED_POINT */
173
174 /* De-normalise the energy to produce the synthesis from the unit-energy bands */
denormalise_bands(const CELTMode * m,const celt_norm * OPUS_RESTRICT X,celt_sig * OPUS_RESTRICT freq,const celt_ener * bandE,int end,int C,int M)175 void denormalise_bands(const CELTMode *m, const celt_norm * OPUS_RESTRICT X, celt_sig * OPUS_RESTRICT freq, const celt_ener *bandE, int end, int C, int M)
176 {
177 int i, c, N;
178 const opus_int16 *eBands = m->eBands;
179 N = M*m->shortMdctSize;
180 celt_assert2(C<=2, "denormalise_bands() not implemented for >2 channels");
181 c=0; do {
182 celt_sig * OPUS_RESTRICT f;
183 const celt_norm * OPUS_RESTRICT x;
184 f = freq+c*N;
185 x = X+c*N;
186 for (i=0;i<end;i++)
187 {
188 int j, band_end;
189 opus_val32 g = SHR32(bandE[i+c*m->nbEBands],1);
190 j=M*eBands[i];
191 band_end = M*eBands[i+1];
192 do {
193 *f++ = SHL32(MULT16_32_Q15(*x, g),2);
194 x++;
195 } while (++j<band_end);
196 }
197 for (i=M*eBands[end];i<N;i++)
198 *f++ = 0;
199 } while (++c<C);
200 }
201
202 /* This prevents energy collapse for transients with multiple short MDCTs */
anti_collapse(const CELTMode * m,celt_norm * X_,unsigned char * collapse_masks,int LM,int C,int size,int start,int end,opus_val16 * logE,opus_val16 * prev1logE,opus_val16 * prev2logE,int * pulses,opus_uint32 seed)203 void anti_collapse(const CELTMode *m, celt_norm *X_, unsigned char *collapse_masks, int LM, int C, int size,
204 int start, int end, opus_val16 *logE, opus_val16 *prev1logE,
205 opus_val16 *prev2logE, int *pulses, opus_uint32 seed)
206 {
207 int c, i, j, k;
208 for (i=start;i<end;i++)
209 {
210 int N0;
211 opus_val16 thresh, sqrt_1;
212 int depth;
213 #ifdef FIXED_POINT
214 int shift;
215 opus_val32 thresh32;
216 #endif
217
218 N0 = m->eBands[i+1]-m->eBands[i];
219 /* depth in 1/8 bits */
220 depth = (1+pulses[i])/((m->eBands[i+1]-m->eBands[i])<<LM);
221
222 #ifdef FIXED_POINT
223 thresh32 = SHR32(celt_exp2(-SHL16(depth, 10-BITRES)),1);
224 thresh = MULT16_32_Q15(QCONST16(0.5f, 15), MIN32(32767,thresh32));
225 {
226 opus_val32 t;
227 t = N0<<LM;
228 shift = celt_ilog2(t)>>1;
229 t = SHL32(t, (7-shift)<<1);
230 sqrt_1 = celt_rsqrt_norm(t);
231 }
232 #else
233 thresh = .5f*celt_exp2(-.125f*depth);
234 sqrt_1 = celt_rsqrt(N0<<LM);
235 #endif
236
237 c=0; do
238 {
239 celt_norm *X;
240 opus_val16 prev1;
241 opus_val16 prev2;
242 opus_val32 Ediff;
243 opus_val16 r;
244 int renormalize=0;
245 prev1 = prev1logE[c*m->nbEBands+i];
246 prev2 = prev2logE[c*m->nbEBands+i];
247 if (C==1)
248 {
249 prev1 = MAX16(prev1,prev1logE[m->nbEBands+i]);
250 prev2 = MAX16(prev2,prev2logE[m->nbEBands+i]);
251 }
252 Ediff = EXTEND32(logE[c*m->nbEBands+i])-EXTEND32(MIN16(prev1,prev2));
253 Ediff = MAX32(0, Ediff);
254
255 #ifdef FIXED_POINT
256 if (Ediff < 16384)
257 {
258 opus_val32 r32 = SHR32(celt_exp2(-EXTRACT16(Ediff)),1);
259 r = 2*MIN16(16383,r32);
260 } else {
261 r = 0;
262 }
263 if (LM==3)
264 r = MULT16_16_Q14(23170, MIN32(23169, r));
265 r = SHR16(MIN16(thresh, r),1);
266 r = SHR32(MULT16_16_Q15(sqrt_1, r),shift);
267 #else
268 /* r needs to be multiplied by 2 or 2*sqrt(2) depending on LM because
269 short blocks don't have the same energy as long */
270 r = 2.f*celt_exp2(-Ediff);
271 if (LM==3)
272 r *= 1.41421356f;
273 r = MIN16(thresh, r);
274 r = r*sqrt_1;
275 #endif
276 X = X_+c*size+(m->eBands[i]<<LM);
277 for (k=0;k<1<<LM;k++)
278 {
279 /* Detect collapse */
280 if (!(collapse_masks[i*C+c]&1<<k))
281 {
282 /* Fill with noise */
283 for (j=0;j<N0;j++)
284 {
285 seed = celt_lcg_rand(seed);
286 X[(j<<LM)+k] = (seed&0x8000 ? r : -r);
287 }
288 renormalize = 1;
289 }
290 }
291 /* We just added some energy, so we need to renormalise */
292 if (renormalize)
293 renormalise_vector(X, N0<<LM, Q15ONE);
294 } while (++c<C);
295 }
296 }
297
298 #ifdef OPUS_ENABLE_ENCODER
intensity_stereo(const CELTMode * m,celt_norm * X,celt_norm * Y,const celt_ener * bandE,int bandID,int N)299 static void intensity_stereo(const CELTMode *m, celt_norm *X, celt_norm *Y, const celt_ener *bandE, int bandID, int N)
300 {
301 int i = bandID;
302 int j;
303 opus_val16 a1, a2;
304 opus_val16 left, right;
305 opus_val16 norm;
306 #ifdef FIXED_POINT
307 int shift = celt_zlog2(MAX32(bandE[i], bandE[i+m->nbEBands]))-13;
308 #endif
309 left = VSHR32(bandE[i],shift);
310 right = VSHR32(bandE[i+m->nbEBands],shift);
311 norm = EPSILON + celt_sqrt(EPSILON+MULT16_16(left,left)+MULT16_16(right,right));
312 a1 = DIV32_16(SHL32(EXTEND32(left),14),norm);
313 a2 = DIV32_16(SHL32(EXTEND32(right),14),norm);
314 for (j=0;j<N;j++)
315 {
316 celt_norm r, l;
317 l = X[j];
318 r = Y[j];
319 X[j] = MULT16_16_Q14(a1,l) + MULT16_16_Q14(a2,r);
320 /* Side is not encoded, no need to calculate */
321 }
322 }
323
stereo_split(celt_norm * X,celt_norm * Y,int N)324 static void stereo_split(celt_norm *X, celt_norm *Y, int N)
325 {
326 int j;
327 for (j=0;j<N;j++)
328 {
329 celt_norm r, l;
330 l = MULT16_16_Q15(QCONST16(.70710678f,15), X[j]);
331 r = MULT16_16_Q15(QCONST16(.70710678f,15), Y[j]);
332 X[j] = l+r;
333 Y[j] = r-l;
334 }
335 }
336 #endif /* OPUS_ENABLE_ENCODER */
337
stereo_merge(celt_norm * X,celt_norm * Y,opus_val16 mid,int N)338 static void stereo_merge(celt_norm *X, celt_norm *Y, opus_val16 mid, int N)
339 {
340 int j;
341 opus_val32 xp=0, side=0;
342 opus_val32 El, Er;
343 opus_val16 mid2;
344 #ifdef FIXED_POINT
345 int kl, kr;
346 #endif
347 opus_val32 t, lgain, rgain;
348
349 /* Compute the norm of X+Y and X-Y as |X|^2 + |Y|^2 +/- sum(xy) */
350 for (j=0;j<N;j++)
351 {
352 xp = MAC16_16(xp, X[j], Y[j]);
353 side = MAC16_16(side, Y[j], Y[j]);
354 }
355 /* Compensating for the mid normalization */
356 xp = MULT16_32_Q15(mid, xp);
357 /* mid and side are in Q15, not Q14 like X and Y */
358 mid2 = SHR32(mid, 1);
359 El = MULT16_16(mid2, mid2) + side - 2*xp;
360 Er = MULT16_16(mid2, mid2) + side + 2*xp;
361 if (Er < QCONST32(6e-4f, 28) || El < QCONST32(6e-4f, 28))
362 {
363 for (j=0;j<N;j++)
364 Y[j] = X[j];
365 return;
366 }
367
368 #ifdef FIXED_POINT
369 kl = celt_ilog2(El)>>1;
370 kr = celt_ilog2(Er)>>1;
371 #endif
372 t = VSHR32(El, (kl-7)<<1);
373 lgain = celt_rsqrt_norm(t);
374 t = VSHR32(Er, (kr-7)<<1);
375 rgain = celt_rsqrt_norm(t);
376
377 #ifdef FIXED_POINT
378 if (kl < 7)
379 kl = 7;
380 if (kr < 7)
381 kr = 7;
382 #endif
383
384 for (j=0;j<N;j++)
385 {
386 celt_norm r, l;
387 /* Apply mid scaling (side is already scaled) */
388 l = MULT16_16_Q15(mid, X[j]);
389 r = Y[j];
390 X[j] = EXTRACT16(PSHR32(MULT16_16(lgain, SUB16(l,r)), kl+1));
391 Y[j] = EXTRACT16(PSHR32(MULT16_16(rgain, ADD16(l,r)), kr+1));
392 }
393 }
394
395 /* Decide whether we should spread the pulses in the current frame */
spreading_decision(const CELTMode * m,celt_norm * X,int * average,int last_decision,int * hf_average,int * tapset_decision,int update_hf,int end,int C,int M)396 int spreading_decision(const CELTMode *m, celt_norm *X, int *average,
397 int last_decision, int *hf_average, int *tapset_decision, int update_hf,
398 int end, int C, int M)
399 {
400 int i, c, N0;
401 int sum = 0, nbBands=0;
402 const opus_int16 * OPUS_RESTRICT eBands = m->eBands;
403 int decision;
404 int hf_sum=0;
405
406 celt_assert(end>0);
407
408 N0 = M*m->shortMdctSize;
409
410 if (M*(eBands[end]-eBands[end-1]) <= 8)
411 return SPREAD_NONE;
412 c=0; do {
413 for (i=0;i<end;i++)
414 {
415 int j, N, tmp=0;
416 int tcount[3] = {0,0,0};
417 celt_norm * OPUS_RESTRICT x = X+M*eBands[i]+c*N0;
418 N = M*(eBands[i+1]-eBands[i]);
419 if (N<=8)
420 continue;
421 /* Compute rough CDF of |x[j]| */
422 for (j=0;j<N;j++)
423 {
424 opus_val32 x2N; /* Q13 */
425
426 x2N = MULT16_16(MULT16_16_Q15(x[j], x[j]), N);
427 if (x2N < QCONST16(0.25f,13))
428 tcount[0]++;
429 if (x2N < QCONST16(0.0625f,13))
430 tcount[1]++;
431 if (x2N < QCONST16(0.015625f,13))
432 tcount[2]++;
433 }
434
435 /* Only include four last bands (8 kHz and up) */
436 if (i>m->nbEBands-4)
437 hf_sum += 32*(tcount[1]+tcount[0])/N;
438 tmp = (2*tcount[2] >= N) + (2*tcount[1] >= N) + (2*tcount[0] >= N);
439 sum += tmp*256;
440 nbBands++;
441 }
442 } while (++c<C);
443
444 if (update_hf)
445 {
446 if (hf_sum)
447 hf_sum /= C*(4-m->nbEBands+end);
448 *hf_average = (*hf_average+hf_sum)>>1;
449 hf_sum = *hf_average;
450 if (*tapset_decision==2)
451 hf_sum += 4;
452 else if (*tapset_decision==0)
453 hf_sum -= 4;
454 if (hf_sum > 22)
455 *tapset_decision=2;
456 else if (hf_sum > 18)
457 *tapset_decision=1;
458 else
459 *tapset_decision=0;
460 }
461 /*printf("%d %d %d\n", hf_sum, *hf_average, *tapset_decision);*/
462 celt_assert(nbBands>0); /*M*(eBands[end]-eBands[end-1]) <= 8 assures this*/
463 sum /= nbBands;
464 /* Recursive averaging */
465 sum = (sum+*average)>>1;
466 *average = sum;
467 /* Hysteresis */
468 sum = (3*sum + (((3-last_decision)<<7) + 64) + 2)>>2;
469 if (sum < 80)
470 {
471 decision = SPREAD_AGGRESSIVE;
472 } else if (sum < 256)
473 {
474 decision = SPREAD_NORMAL;
475 } else if (sum < 384)
476 {
477 decision = SPREAD_LIGHT;
478 } else {
479 decision = SPREAD_NONE;
480 }
481 #ifdef FUZZING
482 decision = rand()&0x3;
483 *tapset_decision=rand()%3;
484 #endif
485 return decision;
486 }
487
488 #ifdef MEASURE_NORM_MSE
489
490 float MSE[30] = {0};
491 int nbMSEBands = 0;
492 int MSECount[30] = {0};
493
dump_norm_mse(void)494 void dump_norm_mse(void)
495 {
496 int i;
497 for (i=0;i<nbMSEBands;i++)
498 {
499 printf ("%g ", MSE[i]/MSECount[i]);
500 }
501 printf ("\n");
502 }
503
measure_norm_mse(const CELTMode * m,float * X,float * X0,float * bandE,float * bandE0,int M,int N,int C)504 void measure_norm_mse(const CELTMode *m, float *X, float *X0, float *bandE, float *bandE0, int M, int N, int C)
505 {
506 static int init = 0;
507 int i;
508 if (!init)
509 {
510 atexit(dump_norm_mse);
511 init = 1;
512 }
513 for (i=0;i<m->nbEBands;i++)
514 {
515 int j;
516 int c;
517 float g;
518 if (bandE0[i]<10 || (C==2 && bandE0[i+m->nbEBands]<1))
519 continue;
520 c=0; do {
521 g = bandE[i+c*m->nbEBands]/(1e-15+bandE0[i+c*m->nbEBands]);
522 for (j=M*m->eBands[i];j<M*m->eBands[i+1];j++)
523 MSE[i] += (g*X[j+c*N]-X0[j+c*N])*(g*X[j+c*N]-X0[j+c*N]);
524 } while (++c<C);
525 MSECount[i]+=C;
526 }
527 nbMSEBands = m->nbEBands;
528 }
529
530 #endif
531
532 /* Indexing table for converting from natural Hadamard to ordery Hadamard
533 This is essentially a bit-reversed Gray, on top of which we've added
534 an inversion of the order because we want the DC at the end rather than
535 the beginning. The lines are for N=2, 4, 8, 16 */
536 static const int ordery_table[] = {
537 1, 0,
538 3, 0, 2, 1,
539 7, 0, 4, 3, 6, 1, 5, 2,
540 15, 0, 8, 7, 12, 3, 11, 4, 14, 1, 9, 6, 13, 2, 10, 5,
541 };
542
deinterleave_hadamard(celt_norm * X,int N0,int stride,int hadamard)543 static void deinterleave_hadamard(celt_norm *X, int N0, int stride, int hadamard)
544 {
545 int i,j;
546 VARDECL(celt_norm, tmp);
547 int N;
548 SAVE_STACK;
549 N = N0*stride;
550 ALLOC(tmp, N, celt_norm);
551 celt_assert(stride>0);
552 if (hadamard)
553 {
554 const int *ordery = ordery_table+stride-2;
555 for (i=0;i<stride;i++)
556 {
557 for (j=0;j<N0;j++)
558 tmp[ordery[i]*N0+j] = X[j*stride+i];
559 }
560 } else {
561 for (i=0;i<stride;i++)
562 for (j=0;j<N0;j++)
563 tmp[i*N0+j] = X[j*stride+i];
564 }
565 for (j=0;j<N;j++)
566 X[j] = tmp[j];
567 RESTORE_STACK;
568 }
569
interleave_hadamard(celt_norm * X,int N0,int stride,int hadamard)570 static void interleave_hadamard(celt_norm *X, int N0, int stride, int hadamard)
571 {
572 int i,j;
573 VARDECL(celt_norm, tmp);
574 int N;
575 SAVE_STACK;
576 N = N0*stride;
577 ALLOC(tmp, N, celt_norm);
578 if (hadamard)
579 {
580 const int *ordery = ordery_table+stride-2;
581 for (i=0;i<stride;i++)
582 for (j=0;j<N0;j++)
583 tmp[j*stride+i] = X[ordery[i]*N0+j];
584 } else {
585 for (i=0;i<stride;i++)
586 for (j=0;j<N0;j++)
587 tmp[j*stride+i] = X[i*N0+j];
588 }
589 for (j=0;j<N;j++)
590 X[j] = tmp[j];
591 RESTORE_STACK;
592 }
593
haar1(celt_norm * X,int N0,int stride)594 void haar1(celt_norm *X, int N0, int stride)
595 {
596 int i, j;
597 N0 >>= 1;
598 for (i=0;i<stride;i++)
599 for (j=0;j<N0;j++)
600 {
601 celt_norm tmp1, tmp2;
602 tmp1 = MULT16_16_Q15(QCONST16(.70710678f,15), X[stride*2*j+i]);
603 tmp2 = MULT16_16_Q15(QCONST16(.70710678f,15), X[stride*(2*j+1)+i]);
604 X[stride*2*j+i] = tmp1 + tmp2;
605 X[stride*(2*j+1)+i] = tmp1 - tmp2;
606 }
607 }
608
compute_qn(int N,int b,int offset,int pulse_cap,int stereo)609 static int compute_qn(int N, int b, int offset, int pulse_cap, int stereo)
610 {
611 static const opus_int16 exp2_table8[8] =
612 {16384, 17866, 19483, 21247, 23170, 25267, 27554, 30048};
613 int qn, qb;
614 int N2 = 2*N-1;
615 if (stereo && N==2)
616 N2--;
617 /* The upper limit ensures that in a stereo split with itheta==16384, we'll
618 always have enough bits left over to code at least one pulse in the
619 side; otherwise it would collapse, since it doesn't get folded. */
620 qb = IMIN(b-pulse_cap-(4<<BITRES), (b+N2*offset)/N2);
621
622 qb = IMIN(8<<BITRES, qb);
623
624 if (qb<(1<<BITRES>>1)) {
625 qn = 1;
626 } else {
627 qn = exp2_table8[qb&0x7]>>(14-(qb>>BITRES));
628 qn = (qn+1)>>1<<1;
629 }
630 celt_assert(qn <= 256);
631 return qn;
632 }
633
634 /* This function is responsible for encoding and decoding a band for both
635 the mono and stereo case. Even in the mono case, it can split the band
636 in two and transmit the energy difference with the two half-bands. It
637 can be called recursively so bands can end up being split in 8 parts. */
quant_band(int encode,const CELTMode * m,int i,celt_norm * X,celt_norm * Y,int N,int b,int spread,int B,int intensity,int tf_change,celt_norm * lowband,ec_ctx * ec,opus_int32 * remaining_bits,int LM,celt_norm * lowband_out,const celt_ener * bandE,int level,opus_uint32 * seed,opus_val16 gain,celt_norm * lowband_scratch,int fill)638 static unsigned quant_band(int encode, const CELTMode *m, int i, celt_norm *X, celt_norm *Y,
639 int N, int b, int spread, int B, int intensity, int tf_change, celt_norm *lowband, ec_ctx *ec,
640 opus_int32 *remaining_bits, int LM, celt_norm *lowband_out, const celt_ener *bandE, int level,
641 opus_uint32 *seed, opus_val16 gain, celt_norm *lowband_scratch, int fill)
642 {
643 const unsigned char *cache;
644 int q;
645 int curr_bits;
646 int stereo, split;
647 int imid=0, iside=0;
648 int N0=N;
649 int N_B=N;
650 int N_B0;
651 int B0=B;
652 int time_divide=0;
653 int recombine=0;
654 int inv = 0;
655 opus_val16 mid=0, side=0;
656 int longBlocks;
657 unsigned cm=0;
658 #ifdef RESYNTH
659 int resynth = 1;
660 #else
661 int resynth = !encode;
662 #endif
663
664 longBlocks = B0==1;
665
666 N_B /= B;
667 N_B0 = N_B;
668
669 split = stereo = Y != NULL;
670
671 /* Special case for one sample */
672 if (N==1)
673 {
674 int c;
675 celt_norm *x = X;
676 c=0; do {
677 int sign=0;
678 if (*remaining_bits>=1<<BITRES)
679 {
680 #ifdef OPUS_ENABLE_ENCODER
681 if (encode)
682 {
683 sign = x[0]<0;
684 ec_enc_bits(ec, sign, 1);
685 } else
686 #endif
687 {
688 sign = ec_dec_bits(ec, 1);
689 }
690 *remaining_bits -= 1<<BITRES;
691 b-=1<<BITRES;
692 }
693 if (resynth)
694 x[0] = sign ? -NORM_SCALING : NORM_SCALING;
695 x = Y;
696 } while (++c<1+stereo);
697 if (lowband_out)
698 lowband_out[0] = SHR16(X[0],4);
699 return 1;
700 }
701
702 if (!stereo && level == 0)
703 {
704 int k;
705 if (tf_change>0)
706 recombine = tf_change;
707 /* Band recombining to increase frequency resolution */
708
709 if (lowband && (recombine || ((N_B&1) == 0 && tf_change<0) || B0>1))
710 {
711 int j;
712 for (j=0;j<N;j++)
713 lowband_scratch[j] = lowband[j];
714 lowband = lowband_scratch;
715 }
716
717 for (k=0;k<recombine;k++)
718 {
719 static const unsigned char bit_interleave_table[16]={
720 0,1,1,1,2,3,3,3,2,3,3,3,2,3,3,3
721 };
722 #ifdef OPUS_ENABLE_ENCODER
723 if (encode)
724 haar1(X, N>>k, 1<<k);
725 #endif
726 if (lowband)
727 haar1(lowband, N>>k, 1<<k);
728 fill = bit_interleave_table[fill&0xF]|bit_interleave_table[fill>>4]<<2;
729 }
730 B>>=recombine;
731 N_B<<=recombine;
732
733 /* Increasing the time resolution */
734 while ((N_B&1) == 0 && tf_change<0)
735 {
736 #ifdef OPUS_ENABLE_ENCODER
737 if (encode)
738 haar1(X, N_B, B);
739 #endif
740 if (lowband)
741 haar1(lowband, N_B, B);
742 fill |= fill<<B;
743 B <<= 1;
744 N_B >>= 1;
745 time_divide++;
746 tf_change++;
747 }
748 B0=B;
749 N_B0 = N_B;
750
751 /* Reorganize the samples in time order instead of frequency order */
752 if (B0>1)
753 {
754 #ifdef OPUS_ENABLE_ENCODER
755 if (encode)
756 deinterleave_hadamard(X, N_B>>recombine, B0<<recombine, longBlocks);
757 #endif
758 if (lowband)
759 deinterleave_hadamard(lowband, N_B>>recombine, B0<<recombine, longBlocks);
760 }
761 }
762
763 /* If we need 1.5 more bit than we can produce, split the band in two. */
764 cache = m->cache.bits + m->cache.index[(LM+1)*m->nbEBands+i];
765 if (!stereo && LM != -1 && b > cache[cache[0]]+12 && N>2)
766 {
767 N >>= 1;
768 Y = X+N;
769 split = 1;
770 LM -= 1;
771 if (B==1)
772 fill = (fill&1)|(fill<<1);
773 B = (B+1)>>1;
774 }
775
776 if (split)
777 {
778 int qn;
779 int itheta=0;
780 int mbits, sbits, delta;
781 int qalloc;
782 int pulse_cap;
783 int offset;
784 int orig_fill;
785 opus_int32 tell;
786
787 /* Decide on the resolution to give to the split parameter theta */
788 pulse_cap = m->logN[i]+LM*(1<<BITRES);
789 offset = (pulse_cap>>1) - (stereo&&N==2 ? QTHETA_OFFSET_TWOPHASE : QTHETA_OFFSET);
790 qn = compute_qn(N, b, offset, pulse_cap, stereo);
791 if (stereo && i>=intensity)
792 qn = 1;
793 #ifdef OPUS_ENABLE_ENCODER
794 if (encode)
795 {
796 /* theta is the atan() of the ratio between the (normalized)
797 side and mid. With just that parameter, we can re-scale both
798 mid and side because we know that 1) they have unit norm and
799 2) they are orthogonal. */
800 itheta = stereo_itheta(X, Y, stereo, N);
801 }
802 #endif
803 tell = ec_tell_frac(ec);
804 if (qn!=1)
805 {
806 #ifdef OPUS_ENABLE_ENCODER
807 if (encode)
808 itheta = (itheta*qn+8192)>>14;
809 #endif
810
811 /* Entropy coding of the angle. We use a uniform pdf for the
812 time split, a step for stereo, and a triangular one for the rest. */
813 if (stereo && N>2)
814 {
815 int p0 = 3;
816 int x = itheta;
817 int x0 = qn/2;
818 int ft = p0*(x0+1) + x0;
819 /* Use a probability of p0 up to itheta=8192 and then use 1 after */
820 #ifdef OPUS_ENABLE_ENCODER
821 if (encode)
822 {
823 ec_encode(ec,x<=x0?p0*x:(x-1-x0)+(x0+1)*p0,x<=x0?p0*(x+1):(x-x0)+(x0+1)*p0,ft);
824 } else
825 #endif
826 {
827 int fs;
828 fs=ec_decode(ec,ft);
829 if (fs<(x0+1)*p0)
830 x=fs/p0;
831 else
832 x=x0+1+(fs-(x0+1)*p0);
833 ec_dec_update(ec,x<=x0?p0*x:(x-1-x0)+(x0+1)*p0,x<=x0?p0*(x+1):(x-x0)+(x0+1)*p0,ft);
834 itheta = x;
835 }
836 } else if (B0>1 || stereo) {
837 /* Uniform pdf */
838 #ifdef OPUS_ENABLE_ENCODER
839 if (encode)
840 ec_enc_uint(ec, itheta, qn+1);
841 else
842 #endif
843 itheta = ec_dec_uint(ec, qn+1);
844 } else {
845 int fs=1, ft;
846 ft = ((qn>>1)+1)*((qn>>1)+1);
847 #ifdef OPUS_ENABLE_ENCODER
848 if (encode)
849 {
850 int fl;
851
852 fs = itheta <= (qn>>1) ? itheta + 1 : qn + 1 - itheta;
853 fl = itheta <= (qn>>1) ? itheta*(itheta + 1)>>1 :
854 ft - ((qn + 1 - itheta)*(qn + 2 - itheta)>>1);
855
856 ec_encode(ec, fl, fl+fs, ft);
857 } else
858 #endif
859 {
860 /* Triangular pdf */
861 int fl=0;
862 int fm;
863 fm = ec_decode(ec, ft);
864
865 if (fm < ((qn>>1)*((qn>>1) + 1)>>1))
866 {
867 itheta = (isqrt32(8*(opus_uint32)fm + 1) - 1)>>1;
868 fs = itheta + 1;
869 fl = itheta*(itheta + 1)>>1;
870 }
871 else
872 {
873 itheta = (2*(qn + 1)
874 - isqrt32(8*(opus_uint32)(ft - fm - 1) + 1))>>1;
875 fs = qn + 1 - itheta;
876 fl = ft - ((qn + 1 - itheta)*(qn + 2 - itheta)>>1);
877 }
878
879 ec_dec_update(ec, fl, fl+fs, ft);
880 }
881 }
882 itheta = (opus_int32)itheta*16384/qn;
883 #ifdef OPUS_ENABLE_ENCODER
884 if (encode && stereo)
885 {
886 if (itheta==0)
887 intensity_stereo(m, X, Y, bandE, i, N);
888 else
889 stereo_split(X, Y, N);
890 }
891 #endif
892 /* NOTE: Renormalising X and Y *may* help fixed-point a bit at very high rate.
893 Let's do that at higher complexity */
894 } else if (stereo) {
895 #ifdef OPUS_ENABLE_ENCODER
896 if (encode)
897 {
898 inv = itheta > 8192;
899 if (inv)
900 {
901 int j;
902 for (j=0;j<N;j++)
903 Y[j] = -Y[j];
904 }
905 intensity_stereo(m, X, Y, bandE, i, N);
906 }
907 #endif
908 if (b>2<<BITRES && *remaining_bits > 2<<BITRES)
909 {
910 #ifdef OPUS_ENABLE_ENCODER
911 if (encode)
912 ec_enc_bit_logp(ec, inv, 2);
913 else
914 #endif
915 inv = ec_dec_bit_logp(ec, 2);
916 } else
917 inv = 0;
918 itheta = 0;
919 }
920 qalloc = ec_tell_frac(ec) - tell;
921 b -= qalloc;
922
923 orig_fill = fill;
924 if (itheta == 0)
925 {
926 imid = 32767;
927 iside = 0;
928 fill &= (1<<B)-1;
929 delta = -16384;
930 } else if (itheta == 16384)
931 {
932 imid = 0;
933 iside = 32767;
934 fill &= ((1<<B)-1)<<B;
935 delta = 16384;
936 } else {
937 imid = bitexact_cos((opus_int16)itheta);
938 iside = bitexact_cos((opus_int16)(16384-itheta));
939 /* This is the mid vs side allocation that minimizes squared error
940 in that band. */
941 delta = FRAC_MUL16((N-1)<<7,bitexact_log2tan(iside,imid));
942 }
943
944 #ifdef FIXED_POINT
945 mid = imid;
946 side = iside;
947 #else
948 mid = (1.f/32768)*imid;
949 side = (1.f/32768)*iside;
950 #endif
951
952 /* This is a special case for N=2 that only works for stereo and takes
953 advantage of the fact that mid and side are orthogonal to encode
954 the side with just one bit. */
955 if (N==2 && stereo)
956 {
957 int c;
958 int sign=0;
959 celt_norm *x2, *y2;
960 mbits = b;
961 sbits = 0;
962 /* Only need one bit for the side */
963 if (itheta != 0 && itheta != 16384)
964 sbits = 1<<BITRES;
965 mbits -= sbits;
966 c = itheta > 8192;
967 *remaining_bits -= qalloc+sbits;
968
969 x2 = c ? Y : X;
970 y2 = c ? X : Y;
971 if (sbits)
972 {
973 #ifdef OPUS_ENABLE_ENCODER
974 if (encode)
975 {
976 /* Here we only need to encode a sign for the side */
977 sign = x2[0]*y2[1] - x2[1]*y2[0] < 0;
978 ec_enc_bits(ec, sign, 1);
979 } else
980 #endif
981 {
982 sign = ec_dec_bits(ec, 1);
983 }
984 }
985 sign = 1-2*sign;
986 /* We use orig_fill here because we want to fold the side, but if
987 itheta==16384, we'll have cleared the low bits of fill. */
988 cm = quant_band(encode, m, i, x2, NULL, N, mbits, spread, B, intensity, tf_change, lowband, ec, remaining_bits, LM, lowband_out, NULL, level, seed, gain, lowband_scratch, orig_fill);
989 /* We don't split N=2 bands, so cm is either 1 or 0 (for a fold-collapse),
990 and there's no need to worry about mixing with the other channel. */
991 y2[0] = -sign*x2[1];
992 y2[1] = sign*x2[0];
993 if (resynth)
994 {
995 celt_norm tmp;
996 X[0] = MULT16_16_Q15(mid, X[0]);
997 X[1] = MULT16_16_Q15(mid, X[1]);
998 Y[0] = MULT16_16_Q15(side, Y[0]);
999 Y[1] = MULT16_16_Q15(side, Y[1]);
1000 tmp = X[0];
1001 X[0] = SUB16(tmp,Y[0]);
1002 Y[0] = ADD16(tmp,Y[0]);
1003 tmp = X[1];
1004 X[1] = SUB16(tmp,Y[1]);
1005 Y[1] = ADD16(tmp,Y[1]);
1006 }
1007 } else {
1008 /* "Normal" split code */
1009 celt_norm *next_lowband2=NULL;
1010 celt_norm *next_lowband_out1=NULL;
1011 int next_level=0;
1012 opus_int32 rebalance;
1013
1014 /* Give more bits to low-energy MDCTs than they would otherwise deserve */
1015 if (B0>1 && !stereo && (itheta&0x3fff))
1016 {
1017 if (itheta > 8192)
1018 /* Rough approximation for pre-echo masking */
1019 delta -= delta>>(4-LM);
1020 else
1021 /* Corresponds to a forward-masking slope of 1.5 dB per 10 ms */
1022 delta = IMIN(0, delta + (N<<BITRES>>(5-LM)));
1023 }
1024 mbits = IMAX(0, IMIN(b, (b-delta)/2));
1025 sbits = b-mbits;
1026 *remaining_bits -= qalloc;
1027
1028 if (lowband && !stereo)
1029 next_lowband2 = lowband+N; /* >32-bit split case */
1030
1031 /* Only stereo needs to pass on lowband_out. Otherwise, it's
1032 handled at the end */
1033 if (stereo)
1034 next_lowband_out1 = lowband_out;
1035 else
1036 next_level = level+1;
1037
1038 rebalance = *remaining_bits;
1039 if (mbits >= sbits)
1040 {
1041 /* In stereo mode, we do not apply a scaling to the mid because we need the normalized
1042 mid for folding later */
1043 cm = quant_band(encode, m, i, X, NULL, N, mbits, spread, B, intensity, tf_change,
1044 lowband, ec, remaining_bits, LM, next_lowband_out1,
1045 NULL, next_level, seed, stereo ? Q15ONE : MULT16_16_P15(gain,mid), lowband_scratch, fill);
1046 rebalance = mbits - (rebalance-*remaining_bits);
1047 if (rebalance > 3<<BITRES && itheta!=0)
1048 sbits += rebalance - (3<<BITRES);
1049
1050 /* For a stereo split, the high bits of fill are always zero, so no
1051 folding will be done to the side. */
1052 cm |= quant_band(encode, m, i, Y, NULL, N, sbits, spread, B, intensity, tf_change,
1053 next_lowband2, ec, remaining_bits, LM, NULL,
1054 NULL, next_level, seed, MULT16_16_P15(gain,side), NULL, fill>>B)<<((B0>>1)&(stereo-1));
1055 } else {
1056 /* For a stereo split, the high bits of fill are always zero, so no
1057 folding will be done to the side. */
1058 cm = quant_band(encode, m, i, Y, NULL, N, sbits, spread, B, intensity, tf_change,
1059 next_lowband2, ec, remaining_bits, LM, NULL,
1060 NULL, next_level, seed, MULT16_16_P15(gain,side), NULL, fill>>B)<<((B0>>1)&(stereo-1));
1061 rebalance = sbits - (rebalance-*remaining_bits);
1062 if (rebalance > 3<<BITRES && itheta!=16384)
1063 mbits += rebalance - (3<<BITRES);
1064 /* In stereo mode, we do not apply a scaling to the mid because we need the normalized
1065 mid for folding later */
1066 cm |= quant_band(encode, m, i, X, NULL, N, mbits, spread, B, intensity, tf_change,
1067 lowband, ec, remaining_bits, LM, next_lowband_out1,
1068 NULL, next_level, seed, stereo ? Q15ONE : MULT16_16_P15(gain,mid), lowband_scratch, fill);
1069 }
1070 }
1071
1072 } else {
1073 /* This is the basic no-split case */
1074 q = bits2pulses(m, i, LM, b);
1075 curr_bits = pulses2bits(m, i, LM, q);
1076 *remaining_bits -= curr_bits;
1077
1078 /* Ensures we can never bust the budget */
1079 while (*remaining_bits < 0 && q > 0)
1080 {
1081 *remaining_bits += curr_bits;
1082 q--;
1083 curr_bits = pulses2bits(m, i, LM, q);
1084 *remaining_bits -= curr_bits;
1085 }
1086
1087 if (q!=0)
1088 {
1089 int K = get_pulses(q);
1090
1091 /* Finally do the actual quantization */
1092 #ifdef OPUS_ENABLE_ENCODER
1093 if (encode)
1094 {
1095 cm = alg_quant(X, N, K, spread, B, ec
1096 #ifdef RESYNTH
1097 , gain
1098 #endif
1099 );
1100 } else
1101 #endif
1102 {
1103 cm = alg_unquant(X, N, K, spread, B, ec, gain);
1104 }
1105 } else {
1106 /* If there's no pulse, fill the band anyway */
1107 int j;
1108 if (resynth)
1109 {
1110 unsigned cm_mask;
1111 /*B can be as large as 16, so this shift might overflow an int on a
1112 16-bit platform; use a long to get defined behavior.*/
1113 cm_mask = (unsigned)(1UL<<B)-1;
1114 fill &= cm_mask;
1115 if (!fill)
1116 {
1117 for (j=0;j<N;j++)
1118 X[j] = 0;
1119 } else {
1120 if (lowband == NULL)
1121 {
1122 /* Noise */
1123 for (j=0;j<N;j++)
1124 {
1125 *seed = celt_lcg_rand(*seed);
1126 X[j] = (celt_norm)((opus_int32)*seed>>20);
1127 }
1128 cm = cm_mask;
1129 } else {
1130 /* Folded spectrum */
1131 for (j=0;j<N;j++)
1132 {
1133 opus_val16 tmp;
1134 *seed = celt_lcg_rand(*seed);
1135 /* About 48 dB below the "normal" folding level */
1136 tmp = QCONST16(1.0f/256, 10);
1137 tmp = (*seed)&0x8000 ? tmp : -tmp;
1138 X[j] = lowband[j]+tmp;
1139 }
1140 cm = fill;
1141 }
1142 renormalise_vector(X, N, gain);
1143 }
1144 }
1145 }
1146 }
1147
1148 /* This code is used by the decoder and by the resynthesis-enabled encoder */
1149 if (resynth)
1150 {
1151 if (stereo)
1152 {
1153 if (N!=2)
1154 stereo_merge(X, Y, mid, N);
1155 if (inv)
1156 {
1157 int j;
1158 for (j=0;j<N;j++)
1159 Y[j] = -Y[j];
1160 }
1161 } else if (level == 0)
1162 {
1163 int k;
1164
1165 /* Undo the sample reorganization going from time order to frequency order */
1166 if (B0>1)
1167 interleave_hadamard(X, N_B>>recombine, B0<<recombine, longBlocks);
1168
1169 /* Undo time-freq changes that we did earlier */
1170 N_B = N_B0;
1171 B = B0;
1172 for (k=0;k<time_divide;k++)
1173 {
1174 B >>= 1;
1175 N_B <<= 1;
1176 cm |= cm>>B;
1177 haar1(X, N_B, B);
1178 }
1179
1180 for (k=0;k<recombine;k++)
1181 {
1182 static const unsigned char bit_deinterleave_table[16]={
1183 0x00,0x03,0x0C,0x0F,0x30,0x33,0x3C,0x3F,
1184 0xC0,0xC3,0xCC,0xCF,0xF0,0xF3,0xFC,0xFF
1185 };
1186 cm = bit_deinterleave_table[cm];
1187 haar1(X, N0>>k, 1<<k);
1188 }
1189 B<<=recombine;
1190
1191 /* Scale output for later folding */
1192 if (lowband_out)
1193 {
1194 int j;
1195 opus_val16 n;
1196 n = celt_sqrt(SHL32(EXTEND32(N0),22));
1197 for (j=0;j<N0;j++)
1198 lowband_out[j] = MULT16_16_Q15(n,X[j]);
1199 }
1200 cm &= (1<<B)-1;
1201 }
1202 }
1203 return cm;
1204 }
1205
quant_all_bands(int encode,const CELTMode * m,int start,int end,celt_norm * X_,celt_norm * Y_,unsigned char * collapse_masks,const celt_ener * bandE,int * pulses,int shortBlocks,int spread,int dual_stereo,int intensity,int * tf_res,opus_int32 total_bits,opus_int32 balance,ec_ctx * ec,int LM,int codedBands,opus_uint32 * seed)1206 void quant_all_bands(int encode, const CELTMode *m, int start, int end,
1207 celt_norm *X_, celt_norm *Y_, unsigned char *collapse_masks, const celt_ener *bandE, int *pulses,
1208 int shortBlocks, int spread, int dual_stereo, int intensity, int *tf_res,
1209 opus_int32 total_bits, opus_int32 balance, ec_ctx *ec, int LM, int codedBands, opus_uint32 *seed)
1210 {
1211 int i;
1212 opus_int32 remaining_bits;
1213 const opus_int16 * OPUS_RESTRICT eBands = m->eBands;
1214 celt_norm * OPUS_RESTRICT norm, * OPUS_RESTRICT norm2;
1215 VARDECL(celt_norm, _norm);
1216 VARDECL(celt_norm, lowband_scratch);
1217 int B;
1218 int M;
1219 int lowband_offset;
1220 int update_lowband = 1;
1221 int C = Y_ != NULL ? 2 : 1;
1222 #ifdef RESYNTH
1223 int resynth = 1;
1224 #else
1225 int resynth = !encode;
1226 #endif
1227 SAVE_STACK;
1228
1229 M = 1<<LM;
1230 B = shortBlocks ? M : 1;
1231 ALLOC(_norm, C*M*eBands[m->nbEBands], celt_norm);
1232 ALLOC(lowband_scratch, M*(eBands[m->nbEBands]-eBands[m->nbEBands-1]), celt_norm);
1233 norm = _norm;
1234 norm2 = norm + M*eBands[m->nbEBands];
1235
1236 lowband_offset = 0;
1237 for (i=start;i<end;i++)
1238 {
1239 opus_int32 tell;
1240 int b;
1241 int N;
1242 opus_int32 curr_balance;
1243 int effective_lowband=-1;
1244 celt_norm * OPUS_RESTRICT X, * OPUS_RESTRICT Y;
1245 int tf_change=0;
1246 unsigned x_cm;
1247 unsigned y_cm;
1248
1249 X = X_+M*eBands[i];
1250 if (Y_!=NULL)
1251 Y = Y_+M*eBands[i];
1252 else
1253 Y = NULL;
1254 N = M*eBands[i+1]-M*eBands[i];
1255 tell = ec_tell_frac(ec);
1256
1257 /* Compute how many bits we want to allocate to this band */
1258 if (i != start)
1259 balance -= tell;
1260 remaining_bits = total_bits-tell-1;
1261 if (i <= codedBands-1)
1262 {
1263 curr_balance = balance / IMIN(3, codedBands-i);
1264 b = IMAX(0, IMIN(16383, IMIN(remaining_bits+1,pulses[i]+curr_balance)));
1265 } else {
1266 b = 0;
1267 }
1268
1269 if (resynth && M*eBands[i]-N >= M*eBands[start] && (update_lowband || lowband_offset==0))
1270 lowband_offset = i;
1271
1272 tf_change = tf_res[i];
1273 if (i>=m->effEBands)
1274 {
1275 X=norm;
1276 if (Y_!=NULL)
1277 Y = norm;
1278 }
1279
1280 /* Get a conservative estimate of the collapse_mask's for the bands we're
1281 going to be folding from. */
1282 if (lowband_offset != 0 && (spread!=SPREAD_AGGRESSIVE || B>1 || tf_change<0))
1283 {
1284 int fold_start;
1285 int fold_end;
1286 int fold_i;
1287 /* This ensures we never repeat spectral content within one band */
1288 effective_lowband = IMAX(M*eBands[start], M*eBands[lowband_offset]-N);
1289 fold_start = lowband_offset;
1290 while(M*eBands[--fold_start] > effective_lowband);
1291 fold_end = lowband_offset-1;
1292 while(M*eBands[++fold_end] < effective_lowband+N);
1293 x_cm = y_cm = 0;
1294 fold_i = fold_start; do {
1295 x_cm |= collapse_masks[fold_i*C+0];
1296 y_cm |= collapse_masks[fold_i*C+C-1];
1297 } while (++fold_i<fold_end);
1298 }
1299 /* Otherwise, we'll be using the LCG to fold, so all blocks will (almost
1300 always) be non-zero.*/
1301 else
1302 x_cm = y_cm = (1<<B)-1;
1303
1304 if (dual_stereo && i==intensity)
1305 {
1306 int j;
1307
1308 /* Switch off dual stereo to do intensity */
1309 dual_stereo = 0;
1310 if (resynth)
1311 for (j=M*eBands[start];j<M*eBands[i];j++)
1312 norm[j] = HALF32(norm[j]+norm2[j]);
1313 }
1314 if (dual_stereo)
1315 {
1316 x_cm = quant_band(encode, m, i, X, NULL, N, b/2, spread, B, intensity, tf_change,
1317 effective_lowband != -1 ? norm+effective_lowband : NULL, ec, &remaining_bits, LM,
1318 norm+M*eBands[i], bandE, 0, seed, Q15ONE, lowband_scratch, x_cm);
1319 y_cm = quant_band(encode, m, i, Y, NULL, N, b/2, spread, B, intensity, tf_change,
1320 effective_lowband != -1 ? norm2+effective_lowband : NULL, ec, &remaining_bits, LM,
1321 norm2+M*eBands[i], bandE, 0, seed, Q15ONE, lowband_scratch, y_cm);
1322 } else {
1323 x_cm = quant_band(encode, m, i, X, Y, N, b, spread, B, intensity, tf_change,
1324 effective_lowband != -1 ? norm+effective_lowband : NULL, ec, &remaining_bits, LM,
1325 norm+M*eBands[i], bandE, 0, seed, Q15ONE, lowband_scratch, x_cm|y_cm);
1326 y_cm = x_cm;
1327 }
1328 collapse_masks[i*C+0] = (unsigned char)x_cm;
1329 collapse_masks[i*C+C-1] = (unsigned char)y_cm;
1330 balance += pulses[i] + tell;
1331
1332 /* Update the folding position only as long as we have 1 bit/sample depth */
1333 update_lowband = b>(N<<BITRES);
1334 }
1335 RESTORE_STACK;
1336 }
1337
1338