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