Ticket #2686: aac-improvements-wip-v3-vbr.patch

libavcodec/aaccoder.c

@@ -57 +57 @@
     run_value_bits_long, run_value_bits_short
 };
 
-
 /**
  * Quantize one coefficient.
  * @return absolute value of the quantized coefficient
@@ -66 +65 @@
 static av_always_inline int quant(float coef, const float Q)
 {
     float a = coef * Q;
-    return sqrtf(a * sqrtf(a)) + 0.4054;
+    return sqrtf(a * sqrtf(a)) + 0.4054f;
 }
 
 static void quantize_bands(int *out, const float *in, const float *scaled,
@@ -76 +75 @@
     double qc;
     for (i = 0; i < size; i++) {
         qc = scaled[i] * Q34;
-        out[i] = (int)FFMIN(qc + 0.4054, (double)maxval);
+        out[i] = (int)FFMIN(qc + 0.4054f, (double)maxval);
         if (is_signed && in[i] < 0.0f) {
             out[i] = -out[i];
         }
@@ -282 +281 @@
     return maxval;
 }
 
+static float find_max_absval(int group_len, int swb_size, const float *scaled) {
+    float maxval = 0.0f;
+    int w2, i;
+    for (w2 = 0; w2 < group_len; w2++) {
+        for (i = 0; i < swb_size; i++) {
+            maxval = FFMAX(maxval, fabs(scaled[w2*128+i]));
+        }
+    }
+    return maxval;
+}
+
 static int find_min_book(float maxval, int sf) {
     float Q = ff_aac_pow2sf_tab[POW_SF2_ZERO - sf + SCALE_ONE_POS - SCALE_DIV_512];
     float Q34 = sqrtf(Q * sqrtf(Q));
@@ -701 +711 @@
                 sce->sf_idx[(w+w2)*16+g] = sce->sf_idx[w*16+g];
 }
 
+#define sclip(x) av_clip(x,60,218)
+
 /**
  * two-loop quantizers search taken from ISO 13818-7 Appendix C
  */
 static void search_for_quantizers_twoloop(AVCodecContext *avctx,
                                           AACEncContext *s,
                                           SingleChannelElement *sce,
-                                          const float lambda)
+                                          float lambda)
 {
     int start = 0, i, w, w2, g;
-    int destbits = avctx->bit_rate * 1024.0 / avctx->sample_rate / avctx->channels * (lambda / 120.f);
-    float dists[128] = { 0 }, uplims[128];
+    int destbits = avctx->bit_rate * 1024.0 / avctx->sample_rate 
+        / ((avctx->flags & CODEC_FLAG_QSCALE) ? 2.0f : avctx->channels) 
+        * (lambda / 120.f);
+    int refbits = destbits;
+    int toomanybits, toofewbits;
+    float dists[128] = { 0 }, uplims[128], energies[128];
     float maxvals[128];
-    int fflag, minscaler;
+    
+    /* rdlambda controls the maximum tolerated distortion. Twoloop
+     * will keep iterating until it fails to lower it or it reaches
+     * ulimit * rdlambda. Keeping it low increases quality on difficult
+     * signals, but lower it too much, and bits will be taken from weak
+     * signals, creating "holes". A balance is necesary. 
+     * rdmax and rdmin specify the relative deviation from rdlambda
+     * allowed for tonality compensation
+     */
+    float rdlambda = av_clipf(2 * 120.f / lambda, 0.0625f, 16.0f);
+    float rdmin = 0.0625f;
+    float rdmax = 1.0f;
+    
+    /* sfoffs controls an offset of optmium allocation that will be
+     * applied based on lambda. Keep it real and modest, the loop
+     * will take care of the rest, this just accelerates convergence
+     */
+    float sfoffs = av_clipf(log2f(120.0f / lambda) * 4.0f, -5, 10);
+    
+    int fflag, minscaler, nminscaler, minrdsf;
     int its  = 0;
+    int maxits = 20;
     int allz = 0;
-    float minthr = INFINITY;
+    int tbits;
+    int cutoff = avctx->sample_rate/2;
+    
+    /* zeroscale controls a multiplier of the threshold, if band energy
+     * is below this, a zero is forced. Keep it lower than 1, unless
+     * low lambda is used, because energy < threshold doesn't mean there's
+     * no audible signal outright, it's just energy. Also make it rise
+     * slower than rdlambda, as rdscale has due compensation with
+     * noisy band depriorization below, whereas zeroing logic is rather dumb
+     */
+    float zeroscale;
+    if (lambda > 120.f)
+        zeroscale = av_clipf(powf(120.f / lambda, 0.25f), 0.0625f, 1.0f);
+    else
+        zeroscale = 1.f;
+    
+    if (s->psy.bitres.alloc >= 0) {
+        // Psy granted us extra bits to use, from the reservoire
+        // adjust for lambda except what psy already did
+        destbits = s->psy.bitres.alloc
+            * (lambda / (avctx->global_quality ? avctx->global_quality : 120)); 
+    }
+    
+    if (avctx->flags & CODEC_FLAG_QSCALE) {
+        // When using a constant Q-scale, be lenient on bit under/overflow
+        toomanybits = FFMAX3(refbits * 2, destbits * 8, 1280);
+        toofewbits = destbits / 2;
+        
+        // search further
+        maxits = 40;
+        
+        // and zero out above cutoff frequency
+        {
+            int i, w;
+            int frame_bit_rate = (destbits * 1.4f * avctx->sample_rate / 1024) * (lambda / 40.0f);
+            
+            /* For high lambdas, don't cut off */
+            int bandwidth = FFMAX(3000, _AAC_CUTOFF(frame_bit_rate, 1, avctx->sample_rate));
+            int wlen = 1024 / sce->ics.num_windows;
+            cutoff = bandwidth * 2 * wlen / avctx->sample_rate;
+        }
+        
+    } else {
+        // When using ABR, be strict
+        toomanybits = destbits + destbits/16;
+        toofewbits = destbits - destbits/16;
+    }
 
     // for values above this the decoder might end up in an endless loop
     // due to always having more bits than what can be encoded.
     destbits = FFMIN(destbits, 5800);
+    toomanybits = FFMIN(toomanybits, 5800);
+    toofewbits = FFMIN(toofewbits, 5800);
     //XXX: some heuristic to determine initial quantizers will reduce search time
     //determine zero bands and upper limits
     for (w = 0; w < sce->ics.num_windows; w += sce->ics.group_len[w]) {
-        for (g = 0;  g < sce->ics.num_swb; g++) {
+        for (g = start = 0;  g < sce->ics.num_swb; start += sce->ics.swb_sizes[g++]) {
             int nz = 0;
-            float uplim = 0.0f;
+            float uplim = INFINITY;
+            float energy = 0.0f;
+            float lpfactor = FFMAX(1.0f, (float)start / cutoff);
+            
             for (w2 = 0; w2 < sce->ics.group_len[w]; w2++) {
                 FFPsyBand *band = &s->psy.ch[s->cur_channel].psy_bands[(w+w2)*16+g];
-                uplim += band->threshold;
-                if (band->energy <= band->threshold || band->threshold == 0.0f) {
+                float th = band->threshold * lpfactor;
+                if (band->energy <= th * zeroscale || band->threshold == 0.0f) {
                     sce->zeroes[(w+w2)*16+g] = 1;
                     continue;
                 }
                 nz = 1;
             }
-            uplims[w*16+g] = uplim *512;
+            if (!nz) {
+                uplim = 0.0f;
+            } else {
+                for (w2 = 0; w2 < sce->ics.group_len[w]; w2++) {
+                    FFPsyBand *band = &s->psy.ch[s->cur_channel].psy_bands[(w+w2)*16+g];
+                    float th = band->threshold * lpfactor;
+                    if (band->energy <= band->threshold * zeroscale || band->threshold == 0.0f)
+                        continue;
+                    if (uplim > th) {
+                        uplim = th;
+                        energy = band->energy;
+                    }
+                }
+            }
+            uplims[w*16+g] = uplim;
+            energies[w*16+g] = energy;
             sce->zeroes[w*16+g] = !nz;
-            if (nz)
-                minthr = FFMIN(minthr, uplim);
             allz |= nz;
         }
     }
+    
+    /* Compute initial scalers */
+    minscaler = 65535;
     for (w = 0; w < sce->ics.num_windows; w += sce->ics.group_len[w]) {
         for (g = 0;  g < sce->ics.num_swb; g++) {
             if (sce->zeroes[w*16+g]) {
                 sce->sf_idx[w*16+g] = SCALE_ONE_POS;
                 continue;
             }
-            sce->sf_idx[w*16+g] = SCALE_ONE_POS + FFMIN(log2f(uplims[w*16+g]/minthr)*4,59);
+            /* log2f-to-distortion ratio is, technically, 2 (1.5db = 4, but it's power vs level so it's 2). 
+             * But, as offsets are applied, low-frequency signals are too sensitive to the induced distortion,
+             * so we make scaling more conservative by choosing a lower log2f-to-distortion ratio, and thus
+             * more robust.
+             */
+            sce->sf_idx[w*16+g] = av_clip(
+                SCALE_ONE_POS 
+                    + 1.75*log2f(FFMAX(0.00125f,uplims[w*16+g]) / sce->ics.swb_sizes[g]) 
+                    + sfoffs,
+                60, SCALE_MAX_POS);
+            //fprintf(stderr, "%02x ", sce->sf_idx[w*16+g]);
+            minscaler = FFMIN(minscaler, sce->sf_idx[w*16+g]);
         }
+        //fprintf(stderr, "|\n");
     }
-
+    //fprintf(stderr, "\n");
+    
+    /* Clip */
+    minscaler = av_clip(minscaler, SCALE_ONE_POS - SCALE_DIV_512, SCALE_MAX_POS - SCALE_DIV_512);
+    for (w = 0; w < sce->ics.num_windows; w += sce->ics.group_len[w]) 
+        for (g = 0;  g < sce->ics.num_swb; g++) 
+            if (!sce->zeroes[w*16+g]) 
+                sce->sf_idx[w*16+g] = av_clip(sce->sf_idx[w*16+g], minscaler, minscaler + SCALE_MAX_DIFF - 1);
+    
     if (!allz)
         return;
     abs_pow34_v(s->scoefs, sce->coeffs, 1024);
@@ -766 +889 @@
         }
     }
 
+    /* Scale uplims to match rate distortion to quality 
+     * and apply noisy band depriorization and tonal band priorization.
+     * Maxval-energy ratio gives us an idea of how noisy/tonal the band is. 
+     * If maxval^2 ~ energy, then that band is mostly noise, and we can relax
+     * rate distortion requirements.
+     */
+    for (w = 0; w < sce->ics.num_windows; w += sce->ics.group_len[w]) {
+        start = w*128;
+        for (g = 0;  g < sce->ics.num_swb; g++) {
+            float max = find_max_absval(sce->ics.group_len[w], sce->ics.swb_sizes[g], sce->coeffs + start);
+            if (max > 0) {
+                float energy2uplim = energies[w*16+g] / (max*max*sce->ics.swb_sizes[g]);
+                energy2uplim = FFMAX(0.0625f, FFMIN(1.0f,energy2uplim));
+                uplims[w*16+g] *= av_clipf(rdlambda * rdlambda * energy2uplim, rdmin, rdmax);
+                start += sce->ics.swb_sizes[g];
+            }
+        }
+    }
+
     //perform two-loop search
     //outer loop - improve quality
     do {
-        int tbits, qstep;
-        minscaler = sce->sf_idx[0];
+        int qstep;
         //inner loop - quantize spectrum to fit into given number of bits
         qstep = its ? 1 : 32;
         do {
@@ -790 +931 @@
                         start += sce->ics.swb_sizes[g];
                         continue;
                     }
-                    minscaler = FFMIN(minscaler, sce->sf_idx[w*16+g]);
                     cb = find_min_book(maxvals[w*16+g], sce->sf_idx[w*16+g]);
                     for (w2 = 0; w2 < sce->ics.group_len[w]; w2++) {
                         int b;
@@ -813 +953 @@
                     prev = sce->sf_idx[w*16+g];
                 }
             }
-            if (tbits > destbits) {
+            if (tbits > toomanybits) {
                 for (i = 0; i < 128; i++)
-                    if (sce->sf_idx[i] < 218 - qstep)
-                        sce->sf_idx[i] += qstep;
-            } else {
+                    if (sce->sf_idx[i] < (SCALE_MAX_POS - SCALE_DIV_512))
+                        sce->sf_idx[i] = FFMIN(SCALE_MAX_POS, sce->sf_idx[i] + qstep);
+            } else if (tbits < toofewbits) {
                 for (i = 0; i < 128; i++)
-                    if (sce->sf_idx[i] > 60 - qstep)
-                        sce->sf_idx[i] -= qstep;
+                    if (sce->sf_idx[i] > SCALE_ONE_POS)
+                        sce->sf_idx[i] = FFMAX(SCALE_ONE_POS, sce->sf_idx[i] - qstep);
             }
             qstep >>= 1;
-            if (!qstep && tbits > destbits*1.02 && sce->sf_idx[0] < 217)
+            if (!qstep && tbits > toomanybits && sce->sf_idx[0] < 217)
                 qstep = 1;
         } while (qstep);
 
+        minscaler = SCALE_MAX_POS;
+        for (w = 0; w < sce->ics.num_windows; w += sce->ics.group_len[w])
+            for (g = 0;  g < sce->ics.num_swb; g++)
+                if (!sce->zeroes[w*16+g])
+                    minscaler = FFMIN(minscaler, sce->sf_idx[w*16+g]);
+
         fflag = 0;
-        minscaler = av_clip(minscaler, 60, 255 - SCALE_MAX_DIFF);
+        minscaler = nminscaler = av_clip(minscaler, SCALE_ONE_POS - SCALE_DIV_512, SCALE_MAX_POS - SCALE_DIV_512);
+        minrdsf = (avctx->flags & CODEC_FLAG_QSCALE) ? 60 : minscaler;
         for (w = 0; w < sce->ics.num_windows; w += sce->ics.group_len[w]) {
             for (g = 0; g < sce->ics.num_swb; g++) {
                 int prevsc = sce->sf_idx[w*16+g];
-                if (dists[w*16+g] > uplims[w*16+g] && sce->sf_idx[w*16+g] > 60) {
+                //float hfprio = 1.0 + its * av_clipf(120.0f / lambda - 1.0f, 0.0f, 8.0f) * g / (sce->ics.num_swb * maxits);
+                if (dists[w*16+g] > (uplims[w*16+g]/* * hfprio*/) && sce->sf_idx[w*16+g] > minrdsf) {
                     if (find_min_book(maxvals[w*16+g], sce->sf_idx[w*16+g]-1))
                         sce->sf_idx[w*16+g]--;
                     else //Try to make sure there is some energy in every band
                         sce->sf_idx[w*16+g]-=2;
                 }
-                sce->sf_idx[w*16+g] = av_clip(sce->sf_idx[w*16+g], minscaler, minscaler + SCALE_MAX_DIFF);
-                sce->sf_idx[w*16+g] = FFMIN(sce->sf_idx[w*16+g], 219);
+                sce->sf_idx[w*16+g] = av_clip(sce->sf_idx[w*16+g], minrdsf, minscaler + SCALE_MAX_DIFF);
+                sce->sf_idx[w*16+g] = FFMIN(sce->sf_idx[w*16+g], SCALE_MAX_POS - SCALE_DIV_512);
                 if (sce->sf_idx[w*16+g] != prevsc)
                     fflag = 1;
+                nminscaler = FFMIN(nminscaler, sce->sf_idx[w*16+g]);
                 sce->band_type[w*16+g] = find_min_book(maxvals[w*16+g], sce->sf_idx[w*16+g]);
             }
         }
+        if (nminscaler < minscaler) {
+            // Drecreased some scalers below minscaler. Must re-clamp.
+            minscaler = nminscaler;
+            for (w = 0; w < sce->ics.num_windows; w += sce->ics.group_len[w]) 
+                for (g = 0; g < sce->ics.num_swb; g++) 
+                    sce->sf_idx[w*16+g] = av_clip(sce->sf_idx[w*16+g], minscaler, minscaler + SCALE_MAX_DIFF);
+        }
         its++;
-    } while (fflag && its < 10);
+    } while (fflag && its < maxits);
+    
+    /* Fill implicit zeroes */
+    for (w = 0; w < sce->ics.num_windows; w += sce->ics.group_len[w]) {
+        /* No lowly codebooks beyond cutoff zone, to clean up noisy coefs */
+        int scutoff = cutoff + cutoff/5;
+        for (g = start = 0; g < sce->ics.num_swb; start += sce->ics.swb_sizes[g++]) {
+            int minbt = (start < scutoff) ? 0 : 7;
+            if (sce->band_type[w*16+g] <= minbt) {
+                sce->zeroes[w*16+g] = 1;
+                sce->band_type[w*16+g] = 0;
+            }
+        }
+    }
+    //fprintf(stderr, "ba:%d br:%d  \t\r", s->psy.bitres.alloc, tbits);
 }
 
 static void search_for_quantizers_faac(AVCodecContext *avctx, AACEncContext *s,
@@ -1021 +1191 @@
                                        SingleChannelElement *sce,
                                        const float lambda)
 {
-    int i, w, w2, g;
-    int minq = 255;
+    int w, w2, g;
+    float lowlambda = av_clipf(120.f / lambda, 0.85f, 1.f);
+    float rlambda = av_clipf(120.f / lambda, 0.75f, 10.f);
+    const int minq = av_clip(2 * log2f(120.f / lambda) + 150, 100, 218 - SCALE_MAX_DIFF);
+    const int maxq = minq + SCALE_MAX_DIFF - 1;
 
     memset(sce->sf_idx, 0, sizeof(sce->sf_idx));
     for (w = 0; w < sce->ics.num_windows; w += sce->ics.group_len[w]) {
         for (g = 0; g < sce->ics.num_swb; g++) {
             for (w2 = 0; w2 < sce->ics.group_len[w]; w2++) {
                 FFPsyBand *band = &s->psy.ch[s->cur_channel].psy_bands[(w+w2)*16+g];
-                if (band->energy <= band->threshold) {
-                    sce->sf_idx[(w+w2)*16+g] = 218;
+                if (band->energy <= 0.05 * lowlambda * band->threshold) {
+                    sce->sf_idx[(w+w2)*16+g] = maxq;
                     sce->zeroes[(w+w2)*16+g] = 1;
                 } else {
-                    sce->sf_idx[(w+w2)*16+g] = av_clip(SCALE_ONE_POS - SCALE_DIV_512 + log2f(band->threshold), 80, 218);
+                    sce->sf_idx[(w+w2)*16+g] = av_clip(SCALE_ONE_POS - SCALE_DIV_512 + 1.414*log2f(band->threshold * rlambda), minq, maxq);
                     sce->zeroes[(w+w2)*16+g] = 0;
                 }
-                minq = FFMIN(minq, sce->sf_idx[(w+w2)*16+g]);
             }
         }
     }
-    for (i = 0; i < 128; i++) {
-        sce->sf_idx[i] = 140;
-        //av_clip(sce->sf_idx[i], minq, minq + SCALE_MAX_DIFF - 1);
-    }
     //set the same quantizers inside window groups
-    for (w = 0; w < sce->ics.num_windows; w += sce->ics.group_len[w])
-        for (g = 0;  g < sce->ics.num_swb; g++)
-            for (w2 = 1; w2 < sce->ics.group_len[w]; w2++)
-                sce->sf_idx[(w+w2)*16+g] = sce->sf_idx[w*16+g];
+    if (sce->ics.num_windows > 1) {
+        for (w = 0; w < sce->ics.num_windows; w += sce->ics.group_len[w]) {
+            for (g = 0;  g < sce->ics.num_swb; g++) {
+                if (sce->ics.group_len[w] > 1) {
+                    int avg_sf_idx = 0;
+                    for (w2 = 0; w2 < sce->ics.group_len[w]; w2++)
+                        avg_sf_idx += sce->sf_idx[w*16+g];
+                    avg_sf_idx /= sce->ics.group_len[w];
+                    for (w2 = 1; w2 < sce->ics.group_len[w]; w2++)
+                        sce->sf_idx[(w+w2)*16+g] = avg_sf_idx;
+                }
+            }
+        }
+    }
+}
+
+static float bval2bmax(float b)
+{
+    /* approximates exp10f(-3.0f*(0.5f + 0.5f * cosf(FFMIN(b,15.5f) / 15.5f))) */
+    return 0.001f + 0.0035f * (b*b*b) / (15.5f*15.5f*15.5f);
 }
 
 static void search_for_ms(AACEncContext *s, ChannelElement *cpe,
@@ -1059 +1243 @@
     float *L34 = s->scoefs, *R34 = s->scoefs + 128, *M34 = s->scoefs + 128*2, *S34 = s->scoefs + 128*3;
     SingleChannelElement *sce0 = &cpe->ch[0];
     SingleChannelElement *sce1 = &cpe->ch[1];
+    
     if (!cpe->common_window)
         return;
+    
     for (w = 0; w < sce0->ics.num_windows; w += sce0->ics.group_len[w]) {
+        int min_sf_idx_mid = SCALE_MAX_POS;
+        int min_sf_idx_side = SCALE_MAX_POS;
+        for (g = 0; g < sce0->ics.num_swb; g++) {
+            if (!sce0->zeroes[w*16+g])
+                min_sf_idx_mid = FFMIN(min_sf_idx_mid, sce0->sf_idx[w*16+g]);
+            if (!sce1->zeroes[w*16+g])
+                min_sf_idx_side = FFMIN(min_sf_idx_side, sce1->sf_idx[w*16+g]);
+        }
+        
         for (g = 0;  g < sce0->ics.num_swb; g++) {
+            float bmax = bval2bmax(g * 17.0f / sce0->ics.num_swb) / 0.0045f;
             if (!cpe->ch[0].zeroes[w*16+g] && !cpe->ch[1].zeroes[w*16+g]) {
                 float dist1 = 0.0f, dist2 = 0.0f;
+                int B0 = 0, B1 = 0;
+                int minidx;
+                int mididx, sididx;
+                float Mmax = 0.0f, Smax = 0.0f;
+                int midcb, sidcb;
+                
+                /* Must compute mid/side SF and book for the whole window group */
+                minidx = FFMIN(sce0->sf_idx[w*16+g], sce1->sf_idx[w*16+g]);
+                for (w2 = 0; w2 < sce0->ics.group_len[w]; w2++) {
+                    for (i = 0; i < sce0->ics.swb_sizes[g]; i++) {
+                        M[i] = (sce0->coeffs[start+w2*128+i]
+                              + sce1->coeffs[start+w2*128+i]) * 0.5;
+                        S[i] =  M[i]
+                              - sce1->coeffs[start+w2*128+i];
+                    }
+                    abs_pow34_v(M34, M, sce0->ics.swb_sizes[g]);
+                    abs_pow34_v(S34, S, sce0->ics.swb_sizes[g]);
+                    for (i = 0; i < sce0->ics.swb_sizes[g]; i++ ) {
+                        Mmax = FFMAX(Mmax, M34[i]);
+                        Smax = FFMAX(Smax, S34[i]);
+                    }
+                }
+                mididx = av_clip(minidx, min_sf_idx_mid, min_sf_idx_mid + SCALE_MAX_DIFF);
+                sididx = av_clip(minidx, min_sf_idx_side, min_sf_idx_side + SCALE_MAX_DIFF);
+                midcb = find_min_book(Mmax, mididx);
+                sidcb = find_min_book(Smax, sididx);
+                
                 for (w2 = 0; w2 < sce0->ics.group_len[w]; w2++) {
                     FFPsyBand *band0 = &s->psy.ch[s->cur_channel+0].psy_bands[(w+w2)*16+g];
                     FFPsyBand *band1 = &s->psy.ch[s->cur_channel+1].psy_bands[(w+w2)*16+g];
                     float minthr = FFMIN(band0->threshold, band1->threshold);
-                    float maxthr = FFMAX(band0->threshold, band1->threshold);
+                    int b1,b2,b3,b4;
                     for (i = 0; i < sce0->ics.swb_sizes[g]; i++) {
                         M[i] = (sce0->coeffs[start+w2*128+i]
                               + sce1->coeffs[start+w2*128+i]) * 0.5;
@@ -1085 +1308 @@
                                                 sce0->ics.swb_sizes[g],
                                                 sce0->sf_idx[(w+w2)*16+g],
                                                 sce0->band_type[(w+w2)*16+g],
-                                                lambda / band0->threshold, INFINITY, NULL);
+                                                lambda / band0->threshold, INFINITY, &b1);
                     dist1 += quantize_band_cost(s, sce1->coeffs + start + w2*128,
                                                 R34,
                                                 sce1->ics.swb_sizes[g],
                                                 sce1->sf_idx[(w+w2)*16+g],
                                                 sce1->band_type[(w+w2)*16+g],
-                                                lambda / band1->threshold, INFINITY, NULL);
+                                                lambda / band1->threshold, INFINITY, &b2);
                     dist2 += quantize_band_cost(s, M,
                                                 M34,
                                                 sce0->ics.swb_sizes[g],
-                                                sce0->sf_idx[(w+w2)*16+g],
-                                                sce0->band_type[(w+w2)*16+g],
-                                                lambda / maxthr, INFINITY, NULL);
+                                                mididx,
+                                                midcb,
+                                                lambda / minthr, INFINITY, &b3);
                     dist2 += quantize_band_cost(s, S,
                                                 S34,
                                                 sce1->ics.swb_sizes[g],
-                                                sce1->sf_idx[(w+w2)*16+g],
-                                                sce1->band_type[(w+w2)*16+g],
-                                                lambda / minthr, INFINITY, NULL);
+                                                sididx,
+                                                sidcb,
+                                                lambda * (lambda / 120.0f) / (minthr * bmax), INFINITY, &b4);
+                    B0 += b1+b2;
+                    B1 += b3+b4;
+                    dist1 -= B0;
+                    dist2 -= B1;
+                }
+                cpe->ms_mask[w*16+g] = dist2 <= dist1 && B1 < B0;
+                if (cpe->ms_mask[w*16+g]) {
+                    for (w2 = 0; w2 < sce0->ics.group_len[w]; w2++) {
+                        sce0->sf_idx[(w+w2)*16+g] = mididx;
+                        sce0->band_type[(w+w2)*16+g] = midcb;
+                        sce1->sf_idx[(w+w2)*16+g] = sididx;
+                        sce1->band_type[(w+w2)*16+g] = sidcb;
+                    }
                 }
-                cpe->ms_mask[w*16+g] = dist2 < dist1;
+            } else {
+                cpe->ms_mask[w*16+g] = 0;
             }
             start += sce0->ics.swb_sizes[g];
         }

