FFT cleanup
[opus.git] / libcelt / bands.c
1 /* Copyright (c) 2007-2008 CSIRO
2    Copyright (c) 2007-2009 Xiph.Org Foundation
3    Copyright (c) 2008-2009 Gregory Maxwell 
4    Written by Jean-Marc Valin and Gregory Maxwell */
5 /*
6    Redistribution and use in source and binary forms, with or without
7    modification, are permitted provided that the following conditions
8    are met:
9    
10    - Redistributions of source code must retain the above copyright
11    notice, this list of conditions and the following disclaimer.
12    
13    - Redistributions in binary form must reproduce the above copyright
14    notice, this list of conditions and the following disclaimer in the
15    documentation and/or other materials provided with the distribution.
16    
17    - Neither the name of the Xiph.org Foundation nor the names of its
18    contributors may be used to endorse or promote products derived from
19    this software without specific prior written permission.
20    
21    THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
22    ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
23    LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
24    A PARTICULAR PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE FOUNDATION OR
25    CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
26    EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
27    PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
28    PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
29    LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
30    NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
31    SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
32 */
33
34 #ifdef HAVE_CONFIG_H
35 #include "config.h"
36 #endif
37
38 #include <math.h>
39 #include "bands.h"
40 #include "modes.h"
41 #include "vq.h"
42 #include "cwrs.h"
43 #include "stack_alloc.h"
44 #include "os_support.h"
45 #include "mathops.h"
46 #include "rate.h"
47
48 /* This is a cos() approximation designed to be bit-exact on any platform. Bit exactness
49    with this approximation is important because it has an impact on the bit allocation */
50 static celt_int16 bitexact_cos(celt_int16 x)
51 {
52    celt_int32 tmp;
53    celt_int16 x2;
54    tmp = (4096+((celt_int32)(x)*(x)))>>13;
55    if (tmp > 32767)
56       tmp = 32767;
57    x2 = tmp;
58    x2 = (32767-x2) + FRAC_MUL16(x2, (-7651 + FRAC_MUL16(x2, (8277 + FRAC_MUL16(-626, x2)))));
59    if (x2 > 32766)
60       x2 = 32766;
61    return 1+x2;
62 }
63
64
65 #ifdef FIXED_POINT
66 /* Compute the amplitude (sqrt energy) in each of the bands */
67 void compute_band_energies(const CELTMode *m, const celt_sig *X, celt_ener *bank, int end, int _C, int M)
68 {
69    int i, c, N;
70    const celt_int16 *eBands = m->eBands;
71    const int C = CHANNELS(_C);
72    N = M*m->shortMdctSize;
73    for (c=0;c<C;c++)
74    {
75       for (i=0;i<end;i++)
76       {
77          int j;
78          celt_word32 maxval=0;
79          celt_word32 sum = 0;
80          
81          j=M*eBands[i]; do {
82             maxval = MAX32(maxval, X[j+c*N]);
83             maxval = MAX32(maxval, -X[j+c*N]);
84          } while (++j<M*eBands[i+1]);
85          
86          if (maxval > 0)
87          {
88             int shift = celt_ilog2(maxval)-10;
89             j=M*eBands[i]; do {
90                sum = MAC16_16(sum, EXTRACT16(VSHR32(X[j+c*N],shift)),
91                                    EXTRACT16(VSHR32(X[j+c*N],shift)));
92             } while (++j<M*eBands[i+1]);
93             /* We're adding one here to make damn sure we never end up with a pitch vector that's
94                larger than unity norm */
95             bank[i+c*m->nbEBands] = EPSILON+VSHR32(EXTEND32(celt_sqrt(sum)),-shift);
96          } else {
97             bank[i+c*m->nbEBands] = EPSILON;
98          }
99          /*printf ("%f ", bank[i+c*m->nbEBands]);*/
100       }
101    }
102    /*printf ("\n");*/
103 }
104
105 /* Normalise each band such that the energy is one. */
106 void normalise_bands(const CELTMode *m, const celt_sig * restrict freq, celt_norm * restrict X, const celt_ener *bank, int end, int _C, int M)
107 {
108    int i, c, N;
109    const celt_int16 *eBands = m->eBands;
110    const int C = CHANNELS(_C);
111    N = M*m->shortMdctSize;
112    for (c=0;c<C;c++)
113    {
114       i=0; do {
115          celt_word16 g;
116          int j,shift;
117          celt_word16 E;
118          shift = celt_zlog2(bank[i+c*m->nbEBands])-13;
119          E = VSHR32(bank[i+c*m->nbEBands], shift);
120          g = EXTRACT16(celt_rcp(SHL32(E,3)));
121          j=M*eBands[i]; do {
122             X[j+c*N] = MULT16_16_Q15(VSHR32(freq[j+c*N],shift-1),g);
123          } while (++j<M*eBands[i+1]);
124       } while (++i<end);
125    }
126 }
127
128 #else /* FIXED_POINT */
129 /* Compute the amplitude (sqrt energy) in each of the bands */
130 void compute_band_energies(const CELTMode *m, const celt_sig *X, celt_ener *bank, int end, int _C, int M)
131 {
132    int i, c, N;
133    const celt_int16 *eBands = m->eBands;
134    const int C = CHANNELS(_C);
135    N = M*m->shortMdctSize;
136    for (c=0;c<C;c++)
137    {
138       for (i=0;i<end;i++)
139       {
140          int j;
141          celt_word32 sum = 1e-10;
142          for (j=M*eBands[i];j<M*eBands[i+1];j++)
143             sum += X[j+c*N]*X[j+c*N];
144          bank[i+c*m->nbEBands] = sqrt(sum);
145          /*printf ("%f ", bank[i+c*m->nbEBands]);*/
146       }
147    }
148    /*printf ("\n");*/
149 }
150
151 /* Normalise each band such that the energy is one. */
152 void normalise_bands(const CELTMode *m, const celt_sig * restrict freq, celt_norm * restrict X, const celt_ener *bank, int end, int _C, int M)
153 {
154    int i, c, N;
155    const celt_int16 *eBands = m->eBands;
156    const int C = CHANNELS(_C);
157    N = M*m->shortMdctSize;
158    for (c=0;c<C;c++)
159    {
160       for (i=0;i<end;i++)
161       {
162          int j;
163          celt_word16 g = 1.f/(1e-10f+bank[i+c*m->nbEBands]);
164          for (j=M*eBands[i];j<M*eBands[i+1];j++)
165             X[j+c*N] = freq[j+c*N]*g;
166       }
167    }
168 }
169
170 #endif /* FIXED_POINT */
171
172 void renormalise_bands(const CELTMode *m, celt_norm * restrict X, int end, int _C, int M)
173 {
174    int i, c;
175    const celt_int16 *eBands = m->eBands;
176    const int C = CHANNELS(_C);
177    for (c=0;c<C;c++)
178    {
179       i=0; do {
180          renormalise_vector(X+M*eBands[i]+c*M*m->shortMdctSize, Q15ONE, M*eBands[i+1]-M*eBands[i], 1);
181       } while (++i<end);
182    }
183 }
184
185 /* De-normalise the energy to produce the synthesis from the unit-energy bands */
186 void denormalise_bands(const CELTMode *m, const celt_norm * restrict X, celt_sig * restrict freq, const celt_ener *bank, int end, int _C, int M)
187 {
188    int i, c, N;
189    const celt_int16 *eBands = m->eBands;
190    const int C = CHANNELS(_C);
191    N = M*m->shortMdctSize;
192    celt_assert2(C<=2, "denormalise_bands() not implemented for >2 channels");
193    for (c=0;c<C;c++)
194    {
195       celt_sig * restrict f;
196       const celt_norm * restrict x;
197       f = freq+c*N;
198       x = X+c*N;
199       for (i=0;i<end;i++)
200       {
201          int j, band_end;
202          celt_word32 g = SHR32(bank[i+c*m->nbEBands],1);
203          j=M*eBands[i];
204          band_end = M*eBands[i+1];
205          do {
206             *f++ = SHL32(MULT16_32_Q15(*x, g),2);
207             x++;
208          } while (++j<band_end);
209       }
210       for (i=M*eBands[m->nbEBands];i<N;i++)
211          *f++ = 0;
212    }
213 }
214
215 static void stereo_band_mix(const CELTMode *m, celt_norm *X, celt_norm *Y, const celt_ener *bank, int stereo_mode, int bandID, int dir, int N)
216 {
217    int i = bandID;
218    int j;
219    celt_word16 a1, a2;
220    if (stereo_mode==0)
221    {
222       /* Do mid-side when not doing intensity stereo */
223       a1 = QCONST16(.70711f,14);
224       a2 = dir*QCONST16(.70711f,14);
225    } else {
226       celt_word16 left, right;
227       celt_word16 norm;
228 #ifdef FIXED_POINT
229       int shift = celt_zlog2(MAX32(bank[i], bank[i+m->nbEBands]))-13;
230 #endif
231       left = VSHR32(bank[i],shift);
232       right = VSHR32(bank[i+m->nbEBands],shift);
233       norm = EPSILON + celt_sqrt(EPSILON+MULT16_16(left,left)+MULT16_16(right,right));
234       a1 = DIV32_16(SHL32(EXTEND32(left),14),norm);
235       a2 = dir*DIV32_16(SHL32(EXTEND32(right),14),norm);
236    }
237    for (j=0;j<N;j++)
238    {
239       celt_norm r, l;
240       l = X[j];
241       r = Y[j];
242       X[j] = MULT16_16_Q14(a1,l) + MULT16_16_Q14(a2,r);
243       Y[j] = MULT16_16_Q14(a1,r) - MULT16_16_Q14(a2,l);
244    }
245 }
246
247 /* Decide whether we should spread the pulses in the current frame */
248 int folding_decision(const CELTMode *m, celt_norm *X, celt_word16 *average, int *last_decision, int end, int _C, int M)
249 {
250    int i, c, N0;
251    int NR=0;
252    celt_word32 ratio = EPSILON;
253    const int C = CHANNELS(_C);
254    const celt_int16 * restrict eBands = m->eBands;
255    
256    N0 = M*m->shortMdctSize;
257
258    for (c=0;c<C;c++)
259    {
260    for (i=0;i<end;i++)
261    {
262       int j, N;
263       int max_i=0;
264       celt_word16 max_val=EPSILON;
265       celt_word32 floor_ener=EPSILON;
266       celt_norm * restrict x = X+M*eBands[i]+c*N0;
267       N = M*eBands[i+1]-M*eBands[i];
268       for (j=0;j<N;j++)
269       {
270          if (ABS16(x[j])>max_val)
271          {
272             max_val = ABS16(x[j]);
273             max_i = j;
274          }
275       }
276 #if 0
277       for (j=0;j<N;j++)
278       {
279          if (abs(j-max_i)>2)
280             floor_ener += x[j]*x[j];
281       }
282 #else
283       floor_ener = QCONST32(1.,28)-MULT16_16(max_val,max_val);
284       if (max_i < N-1)
285          floor_ener -= MULT16_16(x[(max_i+1)], x[(max_i+1)]);
286       if (max_i < N-2)
287          floor_ener -= MULT16_16(x[(max_i+2)], x[(max_i+2)]);
288       if (max_i > 0)
289          floor_ener -= MULT16_16(x[(max_i-1)], x[(max_i-1)]);
290       if (max_i > 1)
291          floor_ener -= MULT16_16(x[(max_i-2)], x[(max_i-2)]);
292       floor_ener = MAX32(floor_ener, EPSILON);
293 #endif
294       if (N>7)
295       {
296          celt_word16 r;
297          celt_word16 den = celt_sqrt(floor_ener);
298          den = MAX32(QCONST16(.02f, 15), den);
299          r = DIV32_16(SHL32(EXTEND32(max_val),8),den);
300          ratio = ADD32(ratio, EXTEND32(r));
301          NR++;
302       }
303    }
304    }
305    if (NR>0)
306       ratio = DIV32_16(ratio, NR);
307    ratio = ADD32(HALF32(ratio), HALF32(*average));
308    if (!*last_decision)
309    {
310       *last_decision = (ratio < QCONST16(1.8f,8));
311    } else {
312       *last_decision = (ratio < QCONST16(3.f,8));
313    }
314    *average = EXTRACT16(ratio);
315    return *last_decision;
316 }
317
318 #ifdef MEASURE_NORM_MSE
319
320 float MSE[30] = {0};
321 int nbMSEBands = 0;
322 int MSECount[30] = {0};
323
324 void dump_norm_mse(void)
325 {
326    int i;
327    for (i=0;i<nbMSEBands;i++)
328    {
329       printf ("%g ", MSE[i]/MSECount[i]);
330    }
331    printf ("\n");
332 }
333
334 void measure_norm_mse(const CELTMode *m, float *X, float *X0, float *bandE, float *bandE0, int M, int N, int C)
335 {
336    static int init = 0;
337    int i;
338    if (!init)
339    {
340       atexit(dump_norm_mse);
341       init = 1;
342    }
343    for (i=0;i<m->nbEBands;i++)
344    {
345       int j;
346       int c;
347       float g;
348       if (bandE0[i]<10 || (C==2 && bandE0[i+m->nbEBands]<1))
349          continue;
350       for (c=0;c<C;c++)
351       {
352          g = bandE[i+c*m->nbEBands]/(1e-15+bandE0[i+c*m->nbEBands]);
353          for (j=M*m->eBands[i];j<M*m->eBands[i+1];j++)
354             MSE[i] += (g*X[j+c*N]-X0[j+c*N])*(g*X[j+c*N]-X0[j+c*N]);
355       }
356       MSECount[i]+=C;
357    }
358    nbMSEBands = m->nbEBands;
359 }
360
361 #endif
362
363 static void interleave_vector(celt_norm *X, int N0, int stride)
364 {
365    int i,j;
366    VARDECL(celt_norm, tmp);
367    int N;
368    SAVE_STACK;
369    N = N0*stride;
370    ALLOC(tmp, N, celt_norm);
371    for (i=0;i<stride;i++)
372       for (j=0;j<N0;j++)
373          tmp[j*stride+i] = X[i*N0+j];
374    for (j=0;j<N;j++)
375       X[j] = tmp[j];
376    RESTORE_STACK;
377 }
378
379 static void deinterleave_vector(celt_norm *X, int N0, int stride)
380 {
381    int i,j;
382    VARDECL(celt_norm, tmp);
383    int N;
384    SAVE_STACK;
385    N = N0*stride;
386    ALLOC(tmp, N, celt_norm);
387    for (i=0;i<stride;i++)
388       for (j=0;j<N0;j++)
389          tmp[i*N0+j] = X[j*stride+i];
390    for (j=0;j<N;j++)
391       X[j] = tmp[j];
392    RESTORE_STACK;
393 }
394
395 static void haar1(celt_norm *X, int N0, int stride)
396 {
397    int i, j;
398    N0 >>= 1;
399    for (i=0;i<stride;i++)
400       for (j=0;j<N0;j++)
401       {
402          celt_norm tmp = X[stride*2*j+i];
403          X[stride*2*j+i] = MULT16_16_Q15(QCONST16(.7070678f,15), X[stride*2*j+i] + X[stride*(2*j+1)+i]);
404          X[stride*(2*j+1)+i] = MULT16_16_Q15(QCONST16(.7070678f,15), tmp - X[stride*(2*j+1)+i]);
405       }
406 }
407
408 static int compute_qn(int N, int b, int offset, int stereo)
409 {
410    static const celt_int16 exp2_table8[8] =
411       {16384, 17866, 19483, 21247, 23170, 25267, 27554, 30048};
412    int qn, qb;
413    int N2 = 2*N-1;
414    if (stereo && N==2)
415       N2--;
416    qb = (b+N2*offset)/N2;
417    if (qb > (b>>1)-(1<<BITRES))
418       qb = (b>>1)-(1<<BITRES);
419
420    if (qb<0)
421        qb = 0;
422    if (qb>14<<BITRES)
423      qb = 14<<BITRES;
424
425    if (qb<(1<<BITRES>>1)) {
426       qn = 1;
427    } else {
428       qn = exp2_table8[qb&0x7]>>(14-(qb>>BITRES));
429       qn = (qn+1)>>1<<1;
430       if (qn>1024)
431          qn = 1024;
432    }
433    return qn;
434 }
435
436
437 /* This function is responsible for encoding and decoding a band for both
438    the mono and stereo case. Even in the mono case, it can split the band
439    in two and transmit the energy difference with the two half-bands. It
440    can be called recursively so bands can end up being split in 8 parts. */
441 static void quant_band(int encode, const CELTMode *m, int i, celt_norm *X, celt_norm *Y,
442       int N, int b, int spread, int B, int tf_change, celt_norm *lowband, int resynth, void *ec,
443       celt_int32 *remaining_bits, int LM, celt_norm *lowband_out, const celt_ener *bandE, int level, celt_int32 *seed)
444 {
445    int q;
446    int curr_bits;
447    int stereo, split;
448    int imid=0, iside=0;
449    int N0=N;
450    int N_B=N;
451    int N_B0;
452    int B0=B;
453    int time_divide=0;
454    int recombine=0;
455
456    N_B /= B;
457    N_B0 = N_B;
458
459    split = stereo = Y != NULL;
460
461    /* Special case for one sample */
462    if (N==1)
463    {
464       int c;
465       celt_norm *x = X;
466       for (c=0;c<1+stereo;c++)
467       {
468          int sign=0;
469          if (b>=1<<BITRES && *remaining_bits>=1<<BITRES)
470          {
471             if (encode)
472             {
473                sign = x[0]<0;
474                ec_enc_bits((ec_enc*)ec, sign, 1);
475             } else {
476                sign = ec_dec_bits((ec_dec*)ec, 1);
477             }
478             *remaining_bits -= 1<<BITRES;
479             b-=1<<BITRES;
480          }
481          if (resynth)
482             x[0] = sign ? -NORM_SCALING : NORM_SCALING;
483          x = Y;
484       }
485       if (lowband_out)
486          lowband_out[0] = X[0];
487       return;
488    }
489
490    /* Band recombining to increase frequency resolution */
491    if (!stereo && B > 1 && level == 0 && tf_change>0)
492    {
493       while (B>1 && tf_change>0)
494       {
495          B>>=1;
496          N_B<<=1;
497          if (encode)
498             haar1(X, N_B, B);
499          if (lowband)
500             haar1(lowband, N_B, B);
501          recombine++;
502          tf_change--;
503       }
504       B0=B;
505       N_B0 = N_B;
506    }
507
508    /* Increasing the time resolution */
509    if (!stereo && level==0)
510    {
511       while ((N_B&1) == 0 && tf_change<0 && B <= (1<<LM))
512       {
513          if (encode)
514             haar1(X, N_B, B);
515          if (lowband)
516             haar1(lowband, N_B, B);
517          B <<= 1;
518          N_B >>= 1;
519          time_divide++;
520          tf_change++;
521       }
522       B0 = B;
523       N_B0 = N_B;
524    }
525
526    /* Reorganize the samples in time order instead of frequency order */
527    if (!stereo && B0>1 && level==0)
528    {
529       if (encode)
530          deinterleave_vector(X, N_B, B0);
531       if (lowband)
532          deinterleave_vector(lowband, N_B, B0);
533    }
534
535    /* If we need more than 32 bits, try splitting the band in two. */
536    if (!stereo && LM != -1 && b > 32<<BITRES && N>2)
537    {
538       if (LM>0 || (N&1)==0)
539       {
540          N >>= 1;
541          Y = X+N;
542          split = 1;
543          LM -= 1;
544          B = (B+1)>>1;
545       }
546    }
547
548    if (split)
549    {
550       int qn;
551       int itheta=0;
552       int mbits, sbits, delta;
553       int qalloc;
554       celt_word16 mid, side;
555       int offset;
556
557       /* Decide on the resolution to give to the split parameter theta */
558       offset = ((m->logN[i]+(LM<<BITRES))>>1) - (stereo ? QTHETA_OFFSET_STEREO : QTHETA_OFFSET);
559       qn = compute_qn(N, b, offset, stereo);
560
561       qalloc = 0;
562       if (qn!=1)
563       {
564          if (encode)
565          {
566             if (stereo)
567                stereo_band_mix(m, X, Y, bandE, 0, i, 1, N);
568
569             mid = renormalise_vector(X, Q15ONE, N, 1);
570             side = renormalise_vector(Y, Q15ONE, N, 1);
571
572             /* theta is the atan() of the ratio between the (normalized)
573                side and mid. With just that parameter, we can re-scale both
574                mid and side because we know that 1) they have unit norm and
575                2) they are orthogonal. */
576    #ifdef FIXED_POINT
577             /* 0.63662 = 2/pi */
578             itheta = MULT16_16_Q15(QCONST16(0.63662f,15),celt_atan2p(side, mid));
579    #else
580             itheta = floor(.5f+16384*0.63662f*atan2(side,mid));
581    #endif
582
583             itheta = (itheta*qn+8192)>>14;
584          }
585
586          /* Entropy coding of the angle. We use a uniform pdf for the
587             first stereo split but a triangular one for the rest. */
588          if (stereo || qn>256 || B>1)
589          {
590             if (encode)
591                ec_enc_uint((ec_enc*)ec, itheta, qn+1);
592             else
593                itheta = ec_dec_uint((ec_dec*)ec, qn+1);
594             qalloc = log2_frac(qn+1,BITRES);
595          } else {
596             int fs=1, ft;
597             ft = ((qn>>1)+1)*((qn>>1)+1);
598             if (encode)
599             {
600                int fl;
601
602                fs = itheta <= (qn>>1) ? itheta + 1 : qn + 1 - itheta;
603                fl = itheta <= (qn>>1) ? itheta*(itheta + 1)>>1 :
604                 ft - ((qn + 1 - itheta)*(qn + 2 - itheta)>>1);
605
606                ec_encode((ec_enc*)ec, fl, fl+fs, ft);
607             } else {
608                int fl=0;
609                int fm;
610                fm = ec_decode((ec_dec*)ec, ft);
611
612                if (fm < ((qn>>1)*((qn>>1) + 1)>>1))
613                {
614                   itheta = (isqrt32(8*(celt_uint32)fm + 1) - 1)>>1;
615                   fs = itheta + 1;
616                   fl = itheta*(itheta + 1)>>1;
617                }
618                else
619                {
620                   itheta = (2*(qn + 1)
621                    - isqrt32(8*(celt_uint32)(ft - fm - 1) + 1))>>1;
622                   fs = qn + 1 - itheta;
623                   fl = ft - ((qn + 1 - itheta)*(qn + 2 - itheta)>>1);
624                }
625
626                ec_dec_update((ec_dec*)ec, fl, fl+fs, ft);
627             }
628             qalloc = log2_frac(ft,BITRES) - log2_frac(fs,BITRES) + 1;
629          }
630          itheta = itheta*16384/qn;
631       } else {
632          if (stereo && encode)
633             stereo_band_mix(m, X, Y, bandE, 1, i, 1, N);
634       }
635
636       if (itheta == 0)
637       {
638          imid = 32767;
639          iside = 0;
640          delta = -10000;
641       } else if (itheta == 16384)
642       {
643          imid = 0;
644          iside = 32767;
645          delta = 10000;
646       } else {
647          imid = bitexact_cos(itheta);
648          iside = bitexact_cos(16384-itheta);
649          /* This is the mid vs side allocation that minimizes squared error
650             in that band. */
651          delta = (N-1)*(log2_frac(iside,BITRES+2)-log2_frac(imid,BITRES+2))>>2;
652       }
653
654       /* This is a special case for N=2 that only works for stereo and takes
655          advantage of the fact that mid and side are orthogonal to encode
656          the side with just one bit. */
657       if (N==2 && stereo)
658       {
659          int c;
660          int sign=1;
661          celt_norm *x2, *y2;
662          mbits = b-qalloc;
663          sbits = 0;
664          /* Only need one bit for the side */
665          if (itheta != 0 && itheta != 16384)
666             sbits = 1<<BITRES;
667          mbits -= sbits;
668          c = itheta > 8192;
669          *remaining_bits -= qalloc+sbits;
670
671          x2 = c ? Y : X;
672          y2 = c ? X : Y;
673          if (sbits)
674          {
675             if (encode)
676             {
677                /* Here we only need to encode a sign for the side */
678                sign = x2[0]*y2[1] - x2[1]*y2[0] > 0;
679                ec_enc_bits((ec_enc*)ec, sign, 1);
680             } else {
681                sign = ec_dec_bits((ec_dec*)ec, 1);
682             }
683          }
684          sign = 2*sign - 1;
685          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);
686          y2[0] = -sign*x2[1];
687          y2[1] = sign*x2[0];
688       } else {
689          /* "Normal" split code */
690          celt_norm *next_lowband2=NULL;
691          celt_norm *next_lowband_out1=NULL;
692          int next_level=0;
693
694          /* Give more bits to low-energy MDCTs than they would otherwise deserve */
695          if (B>1 && !stereo)
696             delta >>= 1;
697
698          mbits = (b-qalloc-delta)/2;
699          if (mbits > b-qalloc)
700             mbits = b-qalloc;
701          if (mbits<0)
702             mbits=0;
703          sbits = b-qalloc-mbits;
704          *remaining_bits -= qalloc;
705
706          if (lowband && !stereo)
707             next_lowband2 = lowband+N; /* >32-bit split case */
708
709          /* Only stereo needs to pass on lowband_out. Otherwise, it's handled at the end */
710          if (stereo)
711             next_lowband_out1 = lowband_out;
712          else
713             next_level = level+1;
714
715          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);
716          quant_band(encode, m, i, Y, NULL, N, sbits, spread, B, tf_change, next_lowband2, resynth, ec, remaining_bits, LM, NULL, NULL, level, seed);
717       }
718
719    } else {
720       /* This is the basic no-split case */
721       q = bits2pulses(m, i, LM, b);
722       curr_bits = pulses2bits(m, i, LM, q);
723       *remaining_bits -= curr_bits;
724
725       /* Ensures we can never bust the budget */
726       while (*remaining_bits < 0 && q > 0)
727       {
728          *remaining_bits += curr_bits;
729          q--;
730          curr_bits = pulses2bits(m, i, LM, q);
731          *remaining_bits -= curr_bits;
732       }
733
734       /* Finally do the actual quantization */
735       if (encode)
736          alg_quant(X, N, q, spread, B, lowband, resynth, (ec_enc*)ec, seed);
737       else
738          alg_unquant(X, N, q, spread, B, lowband, (ec_dec*)ec, seed);
739    }
740
741    /* This code is used by the decoder and by the resynthesis-enabled encoder */
742    if (resynth)
743    {
744       int k;
745
746       if (split)
747       {
748          int j;
749          celt_word16 mid, side;
750 #ifdef FIXED_POINT
751          mid = imid;
752          side = iside;
753 #else
754          mid = (1.f/32768)*imid;
755          side = (1.f/32768)*iside;
756 #endif
757          for (j=0;j<N;j++)
758             X[j] = MULT16_16_Q15(X[j], mid);
759          for (j=0;j<N;j++)
760             Y[j] = MULT16_16_Q15(Y[j], side);
761       }
762
763       /* Undo the sample reorganization going from time order to frequency order */
764       if (!stereo && B0>1 && level==0)
765       {
766          interleave_vector(X, N_B, B0);
767          if (lowband)
768             interleave_vector(lowband, N_B, B0);
769       }
770
771       /* Undo time-freq changes that we did earlier */
772       N_B = N_B0;
773       B = B0;
774       for (k=0;k<time_divide;k++)
775       {
776          B >>= 1;
777          N_B <<= 1;
778          haar1(X, N_B, B);
779          if (lowband)
780             haar1(lowband, N_B, B);
781       }
782
783       for (k=0;k<recombine;k++)
784       {
785          haar1(X, N_B, B);
786          if (lowband)
787             haar1(lowband, N_B, B);
788          N_B>>=1;
789          B <<= 1;
790       }
791
792       /* Scale output for later folding */
793       if (lowband_out && !stereo)
794       {
795          int j;
796          celt_word16 n;
797          n = celt_sqrt(SHL32(EXTEND32(N0),22));
798          for (j=0;j<N0;j++)
799             lowband_out[j] = MULT16_16_Q15(n,X[j]);
800       }
801
802       if (stereo)
803       {
804          stereo_band_mix(m, X, Y, bandE, 0, i, -1, N);
805          /* We only need to renormalize because quantization may not
806             have preserved orthogonality of mid and side */
807          renormalise_vector(X, Q15ONE, N, 1);
808          renormalise_vector(Y, Q15ONE, N, 1);
809       }
810    }
811 }
812
813 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)
814 {
815    int i, balance;
816    celt_int32 remaining_bits;
817    const celt_int16 * restrict eBands = m->eBands;
818    celt_norm * restrict norm;
819    VARDECL(celt_norm, _norm);
820    int B;
821    int M;
822    celt_int32 seed;
823    celt_norm *lowband;
824    int update_lowband = 1;
825    int C = _Y != NULL ? 2 : 1;
826    SAVE_STACK;
827
828    M = 1<<LM;
829    B = shortBlocks ? M : 1;
830    ALLOC(_norm, M*eBands[m->nbEBands], celt_norm);
831    norm = _norm;
832
833    if (encode)
834       seed = ((ec_enc*)ec)->rng;
835    else
836       seed = ((ec_dec*)ec)->rng;
837    balance = 0;
838    lowband = NULL;
839    for (i=start;i<end;i++)
840    {
841       int tell;
842       int b;
843       int N;
844       int curr_balance;
845       celt_norm * restrict X, * restrict Y;
846       int tf_change=0;
847       celt_norm *effective_lowband;
848       
849       X = _X+M*eBands[i];
850       if (_Y!=NULL)
851          Y = _Y+M*eBands[i];
852       else
853          Y = NULL;
854       N = M*eBands[i+1]-M*eBands[i];
855       if (encode)
856          tell = ec_enc_tell((ec_enc*)ec, BITRES);
857       else
858          tell = ec_dec_tell((ec_dec*)ec, BITRES);
859
860       /* Compute how many bits we want to allocate to this band */
861       if (i != start)
862          balance -= tell;
863       remaining_bits = (total_bits<<BITRES)-tell-1;
864       curr_balance = (end-i);
865       if (curr_balance > 3)
866          curr_balance = 3;
867       curr_balance = balance / curr_balance;
868       b = IMIN(remaining_bits+1,pulses[i]+curr_balance);
869       if (b<0)
870          b = 0;
871       /* Prevents ridiculous bit depths */
872       if (b > C*16*N<<BITRES)
873          b = C*16*N<<BITRES;
874
875       if (M*eBands[i]-N >= M*eBands[start] && (update_lowband || lowband==NULL))
876             lowband = norm+M*eBands[i]-N;
877
878       tf_change = tf_res[i];
879       if (i>=m->effEBands)
880       {
881          X=norm;
882          if (_Y!=NULL)
883             Y = norm;
884       }
885
886       if (tf_change==0 && !shortBlocks && fold)
887          effective_lowband = NULL;
888       else
889          effective_lowband = lowband;
890       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);
891
892       balance += pulses[i] + tell;
893
894       /* Update the folding position only as long as we have 2 bit/sample depth */
895       update_lowband = (b>>BITRES)>2*N;
896    }
897    RESTORE_STACK;
898 }
899