7050754fc6ef5e4c6b026ad9a96f73fef0c9384d
[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, 17867, 19484, 21247, 23170, 25268, 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, c2;
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 ? 1 : 0;
670          *remaining_bits -= qalloc+sbits;
671
672          x2 = X;
673          y2 = Y;
674          if (encode)
675          {
676             c2 = 1-c;
677
678             /* v is the largest vector between mid and side. w is the other */
679             if (c==0)
680             {
681                v[0] = x2[0];
682                v[1] = x2[1];
683                w[0] = y2[0];
684                w[1] = y2[1];
685             } else {
686                v[0] = y2[0];
687                v[1] = y2[1];
688                w[0] = x2[0];
689                w[1] = x2[1];
690             }
691             /* Here we only need to encode a sign for the side */
692             if (v[0]*w[1] - v[1]*w[0] > 0)
693                sign = 1;
694             else
695                sign = -1;
696          }
697          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);
698          if (sbits)
699          {
700             if (encode)
701             {
702                ec_enc_bits((ec_enc*)ec, sign==1, 1);
703             } else {
704                sign = 2*ec_dec_bits((ec_dec*)ec, 1)-1;
705             }
706          } else {
707             sign = 1;
708          }
709          w[0] = -sign*v[1];
710          w[1] = sign*v[0];
711          if (c==0)
712          {
713             x2[0] = v[0];
714             x2[1] = v[1];
715             y2[0] = w[0];
716             y2[1] = w[1];
717          } else {
718             x2[0] = w[0];
719             x2[1] = w[1];
720             y2[0] = v[0];
721             y2[1] = v[1];
722          }
723       } else
724       {
725          /* "Normal" split code */
726          celt_norm *next_lowband2=NULL;
727          celt_norm *next_lowband_out1=NULL;
728          int next_level=0;
729
730          /* Give more bits to low-energy MDCTs than they would otherwise deserve */
731          if (B>1 && !stereo)
732             delta >>= 1;
733
734          mbits = (b-qalloc-delta)/2;
735          if (mbits > b-qalloc)
736             mbits = b-qalloc;
737          if (mbits<0)
738             mbits=0;
739          sbits = b-qalloc-mbits;
740          *remaining_bits -= qalloc;
741
742          if (lowband && !stereo)
743             next_lowband2 = lowband+N; /* >32-bit split case */
744
745          /* Only stereo needs to pass on lowband_out. Otherwise, it's handled at the end */
746          if (stereo)
747             next_lowband_out1 = lowband_out;
748          else
749             next_level = level+1;
750
751          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);
752          quant_band(encode, m, i, Y, NULL, N, sbits, spread, B, tf_change, next_lowband2, resynth, ec, remaining_bits, LM, NULL, NULL, level, seed);
753       }
754
755    } else {
756       /* This is the basic no-split case */
757       q = bits2pulses(m, m->bits[LM][i], N, b);
758       curr_bits = pulses2bits(m->bits[LM][i], N, q);
759       *remaining_bits -= curr_bits;
760
761       /* Ensures we can never bust the budget */
762       while (*remaining_bits < 0 && q > 0)
763       {
764          *remaining_bits += curr_bits;
765          q--;
766          curr_bits = pulses2bits(m->bits[LM][i], N, q);
767          *remaining_bits -= curr_bits;
768       }
769
770       /* Finally do the actual quantization */
771       if (encode)
772          alg_quant(X, N, q, spread, B, lowband, resynth, (ec_enc*)ec, seed);
773       else
774          alg_unquant(X, N, q, spread, B, lowband, (ec_dec*)ec, seed);
775    }
776
777    /* This code is used by the decoder and by the resynthesis-enabled encoder */
778    if (resynth)
779    {
780       int k;
781
782       if (split)
783       {
784          int j;
785          celt_word16 mid, side;
786 #ifdef FIXED_POINT
787          mid = imid;
788          side = iside;
789 #else
790          mid = (1.f/32768)*imid;
791          side = (1.f/32768)*iside;
792 #endif
793          for (j=0;j<N;j++)
794             X[j] = MULT16_16_Q15(X[j], mid);
795          for (j=0;j<N;j++)
796             Y[j] = MULT16_16_Q15(Y[j], side);
797       }
798
799       /* Undo the sample reorganization going from time order to frequency order */
800       if (!stereo && B0>1 && level==0)
801       {
802          interleave_vector(X, N_B, B0);
803          if (lowband)
804             interleave_vector(lowband, N_B, B0);
805       }
806
807       /* Undo time-freq changes that we did earlier */
808       N_B = N_B0;
809       B = B0;
810       for (k=0;k<time_divide;k++)
811       {
812          B >>= 1;
813          N_B <<= 1;
814          haar1(X, N_B, B);
815          if (lowband)
816             haar1(lowband, N_B, B);
817       }
818
819       for (k=0;k<recombine;k++)
820       {
821          haar1(X, N_B, B);
822          if (lowband)
823             haar1(lowband, N_B, B);
824          N_B>>=1;
825          B <<= 1;
826       }
827
828       /* Scale output for later folding */
829       if (lowband_out && !stereo)
830       {
831          int j;
832          celt_word16 n;
833          n = celt_sqrt(SHL32(EXTEND32(N0),22));
834          for (j=0;j<N0;j++)
835             lowband_out[j] = MULT16_16_Q15(n,X[j]);
836       }
837
838       if (stereo)
839       {
840          stereo_band_mix(m, X, Y, bandE, 0, i, -1, N);
841          /* We only need to renormalize because quantization may not
842             have preserved orthogonality of mid and side */
843          renormalise_vector(X, Q15ONE, N, 1);
844          renormalise_vector(Y, Q15ONE, N, 1);
845       }
846    }
847 }
848
849 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)
850 {
851    int i, balance;
852    celt_int32 remaining_bits;
853    const celt_int16 * restrict eBands = m->eBands;
854    celt_norm * restrict norm;
855    VARDECL(celt_norm, _norm);
856    int B;
857    int M;
858    int spread;
859    celt_int32 seed;
860    celt_norm *lowband;
861    int update_lowband = 1;
862    int C = _Y != NULL ? 2 : 1;
863    SAVE_STACK;
864
865    M = 1<<LM;
866    B = shortBlocks ? M : 1;
867    spread = fold ? B : 0;
868    ALLOC(_norm, M*eBands[m->nbEBands], celt_norm);
869    norm = _norm;
870
871    if (encode)
872       seed = ((ec_enc*)ec)->rng;
873    else
874       seed = ((ec_dec*)ec)->rng;
875    balance = 0;
876    lowband = NULL;
877    for (i=start;i<end;i++)
878    {
879       int tell;
880       int b;
881       int N;
882       int curr_balance;
883       celt_norm * restrict X, * restrict Y;
884       int tf_change=0;
885       celt_norm *effective_lowband;
886       
887       X = _X+M*eBands[i];
888       if (_Y!=NULL)
889          Y = _Y+M*eBands[i];
890       else
891          Y = NULL;
892       N = M*eBands[i+1]-M*eBands[i];
893       if (encode)
894          tell = ec_enc_tell((ec_enc*)ec, BITRES);
895       else
896          tell = ec_dec_tell((ec_dec*)ec, BITRES);
897
898       /* Compute how many bits we want to allocate to this band */
899       if (i != start)
900          balance -= tell;
901       remaining_bits = (total_bits<<BITRES)-tell-1;
902       curr_balance = (end-i);
903       if (curr_balance > 3)
904          curr_balance = 3;
905       curr_balance = balance / curr_balance;
906       b = IMIN(remaining_bits+1,pulses[i]+curr_balance);
907       if (b<0)
908          b = 0;
909       /* Prevents ridiculous bit depths */
910       if (b > C*16*N<<BITRES)
911          b = C*16*N<<BITRES;
912
913       if (M*eBands[i]-N >= M*eBands[start] && (update_lowband || lowband==NULL))
914             lowband = norm+M*eBands[i]-N;
915
916       tf_change = tf_res[i];
917       if (i>=m->effEBands)
918       {
919          X=norm;
920          if (_Y!=NULL)
921             Y = norm;
922       }
923
924       if (tf_change==0 && !shortBlocks && fold)
925          effective_lowband = NULL;
926       else
927          effective_lowband = lowband;
928       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);
929
930       balance += pulses[i] + tell;
931
932       /* Update the folding position only as long as we have 2 bit/sample depth */
933       update_lowband = (b>>BITRES)>2*N;
934    }
935    RESTORE_STACK;
936 }
937