nothing to see here
[opus.git] / libcelt / bands.c
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    - Neither the name of the Xiph.org Foundation nor the names of its
18    contributors may be used to endorse or promote products derived from
19    this software without specific prior written permission.
20    
21    THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
22    ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
23    LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
24    A PARTICULAR PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE FOUNDATION OR
25    CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
26    EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
27    PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
28    PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
29    LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
30    NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
31    SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
32 */
33
34 #ifdef HAVE_CONFIG_H
35 #include "config.h"
36 #endif
37
38 #include <math.h>
39 #include "bands.h"
40 #include "modes.h"
41 #include "vq.h"
42 #include "cwrs.h"
43 #include "stack_alloc.h"
44 #include "os_support.h"
45 #include "mathops.h"
46 #include "rate.h"
47
48 /* This is a cos() approximation designed to be bit-exact on any platform. Bit exactness
49    with this approximation is important because it has an impact on the bit allocation */
50 static celt_int16 bitexact_cos(celt_int16 x)
51 {
52    celt_int32 tmp;
53    celt_int16 x2;
54    tmp = (4096+((celt_int32)(x)*(x)))>>13;
55    if (tmp > 32767)
56       tmp = 32767;
57    x2 = tmp;
58    x2 = (32767-x2) + FRAC_MUL16(x2, (-7651 + FRAC_MUL16(x2, (8277 + FRAC_MUL16(-626, x2)))));
59    if (x2 > 32766)
60       x2 = 32766;
61    return 1+x2;
62 }
63
64
65 #ifdef FIXED_POINT
66 /* Compute the amplitude (sqrt energy) in each of the bands */
67 void compute_band_energies(const CELTMode *m, const celt_sig *X, celt_ener *bank, int end, int _C, int M)
68 {
69    int i, c, N;
70    const celt_int16 *eBands = m->eBands;
71    const int C = CHANNELS(_C);
72    N = M*m->shortMdctSize;
73    c=0; do {
74       for (i=0;i<end;i++)
75       {
76          int j;
77          celt_word32 maxval=0;
78          celt_word32 sum = 0;
79          
80          j=M*eBands[i]; do {
81             maxval = MAX32(maxval, X[j+c*N]);
82             maxval = MAX32(maxval, -X[j+c*N]);
83          } while (++j<M*eBands[i+1]);
84          
85          if (maxval > 0)
86          {
87             int shift = celt_ilog2(maxval)-10;
88             j=M*eBands[i]; do {
89                sum = MAC16_16(sum, EXTRACT16(VSHR32(X[j+c*N],shift)),
90                                    EXTRACT16(VSHR32(X[j+c*N],shift)));
91             } while (++j<M*eBands[i+1]);
92             /* We're adding one here to make damn sure we never end up with a pitch vector that's
93                larger than unity norm */
94             bank[i+c*m->nbEBands] = EPSILON+VSHR32(EXTEND32(celt_sqrt(sum)),-shift);
95          } else {
96             bank[i+c*m->nbEBands] = EPSILON;
97          }
98          /*printf ("%f ", bank[i+c*m->nbEBands]);*/
99       }
100    } while (++c<C);
101    /*printf ("\n");*/
102 }
103
104 /* Normalise each band such that the energy is one. */
105 void normalise_bands(const CELTMode *m, const celt_sig * restrict freq, celt_norm * restrict X, const celt_ener *bank, int end, int _C, int M)
106 {
107    int i, c, N;
108    const celt_int16 *eBands = m->eBands;
109    const int C = CHANNELS(_C);
110    N = M*m->shortMdctSize;
111    c=0; do {
112       i=0; do {
113          celt_word16 g;
114          int j,shift;
115          celt_word16 E;
116          shift = celt_zlog2(bank[i+c*m->nbEBands])-13;
117          E = VSHR32(bank[i+c*m->nbEBands], shift);
118          g = EXTRACT16(celt_rcp(SHL32(E,3)));
119          j=M*eBands[i]; do {
120             X[j+c*N] = MULT16_16_Q15(VSHR32(freq[j+c*N],shift-1),g);
121          } while (++j<M*eBands[i+1]);
122       } while (++i<end);
123    } while (++c<C);
124 }
125
126 #else /* FIXED_POINT */
127 /* Compute the amplitude (sqrt energy) in each of the bands */
128 void compute_band_energies(const CELTMode *m, const celt_sig *X, celt_ener *bank, int end, int _C, int M)
129 {
130    int i, c, N;
131    const celt_int16 *eBands = m->eBands;
132    const int C = CHANNELS(_C);
133    N = M*m->shortMdctSize;
134    c=0; do {
135       for (i=0;i<end;i++)
136       {
137          int j;
138          celt_word32 sum = 1e-10f;
139          for (j=M*eBands[i];j<M*eBands[i+1];j++)
140             sum += X[j+c*N]*X[j+c*N];
141          bank[i+c*m->nbEBands] = celt_sqrt(sum);
142          /*printf ("%f ", bank[i+c*m->nbEBands]);*/
143       }
144    } while (++c<C);
145    /*printf ("\n");*/
146 }
147
148 /* Normalise each band such that the energy is one. */
149 void normalise_bands(const CELTMode *m, const celt_sig * restrict freq, celt_norm * restrict X, const celt_ener *bank, int end, int _C, int M)
150 {
151    int i, c, N;
152    const celt_int16 *eBands = m->eBands;
153    const int C = CHANNELS(_C);
154    N = M*m->shortMdctSize;
155    c=0; do {
156       for (i=0;i<end;i++)
157       {
158          int j;
159          celt_word16 g = 1.f/(1e-10f+bank[i+c*m->nbEBands]);
160          for (j=M*eBands[i];j<M*eBands[i+1];j++)
161             X[j+c*N] = freq[j+c*N]*g;
162       }
163    } while (++c<C);
164 }
165
166 #endif /* FIXED_POINT */
167
168 /* De-normalise the energy to produce the synthesis from the unit-energy bands */
169 void denormalise_bands(const CELTMode *m, const celt_norm * restrict X, celt_sig * restrict freq, const celt_ener *bank, int end, int _C, int M)
170 {
171    int i, c, N;
172    const celt_int16 *eBands = m->eBands;
173    const int C = CHANNELS(_C);
174    N = M*m->shortMdctSize;
175    celt_assert2(C<=2, "denormalise_bands() not implemented for >2 channels");
176    c=0; do {
177       celt_sig * restrict f;
178       const celt_norm * restrict x;
179       f = freq+c*N;
180       x = X+c*N;
181       for (i=0;i<end;i++)
182       {
183          int j, band_end;
184          celt_word32 g = SHR32(bank[i+c*m->nbEBands],1);
185          j=M*eBands[i];
186          band_end = M*eBands[i+1];
187          do {
188             *f++ = SHL32(MULT16_32_Q15(*x, g),2);
189             x++;
190          } while (++j<band_end);
191       }
192       for (i=M*eBands[m->nbEBands];i<N;i++)
193          *f++ = 0;
194    } while (++c<C);
195 }
196
197 static void intensity_stereo(const CELTMode *m, celt_norm *X, celt_norm *Y, const celt_ener *bank, int bandID, int N)
198 {
199    int i = bandID;
200    int j;
201    celt_word16 a1, a2;
202    celt_word16 left, right;
203    celt_word16 norm;
204 #ifdef FIXED_POINT
205    int shift = celt_zlog2(MAX32(bank[i], bank[i+m->nbEBands]))-13;
206 #endif
207    left = VSHR32(bank[i],shift);
208    right = VSHR32(bank[i+m->nbEBands],shift);
209    norm = EPSILON + celt_sqrt(EPSILON+MULT16_16(left,left)+MULT16_16(right,right));
210    a1 = DIV32_16(SHL32(EXTEND32(left),14),norm);
211    a2 = DIV32_16(SHL32(EXTEND32(right),14),norm);
212    for (j=0;j<N;j++)
213    {
214       celt_norm r, l;
215       l = X[j];
216       r = Y[j];
217       X[j] = MULT16_16_Q14(a1,l) + MULT16_16_Q14(a2,r);
218       /* Side is not encoded, no need to calculate */
219    }
220 }
221
222 static void stereo_split(celt_norm *X, celt_norm *Y, int N)
223 {
224    int j;
225    for (j=0;j<N;j++)
226    {
227       celt_norm r, l;
228       l = MULT16_16_Q15(QCONST16(.70711f,15), X[j]);
229       r = MULT16_16_Q15(QCONST16(.70711f,15), Y[j]);
230       X[j] = l+r;
231       Y[j] = r-l;
232    }
233 }
234
235 static void stereo_merge(celt_norm *X, celt_norm *Y, celt_word16 mid, celt_word16 side, int N)
236 {
237    int j;
238    celt_word32 xp=0;
239    celt_word32 El, Er;
240 #ifdef FIXED_POINT
241    int kl, kr;
242 #endif
243    celt_word32 t, lgain, rgain;
244
245    /* Compute the norm of X+Y and X-Y as |X|^2 + |Y|^2 +/- sum(xy) */
246    for (j=0;j<N;j++)
247       xp = MAC16_16(xp, X[j], Y[j]);
248    /* mid and side are in Q15, not Q14 like X and Y */
249    mid = SHR32(mid, 1);
250    side = SHR32(side, 1);
251    El = MULT16_16(mid, mid) + MULT16_16(side, side) - 2*xp;
252    Er = MULT16_16(mid, mid) + MULT16_16(side, side) + 2*xp;
253    if (Er < EPSILON)
254       Er = EPSILON;
255    if (El < EPSILON)
256       El = EPSILON;
257
258 #ifdef FIXED_POINT
259    kl = celt_ilog2(El)>>1;
260    kr = celt_ilog2(Er)>>1;
261 #endif
262    t = VSHR32(El, (kl-7)<<1);
263    lgain = celt_rsqrt_norm(t);
264    t = VSHR32(Er, (kr-7)<<1);
265    rgain = celt_rsqrt_norm(t);
266
267 #ifdef FIXED_POINT
268    if (kl < 7)
269       kl = 7;
270    if (kr < 7)
271       kr = 7;
272 #endif
273
274    for (j=0;j<N;j++)
275    {
276       celt_norm r, l;
277       l = X[j];
278       r = Y[j];
279       X[j] = EXTRACT16(PSHR32(MULT16_16(lgain, SUB16(l,r)), kl+1));
280       Y[j] = EXTRACT16(PSHR32(MULT16_16(rgain, ADD16(l,r)), kr+1));
281    }
282 }
283
284 /* Decide whether we should spread the pulses in the current frame */
285 int folding_decision(const CELTMode *m, celt_norm *X, int *average, int *last_decision, int end, int _C, int M)
286 {
287    int i, c, N0;
288    int sum = 0, nbBands=0;
289    const int C = CHANNELS(_C);
290    const celt_int16 * restrict eBands = m->eBands;
291    int decision;
292    
293    N0 = M*m->shortMdctSize;
294
295    if (M*(eBands[end]-eBands[end-1]) <= 8)
296       return 0;
297    c=0; do {
298       for (i=0;i<end;i++)
299       {
300          int j, N, tmp=0;
301          int tcount[3] = {0};
302          celt_norm * restrict x = X+M*eBands[i]+c*N0;
303          N = M*(eBands[i+1]-eBands[i]);
304          if (N<=8)
305             continue;
306          /* Compute rough CDF of |x[j]| */
307          for (j=0;j<N;j++)
308          {
309             celt_word32 x2N; /* Q13 */
310
311             x2N = MULT16_16(MULT16_16_Q15(x[j], x[j]), N);
312             if (x2N < QCONST16(0.25f,13))
313                tcount[0]++;
314             if (x2N < QCONST16(0.0625f,13))
315                tcount[1]++;
316             if (x2N < QCONST16(0.015625f,13))
317                tcount[2]++;
318          }
319
320          tmp = (2*tcount[2] >= N) + (2*tcount[1] >= N) + (2*tcount[0] >= N);
321          sum += tmp*256;
322          nbBands++;
323       }
324    } while (++c<C);
325    sum /= nbBands;
326    /* Recursive averaging */
327    sum = (sum+*average)>>1;
328    *average = sum;
329    /* Hysteresis */
330    sum = (3*sum + ((*last_decision<<7) + 64) + 2)>>2;
331    /* decision and last_decision do not use the same encoding */
332    if (sum < 80)
333    {
334       decision = 2;
335       *last_decision = 0;
336    } else if (sum < 256)
337    {
338       decision = 1;
339       *last_decision = 1;
340    } else if (sum < 384)
341    {
342       decision = 3;
343       *last_decision = 2;
344    } else {
345       decision = 0;
346       *last_decision = 3;
347    }
348    return decision;
349 }
350
351 #ifdef MEASURE_NORM_MSE
352
353 float MSE[30] = {0};
354 int nbMSEBands = 0;
355 int MSECount[30] = {0};
356
357 void dump_norm_mse(void)
358 {
359    int i;
360    for (i=0;i<nbMSEBands;i++)
361    {
362       printf ("%g ", MSE[i]/MSECount[i]);
363    }
364    printf ("\n");
365 }
366
367 void measure_norm_mse(const CELTMode *m, float *X, float *X0, float *bandE, float *bandE0, int M, int N, int C)
368 {
369    static int init = 0;
370    int i;
371    if (!init)
372    {
373       atexit(dump_norm_mse);
374       init = 1;
375    }
376    for (i=0;i<m->nbEBands;i++)
377    {
378       int j;
379       int c;
380       float g;
381       if (bandE0[i]<10 || (C==2 && bandE0[i+m->nbEBands]<1))
382          continue;
383       c=0; do {
384          g = bandE[i+c*m->nbEBands]/(1e-15+bandE0[i+c*m->nbEBands]);
385          for (j=M*m->eBands[i];j<M*m->eBands[i+1];j++)
386             MSE[i] += (g*X[j+c*N]-X0[j+c*N])*(g*X[j+c*N]-X0[j+c*N]);
387       } while (++c<C);
388       MSECount[i]+=C;
389    }
390    nbMSEBands = m->nbEBands;
391 }
392
393 #endif
394
395 static void interleave_vector(celt_norm *X, int N0, int stride)
396 {
397    int i,j;
398    VARDECL(celt_norm, tmp);
399    int N;
400    SAVE_STACK;
401    N = N0*stride;
402    ALLOC(tmp, N, celt_norm);
403    for (i=0;i<stride;i++)
404       for (j=0;j<N0;j++)
405          tmp[j*stride+i] = X[i*N0+j];
406    for (j=0;j<N;j++)
407       X[j] = tmp[j];
408    RESTORE_STACK;
409 }
410
411 static void deinterleave_vector(celt_norm *X, int N0, int stride)
412 {
413    int i,j;
414    VARDECL(celt_norm, tmp);
415    int N;
416    SAVE_STACK;
417    N = N0*stride;
418    ALLOC(tmp, N, celt_norm);
419    for (i=0;i<stride;i++)
420       for (j=0;j<N0;j++)
421          tmp[i*N0+j] = X[j*stride+i];
422    for (j=0;j<N;j++)
423       X[j] = tmp[j];
424    RESTORE_STACK;
425 }
426
427 void haar1(celt_norm *X, int N0, int stride)
428 {
429    int i, j;
430    N0 >>= 1;
431    for (i=0;i<stride;i++)
432       for (j=0;j<N0;j++)
433       {
434          celt_norm tmp1, tmp2;
435          tmp1 = MULT16_16_Q15(QCONST16(.7070678f,15), X[stride*2*j+i]);
436          tmp2 = MULT16_16_Q15(QCONST16(.7070678f,15), X[stride*(2*j+1)+i]);
437          X[stride*2*j+i] = tmp1 + tmp2;
438          X[stride*(2*j+1)+i] = tmp1 - tmp2;
439       }
440 }
441
442 static int compute_qn(int N, int b, int offset, int stereo)
443 {
444    static const celt_int16 exp2_table8[8] =
445       {16384, 17866, 19483, 21247, 23170, 25267, 27554, 30048};
446    int qn, qb;
447    int N2 = 2*N-1;
448    if (stereo && N==2)
449       N2--;
450    qb = (b+N2*offset)/N2;
451    if (qb > (b>>1)-(1<<BITRES))
452       qb = (b>>1)-(1<<BITRES);
453
454    if (qb<0)
455        qb = 0;
456    if (qb>8<<BITRES)
457      qb = 8<<BITRES;
458
459    if (qb<(1<<BITRES>>1)) {
460       qn = 1;
461    } else {
462       qn = exp2_table8[qb&0x7]>>(14-(qb>>BITRES));
463       qn = (qn+1)>>1<<1;
464    }
465    celt_assert(qn <= 256);
466    return qn;
467 }
468
469 static celt_uint32 lcg_rand(celt_uint32 seed)
470 {
471    return 1664525 * seed + 1013904223;
472 }
473
474 /* This function is responsible for encoding and decoding a band for both
475    the mono and stereo case. Even in the mono case, it can split the band
476    in two and transmit the energy difference with the two half-bands. It
477    can be called recursively so bands can end up being split in 8 parts. */
478 static void quant_band(int encode, const CELTMode *m, int i, celt_norm *X, celt_norm *Y,
479       int N, int b, int spread, int B, int intensity, int tf_change, celt_norm *lowband, int resynth, void *ec,
480       celt_int32 *remaining_bits, int LM, celt_norm *lowband_out, const celt_ener *bandE, int level,
481       celt_int32 *seed, celt_word16 gain, celt_norm *lowband_scratch)
482 {
483    int q;
484    int curr_bits;
485    int stereo, split;
486    int imid=0, iside=0;
487    int N0=N;
488    int N_B=N;
489    int N_B0;
490    int B0=B;
491    int time_divide=0;
492    int recombine=0;
493    int inv = 0;
494    celt_word16 mid=0, side=0;
495
496    N_B /= B;
497    N_B0 = N_B;
498
499    split = stereo = Y != NULL;
500
501    /* Special case for one sample */
502    if (N==1)
503    {
504       int c;
505       celt_norm *x = X;
506       c=0; do {
507          int sign=0;
508          if (*remaining_bits>=1<<BITRES)
509          {
510             if (encode)
511             {
512                sign = x[0]<0;
513                ec_enc_bits((ec_enc*)ec, sign, 1);
514             } else {
515                sign = ec_dec_bits((ec_dec*)ec, 1);
516             }
517             *remaining_bits -= 1<<BITRES;
518             b-=1<<BITRES;
519          }
520          if (resynth)
521             x[0] = sign ? -NORM_SCALING : NORM_SCALING;
522          x = Y;
523       } while (++c<1+stereo);
524       if (lowband_out)
525          lowband_out[0] = SHR16(X[0],4);
526       return;
527    }
528
529    if (!stereo && level == 0)
530    {
531       int k;
532       if (tf_change>0)
533          recombine = tf_change;
534       /* Band recombining to increase frequency resolution */
535
536       if (lowband && (recombine || ((N_B&1) == 0 && tf_change<0) || B0>1))
537       {
538          int j;
539          for (j=0;j<N;j++)
540             lowband_scratch[j] = lowband[j];
541          lowband = lowband_scratch;
542       }
543
544       for (k=0;k<recombine;k++)
545       {
546          B>>=1;
547          N_B<<=1;
548          if (encode)
549             haar1(X, N_B, B);
550          if (lowband)
551             haar1(lowband, N_B, B);
552       }
553
554       /* Increasing the time resolution */
555       while ((N_B&1) == 0 && tf_change<0)
556       {
557          if (encode)
558             haar1(X, N_B, B);
559          if (lowband)
560             haar1(lowband, N_B, B);
561          B <<= 1;
562          N_B >>= 1;
563          time_divide++;
564          tf_change++;
565       }
566       B0=B;
567       N_B0 = N_B;
568
569       /* Reorganize the samples in time order instead of frequency order */
570       if (B0>1)
571       {
572          if (encode)
573             deinterleave_vector(X, N_B, B0);
574          if (lowband)
575             deinterleave_vector(lowband, N_B, B0);
576       }
577    }
578
579    /* If we need more than 32 bits, try splitting the band in two. */
580    if (!stereo && LM != -1 && b > 32<<BITRES && N>2)
581    {
582       if (LM>0 || (N&1)==0)
583       {
584          N >>= 1;
585          Y = X+N;
586          split = 1;
587          LM -= 1;
588          B = (B+1)>>1;
589       }
590    }
591
592    if (split)
593    {
594       int qn;
595       int itheta=0;
596       int mbits, sbits, delta;
597       int qalloc;
598       int offset;
599
600       /* Decide on the resolution to give to the split parameter theta */
601       offset = ((m->logN[i]+(LM<<BITRES))>>1) - (stereo ? QTHETA_OFFSET_STEREO : QTHETA_OFFSET);
602       qn = compute_qn(N, b, offset, stereo);
603       qalloc = 0;
604       if (stereo && i>=intensity)
605          qn = 1;
606       if (qn!=1)
607       {
608          if (encode)
609          {
610             /* theta is the atan() of the ratio between the (normalized)
611                side and mid. With just that parameter, we can re-scale both
612                mid and side because we know that 1) they have unit norm and
613                2) they are orthogonal. */
614             itheta = stereo_itheta(X, Y, stereo, N);
615
616             itheta = (itheta*qn+8192)>>14;
617          }
618
619          /* Entropy coding of the angle. We use a uniform pdf for the
620             first stereo split but a triangular one for the rest. */
621          if (stereo || B>1)
622          {
623             if (encode)
624                ec_enc_uint((ec_enc*)ec, itheta, qn+1);
625             else
626                itheta = ec_dec_uint((ec_dec*)ec, qn+1);
627             qalloc = log2_frac(qn+1,BITRES);
628          } else {
629             int fs=1, ft;
630             ft = ((qn>>1)+1)*((qn>>1)+1);
631             if (encode)
632             {
633                int fl;
634
635                fs = itheta <= (qn>>1) ? itheta + 1 : qn + 1 - itheta;
636                fl = itheta <= (qn>>1) ? itheta*(itheta + 1)>>1 :
637                 ft - ((qn + 1 - itheta)*(qn + 2 - itheta)>>1);
638
639                ec_encode((ec_enc*)ec, fl, fl+fs, ft);
640             } else {
641                int fl=0;
642                int fm;
643                fm = ec_decode((ec_dec*)ec, ft);
644
645                if (fm < ((qn>>1)*((qn>>1) + 1)>>1))
646                {
647                   itheta = (isqrt32(8*(celt_uint32)fm + 1) - 1)>>1;
648                   fs = itheta + 1;
649                   fl = itheta*(itheta + 1)>>1;
650                }
651                else
652                {
653                   itheta = (2*(qn + 1)
654                    - isqrt32(8*(celt_uint32)(ft - fm - 1) + 1))>>1;
655                   fs = qn + 1 - itheta;
656                   fl = ft - ((qn + 1 - itheta)*(qn + 2 - itheta)>>1);
657                }
658
659                ec_dec_update((ec_dec*)ec, fl, fl+fs, ft);
660             }
661             qalloc = log2_frac(ft,BITRES) - log2_frac(fs,BITRES) + 1;
662          }
663          itheta = (celt_int32)itheta*16384/qn;
664          if (encode && stereo)
665             stereo_split(X, Y, N);
666          /* TODO: Renormalising X and Y *may* help fixed-point a bit at very high rate.
667                   Let's do that at higher complexity */
668       } else if (stereo) {
669          if (encode)
670          {
671             inv = itheta > 8192;
672             if (inv)
673             {
674                int j;
675                for (j=0;j<N;j++)
676                   Y[j] = -Y[j];
677             }
678             intensity_stereo(m, X, Y, bandE, i, N);
679          }
680          if (b>2<<BITRES)
681          {
682             if (encode)
683                ec_enc_bit_prob(ec, inv, 16384);
684             else
685                inv = ec_dec_bit_prob(ec, 16384);
686             qalloc = 1<<BITRES;
687          } else
688             inv = 0;
689       }
690
691       if (itheta == 0)
692       {
693          imid = 32767;
694          iside = 0;
695          delta = -10000;
696       } else if (itheta == 16384)
697       {
698          imid = 0;
699          iside = 32767;
700          delta = 10000;
701       } else {
702          imid = bitexact_cos(itheta);
703          iside = bitexact_cos(16384-itheta);
704          /* This is the mid vs side allocation that minimizes squared error
705             in that band. */
706          delta = (N-1)*(log2_frac(iside,BITRES+2)-log2_frac(imid,BITRES+2))>>2;
707       }
708
709 #ifdef FIXED_POINT
710       mid = imid;
711       side = iside;
712 #else
713       mid = (1.f/32768)*imid;
714       side = (1.f/32768)*iside;
715 #endif
716
717       /* This is a special case for N=2 that only works for stereo and takes
718          advantage of the fact that mid and side are orthogonal to encode
719          the side with just one bit. */
720       if (N==2 && stereo)
721       {
722          int c;
723          int sign=1;
724          celt_norm *x2, *y2;
725          mbits = b-qalloc;
726          sbits = 0;
727          /* Only need one bit for the side */
728          if (itheta != 0 && itheta != 16384)
729             sbits = 1<<BITRES;
730          mbits -= sbits;
731          c = itheta > 8192;
732          *remaining_bits -= qalloc+sbits;
733
734          x2 = c ? Y : X;
735          y2 = c ? X : Y;
736          if (sbits)
737          {
738             if (encode)
739             {
740                /* Here we only need to encode a sign for the side */
741                sign = x2[0]*y2[1] - x2[1]*y2[0] > 0;
742                ec_enc_bits((ec_enc*)ec, sign, 1);
743             } else {
744                sign = ec_dec_bits((ec_dec*)ec, 1);
745             }
746          }
747          sign = 2*sign - 1;
748          quant_band(encode, m, i, x2, NULL, N, mbits, spread, B, intensity, tf_change, lowband, resynth, ec, remaining_bits, LM, lowband_out, NULL, level, seed, gain, lowband_scratch);
749          y2[0] = -sign*x2[1];
750          y2[1] = sign*x2[0];
751          if (resynth)
752          {
753             celt_norm tmp;
754             X[0] = MULT16_16_Q15(mid, X[0]);
755             X[1] = MULT16_16_Q15(mid, X[1]);
756             Y[0] = MULT16_16_Q15(side, Y[0]);
757             Y[1] = MULT16_16_Q15(side, Y[1]);
758             tmp = X[0];
759             X[0] = SUB16(tmp,Y[0]);
760             Y[0] = ADD16(tmp,Y[0]);
761             tmp = X[1];
762             X[1] = SUB16(tmp,Y[1]);
763             Y[1] = ADD16(tmp,Y[1]);
764          }
765       } else {
766          /* "Normal" split code */
767          celt_norm *next_lowband2=NULL;
768          celt_norm *next_lowband_out1=NULL;
769          int next_level=0;
770
771          /* Give more bits to low-energy MDCTs than they would otherwise deserve */
772          if (B>1 && !stereo && itheta > 8192)
773             delta -= delta>>(1+level);
774
775          mbits = (b-qalloc-delta)/2;
776          if (mbits > b-qalloc)
777             mbits = b-qalloc;
778          if (mbits<0)
779             mbits=0;
780          sbits = b-qalloc-mbits;
781          *remaining_bits -= qalloc;
782
783          if (lowband && !stereo)
784             next_lowband2 = lowband+N; /* >32-bit split case */
785
786          /* Only stereo needs to pass on lowband_out. Otherwise, it's
787             handled at the end */
788          if (stereo)
789             next_lowband_out1 = lowband_out;
790          else
791             next_level = level+1;
792
793          quant_band(encode, m, i, X, NULL, N, mbits, spread, B, intensity, tf_change,
794                lowband, resynth, ec, remaining_bits, LM, next_lowband_out1,
795                NULL, next_level, seed, MULT16_16_P15(gain,mid), lowband_scratch);
796          quant_band(encode, m, i, Y, NULL, N, sbits, spread, B, intensity, tf_change,
797                next_lowband2, resynth, ec, remaining_bits, LM, NULL,
798                NULL, next_level, seed, MULT16_16_P15(gain,side), NULL);
799       }
800
801    } else {
802       /* This is the basic no-split case */
803       q = bits2pulses(m, i, LM, b);
804       curr_bits = pulses2bits(m, i, LM, q);
805       *remaining_bits -= curr_bits;
806
807       /* Ensures we can never bust the budget */
808       while (*remaining_bits < 0 && q > 0)
809       {
810          *remaining_bits += curr_bits;
811          q--;
812          curr_bits = pulses2bits(m, i, LM, q);
813          *remaining_bits -= curr_bits;
814       }
815
816       if (q!=0)
817       {
818          int K = get_pulses(q);
819
820          /* Finally do the actual quantization */
821          if (encode)
822             alg_quant(X, N, K, spread, B, lowband, resynth, (ec_enc*)ec, seed, gain);
823          else
824             alg_unquant(X, N, K, spread, B, lowband, (ec_dec*)ec, seed, gain);
825       } else {
826          /* If there's no pulse, fill the band anyway */
827          int j;
828          if (lowband != NULL && resynth)
829          {
830             if (spread==2 && B<=1)
831             {
832                /* Folded spectrum */
833                for (j=0;j<N;j++)
834                {
835                   *seed = lcg_rand(*seed);
836                   X[j] = (int)(*seed)>>20;
837                }
838             } else {
839                /* Noise */
840                for (j=0;j<N;j++)
841                   X[j] = lowband[j];
842             }
843             renormalise_vector(X, N, gain);
844          } else {
845             /* This is important for encoding the side in stereo mode */
846             for (j=0;j<N;j++)
847                X[j] = 0;
848          }
849       }
850    }
851
852    /* This code is used by the decoder and by the resynthesis-enabled encoder */
853    if (resynth)
854    {
855       if (stereo)
856       {
857          if (N!=2)
858             stereo_merge(X, Y, mid, side, N);
859          if (inv)
860          {
861             int j;
862             for (j=0;j<N;j++)
863                Y[j] = -Y[j];
864          }
865       } else if (level == 0)
866       {
867          int k;
868
869          /* Undo the sample reorganization going from time order to frequency order */
870          if (B0>1)
871             interleave_vector(X, N_B, B0);
872
873          /* Undo time-freq changes that we did earlier */
874          N_B = N_B0;
875          B = B0;
876          for (k=0;k<time_divide;k++)
877          {
878             B >>= 1;
879             N_B <<= 1;
880             haar1(X, N_B, B);
881          }
882
883          for (k=0;k<recombine;k++)
884          {
885             haar1(X, N_B, B);
886             N_B>>=1;
887             B <<= 1;
888          }
889
890          /* Scale output for later folding */
891          if (lowband_out)
892          {
893             int j;
894             celt_word16 n;
895             n = celt_sqrt(SHL32(EXTEND32(N0),22));
896             for (j=0;j<N0;j++)
897                lowband_out[j] = MULT16_16_Q15(n,X[j]);
898          }
899       }
900    }
901 }
902
903 void quant_all_bands(int encode, const CELTMode *m, int start, int end,
904       celt_norm *_X, celt_norm *_Y, const celt_ener *bandE, int *pulses,
905       int shortBlocks, int fold, int intensity, int *tf_res, int resynth,
906       int total_bits, void *ec, int LM, int codedBands)
907 {
908    int i, balance;
909    celt_int32 remaining_bits;
910    const celt_int16 * restrict eBands = m->eBands;
911    celt_norm * restrict norm;
912    VARDECL(celt_norm, _norm);
913    VARDECL(celt_norm, lowband_scratch);
914    int B;
915    int M;
916    celt_int32 seed;
917    int lowband_offset;
918    int update_lowband = 1;
919    int C = _Y != NULL ? 2 : 1;
920    SAVE_STACK;
921
922    M = 1<<LM;
923    B = shortBlocks ? M : 1;
924    ALLOC(_norm, M*eBands[m->nbEBands], celt_norm);
925    ALLOC(lowband_scratch, M*(eBands[m->nbEBands]-eBands[m->nbEBands-1]), celt_norm);
926    norm = _norm;
927
928    if (C==2)
929    {
930       int j;
931       int left = 0;
932       for (j=intensity;j<codedBands;j++)
933       {
934          int tmp = pulses[j]/2;
935          left += tmp;
936          pulses[j] -= tmp;
937       }
938       if (codedBands) {
939          int perband;
940          perband = left/(m->eBands[codedBands]-m->eBands[start]);
941          for (j=start;j<codedBands;j++)
942             pulses[j] += perband*(m->eBands[j+1]-m->eBands[j]);
943          left = left-(m->eBands[codedBands]-m->eBands[start])*perband;
944          for (j=start;j<codedBands;j++)
945          {
946             int tmp = IMIN(left, m->eBands[j+1]-m->eBands[j]);
947             pulses[j] += tmp;
948             left -= tmp;
949          }
950       }
951    }
952
953    if (encode)
954       seed = ((ec_enc*)ec)->rng;
955    else
956       seed = ((ec_dec*)ec)->rng;
957    balance = 0;
958    lowband_offset = -1;
959    for (i=start;i<end;i++)
960    {
961       int tell;
962       int b;
963       int N;
964       int curr_balance;
965       int effective_lowband=-1;
966       celt_norm * restrict X, * restrict Y;
967       int tf_change=0;
968       
969       X = _X+M*eBands[i];
970       if (_Y!=NULL)
971          Y = _Y+M*eBands[i];
972       else
973          Y = NULL;
974       N = M*eBands[i+1]-M*eBands[i];
975       if (encode)
976          tell = ec_enc_tell((ec_enc*)ec, BITRES);
977       else
978          tell = ec_dec_tell((ec_dec*)ec, BITRES);
979
980       /* Compute how many bits we want to allocate to this band */
981       if (i != start)
982          balance -= tell;
983       remaining_bits = (total_bits<<BITRES)-tell-1;
984       if (i <= codedBands-1)
985       {
986          curr_balance = (codedBands-i);
987          if (curr_balance > 3)
988             curr_balance = 3;
989          curr_balance = balance / curr_balance;
990          b = IMIN(remaining_bits+1,pulses[i]+curr_balance);
991          if (b<0)
992             b = 0;
993       } else {
994          b = 0;
995       }
996       /* Prevents ridiculous bit depths */
997       if (b > C*16*N<<BITRES)
998          b = C*16*N<<BITRES;
999
1000       if (M*eBands[i]-N >= M*eBands[start] && (update_lowband || lowband_offset==-1))
1001             lowband_offset = M*eBands[i];
1002
1003       tf_change = tf_res[i];
1004       if (i>=m->effEBands)
1005       {
1006          X=norm;
1007          if (_Y!=NULL)
1008             Y = norm;
1009       }
1010
1011       /* This ensures we never repeat spectral content within one band */
1012       if (lowband_offset != -1)
1013       {
1014          effective_lowband = lowband_offset-N;
1015          if (effective_lowband < M*eBands[start])
1016             effective_lowband = M*eBands[start];
1017       }
1018       quant_band(encode, m, i, X, Y, N, b, fold, B, intensity, tf_change,
1019             effective_lowband != -1 ? norm+effective_lowband : NULL, resynth, ec, &remaining_bits, LM,
1020             norm+M*eBands[i], bandE, 0, &seed, Q15ONE, lowband_scratch);
1021
1022       balance += pulses[i] + tell;
1023
1024       /* Update the folding position only as long as we have 1 bit/sample depth */
1025       update_lowband = (b>>BITRES)>N;
1026    }
1027    RESTORE_STACK;
1028 }
1029