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