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