Simplifications to the N=2 stereo split.
[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-10;
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] = 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, Q15ONE, M*eBands[i+1]-M*eBands[i], 1);
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, celt_word16 *average, int *last_decision, int end, int _C, int M)
249 {
250    int i, c, N0;
251    int NR=0;
252    celt_word32 ratio = EPSILON;
253    const int C = CHANNELS(_C);
254    const celt_int16 * restrict eBands = m->eBands;
255    
256    N0 = M*m->shortMdctSize;
257
258    for (c=0;c<C;c++)
259    {
260    for (i=0;i<end;i++)
261    {
262       int j, N;
263       int max_i=0;
264       celt_word16 max_val=EPSILON;
265       celt_word32 floor_ener=EPSILON;
266       celt_norm * restrict x = X+M*eBands[i]+c*N0;
267       N = M*eBands[i+1]-M*eBands[i];
268       for (j=0;j<N;j++)
269       {
270          if (ABS16(x[j])>max_val)
271          {
272             max_val = ABS16(x[j]);
273             max_i = j;
274          }
275       }
276 #if 0
277       for (j=0;j<N;j++)
278       {
279          if (abs(j-max_i)>2)
280             floor_ener += x[j]*x[j];
281       }
282 #else
283       floor_ener = QCONST32(1.,28)-MULT16_16(max_val,max_val);
284       if (max_i < N-1)
285          floor_ener -= MULT16_16(x[(max_i+1)], x[(max_i+1)]);
286       if (max_i < N-2)
287          floor_ener -= MULT16_16(x[(max_i+2)], x[(max_i+2)]);
288       if (max_i > 0)
289          floor_ener -= MULT16_16(x[(max_i-1)], x[(max_i-1)]);
290       if (max_i > 1)
291          floor_ener -= MULT16_16(x[(max_i-2)], x[(max_i-2)]);
292       floor_ener = MAX32(floor_ener, EPSILON);
293 #endif
294       if (N>7)
295       {
296          celt_word16 r;
297          celt_word16 den = celt_sqrt(floor_ener);
298          den = MAX32(QCONST16(.02f, 15), den);
299          r = DIV32_16(SHL32(EXTEND32(max_val),8),den);
300          ratio = ADD32(ratio, EXTEND32(r));
301          NR++;
302       }
303    }
304    }
305    if (NR>0)
306       ratio = DIV32_16(ratio, NR);
307    ratio = ADD32(HALF32(ratio), HALF32(*average));
308    if (!*last_decision)
309    {
310       *last_decision = (ratio < QCONST16(1.8f,8));
311    } else {
312       *last_decision = (ratio < QCONST16(3.f,8));
313    }
314    *average = EXTRACT16(ratio);
315    return *last_decision;
316 }
317
318 #ifdef MEASURE_NORM_MSE
319
320 float MSE[30] = {0};
321 int nbMSEBands = 0;
322 int MSECount[30] = {0};
323
324 void dump_norm_mse(void)
325 {
326    int i;
327    for (i=0;i<nbMSEBands;i++)
328    {
329       printf ("%f ", MSE[i]/MSECount[i]);
330    }
331    printf ("\n");
332 }
333
334 void measure_norm_mse(const CELTMode *m, float *X, float *X0, float *bandE, float *bandE0, int M, int N, int C)
335 {
336    static int init = 0;
337    int i;
338    if (!init)
339    {
340       atexit(dump_norm_mse);
341       init = 1;
342    }
343    for (i=0;i<m->nbEBands;i++)
344    {
345       int j;
346       int c;
347       float g;
348       if (bandE0[i]<1 || (C==2 && bandE0[i+m->nbEBands]<1))
349          continue;
350       for (c=0;c<C;c++)
351       {
352          g = bandE[i+c*m->nbEBands]/(1e-15+bandE0[i+c*m->nbEBands]);
353          for (j=M*m->eBands[i];j<M*m->eBands[i+1];j++)
354             MSE[i] += (g*X[j+c*N]-X0[j+c*N])*(g*X[j+c*N]-X0[j+c*N]);
355       }
356       MSECount[i]+=C;
357    }
358    nbMSEBands = m->nbEBands;
359 }
360
361 #endif
362
363 static void interleave_vector(celt_norm *X, int N0, int stride)
364 {
365    int i,j;
366    VARDECL(celt_norm, tmp);
367    int N;
368    SAVE_STACK;
369    N = N0*stride;
370    ALLOC(tmp, N, celt_norm);
371    for (i=0;i<stride;i++)
372       for (j=0;j<N0;j++)
373          tmp[j*stride+i] = X[i*N0+j];
374    for (j=0;j<N;j++)
375       X[j] = tmp[j];
376    RESTORE_STACK;
377 }
378
379 static void deinterleave_vector(celt_norm *X, int N0, int stride)
380 {
381    int i,j;
382    VARDECL(celt_norm, tmp);
383    int N;
384    SAVE_STACK;
385    N = N0*stride;
386    ALLOC(tmp, N, celt_norm);
387    for (i=0;i<stride;i++)
388       for (j=0;j<N0;j++)
389          tmp[i*N0+j] = X[j*stride+i];
390    for (j=0;j<N;j++)
391       X[j] = tmp[j];
392    RESTORE_STACK;
393 }
394
395 static void haar1(celt_norm *X, int N0, int stride)
396 {
397    int i, j;
398    N0 >>= 1;
399    for (i=0;i<stride;i++)
400       for (j=0;j<N0;j++)
401       {
402          celt_norm tmp = X[stride*2*j+i];
403          X[stride*2*j+i] = MULT16_16_Q15(QCONST16(.7070678f,15), X[stride*2*j+i] + X[stride*(2*j+1)+i]);
404          X[stride*(2*j+1)+i] = MULT16_16_Q15(QCONST16(.7070678f,15), tmp - X[stride*(2*j+1)+i]);
405       }
406 }
407
408 static int compute_qn(int N, int b, int offset, int stereo)
409 {
410    static const celt_int16 exp2_table8[8] =
411       {16384, 17866, 19483, 21247, 23170, 25267, 27554, 30048};
412    int qn, qb;
413    int N2 = 2*N-1;
414    if (stereo && N==2)
415       N2--;
416    qb = (b+N2*offset)/N2;
417    if (qb > (b>>1)-(1<<BITRES))
418       qb = (b>>1)-(1<<BITRES);
419
420    if (qb<0)
421        qb = 0;
422    if (qb>14<<BITRES)
423      qb = 14<<BITRES;
424
425    if (qb<(1<<BITRES>>1)) {
426       qn = 1;
427    } else {
428       qn = exp2_table8[qb&0x7]>>(14-(qb>>BITRES));
429       qn = (qn+1)>>1<<1;
430       if (qn>1024)
431          qn = 1024;
432    }
433    return qn;
434 }
435
436
437 /* This function is responsible for encoding and decoding a band for both
438    the mono and stereo case. Even in the mono case, it can split the band
439    in two and transmit the energy difference with the two half-bands. It
440    can be called recursively so bands can end up being split in 8 parts. */
441 static void quant_band(int encode, const CELTMode *m, int i, celt_norm *X, celt_norm *Y,
442       int N, int b, int spread, int B, int tf_change, celt_norm *lowband, int resynth, void *ec,
443       celt_int32 *remaining_bits, int LM, celt_norm *lowband_out, const celt_ener *bandE, int level, celt_int32 *seed)
444 {
445    int q;
446    int curr_bits;
447    int stereo, split;
448    int imid=0, iside=0;
449    int N0=N;
450    int N_B=N;
451    int N_B0;
452    int B0=B;
453    int time_divide=0;
454    int recombine=0;
455
456    N_B /= B;
457    N_B0 = N_B;
458
459    split = stereo = Y != NULL;
460
461    /* Special case for one sample */
462    if (N==1)
463    {
464       int c;
465       celt_norm *x = X;
466       for (c=0;c<1+stereo;c++)
467       {
468          int sign=0;
469          if (b>=1<<BITRES && *remaining_bits>=1<<BITRES)
470          {
471             if (encode)
472             {
473                sign = x[0]<0;
474                ec_enc_bits((ec_enc*)ec, sign, 1);
475             } else {
476                sign = ec_dec_bits((ec_dec*)ec, 1);
477             }
478             *remaining_bits -= 1<<BITRES;
479             b-=1<<BITRES;
480          }
481          if (resynth)
482             x[0] = sign ? -NORM_SCALING : NORM_SCALING;
483          x = Y;
484       }
485       if (lowband_out)
486          lowband_out[0] = X[0];
487       return;
488    }
489
490    /* Band recombining to increase frequency resolution */
491    if (!stereo && B > 1 && level == 0 && tf_change>0)
492    {
493       while (B>1 && tf_change>0)
494       {
495          B>>=1;
496          N_B<<=1;
497          if (encode)
498             haar1(X, N_B, B);
499          if (lowband)
500             haar1(lowband, N_B, B);
501          recombine++;
502          tf_change--;
503       }
504       B0=B;
505       N_B0 = N_B;
506    }
507
508    /* Increasing the time resolution */
509    if (!stereo && level==0)
510    {
511       while ((N_B&1) == 0 && tf_change<0 && B <= (1<<LM))
512       {
513          if (encode)
514             haar1(X, N_B, B);
515          if (lowband)
516             haar1(lowband, N_B, B);
517          B <<= 1;
518          N_B >>= 1;
519          time_divide++;
520          tf_change++;
521       }
522       B0 = B;
523       N_B0 = N_B;
524    }
525
526    /* Reorganize the samples in time order instead of frequency order */
527    if (!stereo && B0>1 && level==0)
528    {
529       if (encode)
530          deinterleave_vector(X, N_B, B0);
531       if (lowband)
532          deinterleave_vector(lowband, N_B, B0);
533    }
534
535    /* If we need more than 32 bits, try splitting the band in two. */
536    if (!stereo && LM != -1 && b > 32<<BITRES && N>2)
537    {
538       if (LM>0 || (N&1)==0)
539       {
540          N >>= 1;
541          Y = X+N;
542          split = 1;
543          LM -= 1;
544          B = (B+1)>>1;
545       }
546    }
547
548    if (split)
549    {
550       int qn;
551       int itheta=0;
552       int mbits, sbits, delta;
553       int qalloc;
554       celt_word16 mid, side;
555       int offset;
556
557       /* Decide on the resolution to give to the split parameter theta */
558       offset = ((m->logN[i]+(LM<<BITRES))>>1) - (stereo ? QTHETA_OFFSET_STEREO : QTHETA_OFFSET);
559       qn = compute_qn(N, b, offset, stereo);
560
561       qalloc = 0;
562       if (qn!=1)
563       {
564          if (encode)
565          {
566             if (stereo)
567                stereo_band_mix(m, X, Y, bandE, 0, i, 1, N);
568
569             mid = renormalise_vector(X, Q15ONE, N, 1);
570             side = renormalise_vector(Y, Q15ONE, N, 1);
571
572             /* theta is the atan() of the ratio between the (normalized)
573                side and mid. With just that parameter, we can re-scale both
574                mid and side because we know that 1) they have unit norm and
575                2) they are orthogonal. */
576    #ifdef FIXED_POINT
577             /* 0.63662 = 2/pi */
578             itheta = MULT16_16_Q15(QCONST16(0.63662f,15),celt_atan2p(side, mid));
579    #else
580             itheta = floor(.5f+16384*0.63662f*atan2(side,mid));
581    #endif
582
583             itheta = (itheta*qn+8192)>>14;
584          }
585
586          /* Entropy coding of the angle. We use a uniform pdf for the
587             first stereo split but a triangular one for the rest. */
588          if (stereo || qn>256 || B>1)
589          {
590             if (encode)
591                ec_enc_uint((ec_enc*)ec, itheta, qn+1);
592             else
593                itheta = ec_dec_uint((ec_dec*)ec, qn+1);
594             qalloc = log2_frac(qn+1,BITRES);
595          } else {
596             int fs=1, ft;
597             ft = ((qn>>1)+1)*((qn>>1)+1);
598             if (encode)
599             {
600                int fl;
601
602                fs = itheta <= (qn>>1) ? itheta + 1 : qn + 1 - itheta;
603                fl = itheta <= (qn>>1) ? itheta*(itheta + 1)>>1 :
604                 ft - ((qn + 1 - itheta)*(qn + 2 - itheta)>>1);
605
606                ec_encode((ec_enc*)ec, fl, fl+fs, ft);
607             } else {
608                int fl=0;
609                int fm;
610                fm = ec_decode((ec_dec*)ec, ft);
611
612                if (fm < ((qn>>1)*((qn>>1) + 1)>>1))
613                {
614                   itheta = (isqrt32(8*(celt_uint32)fm + 1) - 1)>>1;
615                   fs = itheta + 1;
616                   fl = itheta*(itheta + 1)>>1;
617                }
618                else
619                {
620                   itheta = (2*(qn + 1)
621                    - isqrt32(8*(celt_uint32)(ft - fm - 1) + 1))>>1;
622                   fs = qn + 1 - itheta;
623                   fl = ft - ((qn + 1 - itheta)*(qn + 2 - itheta)>>1);
624                }
625
626                ec_dec_update((ec_dec*)ec, fl, fl+fs, ft);
627             }
628             qalloc = log2_frac(ft,BITRES) - log2_frac(fs,BITRES) + 1;
629          }
630          itheta = itheta*16384/qn;
631       } else {
632          if (stereo && encode)
633             stereo_band_mix(m, X, Y, bandE, 1, i, 1, N);
634       }
635
636       if (itheta == 0)
637       {
638          imid = 32767;
639          iside = 0;
640          delta = -10000;
641       } else if (itheta == 16384)
642       {
643          imid = 0;
644          iside = 32767;
645          delta = 10000;
646       } else {
647          imid = bitexact_cos(itheta);
648          iside = bitexact_cos(16384-itheta);
649          /* This is the mid vs side allocation that minimizes squared error
650             in that band. */
651          delta = (N-1)*(log2_frac(iside,BITRES+2)-log2_frac(imid,BITRES+2))>>2;
652       }
653
654       /* This is a special case for N=2 that only works for stereo and takes
655          advantage of the fact that mid and side are orthogonal to encode
656          the side with just one bit. */
657       if (N==2 && stereo)
658       {
659          int c;
660          int sign=1;
661          celt_norm v[2], w[2];
662          celt_norm *x2, *y2;
663          mbits = b-qalloc;
664          sbits = 0;
665          /* Only need one bit for the side */
666          if (itheta != 0 && itheta != 16384)
667             sbits = 1<<BITRES;
668          mbits -= sbits;
669          c = itheta > 8192;
670          *remaining_bits -= qalloc+sbits;
671
672          x2 = c ? Y : X;
673          y2 = c ? X : Y;
674          if (encode)
675          {
676             /* v is the largest vector between mid and side. w is the other */
677             v[0] = x2[0];
678             v[1] = x2[1];
679             w[0] = y2[0];
680             w[1] = y2[1];
681             /* Here we only need to encode a sign for the side */
682             sign = v[0]*w[1] - v[1]*w[0] > 0;
683          }
684          quant_band(encode, m, i, v, NULL, N, mbits, spread, B, tf_change, lowband, resynth, ec, remaining_bits, LM, lowband_out, NULL, level+1, seed);
685          if (sbits)
686          {
687             if (encode)
688             {
689                ec_enc_bits((ec_enc*)ec, sign, 1);
690             } else {
691                sign = ec_dec_bits((ec_dec*)ec, 1);
692             }
693             sign = 2*sign - 1;
694          } else {
695             sign = 1;
696          }
697          w[0] = -sign*v[1];
698          w[1] = sign*v[0];
699          x2[0] = v[0];
700          x2[1] = v[1];
701          y2[0] = w[0];
702          y2[1] = w[1];
703       } else
704       {
705          /* "Normal" split code */
706          celt_norm *next_lowband2=NULL;
707          celt_norm *next_lowband_out1=NULL;
708          int next_level=0;
709
710          /* Give more bits to low-energy MDCTs than they would otherwise deserve */
711          if (B>1 && !stereo)
712             delta >>= 1;
713
714          mbits = (b-qalloc-delta)/2;
715          if (mbits > b-qalloc)
716             mbits = b-qalloc;
717          if (mbits<0)
718             mbits=0;
719          sbits = b-qalloc-mbits;
720          *remaining_bits -= qalloc;
721
722          if (lowband && !stereo)
723             next_lowband2 = lowband+N; /* >32-bit split case */
724
725          /* Only stereo needs to pass on lowband_out. Otherwise, it's handled at the end */
726          if (stereo)
727             next_lowband_out1 = lowband_out;
728          else
729             next_level = level+1;
730
731          quant_band(encode, m, i, X, NULL, N, mbits, spread, B, tf_change, lowband, resynth, ec, remaining_bits, LM, next_lowband_out1, NULL, next_level, seed);
732          quant_band(encode, m, i, Y, NULL, N, sbits, spread, B, tf_change, next_lowband2, resynth, ec, remaining_bits, LM, NULL, NULL, level, seed);
733       }
734
735    } else {
736       /* This is the basic no-split case */
737       q = bits2pulses(m, m->bits[LM][i], N, b);
738       curr_bits = pulses2bits(m->bits[LM][i], N, q);
739       *remaining_bits -= curr_bits;
740
741       /* Ensures we can never bust the budget */
742       while (*remaining_bits < 0 && q > 0)
743       {
744          *remaining_bits += curr_bits;
745          q--;
746          curr_bits = pulses2bits(m->bits[LM][i], N, q);
747          *remaining_bits -= curr_bits;
748       }
749
750       /* Finally do the actual quantization */
751       if (encode)
752          alg_quant(X, N, q, spread, B, lowband, resynth, (ec_enc*)ec, seed);
753       else
754          alg_unquant(X, N, q, spread, B, lowband, (ec_dec*)ec, seed);
755    }
756
757    /* This code is used by the decoder and by the resynthesis-enabled encoder */
758    if (resynth)
759    {
760       int k;
761
762       if (split)
763       {
764          int j;
765          celt_word16 mid, side;
766 #ifdef FIXED_POINT
767          mid = imid;
768          side = iside;
769 #else
770          mid = (1.f/32768)*imid;
771          side = (1.f/32768)*iside;
772 #endif
773          for (j=0;j<N;j++)
774             X[j] = MULT16_16_Q15(X[j], mid);
775          for (j=0;j<N;j++)
776             Y[j] = MULT16_16_Q15(Y[j], side);
777       }
778
779       /* Undo the sample reorganization going from time order to frequency order */
780       if (!stereo && B0>1 && level==0)
781       {
782          interleave_vector(X, N_B, B0);
783          if (lowband)
784             interleave_vector(lowband, N_B, B0);
785       }
786
787       /* Undo time-freq changes that we did earlier */
788       N_B = N_B0;
789       B = B0;
790       for (k=0;k<time_divide;k++)
791       {
792          B >>= 1;
793          N_B <<= 1;
794          haar1(X, N_B, B);
795          if (lowband)
796             haar1(lowband, N_B, B);
797       }
798
799       for (k=0;k<recombine;k++)
800       {
801          haar1(X, N_B, B);
802          if (lowband)
803             haar1(lowband, N_B, B);
804          N_B>>=1;
805          B <<= 1;
806       }
807
808       /* Scale output for later folding */
809       if (lowband_out && !stereo)
810       {
811          int j;
812          celt_word16 n;
813          n = celt_sqrt(SHL32(EXTEND32(N0),22));
814          for (j=0;j<N0;j++)
815             lowband_out[j] = MULT16_16_Q15(n,X[j]);
816       }
817
818       if (stereo)
819       {
820          stereo_band_mix(m, X, Y, bandE, 0, i, -1, N);
821          /* We only need to renormalize because quantization may not
822             have preserved orthogonality of mid and side */
823          renormalise_vector(X, Q15ONE, N, 1);
824          renormalise_vector(Y, Q15ONE, N, 1);
825       }
826    }
827 }
828
829 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)
830 {
831    int i, balance;
832    celt_int32 remaining_bits;
833    const celt_int16 * restrict eBands = m->eBands;
834    celt_norm * restrict norm;
835    VARDECL(celt_norm, _norm);
836    int B;
837    int M;
838    int spread;
839    celt_int32 seed;
840    celt_norm *lowband;
841    int update_lowband = 1;
842    int C = _Y != NULL ? 2 : 1;
843    SAVE_STACK;
844
845    M = 1<<LM;
846    B = shortBlocks ? M : 1;
847    spread = fold ? B : 0;
848    ALLOC(_norm, M*eBands[m->nbEBands], celt_norm);
849    norm = _norm;
850
851    if (encode)
852       seed = ((ec_enc*)ec)->rng;
853    else
854       seed = ((ec_dec*)ec)->rng;
855    balance = 0;
856    lowband = NULL;
857    for (i=start;i<end;i++)
858    {
859       int tell;
860       int b;
861       int N;
862       int curr_balance;
863       celt_norm * restrict X, * restrict Y;
864       int tf_change=0;
865       celt_norm *effective_lowband;
866       
867       X = _X+M*eBands[i];
868       if (_Y!=NULL)
869          Y = _Y+M*eBands[i];
870       else
871          Y = NULL;
872       N = M*eBands[i+1]-M*eBands[i];
873       if (encode)
874          tell = ec_enc_tell((ec_enc*)ec, BITRES);
875       else
876          tell = ec_dec_tell((ec_dec*)ec, BITRES);
877
878       /* Compute how many bits we want to allocate to this band */
879       if (i != start)
880          balance -= tell;
881       remaining_bits = (total_bits<<BITRES)-tell-1;
882       curr_balance = (end-i);
883       if (curr_balance > 3)
884          curr_balance = 3;
885       curr_balance = balance / curr_balance;
886       b = IMIN(remaining_bits+1,pulses[i]+curr_balance);
887       if (b<0)
888          b = 0;
889       /* Prevents ridiculous bit depths */
890       if (b > C*16*N<<BITRES)
891          b = C*16*N<<BITRES;
892
893       if (M*eBands[i]-N >= M*eBands[start] && (update_lowband || lowband==NULL))
894             lowband = norm+M*eBands[i]-N;
895
896       tf_change = tf_res[i];
897       if (i>=m->effEBands)
898       {
899          X=norm;
900          if (_Y!=NULL)
901             Y = norm;
902       }
903
904       if (tf_change==0 && !shortBlocks && fold)
905          effective_lowband = NULL;
906       else
907          effective_lowband = lowband;
908       quant_band(encode, m, i, X, Y, N, b, fold, B, tf_change, effective_lowband, resynth, ec, &remaining_bits, LM, norm+M*eBands[i], bandE, 0, &seed);
909
910       balance += pulses[i] + tell;
911
912       /* Update the folding position only as long as we have 2 bit/sample depth */
913       update_lowband = (b>>BITRES)>2*N;
914    }
915    RESTORE_STACK;
916 }
917