Adjust the splitting threshold.
[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 celt_uint32 lcg_rand(celt_uint32 seed)
49 {
50    return 1664525 * seed + 1013904223;
51 }
52
53 /* This is a cos() approximation designed to be bit-exact on any platform. Bit exactness
54    with this approximation is important because it has an impact on the bit allocation */
55 static celt_int16 bitexact_cos(celt_int16 x)
56 {
57    celt_int32 tmp;
58    celt_int16 x2;
59    tmp = (4096+((celt_int32)(x)*(x)))>>13;
60    if (tmp > 32767)
61       tmp = 32767;
62    x2 = tmp;
63    x2 = (32767-x2) + FRAC_MUL16(x2, (-7651 + FRAC_MUL16(x2, (8277 + FRAC_MUL16(-626, x2)))));
64    if (x2 > 32766)
65       x2 = 32766;
66    return 1+x2;
67 }
68
69 static int bitexact_log2tan(int isin,int icos)
70 {
71    int lc;
72    int ls;
73    lc=EC_ILOG(icos);
74    ls=EC_ILOG(isin);
75    icos<<=15-lc;
76    isin<<=15-ls;
77    return (ls-lc<<11)
78          +FRAC_MUL16(isin, FRAC_MUL16(isin, -2597) + 7932)
79          -FRAC_MUL16(icos, FRAC_MUL16(icos, -2597) + 7932);
80 }
81
82 #ifdef FIXED_POINT
83 /* Compute the amplitude (sqrt energy) in each of the bands */
84 void compute_band_energies(const CELTMode *m, const celt_sig *X, celt_ener *bank, int end, int _C, int M)
85 {
86    int i, c, N;
87    const celt_int16 *eBands = m->eBands;
88    const int C = CHANNELS(_C);
89    N = M*m->shortMdctSize;
90    c=0; do {
91       for (i=0;i<end;i++)
92       {
93          int j;
94          celt_word32 maxval=0;
95          celt_word32 sum = 0;
96          
97          j=M*eBands[i]; do {
98             maxval = MAX32(maxval, X[j+c*N]);
99             maxval = MAX32(maxval, -X[j+c*N]);
100          } while (++j<M*eBands[i+1]);
101          
102          if (maxval > 0)
103          {
104             int shift = celt_ilog2(maxval)-10;
105             j=M*eBands[i]; do {
106                sum = MAC16_16(sum, EXTRACT16(VSHR32(X[j+c*N],shift)),
107                                    EXTRACT16(VSHR32(X[j+c*N],shift)));
108             } while (++j<M*eBands[i+1]);
109             /* We're adding one here to make damn sure we never end up with a pitch vector that's
110                larger than unity norm */
111             bank[i+c*m->nbEBands] = EPSILON+VSHR32(EXTEND32(celt_sqrt(sum)),-shift);
112          } else {
113             bank[i+c*m->nbEBands] = EPSILON;
114          }
115          /*printf ("%f ", bank[i+c*m->nbEBands]);*/
116       }
117    } while (++c<C);
118    /*printf ("\n");*/
119 }
120
121 /* Normalise each band such that the energy is one. */
122 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)
123 {
124    int i, c, N;
125    const celt_int16 *eBands = m->eBands;
126    const int C = CHANNELS(_C);
127    N = M*m->shortMdctSize;
128    c=0; do {
129       i=0; do {
130          celt_word16 g;
131          int j,shift;
132          celt_word16 E;
133          shift = celt_zlog2(bank[i+c*m->nbEBands])-13;
134          E = VSHR32(bank[i+c*m->nbEBands], shift);
135          g = EXTRACT16(celt_rcp(SHL32(E,3)));
136          j=M*eBands[i]; do {
137             X[j+c*N] = MULT16_16_Q15(VSHR32(freq[j+c*N],shift-1),g);
138          } while (++j<M*eBands[i+1]);
139       } while (++i<end);
140    } while (++c<C);
141 }
142
143 #else /* FIXED_POINT */
144 /* Compute the amplitude (sqrt energy) in each of the bands */
145 void compute_band_energies(const CELTMode *m, const celt_sig *X, celt_ener *bank, int end, int _C, int M)
146 {
147    int i, c, N;
148    const celt_int16 *eBands = m->eBands;
149    const int C = CHANNELS(_C);
150    N = M*m->shortMdctSize;
151    c=0; do {
152       for (i=0;i<end;i++)
153       {
154          int j;
155          celt_word32 sum = 1e-27f;
156          for (j=M*eBands[i];j<M*eBands[i+1];j++)
157             sum += X[j+c*N]*X[j+c*N];
158          bank[i+c*m->nbEBands] = celt_sqrt(sum);
159          /*printf ("%f ", bank[i+c*m->nbEBands]);*/
160       }
161    } while (++c<C);
162    /*printf ("\n");*/
163 }
164
165 /* Normalise each band such that the energy is one. */
166 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)
167 {
168    int i, c, N;
169    const celt_int16 *eBands = m->eBands;
170    const int C = CHANNELS(_C);
171    N = M*m->shortMdctSize;
172    c=0; do {
173       for (i=0;i<end;i++)
174       {
175          int j;
176          celt_word16 g = 1.f/(1e-27f+bank[i+c*m->nbEBands]);
177          for (j=M*eBands[i];j<M*eBands[i+1];j++)
178             X[j+c*N] = freq[j+c*N]*g;
179       }
180    } while (++c<C);
181 }
182
183 #endif /* FIXED_POINT */
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    c=0; do {
194       celt_sig * restrict f;
195       const celt_norm * restrict x;
196       f = freq+c*N;
197       x = X+c*N;
198       for (i=0;i<end;i++)
199       {
200          int j, band_end;
201          celt_word32 g = SHR32(bank[i+c*m->nbEBands],1);
202          j=M*eBands[i];
203          band_end = M*eBands[i+1];
204          do {
205             *f++ = SHL32(MULT16_32_Q15(*x, g),2);
206             x++;
207          } while (++j<band_end);
208       }
209       for (i=M*eBands[m->nbEBands];i<N;i++)
210          *f++ = 0;
211    } while (++c<C);
212 }
213
214 /* This prevents energy collapse for transients with multiple short MDCTs */
215 void anti_collapse(const CELTMode *m, celt_norm *_X, unsigned char *collapse_masks, int LM, int C, int CC, int size,
216       int start, int end, celt_word16 *logE, celt_word16 *prev1logE,
217       celt_word16 *prev2logE, int *pulses, celt_uint32 seed)
218 {
219    int c, i, j, k;
220    for (i=start;i<end;i++)
221    {
222       int N0;
223       celt_word16 thresh, sqrt_1;
224       int depth;
225 #ifdef FIXED_POINT
226       int shift;
227 #endif
228
229       N0 = m->eBands[i+1]-m->eBands[i];
230       /* depth in 1/8 bits */
231       depth = (1+pulses[i])/(m->eBands[i+1]-m->eBands[i]<<LM);
232
233 #ifdef FIXED_POINT
234       thresh = MULT16_32_Q15(QCONST16(0.5f, 15), MIN32(32767,SHR32(celt_exp2(-SHL16(depth, 10-BITRES)),1) ));
235       {
236          celt_word32 t;
237          t = N0<<LM;
238          shift = celt_ilog2(t)>>1;
239          t = SHL32(t, (7-shift)<<1);
240          sqrt_1 = celt_rsqrt_norm(t);
241       }
242 #else
243       thresh = .5f*celt_exp2(-.125f*depth);
244       sqrt_1 = celt_rsqrt(N0<<LM);
245 #endif
246
247       c=0; do
248       {
249          celt_norm *X;
250          celt_word16 prev1;
251          celt_word16 prev2;
252          celt_word16 Ediff;
253          celt_word16 r;
254          prev1 = prev1logE[c*m->nbEBands+i];
255          prev2 = prev2logE[c*m->nbEBands+i];
256          if (C<CC)
257          {
258             prev1 = MAX16(prev1,prev1logE[m->nbEBands+i]);
259             prev2 = MAX16(prev2,prev2logE[m->nbEBands+i]);
260          }
261          Ediff = logE[c*m->nbEBands+i]-MIN16(prev1,prev2);
262          Ediff = MAX16(0, Ediff);
263
264 #ifdef FIXED_POINT
265          if (Ediff < 16384)
266             r = 2*MIN16(16383,SHR32(celt_exp2(-Ediff),1));
267          else
268             r = 0;
269          if (LM==3)
270             r = MULT16_16_Q14(23170, MIN32(23169, r));
271          r = SHR16(MIN16(thresh, r),1);
272          r = SHR32(MULT16_16_Q15(sqrt_1, r),shift);
273 #else
274          /* r needs to be multiplied by 2 or 2*sqrt(2) depending on LM because
275             short blocks don't have the same energy as long */
276          r = 2.f*celt_exp2(-Ediff);
277          if (LM==3)
278             r *= 1.41421356f;
279          r = MIN16(thresh, r);
280          r = r*sqrt_1;
281 #endif
282          X = _X+c*size+(m->eBands[i]<<LM);
283          for (k=0;k<1<<LM;k++)
284          {
285             /* Detect collapse */
286             if (!(collapse_masks[i*C+c]&1<<k))
287             {
288                /* Fill with noise */
289                for (j=0;j<N0;j++)
290                {
291                   seed = lcg_rand(seed);
292                   X[(j<<LM)+k] = (seed&0x8000 ? r : -r);
293                }
294             }
295          }
296          /* We just added some energy, so we need to renormalise */
297          renormalise_vector(X, N0<<LM, Q15ONE);
298       } while (++c<C);
299    }
300
301 }
302
303
304 static void intensity_stereo(const CELTMode *m, celt_norm *X, celt_norm *Y, const celt_ener *bank, int bandID, int N)
305 {
306    int i = bandID;
307    int j;
308    celt_word16 a1, a2;
309    celt_word16 left, right;
310    celt_word16 norm;
311 #ifdef FIXED_POINT
312    int shift = celt_zlog2(MAX32(bank[i], bank[i+m->nbEBands]))-13;
313 #endif
314    left = VSHR32(bank[i],shift);
315    right = VSHR32(bank[i+m->nbEBands],shift);
316    norm = EPSILON + celt_sqrt(EPSILON+MULT16_16(left,left)+MULT16_16(right,right));
317    a1 = DIV32_16(SHL32(EXTEND32(left),14),norm);
318    a2 = DIV32_16(SHL32(EXTEND32(right),14),norm);
319    for (j=0;j<N;j++)
320    {
321       celt_norm r, l;
322       l = X[j];
323       r = Y[j];
324       X[j] = MULT16_16_Q14(a1,l) + MULT16_16_Q14(a2,r);
325       /* Side is not encoded, no need to calculate */
326    }
327 }
328
329 static void stereo_split(celt_norm *X, celt_norm *Y, int N)
330 {
331    int j;
332    for (j=0;j<N;j++)
333    {
334       celt_norm r, l;
335       l = MULT16_16_Q15(QCONST16(.70710678f,15), X[j]);
336       r = MULT16_16_Q15(QCONST16(.70710678f,15), Y[j]);
337       X[j] = l+r;
338       Y[j] = r-l;
339    }
340 }
341
342 static void stereo_merge(celt_norm *X, celt_norm *Y, celt_word16 mid, int N)
343 {
344    int j;
345    celt_word32 xp=0, side=0;
346    celt_word32 El, Er;
347    celt_word16 mid2;
348 #ifdef FIXED_POINT
349    int kl, kr;
350 #endif
351    celt_word32 t, lgain, rgain;
352
353    /* Compute the norm of X+Y and X-Y as |X|^2 + |Y|^2 +/- sum(xy) */
354    for (j=0;j<N;j++)
355    {
356       xp = MAC16_16(xp, X[j], Y[j]);
357       side = MAC16_16(side, Y[j], Y[j]);
358    }
359    /* Compensating for the mid normalization */
360    xp = MULT16_32_Q15(mid, xp);
361    /* mid and side are in Q15, not Q14 like X and Y */
362    mid2 = SHR32(mid, 1);
363    El = MULT16_16(mid2, mid2) + side - 2*xp;
364    Er = MULT16_16(mid2, mid2) + side + 2*xp;
365    if (Er < QCONST32(6e-4f, 28) || El < QCONST32(6e-4f, 28))
366    {
367       for (j=0;j<N;j++)
368          Y[j] = X[j];
369       return;
370    }
371
372 #ifdef FIXED_POINT
373    kl = celt_ilog2(El)>>1;
374    kr = celt_ilog2(Er)>>1;
375 #endif
376    t = VSHR32(El, (kl-7)<<1);
377    lgain = celt_rsqrt_norm(t);
378    t = VSHR32(Er, (kr-7)<<1);
379    rgain = celt_rsqrt_norm(t);
380
381 #ifdef FIXED_POINT
382    if (kl < 7)
383       kl = 7;
384    if (kr < 7)
385       kr = 7;
386 #endif
387
388    for (j=0;j<N;j++)
389    {
390       celt_norm r, l;
391       /* Apply mid scaling (side is already scaled) */
392       l = MULT16_16_Q15(mid, X[j]);
393       r = Y[j];
394       X[j] = EXTRACT16(PSHR32(MULT16_16(lgain, SUB16(l,r)), kl+1));
395       Y[j] = EXTRACT16(PSHR32(MULT16_16(rgain, ADD16(l,r)), kr+1));
396    }
397 }
398
399 /* Decide whether we should spread the pulses in the current frame */
400 int spreading_decision(const CELTMode *m, celt_norm *X, int *average,
401       int last_decision, int *hf_average, int *tapset_decision, int update_hf,
402       int end, int _C, int M)
403 {
404    int i, c, N0;
405    int sum = 0, nbBands=0;
406    const int C = CHANNELS(_C);
407    const celt_int16 * restrict eBands = m->eBands;
408    int decision;
409    int hf_sum=0;
410    
411    N0 = M*m->shortMdctSize;
412
413    if (M*(eBands[end]-eBands[end-1]) <= 8)
414       return SPREAD_NONE;
415    c=0; do {
416       for (i=0;i<end;i++)
417       {
418          int j, N, tmp=0;
419          int tcount[3] = {0};
420          celt_norm * restrict x = X+M*eBands[i]+c*N0;
421          N = M*(eBands[i+1]-eBands[i]);
422          if (N<=8)
423             continue;
424          /* Compute rough CDF of |x[j]| */
425          for (j=0;j<N;j++)
426          {
427             celt_word32 x2N; /* Q13 */
428
429             x2N = MULT16_16(MULT16_16_Q15(x[j], x[j]), N);
430             if (x2N < QCONST16(0.25f,13))
431                tcount[0]++;
432             if (x2N < QCONST16(0.0625f,13))
433                tcount[1]++;
434             if (x2N < QCONST16(0.015625f,13))
435                tcount[2]++;
436          }
437
438          /* Only include four last bands (8 kHz and up) */
439          if (i>m->nbEBands-4)
440             hf_sum += 32*(tcount[1]+tcount[0])/N;
441          tmp = (2*tcount[2] >= N) + (2*tcount[1] >= N) + (2*tcount[0] >= N);
442          sum += tmp*256;
443          nbBands++;
444       }
445    } while (++c<C);
446
447    if (update_hf)
448    {
449       if (hf_sum)
450          hf_sum /= C*(4-m->nbEBands+end);
451       *hf_average = (*hf_average+hf_sum)>>1;
452       hf_sum = *hf_average;
453       if (*tapset_decision==2)
454          hf_sum += 4;
455       else if (*tapset_decision==0)
456          hf_sum -= 4;
457       if (hf_sum > 22)
458          *tapset_decision=2;
459       else if (hf_sum > 18)
460          *tapset_decision=1;
461       else
462          *tapset_decision=0;
463    }
464    /*printf("%d %d %d\n", hf_sum, *hf_average, *tapset_decision);*/
465    sum /= nbBands;
466    /* Recursive averaging */
467    sum = (sum+*average)>>1;
468    *average = sum;
469    /* Hysteresis */
470    sum = (3*sum + (((3-last_decision)<<7) + 64) + 2)>>2;
471    if (sum < 80)
472    {
473       decision = SPREAD_AGGRESSIVE;
474    } else if (sum < 256)
475    {
476       decision = SPREAD_NORMAL;
477    } else if (sum < 384)
478    {
479       decision = SPREAD_LIGHT;
480    } else {
481       decision = SPREAD_NONE;
482    }
483    return decision;
484 }
485
486 #ifdef MEASURE_NORM_MSE
487
488 float MSE[30] = {0};
489 int nbMSEBands = 0;
490 int MSECount[30] = {0};
491
492 void dump_norm_mse(void)
493 {
494    int i;
495    for (i=0;i<nbMSEBands;i++)
496    {
497       printf ("%g ", MSE[i]/MSECount[i]);
498    }
499    printf ("\n");
500 }
501
502 void measure_norm_mse(const CELTMode *m, float *X, float *X0, float *bandE, float *bandE0, int M, int N, int C)
503 {
504    static int init = 0;
505    int i;
506    if (!init)
507    {
508       atexit(dump_norm_mse);
509       init = 1;
510    }
511    for (i=0;i<m->nbEBands;i++)
512    {
513       int j;
514       int c;
515       float g;
516       if (bandE0[i]<10 || (C==2 && bandE0[i+m->nbEBands]<1))
517          continue;
518       c=0; do {
519          g = bandE[i+c*m->nbEBands]/(1e-15+bandE0[i+c*m->nbEBands]);
520          for (j=M*m->eBands[i];j<M*m->eBands[i+1];j++)
521             MSE[i] += (g*X[j+c*N]-X0[j+c*N])*(g*X[j+c*N]-X0[j+c*N]);
522       } while (++c<C);
523       MSECount[i]+=C;
524    }
525    nbMSEBands = m->nbEBands;
526 }
527
528 #endif
529
530 /* Indexing table for converting from natural Hadamard to ordery Hadamard
531    This is essentially a bit-reversed Gray, on top of which we've added
532    an inversion of the order because we want the DC at the end rather than
533    the beginning. The lines are for N=2, 4, 8, 16 */
534 static const int ordery_table[] = {
535        1,  0,
536        3,  0,  2,  1,
537        7,  0,  4,  3,  6,  1,  5,  2,
538       15,  0,  8,  7, 12,  3, 11,  4, 14,  1,  9,  6, 13,  2, 10,  5,
539 };
540
541 static void deinterleave_hadamard(celt_norm *X, int N0, int stride, int hadamard)
542 {
543    int i,j;
544    VARDECL(celt_norm, tmp);
545    int N;
546    SAVE_STACK;
547    N = N0*stride;
548    ALLOC(tmp, N, celt_norm);
549    if (hadamard)
550    {
551       const int *ordery = ordery_table+stride-2;
552       for (i=0;i<stride;i++)
553       {
554          for (j=0;j<N0;j++)
555             tmp[ordery[i]*N0+j] = X[j*stride+i];
556       }
557    } else {
558       for (i=0;i<stride;i++)
559          for (j=0;j<N0;j++)
560             tmp[i*N0+j] = X[j*stride+i];
561    }
562    for (j=0;j<N;j++)
563       X[j] = tmp[j];
564    RESTORE_STACK;
565 }
566
567 static void interleave_hadamard(celt_norm *X, int N0, int stride, int hadamard)
568 {
569    int i,j;
570    VARDECL(celt_norm, tmp);
571    int N;
572    SAVE_STACK;
573    N = N0*stride;
574    ALLOC(tmp, N, celt_norm);
575    if (hadamard)
576    {
577       const int *ordery = ordery_table+stride-2;
578       for (i=0;i<stride;i++)
579          for (j=0;j<N0;j++)
580             tmp[j*stride+i] = X[ordery[i]*N0+j];
581    } else {
582       for (i=0;i<stride;i++)
583          for (j=0;j<N0;j++)
584             tmp[j*stride+i] = X[i*N0+j];
585    }
586    for (j=0;j<N;j++)
587       X[j] = tmp[j];
588    RESTORE_STACK;
589 }
590
591 void haar1(celt_norm *X, int N0, int stride)
592 {
593    int i, j;
594    N0 >>= 1;
595    for (i=0;i<stride;i++)
596       for (j=0;j<N0;j++)
597       {
598          celt_norm tmp1, tmp2;
599          tmp1 = MULT16_16_Q15(QCONST16(.70710678f,15), X[stride*2*j+i]);
600          tmp2 = MULT16_16_Q15(QCONST16(.70710678f,15), X[stride*(2*j+1)+i]);
601          X[stride*2*j+i] = tmp1 + tmp2;
602          X[stride*(2*j+1)+i] = tmp1 - tmp2;
603       }
604 }
605
606 static int compute_qn(int N, int b, int offset, int pulse_cap, int stereo)
607 {
608    static const celt_int16 exp2_table8[8] =
609       {16384, 17866, 19483, 21247, 23170, 25267, 27554, 30048};
610    int qn, qb;
611    int N2 = 2*N-1;
612    if (stereo && N==2)
613       N2--;
614    /* The upper limit ensures that in a stereo split with itheta==16384, we'll
615        always have enough bits left over to code at least one pulse in the
616        side; otherwise it would collapse, since it doesn't get folded. */
617    qb = IMIN(b-pulse_cap-(4<<BITRES), (b+N2*offset)/N2);
618
619    qb = IMIN(8<<BITRES, qb);
620
621    if (qb<(1<<BITRES>>1)) {
622       qn = 1;
623    } else {
624       qn = exp2_table8[qb&0x7]>>(14-(qb>>BITRES));
625       qn = (qn+1)>>1<<1;
626    }
627    celt_assert(qn <= 256);
628    return qn;
629 }
630
631 /* This function is responsible for encoding and decoding a band for both
632    the mono and stereo case. Even in the mono case, it can split the band
633    in two and transmit the energy difference with the two half-bands. It
634    can be called recursively so bands can end up being split in 8 parts. */
635 static unsigned quant_band(int encode, const CELTMode *m, int i, celt_norm *X, celt_norm *Y,
636       int N, int b, int spread, int B, int intensity, int tf_change, celt_norm *lowband, int resynth, void *ec,
637       celt_int32 *remaining_bits, int LM, celt_norm *lowband_out, const celt_ener *bandE, int level,
638       celt_uint32 *seed, celt_word16 gain, celt_norm *lowband_scratch, int fill)
639 {
640    const unsigned char *cache;
641    int q;
642    int curr_bits;
643    int stereo, split;
644    int imid=0, iside=0;
645    int N0=N;
646    int N_B=N;
647    int N_B0;
648    int B0=B;
649    int time_divide=0;
650    int recombine=0;
651    int inv = 0;
652    celt_word16 mid=0, side=0;
653    int longBlocks;
654    unsigned cm=0;
655
656    longBlocks = B0==1;
657
658    N_B /= B;
659    N_B0 = N_B;
660
661    split = stereo = Y != NULL;
662
663    /* Special case for one sample */
664    if (N==1)
665    {
666       int c;
667       celt_norm *x = X;
668       c=0; do {
669          int sign=0;
670          if (*remaining_bits>=1<<BITRES)
671          {
672             if (encode)
673             {
674                sign = x[0]<0;
675                ec_enc_bits((ec_enc*)ec, sign, 1);
676             } else {
677                sign = ec_dec_bits((ec_dec*)ec, 1);
678             }
679             *remaining_bits -= 1<<BITRES;
680             b-=1<<BITRES;
681          }
682          if (resynth)
683             x[0] = sign ? -NORM_SCALING : NORM_SCALING;
684          x = Y;
685       } while (++c<1+stereo);
686       if (lowband_out)
687          lowband_out[0] = SHR16(X[0],4);
688       return 1;
689    }
690
691    if (!stereo && level == 0)
692    {
693       int k;
694       if (tf_change>0)
695          recombine = tf_change;
696       /* Band recombining to increase frequency resolution */
697
698       if (lowband && (recombine || ((N_B&1) == 0 && tf_change<0) || B0>1))
699       {
700          int j;
701          for (j=0;j<N;j++)
702             lowband_scratch[j] = lowband[j];
703          lowband = lowband_scratch;
704       }
705
706       for (k=0;k<recombine;k++)
707       {
708          if (encode)
709             haar1(X, N>>k, 1<<k);
710          if (lowband)
711             haar1(lowband, N>>k, 1<<k);
712          fill |= fill<<(1<<k);
713       }
714       B>>=recombine;
715       N_B<<=recombine;
716
717       /* Increasing the time resolution */
718       while ((N_B&1) == 0 && tf_change<0)
719       {
720          if (encode)
721             haar1(X, N_B, B);
722          if (lowband)
723             haar1(lowband, N_B, B);
724          fill |= fill<<B;
725          B <<= 1;
726          N_B >>= 1;
727          time_divide++;
728          tf_change++;
729       }
730       B0=B;
731       N_B0 = N_B;
732
733       /* Reorganize the samples in time order instead of frequency order */
734       if (B0>1)
735       {
736          if (encode)
737             deinterleave_hadamard(X, N_B>>recombine, B0<<recombine, longBlocks);
738          if (lowband)
739             deinterleave_hadamard(lowband, N_B>>recombine, B0<<recombine, longBlocks);
740       }
741    }
742
743    /* If we need 1 more bit than we can produce, split the band in two. */
744    cache = m->cache.bits + m->cache.index[(LM+1)*m->nbEBands+i];
745    if (!stereo && LM != -1 && b > cache[cache[0]]+9 && N>2)
746    {
747       if (LM>0 || (N&1)==0)
748       {
749          N >>= 1;
750          Y = X+N;
751          split = 1;
752          LM -= 1;
753          if (B==1)
754             fill = fill&1|fill<<1;
755          B = (B+1)>>1;
756       }
757    }
758
759    if (split)
760    {
761       int qn;
762       int itheta=0;
763       int mbits, sbits, delta;
764       int qalloc;
765       int pulse_cap;
766       int offset;
767       int orig_fill;
768       celt_int32 tell;
769
770       /* Decide on the resolution to give to the split parameter theta */
771       pulse_cap = m->logN[i]+(LM<<BITRES);
772       offset = (pulse_cap>>1) - (stereo ? QTHETA_OFFSET_STEREO : QTHETA_OFFSET);
773       qn = compute_qn(N, b, offset, pulse_cap, stereo);
774       if (stereo && i>=intensity)
775          qn = 1;
776       if (encode)
777       {
778          /* theta is the atan() of the ratio between the (normalized)
779             side and mid. With just that parameter, we can re-scale both
780             mid and side because we know that 1) they have unit norm and
781             2) they are orthogonal. */
782          itheta = stereo_itheta(X, Y, stereo, N);
783       }
784       tell = encode ? ec_enc_tell(ec, BITRES) : ec_dec_tell(ec, BITRES);
785       if (qn!=1)
786       {
787          if (encode)
788             itheta = (itheta*qn+8192)>>14;
789
790          /* Entropy coding of the angle. We use a uniform pdf for the
791             time split, a step for stereo, and a triangular one for the rest. */
792          if (stereo && N>2)
793          {
794             int p0 = 3;
795             int x = itheta;
796             int x0 = qn/2;
797             int ft = p0*(x0+1) + x0;
798             /* Use a probability of p0 up to itheta=8192 and then use 1 after */
799             if (encode)
800             {
801                ec_encode((ec_enc*)ec,x<=x0?p0*x:(x-1-x0)+(x0+1)*p0,x<=x0?p0*(x+1):(x-x0)+(x0+1)*p0,ft);
802             } else {
803                int fs;
804                fs=ec_decode(ec,ft);
805                if (fs<(x0+1)*p0)
806                   x=fs/p0;
807                else
808                   x=x0+1+(fs-(x0+1)*p0);
809                ec_dec_update(ec,x<=x0?p0*x:(x-1-x0)+(x0+1)*p0,x<=x0?p0*(x+1):(x-x0)+(x0+1)*p0,ft);
810                itheta = x;
811             }
812          } else if (B0>1 || stereo) {
813             /* Uniform pdf */
814             if (encode)
815                ec_enc_uint((ec_enc*)ec, itheta, qn+1);
816             else
817                itheta = ec_dec_uint((ec_dec*)ec, qn+1);
818          } else {
819             int fs=1, ft;
820             ft = ((qn>>1)+1)*((qn>>1)+1);
821             if (encode)
822             {
823                int fl;
824
825                fs = itheta <= (qn>>1) ? itheta + 1 : qn + 1 - itheta;
826                fl = itheta <= (qn>>1) ? itheta*(itheta + 1)>>1 :
827                 ft - ((qn + 1 - itheta)*(qn + 2 - itheta)>>1);
828
829                ec_encode((ec_enc*)ec, fl, fl+fs, ft);
830             } else {
831                /* Triangular pdf */
832                int fl=0;
833                int fm;
834                fm = ec_decode((ec_dec*)ec, ft);
835
836                if (fm < ((qn>>1)*((qn>>1) + 1)>>1))
837                {
838                   itheta = (isqrt32(8*(celt_uint32)fm + 1) - 1)>>1;
839                   fs = itheta + 1;
840                   fl = itheta*(itheta + 1)>>1;
841                }
842                else
843                {
844                   itheta = (2*(qn + 1)
845                    - isqrt32(8*(celt_uint32)(ft - fm - 1) + 1))>>1;
846                   fs = qn + 1 - itheta;
847                   fl = ft - ((qn + 1 - itheta)*(qn + 2 - itheta)>>1);
848                }
849
850                ec_dec_update((ec_dec*)ec, fl, fl+fs, ft);
851             }
852          }
853          itheta = (celt_int32)itheta*16384/qn;
854          if (encode && stereo)
855          {
856             if (itheta==0)
857                intensity_stereo(m, X, Y, bandE, i, N);
858             else
859                stereo_split(X, Y, N);
860          }
861          /* TODO: Renormalising X and Y *may* help fixed-point a bit at very high rate.
862                   Let's do that at higher complexity */
863       } else if (stereo) {
864          if (encode)
865          {
866             inv = itheta > 8192;
867             if (inv)
868             {
869                int j;
870                for (j=0;j<N;j++)
871                   Y[j] = -Y[j];
872             }
873             intensity_stereo(m, X, Y, bandE, i, N);
874          }
875          if (b>2<<BITRES && *remaining_bits > 2<<BITRES)
876          {
877             if (encode)
878                ec_enc_bit_logp(ec, inv, 2);
879             else
880                inv = ec_dec_bit_logp(ec, 2);
881          } else
882             inv = 0;
883          itheta = 0;
884       }
885       qalloc = (encode ? ec_enc_tell(ec, BITRES) : ec_dec_tell(ec, BITRES))
886                - tell;
887       b -= qalloc;
888
889       orig_fill = fill;
890       if (itheta == 0)
891       {
892          imid = 32767;
893          iside = 0;
894          fill &= (1<<B)-1;
895          delta = -16384;
896       } else if (itheta == 16384)
897       {
898          imid = 0;
899          iside = 32767;
900          fill &= (1<<B)-1<<B;
901          delta = 16384;
902       } else {
903          imid = bitexact_cos(itheta);
904          iside = bitexact_cos(16384-itheta);
905          /* This is the mid vs side allocation that minimizes squared error
906             in that band. */
907          delta = FRAC_MUL16(N-1<<7,bitexact_log2tan(iside,imid));
908       }
909
910 #ifdef FIXED_POINT
911       mid = imid;
912       side = iside;
913 #else
914       mid = (1.f/32768)*imid;
915       side = (1.f/32768)*iside;
916 #endif
917
918       /* This is a special case for N=2 that only works for stereo and takes
919          advantage of the fact that mid and side are orthogonal to encode
920          the side with just one bit. */
921       if (N==2 && stereo)
922       {
923          int c;
924          int sign=0;
925          celt_norm *x2, *y2;
926          mbits = b;
927          sbits = 0;
928          /* Only need one bit for the side */
929          if (itheta != 0 && itheta != 16384)
930             sbits = 1<<BITRES;
931          mbits -= sbits;
932          c = itheta > 8192;
933          *remaining_bits -= qalloc+sbits;
934
935          x2 = c ? Y : X;
936          y2 = c ? X : Y;
937          if (sbits)
938          {
939             if (encode)
940             {
941                /* Here we only need to encode a sign for the side */
942                sign = x2[0]*y2[1] - x2[1]*y2[0] < 0;
943                ec_enc_bits((ec_enc*)ec, sign, 1);
944             } else {
945                sign = ec_dec_bits((ec_dec*)ec, 1);
946             }
947          }
948          sign = 1-2*sign;
949          /* We use orig_fill here because we want to fold the side, but if
950              itheta==16384, we'll have cleared the low bits of fill. */
951          cm = 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, orig_fill);
952          /* We don't split N=2 bands, so cm is either 1 or 0 (for a fold-collapse),
953              and there's no need to worry about mixing with the other channel. */
954          y2[0] = -sign*x2[1];
955          y2[1] = sign*x2[0];
956          if (resynth)
957          {
958             celt_norm tmp;
959             X[0] = MULT16_16_Q15(mid, X[0]);
960             X[1] = MULT16_16_Q15(mid, X[1]);
961             Y[0] = MULT16_16_Q15(side, Y[0]);
962             Y[1] = MULT16_16_Q15(side, Y[1]);
963             tmp = X[0];
964             X[0] = SUB16(tmp,Y[0]);
965             Y[0] = ADD16(tmp,Y[0]);
966             tmp = X[1];
967             X[1] = SUB16(tmp,Y[1]);
968             Y[1] = ADD16(tmp,Y[1]);
969          }
970       } else {
971          /* "Normal" split code */
972          celt_norm *next_lowband2=NULL;
973          celt_norm *next_lowband_out1=NULL;
974          int next_level=0;
975          celt_int32 rebalance;
976
977          /* Give more bits to low-energy MDCTs than they would otherwise deserve */
978          if (B0>1 && !stereo && (itheta&0x3fff))
979          {
980             if (itheta > 8192)
981                /* Rough approximation for pre-echo masking */
982                delta -= delta>>(4-LM);
983             else
984                /* Corresponds to a forward-masking slope of 1.5 dB per 10 ms */
985                delta = IMIN(0, delta + (N<<BITRES>>(5-LM)));
986          }
987          mbits = IMAX(0, IMIN(b, (b-delta)/2));
988          sbits = b-mbits;
989          *remaining_bits -= qalloc;
990
991          if (lowband && !stereo)
992             next_lowband2 = lowband+N; /* >32-bit split case */
993
994          /* Only stereo needs to pass on lowband_out. Otherwise, it's
995             handled at the end */
996          if (stereo)
997             next_lowband_out1 = lowband_out;
998          else
999             next_level = level+1;
1000
1001          rebalance = *remaining_bits;
1002          if (mbits >= sbits)
1003          {
1004             /* In stereo mode, we do not apply a scaling to the mid because we need the normalized
1005                mid for folding later */
1006             cm = quant_band(encode, m, i, X, NULL, N, mbits, spread, B, intensity, tf_change,
1007                   lowband, resynth, ec, remaining_bits, LM, next_lowband_out1,
1008                   NULL, next_level, seed, stereo ? Q15ONE : MULT16_16_P15(gain,mid), lowband_scratch, fill);
1009             rebalance = mbits - (rebalance-*remaining_bits);
1010             if (rebalance > 3<<BITRES && itheta!=0)
1011                sbits += rebalance - (3<<BITRES);
1012
1013             /* For a stereo split, the high bits of fill are always zero, so no
1014                folding will be done to the side. */
1015             cm |= quant_band(encode, m, i, Y, NULL, N, sbits, spread, B, intensity, tf_change,
1016                   next_lowband2, resynth, ec, remaining_bits, LM, NULL,
1017                   NULL, next_level, seed, MULT16_16_P15(gain,side), NULL, fill>>B)<<(B0>>1&stereo-1);
1018          } else {
1019             /* For a stereo split, the high bits of fill are always zero, so no
1020                folding will be done to the side. */
1021             cm = quant_band(encode, m, i, Y, NULL, N, sbits, spread, B, intensity, tf_change,
1022                   next_lowband2, resynth, ec, remaining_bits, LM, NULL,
1023                   NULL, next_level, seed, MULT16_16_P15(gain,side), NULL, fill>>B)<<(B0>>1&stereo-1);
1024             rebalance = sbits - (rebalance-*remaining_bits);
1025             if (rebalance > 3<<BITRES && itheta!=16384)
1026                mbits += rebalance - (3<<BITRES);
1027             /* In stereo mode, we do not apply a scaling to the mid because we need the normalized
1028                mid for folding later */
1029             cm |= quant_band(encode, m, i, X, NULL, N, mbits, spread, B, intensity, tf_change,
1030                   lowband, resynth, ec, remaining_bits, LM, next_lowband_out1,
1031                   NULL, next_level, seed, stereo ? Q15ONE : MULT16_16_P15(gain,mid), lowband_scratch, fill);
1032          }
1033       }
1034
1035    } else {
1036       /* This is the basic no-split case */
1037       q = bits2pulses(m, i, LM, b);
1038       curr_bits = pulses2bits(m, i, LM, q);
1039       *remaining_bits -= curr_bits;
1040
1041       /* Ensures we can never bust the budget */
1042       while (*remaining_bits < 0 && q > 0)
1043       {
1044          *remaining_bits += curr_bits;
1045          q--;
1046          curr_bits = pulses2bits(m, i, LM, q);
1047          *remaining_bits -= curr_bits;
1048       }
1049
1050       if (q!=0)
1051       {
1052          int K = get_pulses(q);
1053
1054          /* Finally do the actual quantization */
1055          if (encode)
1056             cm = alg_quant(X, N, K, spread, B, resynth, (ec_enc*)ec, gain);
1057          else
1058             cm = alg_unquant(X, N, K, spread, B, (ec_dec*)ec, gain);
1059       } else {
1060          /* If there's no pulse, fill the band anyway */
1061          int j;
1062          if (resynth)
1063          {
1064             unsigned cm_mask;
1065             /*B can be as large as 16, so this shift might overflow an int on a
1066                16-bit platform; use a long to get defined behavior.*/
1067             cm_mask = (unsigned)(1UL<<B)-1;
1068             fill &= cm_mask;
1069             if (!fill)
1070             {
1071                for (j=0;j<N;j++)
1072                   X[j] = 0;
1073             } else {
1074                if (lowband == NULL)
1075                {
1076                   /* Noise */
1077                   for (j=0;j<N;j++)
1078                   {
1079                      *seed = lcg_rand(*seed);
1080                      X[j] = (celt_int32)(*seed)>>20;
1081                   }
1082                   cm = cm_mask;
1083                } else {
1084                   /* Folded spectrum */
1085                   for (j=0;j<N;j++)
1086                      X[j] = lowband[j];
1087                   cm = fill;
1088                }
1089                renormalise_vector(X, N, gain);
1090             }
1091          }
1092       }
1093    }
1094
1095    /* This code is used by the decoder and by the resynthesis-enabled encoder */
1096    if (resynth)
1097    {
1098       if (stereo)
1099       {
1100          if (N!=2)
1101             stereo_merge(X, Y, mid, N);
1102          if (inv)
1103          {
1104             int j;
1105             for (j=0;j<N;j++)
1106                Y[j] = -Y[j];
1107          }
1108       } else if (level == 0)
1109       {
1110          int k;
1111
1112          /* Undo the sample reorganization going from time order to frequency order */
1113          if (B0>1)
1114             interleave_hadamard(X, N_B>>recombine, B0<<recombine, longBlocks);
1115
1116          /* Undo time-freq changes that we did earlier */
1117          N_B = N_B0;
1118          B = B0;
1119          for (k=0;k<time_divide;k++)
1120          {
1121             B >>= 1;
1122             N_B <<= 1;
1123             cm |= cm>>B;
1124             haar1(X, N_B, B);
1125          }
1126
1127          for (k=0;k<recombine;k++)
1128          {
1129             cm |= cm<<(1<<k);
1130             haar1(X, N0>>k, 1<<k);
1131          }
1132          B<<=recombine;
1133          N_B>>=recombine;
1134
1135          /* Scale output for later folding */
1136          if (lowband_out)
1137          {
1138             int j;
1139             celt_word16 n;
1140             n = celt_sqrt(SHL32(EXTEND32(N0),22));
1141             for (j=0;j<N0;j++)
1142                lowband_out[j] = MULT16_16_Q15(n,X[j]);
1143          }
1144          cm &= (1<<B)-1;
1145       }
1146    }
1147    return cm;
1148 }
1149
1150 void quant_all_bands(int encode, const CELTMode *m, int start, int end,
1151       celt_norm *_X, celt_norm *_Y, unsigned char *collapse_masks, const celt_ener *bandE, int *pulses,
1152       int shortBlocks, int spread, int dual_stereo, int intensity, int *tf_res, int resynth,
1153       celt_int32 total_bits, celt_int32 balance, void *ec, int LM, int codedBands, ec_uint32 *seed)
1154 {
1155    int i;
1156    celt_int32 remaining_bits;
1157    const celt_int16 * restrict eBands = m->eBands;
1158    celt_norm * restrict norm, * restrict norm2;
1159    VARDECL(celt_norm, _norm);
1160    VARDECL(celt_norm, lowband_scratch);
1161    int B;
1162    int M;
1163    int lowband_offset;
1164    int update_lowband = 1;
1165    int C = _Y != NULL ? 2 : 1;
1166    SAVE_STACK;
1167
1168    M = 1<<LM;
1169    B = shortBlocks ? M : 1;
1170    ALLOC(_norm, C*M*eBands[m->nbEBands], celt_norm);
1171    ALLOC(lowband_scratch, M*(eBands[m->nbEBands]-eBands[m->nbEBands-1]), celt_norm);
1172    norm = _norm;
1173    norm2 = norm + M*eBands[m->nbEBands];
1174
1175    lowband_offset = 0;
1176    for (i=start;i<end;i++)
1177    {
1178       celt_int32 tell;
1179       int b;
1180       int N;
1181       celt_int32 curr_balance;
1182       int effective_lowband=-1;
1183       celt_norm * restrict X, * restrict Y;
1184       int tf_change=0;
1185       unsigned x_cm;
1186       unsigned y_cm;
1187
1188       X = _X+M*eBands[i];
1189       if (_Y!=NULL)
1190          Y = _Y+M*eBands[i];
1191       else
1192          Y = NULL;
1193       N = M*eBands[i+1]-M*eBands[i];
1194       if (encode)
1195          tell = ec_enc_tell((ec_enc*)ec, BITRES);
1196       else
1197          tell = ec_dec_tell((ec_dec*)ec, BITRES);
1198
1199       /* Compute how many bits we want to allocate to this band */
1200       if (i != start)
1201          balance -= tell;
1202       remaining_bits = total_bits-tell-1;
1203       if (i <= codedBands-1)
1204       {
1205          curr_balance = balance / IMIN(3, codedBands-i);
1206          b = IMAX(0, IMIN(16383, IMIN(remaining_bits+1,pulses[i]+curr_balance)));
1207       } else {
1208          b = 0;
1209       }
1210
1211       if (resynth && M*eBands[i]-N >= M*eBands[start] && (update_lowband || lowband_offset==0))
1212             lowband_offset = i;
1213
1214       tf_change = tf_res[i];
1215       if (i>=m->effEBands)
1216       {
1217          X=norm;
1218          if (_Y!=NULL)
1219             Y = norm;
1220       }
1221
1222       /* Get a conservative estimate of the collapse_mask's for the bands we're
1223           going to be folding from. */
1224       if (lowband_offset != 0 && (spread!=SPREAD_AGGRESSIVE || B>1 || tf_change<0))
1225       {
1226          int fold_start;
1227          int fold_end;
1228          int fold_i;
1229          /* This ensures we never repeat spectral content within one band */
1230          effective_lowband = IMAX(M*eBands[start], M*eBands[lowband_offset]-N);
1231          fold_start = lowband_offset;
1232          while(M*eBands[--fold_start] > effective_lowband);
1233          fold_end = lowband_offset-1;
1234          while(M*eBands[++fold_end] < effective_lowband+N);
1235          x_cm = y_cm = 0;
1236          fold_i = fold_start; do {
1237            x_cm |= collapse_masks[fold_i*C+0];
1238            y_cm |= collapse_masks[fold_i*C+C-1];
1239          } while (++fold_i<fold_end);
1240       }
1241       /* Otherwise, we'll be using the LCG to fold, so all blocks will (almost
1242           always) be non-zero.*/
1243       else
1244          x_cm = y_cm = (1<<B)-1;
1245
1246       if (dual_stereo && i==intensity)
1247       {
1248          int j;
1249
1250          /* Switch off dual stereo to do intensity */
1251          dual_stereo = 0;
1252          for (j=M*eBands[start];j<M*eBands[i];j++)
1253             norm[j] = HALF32(norm[j]+norm2[j]);
1254       }
1255       if (dual_stereo)
1256       {
1257          x_cm = quant_band(encode, m, i, X, NULL, N, b/2, spread, B, intensity, tf_change,
1258                effective_lowband != -1 ? norm+effective_lowband : NULL, resynth, ec, &remaining_bits, LM,
1259                norm+M*eBands[i], bandE, 0, seed, Q15ONE, lowband_scratch, x_cm);
1260          y_cm = quant_band(encode, m, i, Y, NULL, N, b/2, spread, B, intensity, tf_change,
1261                effective_lowband != -1 ? norm2+effective_lowband : NULL, resynth, ec, &remaining_bits, LM,
1262                norm2+M*eBands[i], bandE, 0, seed, Q15ONE, lowband_scratch, y_cm);
1263       } else {
1264          x_cm = quant_band(encode, m, i, X, Y, N, b, spread, B, intensity, tf_change,
1265                effective_lowband != -1 ? norm+effective_lowband : NULL, resynth, ec, &remaining_bits, LM,
1266                norm+M*eBands[i], bandE, 0, seed, Q15ONE, lowband_scratch, x_cm|y_cm);
1267          y_cm = x_cm;
1268       }
1269       collapse_masks[i*C+0] = (unsigned char)x_cm;
1270       collapse_masks[i*C+C-1] = (unsigned char)y_cm;
1271       balance += pulses[i] + tell;
1272
1273       /* Update the folding position only as long as we have 1 bit/sample depth */
1274       update_lowband = b>(N<<BITRES);
1275    }
1276    RESTORE_STACK;
1277 }
1278