Bug Summary

File:libs/opus-1.1-p2/celt/bands.c
Location:line 571, column 12
Description:Assigned value is garbage or undefined

Annotated Source Code

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