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