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