Moving code to quant_band()
[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
49 #ifdef FIXED_POINT
50 /* Compute the amplitude (sqrt energy) in each of the bands */
51 void compute_band_energies(const CELTMode *m, const celt_sig *X, celt_ener *bank, int _C, int M)
52 {
53    int i, c, N;
54    const celt_int16 *eBands = m->eBands;
55    const int C = CHANNELS(_C);
56    N = M*m->eBands[m->nbEBands+1];
57    for (c=0;c<C;c++)
58    {
59       for (i=0;i<m->nbEBands;i++)
60       {
61          int j;
62          celt_word32 maxval=0;
63          celt_word32 sum = 0;
64          
65          j=M*eBands[i]; do {
66             maxval = MAX32(maxval, X[j+c*N]);
67             maxval = MAX32(maxval, -X[j+c*N]);
68          } while (++j<M*eBands[i+1]);
69          
70          if (maxval > 0)
71          {
72             int shift = celt_ilog2(maxval)-10;
73             j=M*eBands[i]; do {
74                sum = MAC16_16(sum, EXTRACT16(VSHR32(X[j+c*N],shift)),
75                                    EXTRACT16(VSHR32(X[j+c*N],shift)));
76             } while (++j<M*eBands[i+1]);
77             /* We're adding one here to make damn sure we never end up with a pitch vector that's
78                larger than unity norm */
79             bank[i+c*m->nbEBands] = EPSILON+VSHR32(EXTEND32(celt_sqrt(sum)),-shift);
80          } else {
81             bank[i+c*m->nbEBands] = EPSILON;
82          }
83          /*printf ("%f ", bank[i+c*m->nbEBands]);*/
84       }
85    }
86    /*printf ("\n");*/
87 }
88
89 /* Normalise each band such that the energy is one. */
90 void normalise_bands(const CELTMode *m, const celt_sig * restrict freq, celt_norm * restrict X, const celt_ener *bank, int _C, int M)
91 {
92    int i, c, N;
93    const celt_int16 *eBands = m->eBands;
94    const int C = CHANNELS(_C);
95    N = M*m->eBands[m->nbEBands+1];
96    for (c=0;c<C;c++)
97    {
98       i=0; do {
99          celt_word16 g;
100          int j,shift;
101          celt_word16 E;
102          shift = celt_zlog2(bank[i+c*m->nbEBands])-13;
103          E = VSHR32(bank[i+c*m->nbEBands], shift);
104          g = EXTRACT16(celt_rcp(SHL32(E,3)));
105          j=M*eBands[i]; do {
106             X[j+c*N] = MULT16_16_Q15(VSHR32(freq[j+c*N],shift-1),g);
107          } while (++j<M*eBands[i+1]);
108       } while (++i<m->nbEBands);
109    }
110 }
111
112 #else /* FIXED_POINT */
113 /* Compute the amplitude (sqrt energy) in each of the bands */
114 void compute_band_energies(const CELTMode *m, const celt_sig *X, celt_ener *bank, int _C, int M)
115 {
116    int i, c, N;
117    const celt_int16 *eBands = m->eBands;
118    const int C = CHANNELS(_C);
119    N = M*m->eBands[m->nbEBands+1];
120    for (c=0;c<C;c++)
121    {
122       for (i=0;i<m->nbEBands;i++)
123       {
124          int j;
125          celt_word32 sum = 1e-10;
126          for (j=M*eBands[i];j<M*eBands[i+1];j++)
127             sum += X[j+c*N]*X[j+c*N];
128          bank[i+c*m->nbEBands] = sqrt(sum);
129          /*printf ("%f ", bank[i+c*m->nbEBands]);*/
130       }
131    }
132    /*printf ("\n");*/
133 }
134
135 /* Normalise each band such that the energy is one. */
136 void normalise_bands(const CELTMode *m, const celt_sig * restrict freq, celt_norm * restrict X, const celt_ener *bank, int _C, int M)
137 {
138    int i, c, N;
139    const celt_int16 *eBands = m->eBands;
140    const int C = CHANNELS(_C);
141    N = M*m->eBands[m->nbEBands+1];
142    for (c=0;c<C;c++)
143    {
144       for (i=0;i<m->nbEBands;i++)
145       {
146          int j;
147          celt_word16 g = 1.f/(1e-10f+bank[i+c*m->nbEBands]);
148          for (j=M*eBands[i];j<M*eBands[i+1];j++)
149             X[j+c*N] = freq[j+c*N]*g;
150       }
151    }
152 }
153
154 #endif /* FIXED_POINT */
155
156 void renormalise_bands(const CELTMode *m, celt_norm * restrict X, int _C, int M)
157 {
158    int i, c;
159    const celt_int16 *eBands = m->eBands;
160    const int C = CHANNELS(_C);
161    for (c=0;c<C;c++)
162    {
163       i=0; do {
164          renormalise_vector(X+M*eBands[i]+c*M*eBands[m->nbEBands+1], Q15ONE, M*eBands[i+1]-M*eBands[i], 1);
165       } while (++i<m->nbEBands);
166    }
167 }
168
169 /* De-normalise the energy to produce the synthesis from the unit-energy bands */
170 void denormalise_bands(const CELTMode *m, const celt_norm * restrict X, celt_sig * restrict freq, const celt_ener *bank, int _C, int M)
171 {
172    int i, c, N;
173    const celt_int16 *eBands = m->eBands;
174    const int C = CHANNELS(_C);
175    N = M*m->eBands[m->nbEBands+1];
176    if (C>2)
177       celt_fatal("denormalise_bands() not implemented for >2 channels");
178    for (c=0;c<C;c++)
179    {
180       celt_sig * restrict f;
181       const celt_norm * restrict x;
182       f = freq+c*N;
183       x = X+c*N;
184       for (i=0;i<m->nbEBands;i++)
185       {
186          int j, end;
187          celt_word32 g = SHR32(bank[i+c*m->nbEBands],1);
188          j=M*eBands[i];
189          end = M*eBands[i+1];
190          do {
191             *f++ = SHL32(MULT16_32_Q15(*x, g),2);
192             x++;
193          } while (++j<end);
194       }
195       for (i=M*eBands[m->nbEBands];i<M*eBands[m->nbEBands+1];i++)
196          *f++ = 0;
197    }
198 }
199
200 int compute_pitch_gain(const CELTMode *m, const celt_sig *X, const celt_sig *P, int norm_rate, int *gain_id, int _C, celt_word16 *gain_prod, int M)
201 {
202    int j, c;
203    celt_word16 g;
204    celt_word16 delta;
205    const int C = CHANNELS(_C);
206    celt_word32 Sxy=0, Sxx=0, Syy=0;
207    int len = M*m->pitchEnd;
208    int N = M*m->eBands[m->nbEBands+1];
209 #ifdef FIXED_POINT
210    int shift = 0;
211    celt_word32 maxabs=0;
212
213    for (c=0;c<C;c++)
214    {
215       for (j=0;j<len;j++)
216       {
217          maxabs = MAX32(maxabs, ABS32(X[j+c*N]));
218          maxabs = MAX32(maxabs, ABS32(P[j+c*N]));
219       }
220    }
221    shift = celt_ilog2(maxabs)-12;
222    if (shift<0)
223       shift = 0;
224 #endif
225    delta = PDIV32_16(Q15ONE, len);
226    for (c=0;c<C;c++)
227    {
228       celt_word16 gg = Q15ONE;
229       for (j=0;j<len;j++)
230       {
231          celt_word16 Xj, Pj;
232          Xj = EXTRACT16(SHR32(X[j+c*N], shift));
233          Pj = MULT16_16_P15(gg,EXTRACT16(SHR32(P[j+c*N], shift)));
234          Sxy = MAC16_16(Sxy, Xj, Pj);
235          Sxx = MAC16_16(Sxx, Pj, Pj);
236          Syy = MAC16_16(Syy, Xj, Xj);
237          gg = SUB16(gg, delta);
238       }
239    }
240 #ifdef FIXED_POINT
241    {
242       celt_word32 num, den;
243       celt_word16 fact;
244       fact = MULT16_16(QCONST16(.04f, 14), norm_rate);
245       if (fact < QCONST16(1.f, 14))
246          fact = QCONST16(1.f, 14);
247       num = Sxy;
248       den = EPSILON+Sxx+MULT16_32_Q15(QCONST16(.03f,15),Syy);
249       shift = celt_zlog2(Sxy)-16;
250       if (shift < 0)
251          shift = 0;
252       if (Sxy < MULT16_32_Q15(fact, MULT16_16(celt_sqrt(EPSILON+Sxx),celt_sqrt(EPSILON+Syy))))
253          g = 0;
254       else
255          g = DIV32(SHL32(SHR32(num,shift),14),ADD32(EPSILON,SHR32(den,shift)));
256
257       /* This MUST round down so that we don't over-estimate the gain */
258       *gain_id = EXTRACT16(SHR32(MULT16_16(20,(g-QCONST16(.5f,14))),14));
259    }
260 #else
261    {
262       float fact = .04f*norm_rate;
263       if (fact < 1)
264          fact = 1;
265       g = Sxy/(.1f+Sxx+.03f*Syy);
266       if (Sxy < .5f*fact*celt_sqrt(1+Sxx*Syy))
267          g = 0;
268       /* This MUST round down so that we don't over-estimate the gain */
269       *gain_id = floor(20*(g-.5f));
270    }
271 #endif
272    /* This prevents the pitch gain from being above 1.0 for too long by bounding the 
273       maximum error amplification factor to 2.0 */
274    g = ADD16(QCONST16(.5f,14), MULT16_16_16(QCONST16(.05f,14),*gain_id));
275    *gain_prod = MAX16(QCONST32(1.f, 13), MULT16_16_Q14(*gain_prod,g));
276    if (*gain_prod>QCONST32(2.f, 13))
277    {
278       *gain_id=9;
279       *gain_prod = QCONST32(2.f, 13);
280    }
281
282    if (*gain_id < 0)
283    {
284       *gain_id = 0;
285       return 0;
286    } else {
287       if (*gain_id > 15)
288          *gain_id = 15;
289       return 1;
290    }
291 }
292
293 void apply_pitch(const CELTMode *m, celt_sig *X, const celt_sig *P, int gain_id, int pred, int _C, int M)
294 {
295    int j, c, N;
296    celt_word16 gain;
297    celt_word16 delta;
298    const int C = CHANNELS(_C);
299    int len = M*m->pitchEnd;
300
301    N = M*m->eBands[m->nbEBands+1];
302    gain = ADD16(QCONST16(.5f,14), MULT16_16_16(QCONST16(.05f,14),gain_id));
303    delta = PDIV32_16(gain, len);
304    if (pred)
305       gain = -gain;
306    else
307       delta = -delta;
308    for (c=0;c<C;c++)
309    {
310       celt_word16 gg = gain;
311       for (j=0;j<len;j++)
312       {
313          X[j+c*N] += SHL32(MULT16_32_Q15(gg,P[j+c*N]),1);
314          gg = ADD16(gg, delta);
315       }
316    }
317 }
318
319 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)
320 {
321    int i = bandID;
322    const celt_int16 *eBands = m->eBands;
323    int j;
324    celt_word16 a1, a2;
325    if (stereo_mode==0)
326    {
327       /* Do mid-side when not doing intensity stereo */
328       a1 = QCONST16(.70711f,14);
329       a2 = dir*QCONST16(.70711f,14);
330    } else {
331       celt_word16 left, right;
332       celt_word16 norm;
333 #ifdef FIXED_POINT
334       int shift = celt_zlog2(MAX32(bank[i], bank[i+m->nbEBands]))-13;
335 #endif
336       left = VSHR32(bank[i],shift);
337       right = VSHR32(bank[i+m->nbEBands],shift);
338       norm = EPSILON + celt_sqrt(EPSILON+MULT16_16(left,left)+MULT16_16(right,right));
339       a1 = DIV32_16(SHL32(EXTEND32(left),14),norm);
340       a2 = dir*DIV32_16(SHL32(EXTEND32(right),14),norm);
341    }
342    for (j=0;j<N;j++)
343    {
344       celt_norm r, l;
345       l = X[j];
346       r = Y[j];
347       X[j] = MULT16_16_Q14(a1,l) + MULT16_16_Q14(a2,r);
348       Y[j] = MULT16_16_Q14(a1,r) - MULT16_16_Q14(a2,l);
349    }
350 }
351
352
353 int folding_decision(const CELTMode *m, celt_norm *X, celt_word16 *average, int *last_decision, int _C, int M)
354 {
355    int i, c, N0;
356    int NR=0;
357    celt_word32 ratio = EPSILON;
358    const int C = CHANNELS(_C);
359    const celt_int16 * restrict eBands = m->eBands;
360    
361    N0 = M*m->eBands[m->nbEBands+1];
362
363    for (c=0;c<C;c++)
364    {
365    for (i=0;i<m->nbEBands;i++)
366    {
367       int j, N;
368       int max_i=0;
369       celt_word16 max_val=EPSILON;
370       celt_word32 floor_ener=EPSILON;
371       celt_norm * restrict x = X+M*eBands[i]+c*N0;
372       N = M*eBands[i+1]-M*eBands[i];
373       for (j=0;j<N;j++)
374       {
375          if (ABS16(x[j])>max_val)
376          {
377             max_val = ABS16(x[j]);
378             max_i = j;
379          }
380       }
381 #if 0
382       for (j=0;j<N;j++)
383       {
384          if (abs(j-max_i)>2)
385             floor_ener += x[j]*x[j];
386       }
387 #else
388       floor_ener = QCONST32(1.,28)-MULT16_16(max_val,max_val);
389       if (max_i < N-1)
390          floor_ener -= MULT16_16(x[(max_i+1)], x[(max_i+1)]);
391       if (max_i < N-2)
392          floor_ener -= MULT16_16(x[(max_i+2)], x[(max_i+2)]);
393       if (max_i > 0)
394          floor_ener -= MULT16_16(x[(max_i-1)], x[(max_i-1)]);
395       if (max_i > 1)
396          floor_ener -= MULT16_16(x[(max_i-2)], x[(max_i-2)]);
397       floor_ener = MAX32(floor_ener, EPSILON);
398 #endif
399       if (N>7)
400       {
401          celt_word16 r;
402          celt_word16 den = celt_sqrt(floor_ener);
403          den = MAX32(QCONST16(.02f, 15), den);
404          r = DIV32_16(SHL32(EXTEND32(max_val),8),den);
405          ratio = ADD32(ratio, EXTEND32(r));
406          NR++;
407       }
408    }
409    }
410    if (NR>0)
411       ratio = DIV32_16(ratio, NR);
412    ratio = ADD32(HALF32(ratio), HALF32(*average));
413    if (!*last_decision)
414    {
415       *last_decision = (ratio < QCONST16(1.8f,8));
416    } else {
417       *last_decision = (ratio < QCONST16(3.f,8));
418    }
419    *average = EXTRACT16(ratio);
420    return *last_decision;
421 }
422
423 /* This function is responsible for encoding and decoding a band for both
424    the mono and stereo case. Even in the mono case, it can split the band
425    in two and transmit the energy difference with the two half-bands. It
426    can be called recursively so bands can end up being split in 8 parts. */
427 static void quant_band(int encode, const CELTMode *m, int i, celt_norm *X, celt_norm *Y, int N, int b, int spread, celt_norm *lowband, int resynth, ec_enc *ec, celt_int32 *remaining_bits, int LM, celt_norm *lowband_out, const celt_ener *bandE)
428 {
429    int q;
430    int curr_bits;
431    int stereo, split;
432    int imid=0, iside=0;
433    int N0=N;
434
435    split = stereo = Y != NULL;
436
437    /* If we need more than 32 bits, try splitting the band in two. */
438    if (!stereo && LM != -1 && !fits_in32(N, get_pulses(bits2pulses(m, m->bits[LM][i], N, b))))
439    {
440       if (LM>0 || (N&1)==0)
441       {
442          N >>= 1;
443          Y = X+N;
444          split = 1;
445          LM -= 1;
446       }
447    }
448
449    if (split)
450    {
451       int qb;
452       int itheta;
453       int mbits, sbits, delta;
454       int qalloc;
455       celt_word16 mid, side;
456       if (N>1)
457       {
458          qb = (b-2*(N-1)*(QTHETA_OFFSET-m->logN[i]-(LM<<BITRES)))/(2*(N-1)<<BITRES);
459          if (qb > (b>>(BITRES+1))-1)
460             qb = (b>>(BITRES+1))-1;
461       } else {
462          /* For N==1, allocate two bits for the signs and the rest goes to qb */
463          qb = b-(2<<BITRES);
464       }
465       if (qb<0)
466          qb = 0;
467       if (qb>14)
468          qb = 14;
469
470       if (encode)
471       {
472          if (stereo)
473             stereo_band_mix(m, X, Y, bandE, qb==0, i, 1, N);
474
475          mid = renormalise_vector(X, Q15ONE, N, 1);
476          side = renormalise_vector(Y, Q15ONE, N, 1);
477
478          /* 0.63662 = 2/pi */
479 #ifdef FIXED_POINT
480          itheta = MULT16_16_Q15(QCONST16(0.63662f,15),celt_atan2p(side, mid));
481 #else
482          itheta = floor(.5f+16384*0.63662f*atan2(side,mid));
483 #endif
484       }
485
486       qalloc = 0;
487       if (qb==0)
488       {
489          itheta=0;
490       } else {
491          int shift;
492          shift = 14-qb;
493
494          /* Entropy coding of the angle. We use a uniform pdf for the
495             first stereo split but a triangular one for the rest. */
496          if (encode)
497             itheta = (itheta+(1<<shift>>1))>>shift;
498          if (stereo || qb>9)
499          {
500             if (encode)
501                ec_enc_uint(ec, itheta, (1<<qb)+1);
502             else
503                itheta = ec_dec_uint((ec_dec*)ec, (1<<qb)+1);
504             qalloc = log2_frac((1<<qb)+1,BITRES);
505          } else {
506             int fs=1, ft;
507             ft = ((1<<qb>>1)+1)*((1<<qb>>1)+1);
508             if (encode)
509             {
510                int j;
511                int fl=0;
512                j=0;
513                while(1)
514                {
515                   if (j==itheta)
516                      break;
517                   fl+=fs;
518                   if (j<(1<<qb>>1))
519                      fs++;
520                   else
521                      fs--;
522                   j++;
523                }
524                ec_encode(ec, fl, fl+fs, ft);
525             } else {
526                int fl=0;
527                int j, fm;
528                fm = ec_decode((ec_dec*)ec, ft);
529                j=0;
530                while (1)
531                {
532                   if (fm < fl+fs)
533                      break;
534                   fl+=fs;
535                   if (j<(1<<qb>>1))
536                      fs++;
537                   else
538                      fs--;
539                   j++;
540                }
541                itheta = j;
542                ec_dec_update((ec_dec*)ec, fl, fl+fs, ft);
543             }
544             qalloc = log2_frac(ft,BITRES) - log2_frac(fs,BITRES) + 1;
545          }
546          itheta <<= shift;
547       }
548
549       if (itheta == 0)
550       {
551          imid = 32767;
552          iside = 0;
553          delta = -10000;
554       } else if (itheta == 16384)
555       {
556          imid = 0;
557          iside = 32767;
558          delta = 10000;
559       } else {
560          imid = bitexact_cos(itheta);
561          iside = bitexact_cos(16384-itheta);
562          delta = (N-1)*(log2_frac(iside,BITRES+2)-log2_frac(imid,BITRES+2))>>2;
563       }
564
565       /* This is a special case for N=2 that only works for stereo and takes
566          advantage of the fact that mid and side are orthogonal to encode
567          the side with just one bit. */
568       if (N==2 && stereo)
569       {
570          int c, c2;
571          int sign=1;
572          celt_norm v[2], w[2];
573          celt_norm *x2, *y2;
574          mbits = b-qalloc;
575          sbits = 0;
576          if (itheta != 0 && itheta != 16384)
577             sbits = 1<<BITRES;
578          mbits -= sbits;
579          c = itheta > 8192 ? 1 : 0;
580          *remaining_bits -= qalloc+sbits;
581
582          x2 = X;
583          y2 = Y;
584          if (encode)
585          {
586             c2 = 1-c;
587
588             if (c==0)
589             {
590                v[0] = x2[0];
591                v[1] = x2[1];
592                w[0] = y2[0];
593                w[1] = y2[1];
594             } else {
595                v[0] = y2[0];
596                v[1] = y2[1];
597                w[0] = x2[0];
598                w[1] = x2[1];
599             }
600          }
601          quant_band(encode, m, i, v, NULL, N, mbits, spread, lowband, resynth, ec, remaining_bits, LM, NULL, NULL);
602          if (sbits)
603          {
604             if (encode)
605             {
606                /* Here we only need to encode a sign for the side */
607                if (v[0]*w[1] - v[1]*w[0] > 0)
608                   sign = 1;
609                else
610                   sign = -1;
611                ec_enc_bits(ec, sign==1, 1);
612             } else {
613                sign = 2*ec_dec_bits((ec_dec*)ec, 1)-1;
614             }
615          } else {
616             sign = 1;
617          }
618          w[0] = -sign*v[1];
619          w[1] = sign*v[0];
620          if (c==0)
621          {
622             x2[0] = v[0];
623             x2[1] = v[1];
624             y2[0] = w[0];
625             y2[1] = w[1];
626          } else {
627             x2[0] = w[0];
628             x2[1] = w[1];
629             y2[0] = v[0];
630             y2[1] = v[1];
631          }
632       } else
633       {
634          /* "Normal" split code */
635          mbits = (b-qalloc/2-delta)/2;
636          if (mbits > b-qalloc)
637             mbits = b-qalloc;
638          if (mbits<0)
639             mbits=0;
640          sbits = b-qalloc-mbits;
641          *remaining_bits -= qalloc;
642          quant_band(encode, m, i, X, NULL, N, mbits, spread, lowband, resynth, ec, remaining_bits, LM, NULL, NULL);
643          if (stereo)
644             quant_band(encode, m, i, Y, NULL, N, sbits, spread, NULL, resynth, ec, remaining_bits, LM, NULL, NULL);
645          else
646             quant_band(encode, m, i, Y, NULL, N, sbits, spread, lowband ? lowband+N : NULL, resynth, ec, remaining_bits, LM, NULL, NULL);
647       }
648
649    } else {
650       /* This is the basic no-split case */
651       q = bits2pulses(m, m->bits[LM][i], N, b);
652       curr_bits = pulses2bits(m->bits[LM][i], N, q);
653       *remaining_bits -= curr_bits;
654
655       /* Ensures we can never bust the budget */
656       while (*remaining_bits < 0 && q > 0)
657       {
658          *remaining_bits += curr_bits;
659          q--;
660          curr_bits = pulses2bits(m->bits[LM][i], N, q);
661          *remaining_bits -= curr_bits;
662       }
663
664       /* Making sure we will *never* need more than 32 bits for the PVQ */
665       while (!fits_in32(N, get_pulses(q)))
666          q--;
667
668       if (encode)
669          alg_quant(X, N, q, spread, lowband, resynth, ec);
670       else
671          alg_unquant(X, N, q, spread, lowband, (ec_dec*)ec);
672    }
673
674    if (resynth && lowband_out)
675    {
676       int j;
677       celt_word16 n;
678       n = celt_sqrt(SHL32(EXTEND32(N0),22));
679       for (j=0;j<N0;j++)
680          lowband_out[j] = MULT16_16_Q15(n,X[j]);
681    }
682
683    if (split && resynth)
684    {
685       int j;
686       celt_word16 mid, side;
687 #ifdef FIXED_POINT
688       mid = imid;
689       side = iside;
690 #else
691       mid = (1.f/32768)*imid;
692       side = (1.f/32768)*iside;
693 #endif
694       for (j=0;j<N;j++)
695          X[j] = MULT16_16_Q15(X[j], mid);
696       for (j=0;j<N;j++)
697          Y[j] = MULT16_16_Q15(Y[j], side);
698
699       if (stereo)
700       {
701          stereo_band_mix(m, X, Y, bandE, 0, i, -1, N);
702          renormalise_vector(X, Q15ONE, N, 1);
703          renormalise_vector(Y, Q15ONE, N, 1);
704       }
705    }
706 }
707
708 void quant_all_bands(int encode, const CELTMode *m, int start, celt_norm *_X, celt_norm *_Y, const celt_ener *bandE, int *pulses, int shortBlocks, int fold, int resynth, int total_bits, ec_enc *ec, int LM)
709 {
710    int i, remaining_bits, balance;
711    const celt_int16 * restrict eBands = m->eBands;
712    celt_norm * restrict norm;
713    VARDECL(celt_norm, _norm);
714    int B;
715    int M;
716    int spread;
717    SAVE_STACK;
718
719    M = 1<<LM;
720    B = shortBlocks ? M : 1;
721    spread = fold ? B : 0;
722    ALLOC(_norm, M*eBands[m->nbEBands+1], celt_norm);
723    norm = _norm;
724    /* Just in case the first bands attempts to fold -- not that rare for stereo */
725    for (i=0;i<M;i++)
726       norm[i] = 0;
727
728    balance = 0;
729    for (i=start;i<m->nbEBands;i++)
730    {
731       int tell;
732       int b;
733       int N;
734       int curr_balance;
735       celt_norm * restrict X, * restrict Y;
736       
737       X = _X+M*eBands[i];
738       if (_Y!=NULL)
739          Y = _Y+M*eBands[i];
740       else
741          Y = NULL;
742       N = M*eBands[i+1]-M*eBands[i];
743       if (encode)
744          tell = ec_enc_tell(ec, BITRES);
745       else
746          tell = ec_dec_tell((ec_dec*)ec, BITRES);
747
748       if (i != start)
749          balance -= tell;
750       remaining_bits = (total_bits<<BITRES)-tell-1;
751       curr_balance = (m->nbEBands-i);
752       if (curr_balance > 3)
753          curr_balance = 3;
754       curr_balance = balance / curr_balance;
755       b = IMIN(remaining_bits+1,pulses[i]+curr_balance);
756       if (b<0)
757          b = 0;
758
759       quant_band(encode, m, i, X, Y, N, b, spread, norm+M*eBands[start], resynth, ec, &remaining_bits, LM, norm+M*eBands[i], bandE);
760
761       balance += pulses[i] + tell;
762    }
763    RESTORE_STACK;
764 }
765