libavcodec/aacenc.c

@@ -317 +317 @@
                 if (cpe->common_window && !ch && cpe->ms_mask[w + g]) {
                     for (i = 0; i < ics->swb_sizes[g]; i++) {
                         cpe->ch[0].coeffs[start+i] = (cpe->ch[0].coeffs[start+i] + cpe->ch[1].coeffs[start+i]) / 2.0;
-                        cpe->ch[1].coeffs[start+i] =  cpe->ch[0].coeffs[start+i] - cpe->ch[1].coeffs[start+i];
+                        cpe->ch[1].coeffs[start+i] = cpe->ch[0].coeffs[start+i] - cpe->ch[1].coeffs[start+i];
                     }
                 }
                 start += ics->swb_sizes[g];
@@ -507 +507 @@
     AACEncContext *s = avctx->priv_data;
     float **samples = s->planar_samples, *samples2, *la, *overlap;
     ChannelElement *cpe;
-    int i, ch, w, g, chans, tag, start_ch, ret;
+    int i, ch, w, g, chans, tag, start_ch, ret, frame_bits, its;
     int chan_el_counter[4];
     FFPsyWindowInfo windows[AAC_MAX_CHANNELS];
 
@@ -572 +572 @@
     }
     if ((ret = ff_alloc_packet2(avctx, avpkt, 8192 * s->channels)) < 0)
         return ret;
