Changed all the celt*int*_t types to remove the _t suffix, which is reserved
[opus.git] / libcelt / bands.c
1 /* (C) 2007-2008 Jean-Marc Valin, CSIRO
2    (C) 2008-2009 Gregory Maxwell */
3 /*
4    Redistribution and use in source and binary forms, with or without
5    modification, are permitted provided that the following conditions
6    are met:
7    
8    - Redistributions of source code must retain the above copyright
9    notice, this list of conditions and the following disclaimer.
10    
11    - Redistributions in binary form must reproduce the above copyright
12    notice, this list of conditions and the following disclaimer in the
13    documentation and/or other materials provided with the distribution.
14    
15    - Neither the name of the Xiph.org Foundation nor the names of its
16    contributors may be used to endorse or promote products derived from
17    this software without specific prior written permission.
18    
19    THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
20    ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
21    LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
22    A PARTICULAR PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE FOUNDATION OR
23    CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
24    EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
25    PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
26    PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
27    LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
28    NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
29    SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
30 */
31
32 #ifdef HAVE_CONFIG_H
33 #include "config.h"
34 #endif
35
36 #include <math.h>
37 #include "bands.h"
38 #include "modes.h"
39 #include "vq.h"
40 #include "cwrs.h"
41 #include "stack_alloc.h"
42 #include "os_support.h"
43 #include "mathops.h"
44 #include "rate.h"
45
46
47 #ifdef FIXED_POINT
48 /* Compute the amplitude (sqrt energy) in each of the bands */
49 void compute_band_energies(const CELTMode *m, const celt_sig_t *X, celt_ener_t *bank, int _C)
50 {
51    int i, c, N;
52    const celt_int16 *eBands = m->eBands;
53    const int C = CHANNELS(_C);
54    N = FRAMESIZE(m);
55    for (c=0;c<C;c++)
56    {
57       for (i=0;i<m->nbEBands;i++)
58       {
59          int j;
60          celt_word32_t maxval=0;
61          celt_word32_t sum = 0;
62          
63          j=eBands[i]; do {
64             maxval = MAX32(maxval, X[j+c*N]);
65             maxval = MAX32(maxval, -X[j+c*N]);
66          } while (++j<eBands[i+1]);
67          
68          if (maxval > 0)
69          {
70             int shift = celt_ilog2(maxval)-10;
71             j=eBands[i]; do {
72                sum = MAC16_16(sum, EXTRACT16(VSHR32(X[j+c*N],shift)),
73                                    EXTRACT16(VSHR32(X[j+c*N],shift)));
74             } while (++j<eBands[i+1]);
75             /* We're adding one here to make damn sure we never end up with a pitch vector that's
76                larger than unity norm */
77             bank[i+c*m->nbEBands] = EPSILON+VSHR32(EXTEND32(celt_sqrt(sum)),-shift);
78          } else {
79             bank[i+c*m->nbEBands] = EPSILON;
80          }
81          /*printf ("%f ", bank[i+c*m->nbEBands]);*/
82       }
83    }
84    /*printf ("\n");*/
85 }
86
87 /* Normalise each band such that the energy is one. */
88 void normalise_bands(const CELTMode *m, const celt_sig_t * restrict freq, celt_norm_t * restrict X, const celt_ener_t *bank, int _C)
89 {
90    int i, c, N;
91    const celt_int16 *eBands = m->eBands;
92    const int C = CHANNELS(_C);
93    N = FRAMESIZE(m);
94    for (c=0;c<C;c++)
95    {
96       i=0; do {
97          celt_word16_t g;
98          int j,shift;
99          celt_word16_t E;
100          shift = celt_zlog2(bank[i+c*m->nbEBands])-13;
101          E = VSHR32(bank[i+c*m->nbEBands], shift);
102          g = EXTRACT16(celt_rcp(SHL32(E,3)));
103          j=eBands[i]; do {
104             X[j+c*N] = MULT16_16_Q15(VSHR32(freq[j+c*N],shift-1),g);
105          } while (++j<eBands[i+1]);
106       } while (++i<m->nbEBands);
107    }
108 }
109
110 #else /* FIXED_POINT */
111 /* Compute the amplitude (sqrt energy) in each of the bands */
112 void compute_band_energies(const CELTMode *m, const celt_sig_t *X, celt_ener_t *bank, int _C)
113 {
114    int i, c, N;
115    const celt_int16 *eBands = m->eBands;
116    const int C = CHANNELS(_C);
117    N = FRAMESIZE(m);
118    for (c=0;c<C;c++)
119    {
120       for (i=0;i<m->nbEBands;i++)
121       {
122          int j;
123          celt_word32_t sum = 1e-10;
124          for (j=eBands[i];j<eBands[i+1];j++)
125             sum += X[j+c*N]*X[j+c*N];
126          bank[i+c*m->nbEBands] = sqrt(sum);
127          /*printf ("%f ", bank[i+c*m->nbEBands]);*/
128       }
129    }
130    /*printf ("\n");*/
131 }
132
133 #ifdef EXP_PSY
134 void compute_noise_energies(const CELTMode *m, const celt_sig_t *X, const celt_word16_t *tonality, celt_ener_t *bank, int _C)
135 {
136    int i, c, N;
137    const celt_int16 *eBands = m->eBands;
138    const int C = CHANNELS(_C);
139    N = FRAMESIZE(m);
140    for (c=0;c<C;c++)
141    {
142       for (i=0;i<m->nbEBands;i++)
143       {
144          int j;
145          celt_word32_t sum = 1e-10;
146          for (j=eBands[i];j<eBands[i+1];j++)
147             sum += X[j*C+c]*X[j+c*N]*tonality[j];
148          bank[i+c*m->nbEBands] = sqrt(sum);
149          /*printf ("%f ", bank[i+c*m->nbEBands]);*/
150       }
151    }
152    /*printf ("\n");*/
153 }
154 #endif
155
156 /* Normalise each band such that the energy is one. */
157 void normalise_bands(const CELTMode *m, const celt_sig_t * restrict freq, celt_norm_t * restrict X, const celt_ener_t *bank, int _C)
158 {
159    int i, c, N;
160    const celt_int16 *eBands = m->eBands;
161    const int C = CHANNELS(_C);
162    N = FRAMESIZE(m);
163    for (c=0;c<C;c++)
164    {
165       for (i=0;i<m->nbEBands;i++)
166       {
167          int j;
168          celt_word16_t g = 1.f/(1e-10+bank[i+c*m->nbEBands]);
169          for (j=eBands[i];j<eBands[i+1];j++)
170             X[j+c*N] = freq[j+c*N]*g;
171       }
172    }
173 }
174
175 #endif /* FIXED_POINT */
176
177 void renormalise_bands(const CELTMode *m, celt_norm_t * restrict X, int _C)
178 {
179    int i, c;
180    const celt_int16 *eBands = m->eBands;
181    const int C = CHANNELS(_C);
182    for (c=0;c<C;c++)
183    {
184       i=0; do {
185          renormalise_vector(X+eBands[i]+c*m->eBands[m->nbEBands+1], Q15ONE, eBands[i+1]-eBands[i], 1);
186       } while (++i<m->nbEBands);
187    }
188 }
189
190 /* De-normalise the energy to produce the synthesis from the unit-energy bands */
191 void denormalise_bands(const CELTMode *m, const celt_norm_t * restrict X, celt_sig_t * restrict freq, const celt_ener_t *bank, int _C)
192 {
193    int i, c, N;
194    const celt_int16 *eBands = m->eBands;
195    const int C = CHANNELS(_C);
196    N = FRAMESIZE(m);
197    if (C>2)
198       celt_fatal("denormalise_bands() not implemented for >2 channels");
199    for (c=0;c<C;c++)
200    {
201       for (i=0;i<m->nbEBands;i++)
202       {
203          int j;
204          celt_word32_t g = SHR32(bank[i+c*m->nbEBands],1);
205          j=eBands[i]; do {
206             freq[j+c*N] = SHL32(MULT16_32_Q15(X[j+c*N], g),2);
207          } while (++j<eBands[i+1]);
208       }
209       for (i=eBands[m->nbEBands];i<eBands[m->nbEBands+1];i++)
210          freq[i+c*N] = 0;
211    }
212 }
213
214 int compute_pitch_gain(const CELTMode *m, const celt_sig_t *X, const celt_sig_t *P, int norm_rate, int *gain_id, int _C, celt_word16_t *gain_prod)
215 {
216    int j, c;
217    celt_word16_t g;
218    celt_word16_t delta;
219    const int C = CHANNELS(_C);
220    celt_word32_t Sxy=0, Sxx=0, Syy=0;
221    int len = m->pitchEnd;
222    const int N = FRAMESIZE(m);
223 #ifdef FIXED_POINT
224    int shift = 0;
225    celt_word32_t maxabs=0;
226
227    for (c=0;c<C;c++)
228    {
229       for (j=0;j<len;j++)
230       {
231          maxabs = MAX32(maxabs, ABS32(X[j+c*N]));
232          maxabs = MAX32(maxabs, ABS32(P[j+c*N]));
233       }
234    }
235    shift = celt_ilog2(maxabs)-12;
236    if (shift<0)
237       shift = 0;
238 #endif
239    delta = PDIV32_16(Q15ONE, len);
240    for (c=0;c<C;c++)
241    {
242       celt_word16_t gg = Q15ONE;
243       for (j=0;j<len;j++)
244       {
245          celt_word16_t Xj, Pj;
246          Xj = EXTRACT16(SHR32(X[j+c*N], shift));
247          Pj = MULT16_16_P15(gg,EXTRACT16(SHR32(P[j+c*N], shift)));
248          Sxy = MAC16_16(Sxy, Xj, Pj);
249          Sxx = MAC16_16(Sxx, Pj, Pj);
250          Syy = MAC16_16(Syy, Xj, Xj);
251          gg = SUB16(gg, delta);
252       }
253    }
254 #ifdef FIXED_POINT
255    {
256       celt_word32_t num, den;
257       celt_word16_t fact;
258       fact = MULT16_16(QCONST16(.04, 14), norm_rate);
259       if (fact < QCONST16(1., 14))
260          fact = QCONST16(1., 14);
261       num = Sxy;
262       den = EPSILON+Sxx+MULT16_32_Q15(QCONST16(.03,15),Syy);
263       shift = celt_ilog2(Sxy)-16;
264       if (shift < 0)
265          shift = 0;
266       g = DIV32(SHL32(SHR32(num,shift),14),SHR32(den,shift));
267       if (Sxy < MULT16_32_Q15(fact, MULT16_16(celt_sqrt(EPSILON+Sxx),celt_sqrt(EPSILON+Syy))))
268          g = 0;
269       /* This MUST round down so that we don't over-estimate the gain */
270       *gain_id = EXTRACT16(SHR32(MULT16_16(20,(g-QCONST16(.5,14))),14));
271    }
272 #else
273    {
274       float fact = .04*norm_rate;
275       if (fact < 1)
276          fact = 1;
277       g = Sxy/(.1+Sxx+.03*Syy);
278       if (Sxy < .5*fact*celt_sqrt(1+Sxx*Syy))
279          g = 0;
280       /* This MUST round down so that we don't over-estimate the gain */
281       *gain_id = floor(20*(g-.5));
282    }
283 #endif
284    /* This prevents the pitch gain from being above 1.0 for too long by bounding the 
285       maximum error amplification factor to 2.0 */
286    g = ADD16(QCONST16(.5,14), MULT16_16_16(QCONST16(.05,14),*gain_id));
287    *gain_prod = MAX16(QCONST32(1., 13), MULT16_16_Q14(*gain_prod,g));
288    if (*gain_prod>QCONST32(2., 13))
289    {
290       *gain_id=9;
291       *gain_prod = QCONST32(2., 13);
292    }
293
294    if (*gain_id < 0)
295    {
296       *gain_id = 0;
297       return 0;
298    } else {
299       if (*gain_id > 15)
300          *gain_id = 15;
301       return 1;
302    }
303 }
304
305 void apply_pitch(const CELTMode *m, celt_sig_t *X, const celt_sig_t *P, int gain_id, int pred, int _C)
306 {
307    int j, c, N;
308    celt_word16_t gain;
309    celt_word16_t delta;
310    const int C = CHANNELS(_C);
311    int len = m->pitchEnd;
312
313    N = FRAMESIZE(m);
314    gain = ADD16(QCONST16(.5,14), MULT16_16_16(QCONST16(.05,14),gain_id));
315    delta = PDIV32_16(gain, len);
316    if (pred)
317       gain = -gain;
318    else
319       delta = -delta;
320    for (c=0;c<C;c++)
321    {
322       celt_word16_t gg = gain;
323       for (j=0;j<len;j++)
324       {
325          X[j+c*N] += SHL32(MULT16_32_Q15(gg,P[j+c*N]),1);
326          gg = ADD16(gg, delta);
327       }
328    }
329 }
330
331 #ifndef DISABLE_STEREO
332
333 static void stereo_band_mix(const CELTMode *m, celt_norm_t *X, celt_norm_t *Y, const celt_ener_t *bank, int stereo_mode, int bandID, int dir)
334 {
335    int i = bandID;
336    const celt_int16 *eBands = m->eBands;
337    int j;
338    celt_word16_t a1, a2;
339    if (stereo_mode==0)
340    {
341       /* Do mid-side when not doing intensity stereo */
342       a1 = QCONST16(.70711f,14);
343       a2 = dir*QCONST16(.70711f,14);
344    } else {
345       celt_word16_t left, right;
346       celt_word16_t norm;
347 #ifdef FIXED_POINT
348       int shift = celt_zlog2(MAX32(bank[i], bank[i+m->nbEBands]))-13;
349 #endif
350       left = VSHR32(bank[i],shift);
351       right = VSHR32(bank[i+m->nbEBands],shift);
352       norm = EPSILON + celt_sqrt(EPSILON+MULT16_16(left,left)+MULT16_16(right,right));
353       a1 = DIV32_16(SHL32(EXTEND32(left),14),norm);
354       a2 = dir*DIV32_16(SHL32(EXTEND32(right),14),norm);
355    }
356    for (j=0;j<eBands[i+1]-eBands[i];j++)
357    {
358       celt_norm_t r, l;
359       l = X[j];
360       r = Y[j];
361       X[j] = MULT16_16_Q14(a1,l) + MULT16_16_Q14(a2,r);
362       Y[j] = MULT16_16_Q14(a1,r) - MULT16_16_Q14(a2,l);
363    }
364 }
365
366
367 #endif /* DISABLE_STEREO */
368
369 int folding_decision(const CELTMode *m, celt_norm_t *X, celt_word16_t *average, int *last_decision, int _C)
370 {
371    int i, c, N0;
372    int NR=0;
373    celt_word32_t ratio = EPSILON;
374    const int C = CHANNELS(_C);
375    const celt_int16 * restrict eBands = m->eBands;
376    
377    N0 = FRAMESIZE(m);
378
379    for (c=0;c<C;c++)
380    {
381    for (i=0;i<m->nbEBands;i++)
382    {
383       int j, N;
384       int max_i=0;
385       celt_word16_t max_val=EPSILON;
386       celt_word32_t floor_ener=EPSILON;
387       celt_norm_t * restrict x = X+eBands[i]+c*N0;
388       N = eBands[i+1]-eBands[i];
389       for (j=0;j<N;j++)
390       {
391          if (ABS16(x[j])>max_val)
392          {
393             max_val = ABS16(x[j]);
394             max_i = j;
395          }
396       }
397 #if 0
398       for (j=0;j<N;j++)
399       {
400          if (abs(j-max_i)>2)
401             floor_ener += x[j]*x[j];
402       }
403 #else
404       floor_ener = QCONST32(1.,28)-MULT16_16(max_val,max_val);
405       if (max_i < N-1)
406          floor_ener -= MULT16_16(x[(max_i+1)], x[(max_i+1)]);
407       if (max_i < N-2)
408          floor_ener -= MULT16_16(x[(max_i+2)], x[(max_i+2)]);
409       if (max_i > 0)
410          floor_ener -= MULT16_16(x[(max_i-1)], x[(max_i-1)]);
411       if (max_i > 1)
412          floor_ener -= MULT16_16(x[(max_i-2)], x[(max_i-2)]);
413       floor_ener = MAX32(floor_ener, EPSILON);
414 #endif
415       if (N>7)
416       {
417          celt_word16_t r;
418          celt_word16_t den = celt_sqrt(floor_ener);
419          den = MAX32(QCONST16(.02, 15), den);
420          r = DIV32_16(SHL32(EXTEND32(max_val),8),den);
421          ratio = ADD32(ratio, EXTEND32(r));
422          NR++;
423       }
424    }
425    }
426    if (NR>0)
427       ratio = DIV32_16(ratio, NR);
428    ratio = ADD32(HALF32(ratio), HALF32(*average));
429    if (!*last_decision)
430    {
431       *last_decision = (ratio < QCONST16(1.8,8));
432    } else {
433       *last_decision = (ratio < QCONST16(3.,8));
434    }
435    *average = EXTRACT16(ratio);
436    return *last_decision;
437 }
438
439 /* Quantisation of the residual */
440 void quant_bands(const CELTMode *m, celt_norm_t * restrict X, const celt_ener_t *bandE, int *pulses, int shortBlocks, int fold, int total_bits, int encode, void *enc_dec)
441 {
442    int i, j, remaining_bits, balance;
443    const celt_int16 * restrict eBands = m->eBands;
444    celt_norm_t * restrict norm;
445    VARDECL(celt_norm_t, _norm);
446    int B;
447    SAVE_STACK;
448
449    B = shortBlocks ? m->nbShortMdcts : 1;
450    ALLOC(_norm, eBands[m->nbEBands+1], celt_norm_t);
451    norm = _norm;
452
453    balance = 0;
454    for (i=0;i<m->nbEBands;i++)
455    {
456       int tell;
457       int N;
458       int q;
459       celt_word16_t n;
460       const celt_int16 * const *BPbits;
461       
462       int curr_balance, curr_bits;
463       
464       N = eBands[i+1]-eBands[i];
465       BPbits = m->bits;
466
467       if (encode)
468          tell = ec_enc_tell(enc_dec, BITRES);
469       else
470          tell = ec_dec_tell(enc_dec, BITRES);
471       if (i != 0)
472          balance -= tell;
473       remaining_bits = (total_bits<<BITRES)-tell-1;
474       curr_balance = (m->nbEBands-i);
475       if (curr_balance > 3)
476          curr_balance = 3;
477       curr_balance = balance / curr_balance;
478       q = bits2pulses(m, BPbits[i], N, pulses[i]+curr_balance);
479       curr_bits = pulses2bits(BPbits[i], N, q);
480       remaining_bits -= curr_bits;
481       while (remaining_bits < 0 && q > 0)
482       {
483          remaining_bits += curr_bits;
484          q--;
485          curr_bits = pulses2bits(BPbits[i], N, q);
486          remaining_bits -= curr_bits;
487       }
488       balance += pulses[i] + tell;
489       
490       n = SHL16(celt_sqrt(eBands[i+1]-eBands[i]),11);
491
492       if (q > 0)
493       {
494          int spread = fold ? B : 0;
495          if (encode)
496             alg_quant(X+eBands[i], eBands[i+1]-eBands[i], q, spread, enc_dec);
497          else
498             alg_unquant(X+eBands[i], eBands[i+1]-eBands[i], q, spread, enc_dec);
499       } else {
500          intra_fold(m, eBands[i+1]-eBands[i], norm, X+eBands[i], eBands[i], B);
501       }
502       for (j=eBands[i];j<eBands[i+1];j++)
503          norm[j] = MULT16_16_Q15(n,X[j]);
504    }
505    RESTORE_STACK;
506 }
507
508 #ifndef DISABLE_STEREO
509
510 void quant_bands_stereo(const CELTMode *m, celt_norm_t *_X, const celt_ener_t *bandE, int *pulses, int shortBlocks, int fold, int total_bits, ec_enc *enc)
511 {
512    int i, j, remaining_bits, balance;
513    const celt_int16 * restrict eBands = m->eBands;
514    celt_norm_t * restrict norm;
515    VARDECL(celt_norm_t, _norm);
516    int B;
517    celt_word16_t mid, side;
518    SAVE_STACK;
519
520    B = shortBlocks ? m->nbShortMdcts : 1;
521    ALLOC(_norm, eBands[m->nbEBands+1], celt_norm_t);
522    norm = _norm;
523
524    balance = 0;
525    for (i=0;i<m->nbEBands;i++)
526    {
527       int c;
528       int tell;
529       int q1, q2;
530       celt_word16_t n;
531       const celt_int16 * const *BPbits;
532       int b, qb;
533       int N;
534       int curr_balance, curr_bits;
535       int imid, iside, itheta;
536       int mbits, sbits, delta;
537       int qalloc;
538       celt_norm_t * restrict X, * restrict Y;
539       
540       X = _X+eBands[i];
541       Y = X+eBands[m->nbEBands+1];
542       BPbits = m->bits;
543
544       N = eBands[i+1]-eBands[i];
545       tell = ec_enc_tell(enc, BITRES);
546       if (i != 0)
547          balance -= tell;
548       remaining_bits = (total_bits<<BITRES)-tell-1;
549       curr_balance = (m->nbEBands-i);
550       if (curr_balance > 3)
551          curr_balance = 3;
552       curr_balance = balance / curr_balance;
553       b = IMIN(remaining_bits+1,pulses[i]+curr_balance);
554       if (b<0)
555          b = 0;
556
557       qb = (b-2*(N-1)*(QTHETA_OFFSET-log2_frac(N,BITRES)))/(32*(N-1));
558       if (qb > (b>>BITRES)-1)
559          qb = (b>>BITRES)-1;
560       if (qb<0)
561          qb = 0;
562       if (qb>14)
563          qb = 14;
564
565       stereo_band_mix(m, X, Y, bandE, qb==0, i, 1);
566
567       mid = renormalise_vector(X, Q15ONE, N, 1);
568       side = renormalise_vector(Y, Q15ONE, N, 1);
569 #ifdef FIXED_POINT
570       itheta = MULT16_16_Q15(QCONST16(0.63662,15),celt_atan2p(side, mid));
571 #else
572       itheta = floor(.5+16384*0.63662*atan2(side,mid));
573 #endif
574       qalloc = log2_frac((1<<qb)+1,BITRES);
575       if (qb==0)
576       {
577          itheta=0;
578       } else {
579          int shift;
580          shift = 14-qb;
581          itheta = (itheta+(1<<shift>>1))>>shift;
582          ec_enc_uint(enc, itheta, (1<<qb)+1);
583          itheta <<= shift;
584       }
585       if (itheta == 0)
586       {
587          imid = 32767;
588          iside = 0;
589          delta = -10000;
590       } else if (itheta == 16384)
591       {
592          imid = 0;
593          iside = 32767;
594          delta = 10000;
595       } else {
596          imid = bitexact_cos(itheta);
597          iside = bitexact_cos(16384-itheta);
598          delta = (N-1)*(log2_frac(iside,BITRES+2)-log2_frac(imid,BITRES+2))>>2;
599       }
600       n = SHL16(celt_sqrt((eBands[i+1]-eBands[i])),11);
601
602       if (N==2)
603       {
604          int c2;
605          int sign=1;
606          celt_norm_t v[2], w[2];
607          celt_norm_t *x2, *y2;
608          mbits = b-qalloc;
609          sbits = 0;
610          if (itheta != 0 && itheta != 16384)
611             sbits = 1<<BITRES;
612          mbits -= sbits;
613          c = itheta > 8192 ? 1 : 0;
614          c2 = 1-c;
615
616          x2 = X;
617          y2 = Y;
618          if (c==0)
619          {
620             v[0] = x2[0];
621             v[1] = x2[1];
622             w[0] = y2[0];
623             w[1] = y2[1];
624          } else {
625             v[0] = y2[0];
626             v[1] = y2[1];
627             w[0] = x2[0];
628             w[1] = x2[1];
629          }
630          q1 = bits2pulses(m, BPbits[i], N, mbits);
631          curr_bits = pulses2bits(BPbits[i], N, q1)+qalloc+sbits;
632          remaining_bits -= curr_bits;
633          while (remaining_bits < 0 && q1 > 0)
634          {
635             remaining_bits += curr_bits;
636             q1--;
637             curr_bits = pulses2bits(BPbits[i], N, q1)+qalloc;
638             remaining_bits -= curr_bits;
639          }
640
641          if (q1 > 0)
642          {
643             int spread = fold ? B : 0;
644             alg_quant(v, N, q1, spread, enc);
645          } else {
646             v[0] = QCONST16(1.f, 14);
647             v[1] = 0;
648          }
649          if (sbits)
650          {
651             if (v[0]*w[1] - v[1]*w[0] > 0)
652                sign = 1;
653             else
654                sign = -1;
655             ec_enc_bits(enc, sign==1, 1);
656          } else {
657             sign = 1;
658          }
659          w[0] = -sign*v[1];
660          w[1] = sign*v[0];
661          if (c==0)
662          {
663             x2[0] = v[0];
664             x2[1] = v[1];
665             y2[0] = w[0];
666             y2[1] = w[1];
667          } else {
668             x2[0] = w[0];
669             x2[1] = w[1];
670             y2[0] = v[0];
671             y2[1] = v[1];
672          }
673       } else {
674          
675          mbits = (b-qalloc/2-delta)/2;
676          if (mbits > b-qalloc)
677             mbits = b-qalloc;
678          if (mbits<0)
679             mbits=0;
680          sbits = b-qalloc-mbits;
681          q1 = bits2pulses(m, BPbits[i], N, mbits);
682          q2 = bits2pulses(m, BPbits[i], N, sbits);
683          curr_bits = pulses2bits(BPbits[i], N, q1)+pulses2bits(BPbits[i], N, q2)+qalloc;
684          remaining_bits -= curr_bits;
685          while (remaining_bits < 0 && (q1 > 0 || q2 > 0))
686          {
687             remaining_bits += curr_bits;
688             if (q1>q2)
689             {
690                q1--;
691                curr_bits = pulses2bits(BPbits[i], N, q1)+pulses2bits(BPbits[i], N, q2)+qalloc;
692             } else {
693                q2--;
694                curr_bits = pulses2bits(BPbits[i], N, q1)+pulses2bits(BPbits[i], N, q2)+qalloc;
695             }
696             remaining_bits -= curr_bits;
697          }
698
699          if (q1 > 0) {
700             int spread = fold ? B : 0;
701             alg_quant(X, N, q1, spread, enc);
702          } else {
703             intra_fold(m, eBands[i+1]-eBands[i], norm, X, eBands[i], B);
704          }
705          if (q2 > 0) {
706             int spread = fold ? B : 0;
707             alg_quant(Y, N, q2, spread, enc);
708          } else
709             for (j=0;j<N;j++)
710                Y[j] = 0;
711       }
712       
713       balance += pulses[i] + tell;
714
715 #ifdef FIXED_POINT
716       mid = imid;
717       side = iside;
718 #else
719       mid = (1./32768)*imid;
720       side = (1./32768)*iside;
721 #endif
722       for (j=0;j<N;j++)
723          norm[eBands[i]+j] = MULT16_16_Q15(n,X[j]);
724
725       for (j=0;j<N;j++)
726          X[j] = MULT16_16_Q15(X[j], mid);
727       for (j=0;j<N;j++)
728          Y[j] = MULT16_16_Q15(Y[j], side);
729
730       stereo_band_mix(m, X, Y, bandE, 0, i, -1);
731       renormalise_vector(X, Q15ONE, N, 1);
732       renormalise_vector(Y, Q15ONE, N, 1);
733    }
734    RESTORE_STACK;
735 }
736 #endif /* DISABLE_STEREO */
737
738
739 #ifndef DISABLE_STEREO
740
741 void unquant_bands_stereo(const CELTMode *m, celt_norm_t *_X, const celt_ener_t *bandE, int *pulses, int shortBlocks, int fold, int total_bits, ec_dec *dec)
742 {
743    int i, j, remaining_bits, balance;
744    const celt_int16 * restrict eBands = m->eBands;
745    celt_norm_t * restrict norm;
746    VARDECL(celt_norm_t, _norm);
747    int B;
748    celt_word16_t mid, side;
749    SAVE_STACK;
750
751    B = shortBlocks ? m->nbShortMdcts : 1;
752    ALLOC(_norm, eBands[m->nbEBands+1], celt_norm_t);
753    norm = _norm;
754
755    balance = 0;
756    for (i=0;i<m->nbEBands;i++)
757    {
758       int c;
759       int tell;
760       int q1, q2;
761       celt_word16_t n;
762       const celt_int16 * const *BPbits;
763       int b, qb;
764       int N;
765       int curr_balance, curr_bits;
766       int imid, iside, itheta;
767       int mbits, sbits, delta;
768       int qalloc;
769       celt_norm_t * restrict X, * restrict Y;
770       
771       X = _X+eBands[i];
772       Y = X+eBands[m->nbEBands+1];
773       BPbits = m->bits;
774
775       N = eBands[i+1]-eBands[i];
776       tell = ec_dec_tell(dec, BITRES);
777       if (i != 0)
778          balance -= tell;
779       remaining_bits = (total_bits<<BITRES)-tell-1;
780       curr_balance = (m->nbEBands-i);
781       if (curr_balance > 3)
782          curr_balance = 3;
783       curr_balance = balance / curr_balance;
784       b = IMIN(remaining_bits+1,pulses[i]+curr_balance);
785       if (b<0)
786          b = 0;
787
788       qb = (b-2*(N-1)*(QTHETA_OFFSET-log2_frac(N,BITRES)))/(32*(N-1));
789       if (qb > (b>>BITRES)-1)
790          qb = (b>>BITRES)-1;
791       if (qb>14)
792          qb = 14;
793       if (qb<0)
794          qb = 0;
795       qalloc = log2_frac((1<<qb)+1,BITRES);
796       if (qb==0)
797       {
798          itheta=0;
799       } else {
800          int shift;
801          shift = 14-qb;
802          itheta = ec_dec_uint(dec, (1<<qb)+1);
803          itheta <<= shift;
804       }
805       if (itheta == 0)
806       {
807          imid = 32767;
808          iside = 0;
809          delta = -10000;
810       } else if (itheta == 16384)
811       {
812          imid = 0;
813          iside = 32767;
814          delta = 10000;
815       } else {
816          imid = bitexact_cos(itheta);
817          iside = bitexact_cos(16384-itheta);
818          delta = (N-1)*(log2_frac(iside,BITRES+2)-log2_frac(imid,BITRES+2))>>2;
819       }
820       n = SHL16(celt_sqrt((eBands[i+1]-eBands[i])),11);
821
822       if (N==2)
823       {
824          int c2;
825          int sign=1;
826          celt_norm_t v[2], w[2];
827          celt_norm_t *x2, *y2;
828          mbits = b-qalloc;
829          sbits = 0;
830          if (itheta != 0 && itheta != 16384)
831             sbits = 1<<BITRES;
832          mbits -= sbits;
833          c = itheta > 8192 ? 1 : 0;
834          c2 = 1-c;
835
836          x2 = X;
837          y2 = Y;
838          v[0] = x2[c];
839          v[1] = y2[c];
840          w[0] = x2[c2];
841          w[1] = y2[c2];
842          q1 = bits2pulses(m, BPbits[i], N, mbits);
843          curr_bits = pulses2bits(BPbits[i], N, q1)+qalloc+sbits;
844          remaining_bits -= curr_bits;
845          while (remaining_bits < 0 && q1 > 0)
846          {
847             remaining_bits += curr_bits;
848             q1--;
849             curr_bits = pulses2bits(BPbits[i], N, q1)+qalloc;
850             remaining_bits -= curr_bits;
851          }
852
853          if (q1 > 0)
854          {
855             int spread = fold ? B : 0;
856             alg_unquant(v, N, q1, spread, dec);
857          } else {
858             v[0] = QCONST16(1.f, 14);
859             v[1] = 0;
860          }
861          if (sbits)
862             sign = 2*ec_dec_bits(dec, 1)-1;
863          else
864             sign = 1;
865          w[0] = -sign*v[1];
866          w[1] = sign*v[0];
867          if (c==0)
868          {
869             x2[0] = v[0];
870             x2[1] = v[1];
871             y2[0] = w[0];
872             y2[1] = w[1];
873          } else {
874             x2[0] = w[0];
875             x2[1] = w[1];
876             y2[0] = v[0];
877             y2[1] = v[1];
878          }
879       } else {
880          mbits = (b-qalloc/2-delta)/2;
881          if (mbits > b-qalloc)
882             mbits = b-qalloc;
883          if (mbits<0)
884             mbits=0;
885          sbits = b-qalloc-mbits;
886          q1 = bits2pulses(m, BPbits[i], N, mbits);
887          q2 = bits2pulses(m, BPbits[i], N, sbits);
888          curr_bits = pulses2bits(BPbits[i], N, q1)+pulses2bits(BPbits[i], N, q2)+qalloc;
889          remaining_bits -= curr_bits;
890          while (remaining_bits < 0 && (q1 > 0 || q2 > 0))
891          {
892             remaining_bits += curr_bits;
893             if (q1>q2)
894             {
895                q1--;
896                curr_bits = pulses2bits(BPbits[i], N, q1)+pulses2bits(BPbits[i], N, q2)+qalloc;
897             } else {
898                q2--;
899                curr_bits = pulses2bits(BPbits[i], N, q1)+pulses2bits(BPbits[i], N, q2)+qalloc;
900             }
901             remaining_bits -= curr_bits;
902          }
903          
904          if (q1 > 0)
905          {
906             int spread = fold ? B : 0;
907             alg_unquant(X, N, q1, spread, dec);
908          } else
909             intra_fold(m, eBands[i+1]-eBands[i], norm, X, eBands[i], B);
910          if (q2 > 0)
911          {
912             int spread = fold ? B : 0;
913             alg_unquant(Y, N, q2, spread, dec);
914          } else
915             for (j=0;j<N;j++)
916                Y[j] = 0;
917             /*orthogonalize(X+C*eBands[i], X+C*eBands[i]+N, N);*/
918       }
919       balance += pulses[i] + tell;
920       
921 #ifdef FIXED_POINT
922       mid = imid;
923       side = iside;
924 #else
925       mid = (1./32768)*imid;
926       side = (1./32768)*iside;
927 #endif
928       for (j=0;j<N;j++)
929          norm[eBands[i]+j] = MULT16_16_Q15(n,X[j]);
930
931       for (j=0;j<N;j++)
932          X[j] = MULT16_16_Q15(X[j], mid);
933       for (j=0;j<N;j++)
934          Y[j] = MULT16_16_Q15(Y[j], side);
935
936       stereo_band_mix(m, X, Y, bandE, 0, i, -1);
937       renormalise_vector(X, Q15ONE, N, 1);
938       renormalise_vector(Y, Q15ONE, N, 1);
939    }
940    RESTORE_STACK;
941 }
942
943 #endif /* DISABLE_STEREO */