920a9e961a812a5d03e8b545e908dcc0518cab5b
[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    c=0; do {
74       for (i=0;i<end;i++)
75       {
76          int j;
77          celt_word32 maxval=0;
78          celt_word32 sum = 0;
79          
80          j=M*eBands[i]; do {
81             maxval = MAX32(maxval, X[j+c*N]);
82             maxval = MAX32(maxval, -X[j+c*N]);
83          } while (++j<M*eBands[i+1]);
84          
85          if (maxval > 0)
86          {
87             int shift = celt_ilog2(maxval)-10;
88             j=M*eBands[i]; do {
89                sum = MAC16_16(sum, EXTRACT16(VSHR32(X[j+c*N],shift)),
90                                    EXTRACT16(VSHR32(X[j+c*N],shift)));
91             } while (++j<M*eBands[i+1]);
92             /* We're adding one here to make damn sure we never end up with a pitch vector that's
93                larger than unity norm */
94             bank[i+c*m->nbEBands] = EPSILON+VSHR32(EXTEND32(celt_sqrt(sum)),-shift);
95          } else {
96             bank[i+c*m->nbEBands] = EPSILON;
97          }
98          /*printf ("%f ", bank[i+c*m->nbEBands]);*/
99       }
100    } while (++c<C);
101    /*printf ("\n");*/
102 }
103
104 /* Normalise each band such that the energy is one. */
105 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)
106 {
107    int i, c, N;
108    const celt_int16 *eBands = m->eBands;
109    const int C = CHANNELS(_C);
110    N = M*m->shortMdctSize;
111    c=0; do {
112       i=0; do {
113          celt_word16 g;
114          int j,shift;
115          celt_word16 E;
116          shift = celt_zlog2(bank[i+c*m->nbEBands])-13;
117          E = VSHR32(bank[i+c*m->nbEBands], shift);
118          g = EXTRACT16(celt_rcp(SHL32(E,3)));
119          j=M*eBands[i]; do {
120             X[j+c*N] = MULT16_16_Q15(VSHR32(freq[j+c*N],shift-1),g);
121          } while (++j<M*eBands[i+1]);
122       } while (++i<end);
123    } while (++c<C);
124 }
125
126 #else /* FIXED_POINT */
127 /* Compute the amplitude (sqrt energy) in each of the bands */
128 void compute_band_energies(const CELTMode *m, const celt_sig *X, celt_ener *bank, int end, int _C, int M)
129 {
130    int i, c, N;
131    const celt_int16 *eBands = m->eBands;
132    const int C = CHANNELS(_C);
133    N = M*m->shortMdctSize;
134    c=0; do {
135       for (i=0;i<end;i++)
136       {
137          int j;
138          celt_word32 sum = 1e-10f;
139          for (j=M*eBands[i];j<M*eBands[i+1];j++)
140             sum += X[j+c*N]*X[j+c*N];
141          bank[i+c*m->nbEBands] = celt_sqrt(sum);
142          /*printf ("%f ", bank[i+c*m->nbEBands]);*/
143       }
144    } while (++c<C);
145    /*printf ("\n");*/
146 }
147
148 /* Normalise each band such that the energy is one. */
149 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)
150 {
151    int i, c, N;
152    const celt_int16 *eBands = m->eBands;
153    const int C = CHANNELS(_C);
154    N = M*m->shortMdctSize;
155    c=0; do {
156       for (i=0;i<end;i++)
157       {
158          int j;
159          celt_word16 g = 1.f/(1e-10f+bank[i+c*m->nbEBands]);
160          for (j=M*eBands[i];j<M*eBands[i+1];j++)
161             X[j+c*N] = freq[j+c*N]*g;
162       }
163    } while (++c<C);
164 }
165
166 #endif /* FIXED_POINT */
167
168 /* De-normalise the energy to produce the synthesis from the unit-energy bands */
169 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)
170 {
171    int i, c, N;
172    const celt_int16 *eBands = m->eBands;
173    const int C = CHANNELS(_C);
174    N = M*m->shortMdctSize;
175    celt_assert2(C<=2, "denormalise_bands() not implemented for >2 channels");
176    c=0; do {
177       celt_sig * restrict f;
178       const celt_norm * restrict x;
179       f = freq+c*N;
180       x = X+c*N;
181       for (i=0;i<end;i++)
182       {
183          int j, band_end;
184          celt_word32 g = SHR32(bank[i+c*m->nbEBands],1);
185          j=M*eBands[i];
186          band_end = M*eBands[i+1];
187          do {
188             *f++ = SHL32(MULT16_32_Q15(*x, g),2);
189             x++;
190          } while (++j<band_end);
191       }
192       for (i=M*eBands[m->nbEBands];i<N;i++)
193          *f++ = 0;
194    } while (++c<C);
195 }
196
197 static void intensity_stereo(const CELTMode *m, celt_norm *X, celt_norm *Y, const celt_ener *bank, int bandID, int N)
198 {
199    int i = bandID;
200    int j;
201    celt_word16 a1, a2;
202    celt_word16 left, right;
203    celt_word16 norm;
204 #ifdef FIXED_POINT
205    int shift = celt_zlog2(MAX32(bank[i], bank[i+m->nbEBands]))-13;
206 #endif
207    left = VSHR32(bank[i],shift);
208    right = VSHR32(bank[i+m->nbEBands],shift);
209    norm = EPSILON + celt_sqrt(EPSILON+MULT16_16(left,left)+MULT16_16(right,right));
210    a1 = DIV32_16(SHL32(EXTEND32(left),14),norm);
211    a2 = DIV32_16(SHL32(EXTEND32(right),14),norm);
212    for (j=0;j<N;j++)
213    {
214       celt_norm r, l;
215       l = X[j];
216       r = Y[j];
217       X[j] = MULT16_16_Q14(a1,l) + MULT16_16_Q14(a2,r);
218       /* Side is not encoded, no need to calculate */
219    }
220 }
221
222 static void stereo_split(celt_norm *X, celt_norm *Y, int N)
223 {
224    int j;
225    for (j=0;j<N;j++)
226    {
227       celt_norm r, l;
228       l = MULT16_16_Q15(QCONST16(.70711f,15), X[j]);
229       r = MULT16_16_Q15(QCONST16(.70711f,15), Y[j]);
230       X[j] = l+r;
231       Y[j] = r-l;
232    }
233 }
234
235 static void stereo_merge(celt_norm *X, celt_norm *Y, celt_word16 mid, celt_word16 side, int N)
236 {
237    int j;
238    celt_word32 xp=0;
239    celt_word32 El, Er;
240 #ifdef FIXED_POINT
241    int kl, kr;
242 #endif
243    celt_word32 t, lgain, rgain;
244
245    /* Compute the norm of X+Y and X-Y as |X|^2 + |Y|^2 +/- sum(xy) */
246    for (j=0;j<N;j++)
247       xp = MAC16_16(xp, X[j], Y[j]);
248    /* mid and side are in Q15, not Q14 like X and Y */
249    mid = SHR32(mid, 1);
250    side = SHR32(side, 1);
251    El = MULT16_16(mid, mid) + MULT16_16(side, side) - 2*xp;
252    Er = MULT16_16(mid, mid) + MULT16_16(side, side) + 2*xp;
253    if (Er < EPSILON)
254       Er = EPSILON;
255    if (El < EPSILON)
256       El = EPSILON;
257
258 #ifdef FIXED_POINT
259    kl = celt_ilog2(El)>>1;
260    kr = celt_ilog2(Er)>>1;
261 #endif
262    t = VSHR32(El, (kl-7)<<1);
263    lgain = celt_rsqrt_norm(t);
264    t = VSHR32(Er, (kr-7)<<1);
265    rgain = celt_rsqrt_norm(t);
266
267 #ifdef FIXED_POINT
268    if (kl < 7)
269       kl = 7;
270    if (kr < 7)
271       kr = 7;
272 #endif
273
274    for (j=0;j<N;j++)
275    {
276       celt_norm r, l;
277       l = X[j];
278       r = Y[j];
279       X[j] = EXTRACT16(PSHR32(MULT16_16(lgain, SUB16(l,r)), kl+1));
280       Y[j] = EXTRACT16(PSHR32(MULT16_16(rgain, ADD16(l,r)), kr+1));
281    }
282 }
283
284 /* Decide whether we should spread the pulses in the current frame */
285 int folding_decision(const CELTMode *m, celt_norm *X, int *average, int *last_decision, int end, int _C, int M)
286 {
287    int i, c, N0;
288    int sum = 0, nbBands=0;
289    const int C = CHANNELS(_C);
290    const celt_int16 * restrict eBands = m->eBands;
291    int decision;
292    
293    N0 = M*m->shortMdctSize;
294
295    if (M*(eBands[end]-eBands[end-1]) <= 8)
296       return 0;
297    c=0; do {
298       for (i=0;i<end;i++)
299       {
300          int j, N, tmp=0;
301          int tcount[3] = {0};
302          celt_norm * restrict x = X+M*eBands[i]+c*N0;
303          N = M*(eBands[i+1]-eBands[i]);
304          if (N<=8)
305             continue;
306          /* Compute rough CDF of |x[j]| */
307          for (j=0;j<N;j++)
308          {
309             celt_word32 x2N; /* Q13 */
310
311             x2N = MULT16_16(MULT16_16_Q15(x[j], x[j]), N);
312             if (x2N < QCONST16(0.25f,13))
313                tcount[0]++;
314             if (x2N < QCONST16(0.0625f,13))
315                tcount[1]++;
316             if (x2N < QCONST16(0.015625f,13))
317                tcount[2]++;
318          }
319
320          tmp = (2*tcount[2] >= N) + (2*tcount[1] >= N) + (2*tcount[0] >= N);
321          sum += tmp*256;
322          nbBands++;
323       }
324    } while (++c<C);
325    sum /= nbBands;
326    /* Recursive averaging */
327    sum = (sum+*average)>>1;
328    *average = sum;
329    /* Hysteresis */
330    sum = (3*sum + ((*last_decision<<7) + 64) + 2)>>2;
331    /* decision and last_decision do not use the same encoding */
332    if (sum < 80)
333    {
334       decision = 2;
335       *last_decision = 0;
336    } else if (sum < 256)
337    {
338       decision = 1;
339       *last_decision = 1;
340    } else if (sum < 384)
341    {
342       decision = 3;
343       *last_decision = 2;
344    } else {
345       decision = 0;
346       *last_decision = 3;
347    }
348    return decision;
349 }
350
351 #ifdef MEASURE_NORM_MSE
352
353 float MSE[30] = {0};
354 int nbMSEBands = 0;
355 int MSECount[30] = {0};
356
357 void dump_norm_mse(void)
358 {
359    int i;
360    for (i=0;i<nbMSEBands;i++)
361    {
362       printf ("%g ", MSE[i]/MSECount[i]);
363    }
364    printf ("\n");
365 }
366
367 void measure_norm_mse(const CELTMode *m, float *X, float *X0, float *bandE, float *bandE0, int M, int N, int C)
368 {
369    static int init = 0;
370    int i;
371    if (!init)
372    {
373       atexit(dump_norm_mse);
374       init = 1;
375    }
376    for (i=0;i<m->nbEBands;i++)
377    {
378       int j;
379       int c;
380       float g;
381       if (bandE0[i]<10 || (C==2 && bandE0[i+m->nbEBands]<1))
382          continue;
383       c=0; do {
384          g = bandE[i+c*m->nbEBands]/(1e-15+bandE0[i+c*m->nbEBands]);
385          for (j=M*m->eBands[i];j<M*m->eBands[i+1];j++)
386             MSE[i] += (g*X[j+c*N]-X0[j+c*N])*(g*X[j+c*N]-X0[j+c*N]);
387       } while (++c<C);
388       MSECount[i]+=C;
389    }
390    nbMSEBands = m->nbEBands;
391 }
392
393 #endif
394
395 static void interleave_vector(celt_norm *X, int N0, int stride)
396 {
397    int i,j;
398    VARDECL(celt_norm, tmp);
399    int N;
400    SAVE_STACK;
401    N = N0*stride;
402    ALLOC(tmp, N, celt_norm);
403    for (i=0;i<stride;i++)
404       for (j=0;j<N0;j++)
405          tmp[j*stride+i] = X[i*N0+j];
406    for (j=0;j<N;j++)
407       X[j] = tmp[j];
408    RESTORE_STACK;
409 }
410
411 static void deinterleave_vector(celt_norm *X, int N0, int stride)
412 {
413    int i,j;
414    VARDECL(celt_norm, tmp);
415    int N;
416    SAVE_STACK;
417    N = N0*stride;
418    ALLOC(tmp, N, celt_norm);
419    for (i=0;i<stride;i++)
420       for (j=0;j<N0;j++)
421          tmp[i*N0+j] = X[j*stride+i];
422    for (j=0;j<N;j++)
423       X[j] = tmp[j];
424    RESTORE_STACK;
425 }
426
427 void haar1(celt_norm *X, int N0, int stride)
428 {
429    int i, j;
430    N0 >>= 1;
431    for (i=0;i<stride;i++)
432       for (j=0;j<N0;j++)
433       {
434          celt_norm tmp1, tmp2;
435          tmp1 = MULT16_16_Q15(QCONST16(.7070678f,15), X[stride*2*j+i]);
436          tmp2 = MULT16_16_Q15(QCONST16(.7070678f,15), X[stride*(2*j+1)+i]);
437          X[stride*2*j+i] = tmp1 + tmp2;
438          X[stride*(2*j+1)+i] = tmp1 - tmp2;
439       }
440 }
441
442 static int compute_qn(int N, int b, int offset, int stereo)
443 {
444    static const celt_int16 exp2_table8[8] =
445       {16384, 17866, 19483, 21247, 23170, 25267, 27554, 30048};
446    int qn, qb;
447    int N2 = 2*N-1;
448    if (stereo && N==2)
449       N2--;
450    qb = (b+N2*offset)/N2;
451    if (qb > (b>>1)-(1<<BITRES))
452       qb = (b>>1)-(1<<BITRES);
453
454    if (qb<0)
455        qb = 0;
456    if (qb>8<<BITRES)
457      qb = 8<<BITRES;
458
459    if (qb<(1<<BITRES>>1)) {
460       qn = 1;
461    } else {
462       qn = exp2_table8[qb&0x7]>>(14-(qb>>BITRES));
463       qn = (qn+1)>>1<<1;
464    }
465    celt_assert(qn <= 256);
466    return qn;
467 }
468
469 static celt_uint32 lcg_rand(celt_uint32 seed)
470 {
471    return 1664525 * seed + 1013904223;
472 }
473
474 /* This function is responsible for encoding and decoding a band for both
475    the mono and stereo case. Even in the mono case, it can split the band
476    in two and transmit the energy difference with the two half-bands. It
477    can be called recursively so bands can end up being split in 8 parts. */
478 static void quant_band(int encode, const CELTMode *m, int i, celt_norm *X, celt_norm *Y,
479       int N, int b, int spread, int B, int intensity, int tf_change, celt_norm *lowband, int resynth, void *ec,
480       celt_int32 *remaining_bits, int LM, celt_norm *lowband_out, const celt_ener *bandE, int level,
481       celt_int32 *seed, celt_word16 gain, celt_norm *lowband_scratch)
482 {
483    int q;
484    int curr_bits;
485    int stereo, split;
486    int imid=0, iside=0;
487    int N0=N;
488    int N_B=N;
489    int N_B0;
490    int B0=B;
491    int time_divide=0;
492    int recombine=0;
493    int inv = 0;
494    celt_word16 mid=0, side=0;
495
496    N_B /= B;
497    N_B0 = N_B;
498
499    split = stereo = Y != NULL;
500
501    /* Special case for one sample */
502    if (N==1)
503    {
504       int c;
505       celt_norm *x = X;
506       c=0; do {
507          int sign=0;
508          if (*remaining_bits>=1<<BITRES)
509          {
510             if (encode)
511             {
512                sign = x[0]<0;
513                ec_enc_bits((ec_enc*)ec, sign, 1);
514             } else {
515                sign = ec_dec_bits((ec_dec*)ec, 1);
516             }
517             *remaining_bits -= 1<<BITRES;
518             b-=1<<BITRES;
519          }
520          if (resynth)
521             x[0] = sign ? -NORM_SCALING : NORM_SCALING;
522          x = Y;
523       } while (++c<1+stereo);
524       if (lowband_out)
525          lowband_out[0] = SHR16(X[0],4);
526       return;
527    }
528
529    if (!stereo && level == 0)
530    {
531       int k;
532       if (tf_change>0)
533          recombine = tf_change;
534       /* Band recombining to increase frequency resolution */
535
536       if (lowband && (recombine || ((N_B&1) == 0 && tf_change<0) || B0>1))
537       {
538          int j;
539          for (j=0;j<N;j++)
540             lowband_scratch[j] = lowband[j];
541          lowband = lowband_scratch;
542       }
543
544       for (k=0;k<recombine;k++)
545       {
546          B>>=1;
547          N_B<<=1;
548          if (encode)
549             haar1(X, N_B, B);
550          if (lowband)
551             haar1(lowband, N_B, B);
552       }
553
554       /* Increasing the time resolution */
555       while ((N_B&1) == 0 && tf_change<0)
556       {
557          if (encode)
558             haar1(X, N_B, B);
559          if (lowband)
560             haar1(lowband, N_B, B);
561          B <<= 1;
562          N_B >>= 1;
563          time_divide++;
564          tf_change++;
565       }
566       B0=B;
567       N_B0 = N_B;
568
569       /* Reorganize the samples in time order instead of frequency order */
570       if (B0>1)
571       {
572          if (encode)
573             deinterleave_vector(X, N_B, B0);
574          if (lowband)
575             deinterleave_vector(lowband, N_B, B0);
576       }
577    }
578
579    /* If we need more than 32 bits, try splitting the band in two. */
580    if (!stereo && LM != -1 && b > 32<<BITRES && N>2)
581    {
582       if (LM>0 || (N&1)==0)
583       {
584          N >>= 1;
585          Y = X+N;
586          split = 1;
587          LM -= 1;
588          B = (B+1)>>1;
589       }
590    }
591
592    if (split)
593    {
594       int qn;
595       int itheta=0;
596       int mbits, sbits, delta;
597       int qalloc;
598       int offset;
599
600       /* Decide on the resolution to give to the split parameter theta */
601       offset = ((m->logN[i]+(LM<<BITRES))>>1) - (stereo ? QTHETA_OFFSET_STEREO : QTHETA_OFFSET);
602       qn = compute_qn(N, b, offset, stereo);
603       qalloc = 0;
604       if (stereo && i>=intensity)
605          qn = 1;
606       if (encode)
607       {
608          /* theta is the atan() of the ratio between the (normalized)
609             side and mid. With just that parameter, we can re-scale both
610             mid and side because we know that 1) they have unit norm and
611             2) they are orthogonal. */
612          itheta = stereo_itheta(X, Y, stereo, N);
613       }
614       if (qn!=1)
615       {
616          if (encode)
617             itheta = (itheta*qn+8192)>>14;
618
619          /* Entropy coding of the angle. We use a uniform pdf for the
620             first stereo split but a triangular one for the rest. */
621          if (stereo || B>1)
622          {
623             if (encode)
624                ec_enc_uint((ec_enc*)ec, itheta, qn+1);
625             else
626                itheta = ec_dec_uint((ec_dec*)ec, qn+1);
627             qalloc = log2_frac(qn+1,BITRES);
628          } else {
629             int fs=1, ft;
630             ft = ((qn>>1)+1)*((qn>>1)+1);
631             if (encode)
632             {
633                int fl;
634
635                fs = itheta <= (qn>>1) ? itheta + 1 : qn + 1 - itheta;
636                fl = itheta <= (qn>>1) ? itheta*(itheta + 1)>>1 :
637                 ft - ((qn + 1 - itheta)*(qn + 2 - itheta)>>1);
638
639                ec_encode((ec_enc*)ec, fl, fl+fs, ft);
640             } else {
641                int fl=0;
642                int fm;
643                fm = ec_decode((ec_dec*)ec, ft);
644
645                if (fm < ((qn>>1)*((qn>>1) + 1)>>1))
646                {
647                   itheta = (isqrt32(8*(celt_uint32)fm + 1) - 1)>>1;
648                   fs = itheta + 1;
649                   fl = itheta*(itheta + 1)>>1;
650                }
651                else
652                {
653                   itheta = (2*(qn + 1)
654                    - isqrt32(8*(celt_uint32)(ft - fm - 1) + 1))>>1;
655                   fs = qn + 1 - itheta;
656                   fl = ft - ((qn + 1 - itheta)*(qn + 2 - itheta)>>1);
657                }
658
659                ec_dec_update((ec_dec*)ec, fl, fl+fs, ft);
660             }
661             qalloc = log2_frac(ft,BITRES) - log2_frac(fs,BITRES) + 1;
662          }
663          itheta = (celt_int32)itheta*16384/qn;
664          if (encode && stereo)
665             stereo_split(X, Y, N);
666          /* TODO: Renormalising X and Y *may* help fixed-point a bit at very high rate.
667                   Let's do that at higher complexity */
668       } else if (stereo) {
669          if (encode)
670          {
671             inv = itheta > 8192;
672             if (inv)
673             {
674                int j;
675                for (j=0;j<N;j++)
676                   Y[j] = -Y[j];
677             }
678             intensity_stereo(m, X, Y, bandE, i, N);
679          }
680          if (b>2<<BITRES && *remaining_bits > 2<<BITRES)
681          {
682             if (encode)
683                ec_enc_bit_prob(ec, inv, 16384);
684             else
685                inv = ec_dec_bit_prob(ec, 16384);
686             qalloc = inv ? 16 : 4;
687          } else
688             inv = 0;
689          itheta = 0;
690       }
691
692       if (itheta == 0)
693       {
694          imid = 32767;
695          iside = 0;
696          delta = -10000;
697       } else if (itheta == 16384)
698       {
699          imid = 0;
700          iside = 32767;
701          delta = 10000;
702       } else {
703          imid = bitexact_cos(itheta);
704          iside = bitexact_cos(16384-itheta);
705          /* This is the mid vs side allocation that minimizes squared error
706             in that band. */
707          delta = (N-1)*(log2_frac(iside,BITRES+2)-log2_frac(imid,BITRES+2))>>2;
708       }
709
710 #ifdef FIXED_POINT
711       mid = imid;
712       side = iside;
713 #else
714       mid = (1.f/32768)*imid;
715       side = (1.f/32768)*iside;
716 #endif
717
718       /* This is a special case for N=2 that only works for stereo and takes
719          advantage of the fact that mid and side are orthogonal to encode
720          the side with just one bit. */
721       if (N==2 && stereo)
722       {
723          int c;
724          int sign=1;
725          celt_norm *x2, *y2;
726          mbits = b-qalloc;
727          sbits = 0;
728          /* Only need one bit for the side */
729          if (itheta != 0 && itheta != 16384)
730             sbits = 1<<BITRES;
731          mbits -= sbits;
732          c = itheta > 8192;
733          *remaining_bits -= qalloc+sbits;
734
735          x2 = c ? Y : X;
736          y2 = c ? X : Y;
737          if (sbits)
738          {
739             if (encode)
740             {
741                /* Here we only need to encode a sign for the side */
742                sign = x2[0]*y2[1] - x2[1]*y2[0] > 0;
743                ec_enc_bits((ec_enc*)ec, sign, 1);
744             } else {
745                sign = ec_dec_bits((ec_dec*)ec, 1);
746             }
747          }
748          sign = 2*sign - 1;
749          quant_band(encode, m, i, x2, NULL, N, mbits, spread, B, intensity, tf_change, lowband, resynth, ec, remaining_bits, LM, lowband_out, NULL, level, seed, gain, lowband_scratch);
750          y2[0] = -sign*x2[1];
751          y2[1] = sign*x2[0];
752          if (resynth)
753          {
754             celt_norm tmp;
755             X[0] = MULT16_16_Q15(mid, X[0]);
756             X[1] = MULT16_16_Q15(mid, X[1]);
757             Y[0] = MULT16_16_Q15(side, Y[0]);
758             Y[1] = MULT16_16_Q15(side, Y[1]);
759             tmp = X[0];
760             X[0] = SUB16(tmp,Y[0]);
761             Y[0] = ADD16(tmp,Y[0]);
762             tmp = X[1];
763             X[1] = SUB16(tmp,Y[1]);
764             Y[1] = ADD16(tmp,Y[1]);
765          }
766       } else {
767          /* "Normal" split code */
768          celt_norm *next_lowband2=NULL;
769          celt_norm *next_lowband_out1=NULL;
770          int next_level=0;
771
772          /* Give more bits to low-energy MDCTs than they would otherwise deserve */
773          if (B>1 && !stereo && itheta > 8192)
774             delta -= delta>>(1+level);
775
776          mbits = (b-qalloc-delta)/2;
777          if (mbits > b-qalloc)
778             mbits = b-qalloc;
779          if (mbits<0)
780             mbits=0;
781          sbits = b-qalloc-mbits;
782          *remaining_bits -= qalloc;
783
784          if (lowband && !stereo)
785             next_lowband2 = lowband+N; /* >32-bit split case */
786
787          /* Only stereo needs to pass on lowband_out. Otherwise, it's
788             handled at the end */
789          if (stereo)
790             next_lowband_out1 = lowband_out;
791          else
792             next_level = level+1;
793
794          quant_band(encode, m, i, X, NULL, N, mbits, spread, B, intensity, tf_change,
795                lowband, resynth, ec, remaining_bits, LM, next_lowband_out1,
796                NULL, next_level, seed, MULT16_16_P15(gain,mid), lowband_scratch);
797          quant_band(encode, m, i, Y, NULL, N, sbits, spread, B, intensity, tf_change,
798                next_lowband2, resynth, ec, remaining_bits, LM, NULL,
799                NULL, next_level, seed, MULT16_16_P15(gain,side), NULL);
800       }
801
802    } else {
803       /* This is the basic no-split case */
804       q = bits2pulses(m, i, LM, b);
805       curr_bits = pulses2bits(m, i, LM, q);
806       *remaining_bits -= curr_bits;
807
808       /* Ensures we can never bust the budget */
809       while (*remaining_bits < 0 && q > 0)
810       {
811          *remaining_bits += curr_bits;
812          q--;
813          curr_bits = pulses2bits(m, i, LM, q);
814          *remaining_bits -= curr_bits;
815       }
816
817       if (q!=0)
818       {
819          int K = get_pulses(q);
820
821          /* Finally do the actual quantization */
822          if (encode)
823             alg_quant(X, N, K, spread, B, lowband, resynth, (ec_enc*)ec, seed, gain);
824          else
825             alg_unquant(X, N, K, spread, B, lowband, (ec_dec*)ec, seed, gain);
826       } else {
827          /* If there's no pulse, fill the band anyway */
828          int j;
829          if (lowband != NULL && resynth)
830          {
831             if (spread==2 && B<=1)
832             {
833                /* Noise */
834                for (j=0;j<N;j++)
835                {
836                   *seed = lcg_rand(*seed);
837                   X[j] = (int)(*seed)>>20;
838                }
839             } else {
840                /* Folded spectrum */
841                for (j=0;j<N;j++)
842                   X[j] = lowband[j];
843             }
844             renormalise_vector(X, N, gain);
845          } else {
846             /* This is important for encoding the side in stereo mode */
847             for (j=0;j<N;j++)
848                X[j] = 0;
849          }
850       }
851    }
852
853    /* This code is used by the decoder and by the resynthesis-enabled encoder */
854    if (resynth)
855    {
856       if (stereo)
857       {
858          if (N!=2)
859             stereo_merge(X, Y, mid, side, N);
860          if (inv)
861          {
862             int j;
863             for (j=0;j<N;j++)
864                Y[j] = -Y[j];
865          }
866       } else if (level == 0)
867       {
868          int k;
869
870          /* Undo the sample reorganization going from time order to frequency order */
871          if (B0>1)
872             interleave_vector(X, N_B, B0);
873
874          /* Undo time-freq changes that we did earlier */
875          N_B = N_B0;
876          B = B0;
877          for (k=0;k<time_divide;k++)
878          {
879             B >>= 1;
880             N_B <<= 1;
881             haar1(X, N_B, B);
882          }
883
884          for (k=0;k<recombine;k++)
885          {
886             haar1(X, N_B, B);
887             N_B>>=1;
888             B <<= 1;
889          }
890
891          /* Scale output for later folding */
892          if (lowband_out)
893          {
894             int j;
895             celt_word16 n;
896             n = celt_sqrt(SHL32(EXTEND32(N0),22));
897             for (j=0;j<N0;j++)
898                lowband_out[j] = MULT16_16_Q15(n,X[j]);
899          }
900       }
901    }
902 }
903
904 void quant_all_bands(int encode, const CELTMode *m, int start, int end,
905       celt_norm *_X, celt_norm *_Y, const celt_ener *bandE, int *pulses,
906       int shortBlocks, int fold, int dual_stereo, int intensity, int *tf_res, int resynth,
907       int total_bits, void *ec, int LM, int codedBands)
908 {
909    int i, balance;
910    celt_int32 remaining_bits;
911    const celt_int16 * restrict eBands = m->eBands;
912    celt_norm * restrict norm, * restrict norm2;
913    VARDECL(celt_norm, _norm);
914    VARDECL(celt_norm, lowband_scratch);
915    int B;
916    int M;
917    celt_int32 seed;
918    int lowband_offset;
919    int update_lowband = 1;
920    int C = _Y != NULL ? 2 : 1;
921    SAVE_STACK;
922
923    M = 1<<LM;
924    B = shortBlocks ? M : 1;
925    ALLOC(_norm, C*M*eBands[m->nbEBands], celt_norm);
926    ALLOC(lowband_scratch, M*(eBands[m->nbEBands]-eBands[m->nbEBands-1]), celt_norm);
927    norm = _norm;
928    norm2 = norm + M*eBands[m->nbEBands];
929 #if 0
930    if (C==2)
931    {
932       int j;
933       int left = 0;
934       for (j=intensity;j<codedBands;j++)
935       {
936          int tmp = pulses[j]/2;
937          left += tmp;
938          pulses[j] -= tmp;
939       }
940       if (codedBands) {
941          int perband;
942          perband = left/(m->eBands[codedBands]-m->eBands[start]);
943          for (j=start;j<codedBands;j++)
944             pulses[j] += perband*(m->eBands[j+1]-m->eBands[j]);
945          left = left-(m->eBands[codedBands]-m->eBands[start])*perband;
946          for (j=start;j<codedBands;j++)
947          {
948             int tmp = IMIN(left, m->eBands[j+1]-m->eBands[j]);
949             pulses[j] += tmp;
950             left -= tmp;
951          }
952       }
953    }
954 #endif
955    if (encode)
956       seed = ((ec_enc*)ec)->rng;
957    else
958       seed = ((ec_dec*)ec)->rng;
959    balance = 0;
960    lowband_offset = -1;
961    for (i=start;i<end;i++)
962    {
963       int tell;
964       int b;
965       int N;
966       int curr_balance;
967       int effective_lowband=-1;
968       celt_norm * restrict X, * restrict Y;
969       int tf_change=0;
970       
971       X = _X+M*eBands[i];
972       if (_Y!=NULL)
973          Y = _Y+M*eBands[i];
974       else
975          Y = NULL;
976       N = M*eBands[i+1]-M*eBands[i];
977       if (encode)
978          tell = ec_enc_tell((ec_enc*)ec, BITRES);
979       else
980          tell = ec_dec_tell((ec_dec*)ec, BITRES);
981
982       /* Compute how many bits we want to allocate to this band */
983       if (i != start)
984          balance -= tell;
985       remaining_bits = (total_bits<<BITRES)-tell-1;
986       if (i <= codedBands-1)
987       {
988          curr_balance = (codedBands-i);
989          if (curr_balance > 3)
990             curr_balance = 3;
991          curr_balance = balance / curr_balance;
992          b = IMIN(remaining_bits+1,pulses[i]+curr_balance);
993          if (b<0)
994             b = 0;
995       } else {
996          b = 0;
997       }
998       /* Prevents ridiculous bit depths */
999       if (b > C*16*N<<BITRES)
1000          b = C*16*N<<BITRES;
1001
1002       if (M*eBands[i]-N >= M*eBands[start] && (update_lowband || lowband_offset==-1))
1003             lowband_offset = M*eBands[i];
1004
1005       tf_change = tf_res[i];
1006       if (i>=m->effEBands)
1007       {
1008          X=norm;
1009          if (_Y!=NULL)
1010             Y = norm;
1011       }
1012
1013       /* This ensures we never repeat spectral content within one band */
1014       if (lowband_offset != -1)
1015       {
1016          effective_lowband = lowband_offset-N;
1017          if (effective_lowband < M*eBands[start])
1018             effective_lowband = M*eBands[start];
1019       }
1020       if (dual_stereo && i==intensity)
1021       {
1022          int j;
1023
1024          /* Switch off dual stereo to do intensity */
1025          dual_stereo = 0;
1026          for (j=0;j<M*eBands[i];j++)
1027             norm[j] = HALF32(norm[j]+norm2[j]);
1028       }
1029       if (dual_stereo)
1030       {
1031          quant_band(encode, m, i, X, NULL, N, b/2, fold, B, intensity, tf_change,
1032                effective_lowband != -1 ? norm+effective_lowband : NULL, resynth, ec, &remaining_bits, LM,
1033                norm+M*eBands[i], bandE, 0, &seed, Q15ONE, lowband_scratch);
1034          quant_band(encode, m, i, Y, NULL, N, b/2, fold, B, intensity, tf_change,
1035                effective_lowband != -1 ? norm2+effective_lowband : NULL, resynth, ec, &remaining_bits, LM,
1036                norm2+M*eBands[i], bandE, 0, &seed, Q15ONE, lowband_scratch);
1037       } else {
1038          quant_band(encode, m, i, X, Y, N, b, fold, B, intensity, tf_change,
1039                effective_lowband != -1 ? norm+effective_lowband : NULL, resynth, ec, &remaining_bits, LM,
1040                norm+M*eBands[i], bandE, 0, &seed, Q15ONE, lowband_scratch);
1041       }
1042       balance += pulses[i] + tell;
1043
1044       /* Update the folding position only as long as we have 1 bit/sample depth */
1045       update_lowband = (b>>BITRES)>N;
1046    }
1047    RESTORE_STACK;
1048 }
1049