+    
+    frame_bits = 0;
+    its = 0;
     do {
-        int frame_bits;
-
+        int target_bits, too_many_bits, too_few_bits;
+        
         init_put_bits(&s->pb, avpkt->data, avpkt->size);
 
         if ((avctx->frame_number & 0xFF)==1 && !(avctx->flags & CODEC_FLAG_BITEXACT))
             put_bitstream_info(s, LIBAVCODEC_IDENT);
         start_ch = 0;
+        target_bits = 0;
         memset(chan_el_counter, 0, sizeof(chan_el_counter));
         for (i = 0; i < s->chan_map[0]; i++) {
             FFPsyWindowInfo* wi = windows + start_ch;
@@ -591 +595 @@
             put_bits(&s->pb, 4, chan_el_counter[tag]++);
             for (ch = 0; ch < chans; ch++)
                 coeffs[ch] = cpe->ch[ch].coeffs;
+            s->psy.bitres.alloc = -1;
             s->psy.model->analyze(&s->psy, start_ch, coeffs, wi);
+            target_bits += s->psy.bitres.alloc * (120.f / s->lambda);
             for (ch = 0; ch < chans; ch++) {
                 s->cur_channel = start_ch + ch;
                 s->coder->search_for_quantizers(avctx, s, &cpe->ch[ch], s->lambda);
@@ -635 +641 @@
             start_ch += chans;
         }
 
+        if (avctx->flags & CODEC_FLAG_QSCALE) {
+            // When using a constant Q-scale, don't mess with lambda
+            break;
+        }
+        
+        // rate control stuff
+        // target either the nominal bitrate, or what psy's bit reservoir says to target
+        // whichever is greatest
         frame_bits = put_bits_count(&s->pb);
-        if (frame_bits <= 6144 * s->channels - 3) {
-            s->psy.bitres.bits = frame_bits / s->channels;
+        target_bits = FFMAX(target_bits, avctx->bit_rate * 1024 / avctx->sample_rate);
+        target_bits = FFMIN(target_bits, 6144 * s->channels - 3);
+        
+        // When using ABR, be strict (but only for increasing)
+        too_many_bits = target_bits + target_bits/2;
+        too_few_bits = target_bits - target_bits/8;
+        
+        if (   its == 0 /* for steady-state Q-scale tracking */
+            || (its < 5 && (frame_bits < too_few_bits || frame_bits > too_many_bits))
+            || frame_bits >= 6144 * s->channels - 3  ) 
+        {
+            float prev_lambda = s->lambda;
+            float ratio = ((float)target_bits) / frame_bits;
+            s->lambda = FFMIN(s->lambda * ratio, 65536.f);
+            
+            if (prev_lambda == s->lambda)
+                break;
+            
+            // Keep iterating if we must reduce and lambda is in the sky
+            if (ratio > 0.9f || s->lambda <= 300.f)
+                its++;
+            
+            if (frame_bits >= too_few_bits && frame_bits <= too_many_bits) {
+                /*
+                This path is for steady-state Q-scale tracking
+                When frame bits fall within the stable range, we still need to adjust
+                lambda to maintain it like so in a stable fashion (large jumps in lambda
+                create artifacts and shoulda be avoided)
+                */
+                break;
+            }
+        } else {
             break;
         }
-
-        s->lambda *= avctx->bit_rate * 1024.0f / avctx->sample_rate / frame_bits;
-
     } while (1);
 
     put_bits(&s->pb, 3, TYPE_END);
     flush_put_bits(&s->pb);
     avctx->frame_bits = put_bits_count(&s->pb);
-
-    // rate control stuff
-    if (!(avctx->flags & CODEC_FLAG_QSCALE)) {
-        float ratio = avctx->bit_rate * 1024.0f / avctx->sample_rate / avctx->frame_bits;
-        s->lambda *= ratio;
-        s->lambda = FFMIN(s->lambda, 65536.f);
-    }
+    
+    if (avctx->flags & CODEC_FLAG_QSCALE)
+        // Lie, we probably used more bits due to RD adjustment, but we must not try to save those bits later
+        s->psy.bitres.bits = s->psy.bitres.alloc ? s->psy.bitres.alloc : (avctx->frame_bits / s->channels);
+    else
+        s->psy.bitres.bits = avctx->frame_bits / s->channels;
 
     if (!frame)
         s->last_frame++;

libavcodec/aacpsy.c

@@ -24 +24 @@
  * AAC encoder psychoacoustic model
  */
 
-#include "libavutil/attributes.h"
 #include "libavutil/libm.h"
 
 #include "avcodec.h"
@@ -255 +254 @@
 /**
  * LAME psy model specific initialization
  */
-static av_cold void lame_window_init(AacPsyContext *ctx, AVCodecContext *avctx)
-{
+static void lame_window_init(AacPsyContext *ctx, AVCodecContext *avctx) {
     int i, j;
 
     for (i = 0; i < avctx->channels; i++) {
@@ -299 +297 @@
     float bark;
     int i, j, g, start;
     float prev, minscale, minath, minsnr, pe_min;
-    const int chan_bitrate = ctx->avctx->bit_rate / ctx->avctx->channels;
+    int chan_bitrate = ctx->avctx->bit_rate / ((ctx->avctx->flags & CODEC_FLAG_QSCALE) ? 2.0f : ctx->avctx->channels);
     const int bandwidth    = ctx->avctx->cutoff ? ctx->avctx->cutoff : AAC_CUTOFF(ctx->avctx);
     const float num_bark   = calc_bark((float)bandwidth);
 
     ctx->model_priv_data = av_mallocz(sizeof(AacPsyContext));
     pctx = (AacPsyContext*) ctx->model_priv_data;
 
+    if (ctx->avctx->flags & CODEC_FLAG_QSCALE) {
+        /* Use the target average bitrate to compute spread parameters */
+        chan_bitrate = (int)(chan_bitrate / 120.0 * (ctx->avctx->global_quality ? ctx->avctx->global_quality : 120));
+    }
+    
     pctx->chan_bitrate = chan_bitrate;
-    pctx->frame_bits   = chan_bitrate * AAC_BLOCK_SIZE_LONG / ctx->avctx->sample_rate;
+    pctx->frame_bits   = FFMIN(2560, chan_bitrate * AAC_BLOCK_SIZE_LONG / ctx->avctx->sample_rate);
     pctx->pe.min       =  8.0f * AAC_BLOCK_SIZE_LONG * bandwidth / (ctx->avctx->sample_rate * 2.0f);
     pctx->pe.max       = 12.0f * AAC_BLOCK_SIZE_LONG * bandwidth / (ctx->avctx->sample_rate * 2.0f);
     ctx->bitres.size   = 6144 - pctx->frame_bits;
     ctx->bitres.size  -= ctx->bitres.size % 8;
     pctx->fill_level   = ctx->bitres.size;
     minath = ath(3410, ATH_ADD);
+    
     for (j = 0; j < 2; j++) {
         AacPsyCoeffs *coeffs = pctx->psy_coef[j];
         const uint8_t *band_sizes = ctx->bands[j];
@@ -391 +395 @@
                                                  int channel, int prev_type)
 {
     int i, j;
-    int br               = ctx->avctx->bit_rate / ctx->avctx->channels;
+    int br               = ((AacPsyContext*)ctx->model_priv_data)->chan_bitrate;
     int attack_ratio     = br <= 16000 ? 18 : 10;
     AacPsyContext *pctx = (AacPsyContext*) ctx->model_priv_data;
     AacPsyChannel *pch  = &pctx->ch[channel];
@@ -628 +632 @@
     const uint8_t *band_sizes  = ctx->bands[wi->num_windows == 8];
     AacPsyCoeffs  *coeffs      = pctx->psy_coef[wi->num_windows == 8];
     const float avoid_hole_thr = wi->num_windows == 8 ? PSY_3GPP_AH_THR_SHORT : PSY_3GPP_AH_THR_LONG;
-
+    
     //calculate energies, initial thresholds and related values - 5.4.2 "Threshold Calculation"
     calc_thr_3gpp(wi, num_bands, pch, band_sizes, coefs);
 
@@ -681 +685 @@
         desired_pe *= av_clipf(pctx->pe.previous / PSY_3GPP_BITS_TO_PE(ctx->bitres.bits),
                                0.85f, 1.15f);
     pctx->pe.previous = PSY_3GPP_BITS_TO_PE(desired_bits);
+    ctx->bitres.alloc = desired_bits;
 
     if (desired_pe < pe) {
         /* 5.6.1.3.4 "First Estimation of the reduction value" */

libavcodec/psymodel.c

@@ -75 +75 @@
 
 av_cold void ff_psy_end(FFPsyContext *ctx)
 {
-    if (ctx->model && ctx->model->end)
+    if (ctx->model->end)
         ctx->model->end(ctx);
     av_freep(&ctx->bands);
     av_freep(&ctx->num_bands);
@@ -106 +106 @@
 
     if (!cutoff_coeff && avctx->codec_id == AV_CODEC_ID_AAC)
         cutoff_coeff = 2.0 * AAC_CUTOFF(avctx) / avctx->sample_rate;
-
+    
     if (cutoff_coeff && cutoff_coeff < 0.98)
     ctx->fcoeffs = ff_iir_filter_init_coeffs(avctx, FF_FILTER_TYPE_BUTTERWORTH,
                                              FF_FILTER_MODE_LOWPASS, FILT_ORDER,

libavcodec/psymodel.h

@@ -27 +27 @@
 /** maximum possible number of bands */
 #define PSY_MAX_BANDS 128
 /** maximum number of channels */
-#define PSY_MAX_CHANS 20
+#define PSY_MAX_CHANS 24
 
-#define AAC_CUTOFF(s) (s->bit_rate ? FFMIN3(4000 + s->bit_rate/8, 12000 + s->bit_rate/32, s->sample_rate / 2) : (s->sample_rate / 2))
+/* cutoff for VBR is purposedly increased, since LP filtering actually
+ * hinders VBR performance rather than the opposite
+ */
+#define _AAC_CUTOFF(bit_rate,channels,sample_rate) (bit_rate ? FFMIN3(FFMIN3( \
+    bit_rate/channels/2, \
+    3000 + bit_rate/channels/4, \
+    12000 + bit_rate/channels/16), \
+    20000, \
+    sample_rate / 2): (sample_rate / 2))
+#define AAC_CUTOFF(s) ( \
+    (s->flags & CODEC_FLAG_QSCALE) \
+    ? /*_AAC_CUTOFF(((int)(480000.0f*(s->global_quality ? s->global_quality/120.0f : 1.0f))), 1, s->sample_rate)*/s->sample_rate / 2 \
+    : _AAC_CUTOFF(s->bit_rate, s->channels, s->sample_rate) \
+)
 
 /**
  * single band psychoacoustic information
@@ -88 +101 @@
     struct {
         int size;                     ///< size of the bitresevoir in bits
         int bits;                     ///< number of bits used in the bitresevoir
+        int alloc;                    ///< number of bits allocated by the psy, or -1 if no allocation was done
     } bitres;
 
     void* model_priv_data;            ///< psychoacoustic model implementation private data