Moves the bit-side gain application to the quantizer
[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    for (c=0;c<C;c++)
74    {
75       for (i=0;i<end;i++)
76       {
77          int j;
78          celt_word32 maxval=0;
79          celt_word32 sum = 0;
80          
81          j=M*eBands[i]; do {
82             maxval = MAX32(maxval, X[j+c*N]);
83             maxval = MAX32(maxval, -X[j+c*N]);
84          } while (++j<M*eBands[i+1]);
85          
86          if (maxval > 0)
87          {
88             int shift = celt_ilog2(maxval)-10;
89             j=M*eBands[i]; do {
90                sum = MAC16_16(sum, EXTRACT16(VSHR32(X[j+c*N],shift)),
91                                    EXTRACT16(VSHR32(X[j+c*N],shift)));
92             } while (++j<M*eBands[i+1]);
93             /* We're adding one here to make damn sure we never end up with a pitch vector that's
94                larger than unity norm */
95             bank[i+c*m->nbEBands] = EPSILON+VSHR32(EXTEND32(celt_sqrt(sum)),-shift);
96          } else {
97             bank[i+c*m->nbEBands] = EPSILON;
98          }
99          /*printf ("%f ", bank[i+c*m->nbEBands]);*/
100       }
101    }
102    /*printf ("\n");*/
103 }
104
105 /* Normalise each band such that the energy is one. */
106 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)
107 {
108    int i, c, N;
109    const celt_int16 *eBands = m->eBands;
110    const int C = CHANNELS(_C);
111    N = M*m->shortMdctSize;
112    for (c=0;c<C;c++)
113    {
114       i=0; do {
115          celt_word16 g;
116          int j,shift;
117          celt_word16 E;
118          shift = celt_zlog2(bank[i+c*m->nbEBands])-13;
119          E = VSHR32(bank[i+c*m->nbEBands], shift);
120          g = EXTRACT16(celt_rcp(SHL32(E,3)));
121          j=M*eBands[i]; do {
122             X[j+c*N] = MULT16_16_Q15(VSHR32(freq[j+c*N],shift-1),g);
123          } while (++j<M*eBands[i+1]);
124       } while (++i<end);
125    }
126 }
127
128 #else /* FIXED_POINT */
129 /* Compute the amplitude (sqrt energy) in each of the bands */
130 void compute_band_energies(const CELTMode *m, const celt_sig *X, celt_ener *bank, int end, int _C, int M)
131 {
132    int i, c, N;
133    const celt_int16 *eBands = m->eBands;
134    const int C = CHANNELS(_C);
135    N = M*m->shortMdctSize;
136    for (c=0;c<C;c++)
137    {
138       for (i=0;i<end;i++)
139       {
140          int j;
141          celt_word32 sum = 1e-10f;
142          for (j=M*eBands[i];j<M*eBands[i+1];j++)
143             sum += X[j+c*N]*X[j+c*N];
144          bank[i+c*m->nbEBands] = celt_sqrt(sum);
145          /*printf ("%f ", bank[i+c*m->nbEBands]);*/
146       }
147    }
148    /*printf ("\n");*/
149 }
150
151 /* Normalise each band such that the energy is one. */
152 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)
153 {
154    int i, c, N;
155    const celt_int16 *eBands = m->eBands;
156    const int C = CHANNELS(_C);
157    N = M*m->shortMdctSize;
158    for (c=0;c<C;c++)
159    {
160       for (i=0;i<end;i++)
161       {
162          int j;
163          celt_word16 g = 1.f/(1e-10f+bank[i+c*m->nbEBands]);
164          for (j=M*eBands[i];j<M*eBands[i+1];j++)
165             X[j+c*N] = freq[j+c*N]*g;
166       }
167    }
168 }
169
170 #endif /* FIXED_POINT */
171
172 void renormalise_bands(const CELTMode *m, celt_norm * restrict X, int end, int _C, int M)
173 {
174    int i, c;
175    const celt_int16 *eBands = m->eBands;
176    const int C = CHANNELS(_C);
177    for (c=0;c<C;c++)
178    {
179       i=0; do {
180          renormalise_vector(X+M*eBands[i]+c*M*m->shortMdctSize, M*eBands[i+1]-M*eBands[i], Q15ONE);
181       } while (++i<end);
182    }
183 }
184
185 /* De-normalise the energy to produce the synthesis from the unit-energy bands */
186 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)
187 {
188    int i, c, N;
189    const celt_int16 *eBands = m->eBands;
190    const int C = CHANNELS(_C);
191    N = M*m->shortMdctSize;
192    celt_assert2(C<=2, "denormalise_bands() not implemented for >2 channels");
193    for (c=0;c<C;c++)
194    {
195       celt_sig * restrict f;
196       const celt_norm * restrict x;
197       f = freq+c*N;
198       x = X+c*N;
199       for (i=0;i<end;i++)
200       {
201          int j, band_end;
202          celt_word32 g = SHR32(bank[i+c*m->nbEBands],1);
203          j=M*eBands[i];
204          band_end = M*eBands[i+1];
205          do {
206             *f++ = SHL32(MULT16_32_Q15(*x, g),2);
207             x++;
208          } while (++j<band_end);
209       }
210       for (i=M*eBands[m->nbEBands];i<N;i++)
211          *f++ = 0;
212    }
213 }
214
215 static void stereo_band_mix(const CELTMode *m, celt_norm *X, celt_norm *Y, const celt_ener *bank, int stereo_mode, int bandID, int dir, int N)
216 {
217    int i = bandID;
218    int j;
219    celt_word16 a1, a2;
220    if (stereo_mode==0)
221    {
222       /* Do mid-side when not doing intensity stereo */
223       a1 = QCONST16(.70711f,14);
224       a2 = dir*QCONST16(.70711f,14);
225    } else {
226       celt_word16 left, right;
227       celt_word16 norm;
228 #ifdef FIXED_POINT
229       int shift = celt_zlog2(MAX32(bank[i], bank[i+m->nbEBands]))-13;
230 #endif
231       left = VSHR32(bank[i],shift);
232       right = VSHR32(bank[i+m->nbEBands],shift);
233       norm = EPSILON + celt_sqrt(EPSILON+MULT16_16(left,left)+MULT16_16(right,right));
234       a1 = DIV32_16(SHL32(EXTEND32(left),14),norm);
235       a2 = dir*DIV32_16(SHL32(EXTEND32(right),14),norm);
236    }
237    for (j=0;j<N;j++)
238    {
239       celt_norm r, l;
240       l = X[j];
241       r = Y[j];
242       X[j] = MULT16_16_Q14(a1,l) + MULT16_16_Q14(a2,r);
243       Y[j] = MULT16_16_Q14(a1,r) - MULT16_16_Q14(a2,l);
244    }
245 }
246
247 /* Decide whether we should spread the pulses in the current frame */
248 int folding_decision(const CELTMode *m, celt_norm *X, int *average, int *last_decision, int end, int _C, int M)
249 {
250    int i, c, N0;
251    int sum = 0, nbBands=0;
252    const int C = CHANNELS(_C);
253    const celt_int16 * restrict eBands = m->eBands;
254    int decision;
255    
256    N0 = M*m->shortMdctSize;
257
258    if (M*(eBands[end]-eBands[end-1]) <= 8)
259       return 0;
260    for (c=0;c<C;c++)
261    {
262       for (i=0;i<end;i++)
263       {
264          int j, N, tmp=0;
265          int tcount[3] = {0};
266          celt_norm * restrict x = X+M*eBands[i]+c*N0;
267          N = M*(eBands[i+1]-eBands[i]);
268          if (N<=8)
269             continue;
270          /* Compute rough CDF of |x[j]| */
271          for (j=0;j<N;j++)
272          {
273             celt_word32 x2N; /* Q13 */
274
275             x2N = MULT16_16(MULT16_16_Q15(x[j], x[j]), N);
276             if (x2N < QCONST16(0.25f,13))
277                tcount[0]++;
278             if (x2N < QCONST16(0.0625f,13))
279                tcount[1]++;
280             if (x2N < QCONST16(0.015625f,13))
281                tcount[2]++;
282          }
283
284          tmp = (2*tcount[2] >= N) + (2*tcount[1] >= N) + (2*tcount[0] >= N);
285          sum += tmp*256;
286          nbBands++;
287       }
288    }
289    sum /= nbBands;
290    /* Recursive averaging */
291    sum = (sum+*average)>>1;
292    *average = sum;
293    /* Hysteresis */
294    sum = (3*sum + ((*last_decision<<7) + 64) + 2)>>2;
295    /* decision and last_decision do not use the same encoding */
296    if (sum < 128)
297    {
298       decision = 2;
299       *last_decision = 0;
300    } else if (sum < 256)
301    {
302       decision = 1;
303       *last_decision = 1;
304    } else if (sum < 384)
305    {
306       decision = 3;
307       *last_decision = 2;
308    } else {
309       decision = 0;
310       *last_decision = 3;
311    }
312    return decision;
313 }
314
315 #ifdef MEASURE_NORM_MSE
316
317 float MSE[30] = {0};
318 int nbMSEBands = 0;
319 int MSECount[30] = {0};
320
321 void dump_norm_mse(void)
322 {
323    int i;
324    for (i=0;i<nbMSEBands;i++)
325    {
326       printf ("%g ", MSE[i]/MSECount[i]);
327    }
328    printf ("\n");
329 }
330
331 void measure_norm_mse(const CELTMode *m, float *X, float *X0, float *bandE, float *bandE0, int M, int N, int C)
332 {
333    static int init = 0;
334    int i;
335    if (!init)
336    {
337       atexit(dump_norm_mse);
338       init = 1;
339    }
340    for (i=0;i<m->nbEBands;i++)
341    {
342       int j;
343       int c;
344       float g;
345       if (bandE0[i]<10 || (C==2 && bandE0[i+m->nbEBands]<1))
346          continue;
347       for (c=0;c<C;c++)
348       {
349          g = bandE[i+c*m->nbEBands]/(1e-15+bandE0[i+c*m->nbEBands]);
350          for (j=M*m->eBands[i];j<M*m->eBands[i+1];j++)
351             MSE[i] += (g*X[j+c*N]-X0[j+c*N])*(g*X[j+c*N]-X0[j+c*N]);
352       }
353       MSECount[i]+=C;
354    }
355    nbMSEBands = m->nbEBands;
356 }
357
358 #endif
359
360 static void interleave_vector(celt_norm *X, int N0, int stride)
361 {
362    int i,j;
363    VARDECL(celt_norm, tmp);
364    int N;
365    SAVE_STACK;
366    N = N0*stride;
367    ALLOC(tmp, N, celt_norm);
368    for (i=0;i<stride;i++)
369       for (j=0;j<N0;j++)
370          tmp[j*stride+i] = X[i*N0+j];
371    for (j=0;j<N;j++)
372       X[j] = tmp[j];
373    RESTORE_STACK;
374 }
375
376 static void deinterleave_vector(celt_norm *X, int N0, int stride)
377 {
378    int i,j;
379    VARDECL(celt_norm, tmp);
380    int N;
381    SAVE_STACK;
382    N = N0*stride;
383    ALLOC(tmp, N, celt_norm);
384    for (i=0;i<stride;i++)
385       for (j=0;j<N0;j++)
386          tmp[i*N0+j] = X[j*stride+i];
387    for (j=0;j<N;j++)
388       X[j] = tmp[j];
389    RESTORE_STACK;
390 }
391
392 static void haar1(celt_norm *X, int N0, int stride)
393 {
394    int i, j;
395    N0 >>= 1;
396    for (i=0;i<stride;i++)
397       for (j=0;j<N0;j++)
398       {
399          celt_norm tmp = X[stride*2*j+i];
400          X[stride*2*j+i] = MULT16_16_Q15(QCONST16(.7070678f,15), X[stride*2*j+i] + X[stride*(2*j+1)+i]);
401          X[stride*(2*j+1)+i] = MULT16_16_Q15(QCONST16(.7070678f,15), tmp - X[stride*(2*j+1)+i]);
402       }
403 }
404
405 static int compute_qn(int N, int b, int offset, int stereo)
406 {
407    static const celt_int16 exp2_table8[8] =
408       {16384, 17866, 19483, 21247, 23170, 25267, 27554, 30048};
409    int qn, qb;
410    int N2 = 2*N-1;
411    if (stereo && N==2)
412       N2--;
413    qb = (b+N2*offset)/N2;
414    if (qb > (b>>1)-(1<<BITRES))
415       qb = (b>>1)-(1<<BITRES);
416
417    if (qb<0)
418        qb = 0;
419    if (qb>14<<BITRES)
420      qb = 14<<BITRES;
421
422    if (qb<(1<<BITRES>>1)) {
423       qn = 1;
424    } else {
425       qn = exp2_table8[qb&0x7]>>(14-(qb>>BITRES));
426       qn = (qn+1)>>1<<1;
427       if (qn>1024)
428          qn = 1024;
429    }
430    return qn;
431 }
432
433
434 /* This function is responsible for encoding and decoding a band for both
435    the mono and stereo case. Even in the mono case, it can split the band
436    in two and transmit the energy difference with the two half-bands. It
437    can be called recursively so bands can end up being split in 8 parts. */
438 static void quant_band(int encode, const CELTMode *m, int i, celt_norm *X, celt_norm *Y,
439       int N, int b, int spread, int B, int tf_change, celt_norm *lowband, int resynth, void *ec,
440       celt_int32 *remaining_bits, int LM, celt_norm *lowband_out, const celt_ener *bandE, int level,
441       celt_int32 *seed, celt_word16 gain)
442 {
443    int q;
444    int curr_bits;
445    int stereo, split;
446    int imid=0, iside=0;
447    int N0=N;
448    int N_B=N;
449    int N_B0;
450    int B0=B;
451    int time_divide=0;
452    int recombine=0;
453
454    N_B /= B;
455    N_B0 = N_B;
456
457    split = stereo = Y != NULL;
458
459    /* Special case for one sample */
460    if (N==1)
461    {
462       int c;
463       celt_norm *x = X;
464       for (c=0;c<1+stereo;c++)
465       {
466          int sign=0;
467          if (b>=1<<BITRES && *remaining_bits>=1<<BITRES)
468          {
469             if (encode)
470             {
471                sign = x[0]<0;
472                ec_enc_bits((ec_enc*)ec, sign, 1);
473             } else {
474                sign = ec_dec_bits((ec_dec*)ec, 1);
475             }
476             *remaining_bits -= 1<<BITRES;
477             b-=1<<BITRES;
478          }
479          if (resynth)
480             x[0] = sign ? -NORM_SCALING : NORM_SCALING;
481          x = Y;
482       }
483       if (lowband_out)
484          lowband_out[0] = X[0];
485       return;
486    }
487
488    /* Band recombining to increase frequency resolution */
489    if (!stereo && B > 1 && level == 0 && tf_change>0)
490    {
491       while (B>1 && tf_change>0)
492       {
493          B>>=1;
494          N_B<<=1;
495          if (encode)
496             haar1(X, N_B, B);
497          if (lowband)
498             haar1(lowband, N_B, B);
499          recombine++;
500          tf_change--;
501       }
502       B0=B;
503       N_B0 = N_B;
504    }
505
506    /* Increasing the time resolution */
507    if (!stereo && level==0)
508    {
509       while ((N_B&1) == 0 && tf_change<0 && B <= (1<<LM))
510       {
511          if (encode)
512             haar1(X, N_B, B);
513          if (lowband)
514             haar1(lowband, N_B, B);
515          B <<= 1;
516          N_B >>= 1;
517          time_divide++;
518          tf_change++;
519       }
520       B0 = B;
521       N_B0 = N_B;
522    }
523
524    /* Reorganize the samples in time order instead of frequency order */
525    if (!stereo && B0>1 && level==0)
526    {
527       if (encode)
528          deinterleave_vector(X, N_B, B0);
529       if (lowband)
530          deinterleave_vector(lowband, N_B, B0);
531    }
532
533    /* If we need more than 32 bits, try splitting the band in two. */
534    if (!stereo && LM != -1 && b > 32<<BITRES && N>2)
535    {
536       if (LM>0 || (N&1)==0)
537       {
538          N >>= 1;
539          Y = X+N;
540          split = 1;
541          LM -= 1;
542          B = (B+1)>>1;
543       }
544    }
545
546    if (split)
547    {
548       int qn;
549       int itheta=0;
550       int mbits, sbits, delta;
551       int qalloc;
552       celt_word16 mid, side;
553       int offset;
554
555       /* Decide on the resolution to give to the split parameter theta */
556       offset = ((m->logN[i]+(LM<<BITRES))>>1) - (stereo ? QTHETA_OFFSET_STEREO : QTHETA_OFFSET);
557       qn = compute_qn(N, b, offset, stereo);
558
559       qalloc = 0;
560       if (qn!=1)
561       {
562          if (encode)
563          {
564             if (stereo)
565                stereo_band_mix(m, X, Y, bandE, 0, i, 1, N);
566
567             mid = vector_norm(X, N);
568             side = vector_norm(Y, N);
569             /* TODO: Renormalising X and Y *may* help fixed-point a bit at very high rate.
570                      Let's do that at higher complexity */
571             /*mid = renormalise_vector(X, Q15ONE, N, 1);
572             side = renormalise_vector(Y, Q15ONE, N, 1);*/
573
574             /* theta is the atan() of the ratio between the (normalized)
575                side and mid. With just that parameter, we can re-scale both
576                mid and side because we know that 1) they have unit norm and
577                2) they are orthogonal. */
578    #ifdef FIXED_POINT
579             /* 0.63662 = 2/pi */
580             itheta = MULT16_16_Q15(QCONST16(0.63662f,15),celt_atan2p(side, mid));
581    #else
582             itheta = (int)floor(.5f+16384*0.63662f*atan2(side,mid));
583    #endif
584
585             itheta = (itheta*qn+8192)>>14;
586          }
587
588          /* Entropy coding of the angle. We use a uniform pdf for the
589             first stereo split but a triangular one for the rest. */
590          if (stereo || qn>256 || B>1)
591          {
592             if (encode)
593                ec_enc_uint((ec_enc*)ec, itheta, qn+1);
594             else
595                itheta = ec_dec_uint((ec_dec*)ec, qn+1);
596             qalloc = log2_frac(qn+1,BITRES);
597          } else {
598             int fs=1, ft;
599             ft = ((qn>>1)+1)*((qn>>1)+1);
600             if (encode)
601             {
602                int fl;
603
604                fs = itheta <= (qn>>1) ? itheta + 1 : qn + 1 - itheta;
605                fl = itheta <= (qn>>1) ? itheta*(itheta + 1)>>1 :
606                 ft - ((qn + 1 - itheta)*(qn + 2 - itheta)>>1);
607
608                ec_encode((ec_enc*)ec, fl, fl+fs, ft);
609             } else {
610                int fl=0;
611                int fm;
612                fm = ec_decode((ec_dec*)ec, ft);
613
614                if (fm < ((qn>>1)*((qn>>1) + 1)>>1))
615                {
616                   itheta = (isqrt32(8*(celt_uint32)fm + 1) - 1)>>1;
617                   fs = itheta + 1;
618                   fl = itheta*(itheta + 1)>>1;
619                }
620                else
621                {
622                   itheta = (2*(qn + 1)
623                    - isqrt32(8*(celt_uint32)(ft - fm - 1) + 1))>>1;
624                   fs = qn + 1 - itheta;
625                   fl = ft - ((qn + 1 - itheta)*(qn + 2 - itheta)>>1);
626                }
627
628                ec_dec_update((ec_dec*)ec, fl, fl+fs, ft);
629             }
630             qalloc = log2_frac(ft,BITRES) - log2_frac(fs,BITRES) + 1;
631          }
632          itheta = (celt_int32)itheta*16384/qn;
633       } else {
634          if (stereo && encode)
635             stereo_band_mix(m, X, Y, bandE, 1, i, 1, N);
636       }
637
638       if (itheta == 0)
639       {
640          imid = 32767;
641          iside = 0;
642          delta = -10000;
643       } else if (itheta == 16384)
644       {
645          imid = 0;
646          iside = 32767;
647          delta = 10000;
648       } else {
649          imid = bitexact_cos(itheta);
650          iside = bitexact_cos(16384-itheta);
651          /* This is the mid vs side allocation that minimizes squared error
652             in that band. */
653          delta = (N-1)*(log2_frac(iside,BITRES+2)-log2_frac(imid,BITRES+2))>>2;
654       }
655
656       /* This is a special case for N=2 that only works for stereo and takes
657          advantage of the fact that mid and side are orthogonal to encode
658          the side with just one bit. */
659       if (N==2 && stereo)
660       {
661          int c;
662          int sign=1;
663          celt_norm *x2, *y2;
664          mbits = b-qalloc;
665          sbits = 0;
666          /* Only need one bit for the side */
667          if (itheta != 0 && itheta != 16384)
668             sbits = 1<<BITRES;
669          mbits -= sbits;
670          c = itheta > 8192;
671          *remaining_bits -= qalloc+sbits;
672
673          x2 = c ? Y : X;
674          y2 = c ? X : Y;
675          if (sbits)
676          {
677             if (encode)
678             {
679                /* Here we only need to encode a sign for the side */
680                sign = x2[0]*y2[1] - x2[1]*y2[0] > 0;
681                ec_enc_bits((ec_enc*)ec, sign, 1);
682             } else {
683                sign = ec_dec_bits((ec_dec*)ec, 1);
684             }
685          }
686          sign = 2*sign - 1;
687          quant_band(encode, m, i, x2, NULL, N, mbits, spread, B, tf_change, lowband, resynth, ec, remaining_bits, LM, lowband_out, NULL, level+1, seed, gain);
688          y2[0] = -sign*x2[1];
689          y2[1] = sign*x2[0];
690       } else {
691          /* "Normal" split code */
692          celt_norm *next_lowband2=NULL;
693          celt_norm *next_lowband_out1=NULL;
694          int next_level=0;
695
696 #ifdef FIXED_POINT
697          mid = imid;
698          side = iside;
699 #else
700          mid = (1.f/32768)*imid;
701          side = (1.f/32768)*iside;
702 #endif
703
704          /* Give more bits to low-energy MDCTs than they would otherwise deserve */
705          if (B>1 && !stereo)
706             delta >>= 1;
707
708          mbits = (b-qalloc-delta)/2;
709          if (mbits > b-qalloc)
710             mbits = b-qalloc;
711          if (mbits<0)
712             mbits=0;
713          sbits = b-qalloc-mbits;
714          *remaining_bits -= qalloc;
715
716          if (lowband && !stereo)
717             next_lowband2 = lowband+N; /* >32-bit split case */
718
719          /* Only stereo needs to pass on lowband_out. Otherwise, it's
720             handled at the end */
721          if (stereo)
722             next_lowband_out1 = lowband_out;
723          else
724             next_level = level+1;
725
726          quant_band(encode, m, i, X, NULL, N, mbits, spread, B, tf_change,
727                lowband, resynth, ec, remaining_bits, LM, next_lowband_out1,
728                NULL, next_level, seed, MULT16_16_P15(gain,mid));
729          quant_band(encode, m, i, Y, NULL, N, sbits, spread, B, tf_change,
730                next_lowband2, resynth, ec, remaining_bits, LM, NULL,
731                NULL, next_level, seed, MULT16_16_P15(gain,side));
732       }
733
734    } else {
735       /* This is the basic no-split case */
736       q = bits2pulses(m, i, LM, b);
737       curr_bits = pulses2bits(m, i, LM, q);
738       *remaining_bits -= curr_bits;
739
740       /* Ensures we can never bust the budget */
741       while (*remaining_bits < 0 && q > 0)
742       {
743          *remaining_bits += curr_bits;
744          q--;
745          curr_bits = pulses2bits(m, i, LM, q);
746          *remaining_bits -= curr_bits;
747       }
748
749       /* Finally do the actual quantization */
750       if (encode)
751          alg_quant(X, N, q, spread, B, lowband, resynth, (ec_enc*)ec, seed, gain);
752       else
753          alg_unquant(X, N, q, spread, B, lowband, (ec_dec*)ec, seed, gain);
754    }
755
756    /* This code is used by the decoder and by the resynthesis-enabled encoder */
757    if (resynth)
758    {
759       int k;
760
761       /* Undo the sample reorganization going from time order to frequency order */
762       if (!stereo && B0>1 && level==0)
763       {
764          interleave_vector(X, N_B, B0);
765          if (lowband)
766             interleave_vector(lowband, N_B, B0);
767       }
768
769       /* Undo time-freq changes that we did earlier */
770       N_B = N_B0;
771       B = B0;
772       for (k=0;k<time_divide;k++)
773       {
774          B >>= 1;
775          N_B <<= 1;
776          haar1(X, N_B, B);
777          if (lowband)
778             haar1(lowband, N_B, B);
779       }
780
781       for (k=0;k<recombine;k++)
782       {
783          haar1(X, N_B, B);
784          if (lowband)
785             haar1(lowband, N_B, B);
786          N_B>>=1;
787          B <<= 1;
788       }
789
790       /* Scale output for later folding */
791       if (lowband_out && !stereo)
792       {
793          int j;
794          celt_word16 n;
795          n = celt_sqrt(SHL32(EXTEND32(N0),22));
796          for (j=0;j<N0;j++)
797             lowband_out[j] = MULT16_16_Q15(n,X[j]);
798       }
799
800       if (stereo)
801       {
802          stereo_band_mix(m, X, Y, bandE, 0, i, -1, N);
803          /* We only need to renormalize because quantization may not
804             have preserved orthogonality of mid and side */
805          renormalise_vector(X, N, Q15ONE);
806          renormalise_vector(Y, N, Q15ONE);
807       }
808    }
809 }
810
811 void quant_all_bands(int encode, const CELTMode *m, int start, int end, celt_norm *_X, celt_norm *_Y, const celt_ener *bandE, int *pulses, int shortBlocks, int fold, int *tf_res, int resynth, int total_bits, void *ec, int LM)
812 {
813    int i, balance;
814    celt_int32 remaining_bits;
815    const celt_int16 * restrict eBands = m->eBands;
816    celt_norm * restrict norm;
817    VARDECL(celt_norm, _norm);
818    int B;
819    int M;
820    celt_int32 seed;
821    celt_norm *lowband;
822    int update_lowband = 1;
823    int C = _Y != NULL ? 2 : 1;
824    SAVE_STACK;
825
826    M = 1<<LM;
827    B = shortBlocks ? M : 1;
828    ALLOC(_norm, M*eBands[m->nbEBands], celt_norm);
829    norm = _norm;
830
831    if (encode)
832       seed = ((ec_enc*)ec)->rng;
833    else
834       seed = ((ec_dec*)ec)->rng;
835    balance = 0;
836    lowband = NULL;
837    for (i=start;i<end;i++)
838    {
839       int tell;
840       int b;
841       int N;
842       int curr_balance;
843       celt_norm * restrict X, * restrict Y;
844       int tf_change=0;
845       celt_norm *effective_lowband;
846       
847       X = _X+M*eBands[i];
848       if (_Y!=NULL)
849          Y = _Y+M*eBands[i];
850       else
851          Y = NULL;
852       N = M*eBands[i+1]-M*eBands[i];
853       if (encode)
854          tell = ec_enc_tell((ec_enc*)ec, BITRES);
855       else
856          tell = ec_dec_tell((ec_dec*)ec, BITRES);
857
858       /* Compute how many bits we want to allocate to this band */
859       if (i != start)
860          balance -= tell;
861       remaining_bits = (total_bits<<BITRES)-tell-1;
862       curr_balance = (end-i);
863       if (curr_balance > 3)
864          curr_balance = 3;
865       curr_balance = balance / curr_balance;
866       b = IMIN(remaining_bits+1,pulses[i]+curr_balance);
867       if (b<0)
868          b = 0;
869       /* Prevents ridiculous bit depths */
870       if (b > C*16*N<<BITRES)
871          b = C*16*N<<BITRES;
872
873       if (M*eBands[i]-N >= M*eBands[start] && (update_lowband || lowband==NULL))
874             lowband = norm+M*eBands[i]-N;
875
876       tf_change = tf_res[i];
877       if (i>=m->effEBands)
878       {
879          X=norm;
880          if (_Y!=NULL)
881             Y = norm;
882       }
883
884       if (tf_change==0 && !shortBlocks && fold)
885          effective_lowband = NULL;
886       else
887          effective_lowband = lowband;
888       quant_band(encode, m, i, X, Y, N, b, fold, B, tf_change,
889             effective_lowband, resynth, ec, &remaining_bits, LM,
890             norm+M*eBands[i], bandE, 0, &seed, Q15ONE);
891
892       balance += pulses[i] + tell;
893
894       /* Update the folding position only as long as we have 2 bit/sample depth */
895       update_lowband = (b>>BITRES)>2*N;
896    }
897    RESTORE_STACK;
898 }
899