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