Ticket #2686: aac-improvements-wip-v8g.patch

libavcodec/aac.h

@@ -157 +157 @@
 typedef struct IndividualChannelStream {
     uint8_t max_sfb;            ///< number of scalefactor bands per group
     enum WindowSequence window_sequence[2];
-    uint8_t use_kb_window[2];   ///< If set, use Kaiser-Bessel window, otherwise use a sine window.
+    uint8_t use_kb_window[2];   ///< If set, use Kaiser-Bessel window, otherwise use a sinus window.
     int num_window_groups;
     uint8_t group_len[8];
     LongTermPrediction ltp;
@@ -170 +170 @@
     int predictor_initialized;
     int predictor_reset_group;
     uint8_t prediction_used[41];
+    uint8_t window_clipping[8]; ///< set if a certain window is near clipping
 } IndividualChannelStream;
 
 /**
@@ -233 +234 @@
     float sf[120];                                  ///< scalefactors
     int sf_idx[128];                                ///< scalefactor indices (used by encoder)
     uint8_t zeroes[128];                            ///< band is not coded (used by encoder)
-    DECLARE_ALIGNED(32, float,   coeffs)[1024];     ///< coefficients for IMDCT
-    DECLARE_ALIGNED(32, float,   saved)[1536];      ///< overlap
+    DECLARE_ALIGNED(32, float,   pcoeffs)[1024];    ///< coefficients for IMDCT, pristine
+    DECLARE_ALIGNED(32, float,   coeffs)[1024];     ///< coefficients for IMDCT, maybe processed
+    DECLARE_ALIGNED(32, float,   saved)[1024];      ///< overlap
     DECLARE_ALIGNED(32, float,   ret_buf)[2048];    ///< PCM output buffer
     DECLARE_ALIGNED(16, float,   ltp_state)[3072];  ///< time signal for LTP
     PredictorState predictor_state[MAX_PREDICTORS];

libavcodec/aaccoder.c

@@ -57 +57 @@
     run_value_bits_long, run_value_bits_short
 };
 
+#define ROUND_STANDARD 0.4054f
+#define ROUND_TO_ZERO 0.1054f
 
 /**
  * Quantize one coefficient.
  * @return absolute value of the quantized coefficient
  * @see 3GPP TS26.403 5.6.2 "Scalefactor determination"
  */
-static av_always_inline int quant(float coef, const float Q)
+static av_always_inline int quant(float coef, const float Q, const float rounding)
 {
     float a = coef * Q;
-    return sqrtf(a * sqrtf(a)) + 0.4054;
+    return sqrtf(a * sqrtf(a)) + rounding;
 }
 
 static void quantize_bands(int *out, const float *in, const float *scaled,
-                           int size, float Q34, int is_signed, int maxval)
+                           int size, float Q34, int is_signed, int maxval, const float rounding)
 {
     int i;
     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 + rounding, (double)maxval);
         if (is_signed && in[i] < 0.0f) {
             out[i] = -out[i];
         }
@@ -108 +110 @@
                                 const float *scaled, int size, int scale_idx,
                                 int cb, const float lambda, const float uplim,
                                 int *bits, int BT_ZERO, int BT_UNSIGNED,
-                                int BT_PAIR, int BT_ESC)
+                                int BT_PAIR, int BT_ESC, const float ROUNDING)
 {
     const int q_idx = POW_SF2_ZERO - scale_idx + SCALE_ONE_POS - SCALE_DIV_512;
     const float Q   = ff_aac_pow2sf_tab [q_idx];
@@ -134 +136 @@
         abs_pow34_v(s->scoefs, in, size);
         scaled = s->scoefs;
     }
-    quantize_bands(s->qcoefs, in, scaled, size, Q34, !BT_UNSIGNED, maxval);
+    quantize_bands(s->qcoefs, in, scaled, size, Q34, !BT_UNSIGNED, maxval, ROUNDING);
     if (BT_UNSIGNED) {
         off = 0;
     } else {
@@ -161 +163 @@
                         di = t - CLIPPED_ESCAPE;
                         curbits += 21;
                     } else {
-                        int c = av_clip(quant(t, Q), 0, 8191);
+                        int c = av_clip(quant(t, Q, ROUNDING), 0, 8191);
                         di = t - c*cbrtf(c)*IQ;
                         curbits += av_log2(c)*2 - 4 + 1;
                     }
@@ -191 +193 @@
             if (BT_ESC) {
                 for (j = 0; j < 2; j++) {
                     if (ff_aac_codebook_vectors[cb-1][curidx*2+j] == 64.0f) {
-                        int coef = av_clip(quant(fabsf(in[i+j]), Q), 0, 8191);
+                        int coef = av_clip(quant(fabsf(in[i+j]), Q, ROUNDING), 0, 8191);
                         int len = av_log2(coef);
 
                         put_bits(pb, len - 4 + 1, (1 << (len - 4 + 1)) - 2);
@@ -207 +209 @@
     return cost;
 }
 
-#define QUANTIZE_AND_ENCODE_BAND_COST_FUNC(NAME, BT_ZERO, BT_UNSIGNED, BT_PAIR, BT_ESC) \
+#define QUANTIZE_AND_ENCODE_BAND_COST_FUNC(NAME, BT_ZERO, BT_UNSIGNED, BT_PAIR, BT_ESC, ROUNDING) \
 static float quantize_and_encode_band_cost_ ## NAME(                                        \
                                 struct AACEncContext *s,                                \
                                 PutBitContext *pb, const float *in,                     \
@@ -217 +219 @@
     return quantize_and_encode_band_cost_template(                                      \
                                 s, pb, in, scaled, size, scale_idx,                     \
                                 BT_ESC ? ESC_BT : cb, lambda, uplim, bits,              \
-                                BT_ZERO, BT_UNSIGNED, BT_PAIR, BT_ESC);                 \
+                                BT_ZERO, BT_UNSIGNED, BT_PAIR, BT_ESC, ROUNDING);       \
 }
 
-QUANTIZE_AND_ENCODE_BAND_COST_FUNC(ZERO,  1, 0, 0, 0)
-QUANTIZE_AND_ENCODE_BAND_COST_FUNC(SQUAD, 0, 0, 0, 0)
-QUANTIZE_AND_ENCODE_BAND_COST_FUNC(UQUAD, 0, 1, 0, 0)
-QUANTIZE_AND_ENCODE_BAND_COST_FUNC(SPAIR, 0, 0, 1, 0)
-QUANTIZE_AND_ENCODE_BAND_COST_FUNC(UPAIR, 0, 1, 1, 0)
-QUANTIZE_AND_ENCODE_BAND_COST_FUNC(ESC,   0, 1, 1, 1)
+QUANTIZE_AND_ENCODE_BAND_COST_FUNC(ZERO,  1, 0, 0, 0, ROUND_STANDARD)
+QUANTIZE_AND_ENCODE_BAND_COST_FUNC(SQUAD, 0, 0, 0, 0, ROUND_STANDARD)
+QUANTIZE_AND_ENCODE_BAND_COST_FUNC(UQUAD, 0, 1, 0, 0, ROUND_STANDARD)
+QUANTIZE_AND_ENCODE_BAND_COST_FUNC(SPAIR, 0, 0, 1, 0, ROUND_STANDARD)
+QUANTIZE_AND_ENCODE_BAND_COST_FUNC(UPAIR, 0, 1, 1, 0, ROUND_STANDARD)
+QUANTIZE_AND_ENCODE_BAND_COST_FUNC(ESC,   0, 1, 1, 1, ROUND_STANDARD)
+
+/*QUANTIZE_AND_ENCODE_BAND_COST_FUNC(ZERO_RTZ,  1, 0, 0, 0, ROUND_TO_ZERO)
+QUANTIZE_AND_ENCODE_BAND_COST_FUNC(SQUAD_RTZ, 0, 0, 0, 0, ROUND_TO_ZERO)
+QUANTIZE_AND_ENCODE_BAND_COST_FUNC(UQUAD_RTZ, 0, 1, 0, 0, ROUND_TO_ZERO)
+QUANTIZE_AND_ENCODE_BAND_COST_FUNC(SPAIR_RTZ, 0, 0, 1, 0, ROUND_TO_ZERO)
+QUANTIZE_AND_ENCODE_BAND_COST_FUNC(UPAIR_RTZ, 0, 1, 1, 0, ROUND_TO_ZERO)*/
+QUANTIZE_AND_ENCODE_BAND_COST_FUNC(ESC_RTZ,   0, 1, 1, 1, ROUND_TO_ZERO)
 
 static float (*const quantize_and_encode_band_cost_arr[])(
                                 struct AACEncContext *s,
@@ -247 +256 @@
     quantize_and_encode_band_cost_ESC,
 };
 
+static float (*const quantize_and_encode_band_cost_rtz_arr[])(
+                                struct AACEncContext *s,
+                                PutBitContext *pb, const float *in,
+                                const float *scaled, int size, int scale_idx,
+                                int cb, const float lambda, const float uplim,
+                                int *bits) = {
+    quantize_and_encode_band_cost_ZERO,
+    quantize_and_encode_band_cost_SQUAD,
+    quantize_and_encode_band_cost_SQUAD,
+    quantize_and_encode_band_cost_UQUAD,
+    quantize_and_encode_band_cost_UQUAD,
+    quantize_and_encode_band_cost_SPAIR,
+    quantize_and_encode_band_cost_SPAIR,
+    quantize_and_encode_band_cost_UPAIR,
+    quantize_and_encode_band_cost_UPAIR,
+    quantize_and_encode_band_cost_UPAIR,
+    quantize_and_encode_band_cost_UPAIR,
+    quantize_and_encode_band_cost_ESC_RTZ,
+};
+
 #define quantize_and_encode_band_cost(                                  \
                                 s, pb, in, scaled, size, scale_idx, cb, \
-                                lambda, uplim, bits)                    \
-    quantize_and_encode_band_cost_arr[cb](                              \
+                                lambda, uplim, bits, rtz)               \
+    ((rtz) ? quantize_and_encode_band_cost_rtz_arr : quantize_and_encode_band_cost_arr)[cb]( \
                                 s, pb, in, scaled, size, scale_idx, cb, \
                                 lambda, uplim, bits)
 
 static float quantize_band_cost(struct AACEncContext *s, const float *in,
                                 const float *scaled, int size, int scale_idx,
                                 int cb, const float lambda, const float uplim,
-                                int *bits)
+                                int *bits, int rtz)
 {
     return quantize_and_encode_band_cost(s, NULL, in, scaled, size, scale_idx,
-                                         cb, lambda, uplim, bits);
+                                         cb, lambda, uplim, bits, rtz);
 }
 
 static void quantize_and_encode_band(struct AACEncContext *s, PutBitContext *pb,
                                      const float *in, int size, int scale_idx,
-                                     int cb, const float lambda)
+                                     int cb, const float lambda, int rtz)
 {
     quantize_and_encode_band_cost(s, pb, in, NULL, size, scale_idx, cb, lambda,
-                                  INFINITY, NULL);
+                                  INFINITY, NULL, rtz);
 }
 
 static float find_max_val(int group_len, int swb_size, const float *scaled) {
@@ -282 +311 @@
     return maxval;
 }
 
+static av_const float sqrf(float f) {
+    return f*f;
+}
+
+static float find_form_factor(int group_len, int swb_size, float thresh, const float *scaled, float cleanup_factor) {
+    const float iswb_size = 1.0f / swb_size;
+    const float iswb_sizem1 = 1.0f / (swb_size - 1);
+    const float ethresh = thresh;
+    float form = 0.0f, weight = 0.0f;
+    int w2, i;
+    for (w2 = 0; w2 < group_len; w2++) {
+        float e = 0.0f, e2 = 0.0f, var = 0.0f, maxval = 0.0f;
+        float nzl = 0;
+        for (i = 0; i < swb_size; i++) {
+            float s = fabsf(scaled[w2*128+i]);
+            maxval = FFMAX(maxval, s);
+            e += s;
+            e2 += s *= s;
+            /* We really don't want a hard non-zero-line count, since
+             * even below-threshold lines do add up towards band spectral power.
+             * So, fall steeply towards zero, but smoothly
+             */
+            if (s >= ethresh)
+                nzl += 1.0f;
+            else if (cleanup_factor > 1.0f)
+                nzl += powf(sqrf(s / ethresh), sqrf(FFMIN(4.0f,cleanup_factor)));
+            else
+                nzl += sqrf(s / ethresh);
+        }
+        if (e2 > thresh) {
+            float frm;
+            e *= iswb_size;
+            
+            // compute variance
+            for (i = 0; i < swb_size; i++) {
+                float d = fabsf(scaled[w2*128+i]) - e;
+                var += d*d;
+            }
+            var = sqrtf(var * iswb_sizem1);
+            
+            e2 *= iswb_size;
+            frm = e / FFMIN(e+4*var,maxval);
+            form += e2 * sqrtf(frm) / FFMAX(0.5f,nzl);
+            weight += e2;
+        }
+    }
+    if (weight > 0)
+        return form / weight;
+    else
+        return 1.0f;
+}
+
+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 unsigned char aac_maxval_cb[] = {
+    0, 1, 3, 5, 5, 7, 7, 7, 9, 9, 9, 9, 9, 11
+};
+
 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));
     int qmaxval, cb;
     qmaxval = maxval * Q34 + 0.4054f;
-    if      (qmaxval ==  0) cb = 0;
-    else if (qmaxval ==  1) cb = 1;
-    else if (qmaxval ==  2) cb = 3;
-    else if (qmaxval <=  4) cb = 5;
-    else if (qmaxval <=  7) cb = 7;
-    else if (qmaxval <= 12) cb = 9;
-    else                    cb = 11;
+    if (qmaxval >= (sizeof(aac_maxval_cb) / sizeof(aac_maxval_cb[0])))
+        cb = 11;
+    else
+        cb = aac_maxval_cb[qmaxval];
     return cb;
 }
 
@@ -351 +444 @@
                     rd += quantize_band_cost(s, sce->coeffs + start + w*128,
                                              s->scoefs + start + w*128, size,
                                              sce->sf_idx[(win+w)*16+swb], cb,
-                                             lambda / band->threshold, INFINITY, NULL);
+                                             lambda / band->threshold, INFINITY, NULL, 0);
                 }
                 cost_stay_here = path[swb][cb].cost + rd;
                 cost_get_here  = minrd              + rd + run_bits + 4;
@@ -473 +566 @@
                     bits += quantize_band_cost(s, sce->coeffs + start + w*128,
                                                s->scoefs + start + w*128, size,
                                                sce->sf_idx[(win+w)*16+swb], cb,
-                                               0, INFINITY, NULL);
+                                               0, INFINITY, NULL, 0);
                 }
                 cost_stay_here = path[swb][cb].cost + bits;
                 cost_get_here  = minbits            + bits + run_bits + 4;
@@ -655 +748 @@
                     for (w2 = 0; w2 < sce->ics.group_len[w]; w2++) {
                         FFPsyBand *band = &s->psy.ch[s->cur_channel].psy_bands[(w+w2)*16+g];
                         dist += quantize_band_cost(s, coefs + w2*128, s->scoefs + start + w2*128, sce->ics.swb_sizes[g],
-                                                   q + q0, cb, lambda / band->threshold, INFINITY, NULL);
+                                                   q + q0, cb, lambda / band->threshold, INFINITY, NULL, 0);
                     }
                     minrd = FFMIN(minrd, dist);
 
@@ -701 +794 @@
                 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;
+    char nzs[128];
+    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.0f * 120.f / lambda, 0.0625f, 16.0f);
+    float rdmin = 0.03125f;
+    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, maxscaler, nminscaler, minrdsf;
     int its  = 0;
+    int maxits = 20;
     int allz = 0;
-    float minthr = INFINITY;
+    int tbits;
+    int cutoff = 1024;
+    int cleanoff = 1024;
+    
+    /* 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) {
+        // Constant Q-scale doesn't compensat MS coding on its own
+        // No need to be overly precise, this only controls RD 
+        // adjustment CB limits when going overboard
+        if (s->options.stereo_mode && s->cur_type == TYPE_CPE)
+            destbits *= 2;
+        
+        // When using a constant Q-scale, don't adjust bits, just use RD
+        // Don't let it go overboard, though... 8x psy target is enough
+        toomanybits = 5800; //av_clip(FFMAX(refbits*3, destbits * 8), 768, 5800);
+        toofewbits = 0;
+        
+        // Don't offset scalers, just RD
+        sfoffs = sce->ics.num_windows - 1;
+        rdlambda = sqrtf(rdlambda);
+        
+        // search further
+        maxits = 40;
+    } else {
+        // When using ABR, be strict
+        toomanybits = destbits + destbits/16;
+        toofewbits = destbits - destbits/16;
+        
+        sfoffs = 0;
+        rdlambda = sqrtf(rdlambda);
+    }
+
+    // and zero out above cutoff frequency
+    {
+        int wlen = 1024 / sce->ics.num_windows;
+        int bandwidth;
+        int cleanwidth;
+
+        /* Scale, psy gives us constant quality, this LP only scales
+         * bitrate by lambda, so we save bits on subjectively unimportant HF
+         * rather than increase quantization noise. Adjust nominal bitrate
+         * to effective bitrate according to encoding parameters, _AAC_CUTOFF
+         * is calibrated for effective bitrate.
+         */
+        float rate_bandwidth_multiplier = 1.3f;
+        int frame_bit_rate = (avctx->flags & CODEC_FLAG_QSCALE) 
+            ? (refbits * rate_bandwidth_multiplier * avctx->sample_rate / 1024)
+            : (avctx->bit_rate / avctx->channels);
+        
+        if (avctx->cutoff > 0) 
+            bandwidth = avctx->cutoff;
+        else
+            bandwidth = FFMAX(3000, _AAC_CUTOFF(frame_bit_rate, 1, avctx->sample_rate));
+        cleanwidth = FFMAX(3000, _AAC_CUTOFF(frame_bit_rate * 0.8f, 1, avctx->sample_rate));
+
+        cutoff = bandwidth * 2 * wlen / avctx->sample_rate;
+        cleanoff = cleanwidth * 2 * wlen / avctx->sample_rate;
+    }
 
     // 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 energy = 0.0f;
+            
             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) {
+                if (start >= cutoff || band->energy <= (band->threshold * zeroscale) || band->threshold == 0.0) {
                     sce->zeroes[(w+w2)*16+g] = 1;
                     continue;
                 }
                 nz = 1;
             }
-            uplims[w*16+g] = uplim *512;
+            if (!nz) {
+                uplim = 0.0f;
+            } else {
+                nz = 0;
+                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 * zeroscale) || band->threshold == 0.0f)
+                        continue;
+                    uplim += band->threshold;
+                    energy += band->energy;
+                    nz++;
+                }
+            }
+            uplims[w*16+g] = uplim;
+            energies[w*16+g] = energy;
+            nzs[w*16+g] = nz;
             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 +1001 @@
         }
     }
 
+    /* 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]) {
+        /* psy already priorizes transients to some extent */
+        float de_psy_factor = (sce->ics.num_windows > 1) ? 8.0f / sce->ics.group_len[w] : 1.0f; 
+        start = w*128;
+        for (g = 0;  g < sce->ics.num_swb; g++) {
+            if (nzs[g] > 0) {
+                float energy2uplim = find_form_factor(
+                    sce->ics.group_len[w], sce->ics.swb_sizes[g], 
+                    uplims[w*16+g] / (nzs[g] * sce->ics.swb_sizes[w]), 
+                    sce->coeffs + start,
+                    start / (float)cleanoff);
+                energy2uplim *= de_psy_factor;
+                energy2uplim = FFMAX(0.015625f, FFMIN(1.0f,energy2uplim));
+                //fprintf(stderr, "%.3f ", energy2uplim);
+                uplims[w*16+g] *= av_clipf(rdlambda * energy2uplim, rdmin, rdmax)
+                                  * sce->ics.group_len[w];
+            }
+            start += sce->ics.swb_sizes[g];
+        }
+        //fprintf(stderr, "|");
+    }
+    //fprintf(stderr, "\n");
+
     //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 +1053 @@
                         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;
@@ -801 +1063 @@
                                                    cb,
                                                    1.0f,
                                                    INFINITY,
-                                                   &b);
+                                                   &b,
+                                                   0);
                         bits += b;
                     }
                     dists[w*16+g] = dist - bits;
@@ -813 +1076 @@
                     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;
+        maxscaler = 0;
+        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]);
+                    maxscaler = FFMAX(maxscaler, sce->sf_idx[w*16+g]);
+                }
+            }
+        }
+
+        //NOTE: Emulate v4 patch with n = 0, uncommenting if (mb >= ESC_BT) break,
+        //  and replacing minscaler-1 with minscaler
         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 = FFMAX3(60, minscaler-1, maxscaler - SCALE_MAX_DIFF);
         for (w = 0; w < sce->ics.num_windows; w += sce->ics.group_len[w]) {
+            int depth = (avctx->flags & CODEC_FLAG_QSCALE) ? 4 : 2;
+            int mblim;
+            float uplmax = FFMIN(4.0f, FFMAX(1.0f, tbits / (float)FFMAX(1,destbits)));
+            if (tbits > (destbits*6))
+                mblim = 9;
+            else if (tbits < (destbits/3))
+                mblim = (ESC_BT + 1);
+            else 
+                mblim = ESC_BT;
             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) {
-                    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;
+                if (!sce->zeroes[w*16+g]) {
+                    if (dists[w*16+g] > uplims[w*16+g] && sce->sf_idx[w*16+g] > minrdsf) {
+                        //Try to make sure there is some energy in every nonzero band
+                        //NOTE: This algorithm must be forcibly imbalanced, pushing harder
+                        //  on holes or more distorted bands at first, otherwise there's 
+                        //  no net gain (since the next iteration will offset all bands
+                        //  on the opposite direction to compensate for extra bits)
+                        int n = (its > 4 && g > sce->ics.num_swb/4) ? 0 : (depth/2);
+                        for (i = 0; i < depth; ++i) {
+                            int mb = find_min_book(maxvals[w*16+g], sce->sf_idx[w*16+g]-1);
+                            int cmb = find_min_book(maxvals[w*16+g], sce->sf_idx[w*16+g]);
+                            if ((mb > cmb || cmb > mblim) && mb >= mblim && dists[w*16+g] <= uplims[w*16+g]*uplmax && sce->sf_idx[w*16+g] <= (minrdsf+1)) break;
+                            sce->sf_idx[w*16+g]--;
+                            dists[w*16+g] *= 0.5f;
+                            if (mb && (i >= n || dists[w*16+g] < uplims[w*16+g])) break;
+                        }
+                    } else if (tbits > destbits && dists[w*16+g] < (uplims[w*16+g]*0.5f) && sce->sf_idx[w*16+g] < maxscaler) {
+                        // Um... over target
+                        sce->sf_idx[w*16+g]++;
+                        dists[w*16+g] *= 2.0f;
+                    }
                 }
-                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 */
+    //fprintf(stderr, "its:%d %d/%d: ", its, tbits, destbits);
+    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 : 3;
+            //fprintf(stderr, "%02x ", sce->sf_idx[w*16+g]);
+            if (sce->band_type[w*16+g] <= minbt) {
+                sce->zeroes[w*16+g] = 1;
+                sce->band_type[w*16+g] = 0;
+            }
+        }
+        //fprintf(stderr, "|");
+    }
+    //fprintf(stderr, "\n");
+    //fprintf(stderr, "ba:%d br:%d  \t\r", s->psy.bitres.alloc, tbits);
 }
 
 static void search_for_quantizers_faac(AVCodecContext *avctx, AACEncContext *s,
@@ -967 +1295 @@
                                                ESC_BT,
                                                lambda,
                                                INFINITY,
-                                               &b);
+                                               &b,
+                                               0);
                     dist -= b;
                 }
                 dist *= 1.0f / 512.0f / lambda;
-                quant_max = quant(maxq[w*16+g], ff_aac_pow2sf_tab[POW_SF2_ZERO - scf + SCALE_ONE_POS - SCALE_DIV_512]);
+                quant_max = quant(maxq[w*16+g], ff_aac_pow2sf_tab[POW_SF2_ZERO - scf + SCALE_ONE_POS - SCALE_DIV_512], ROUND_STANDARD);
                 if (quant_max >= 8191) { // too much, return to the previous quantizer
                     sce->sf_idx[w*16+g] = prev_scf;
                     break;
@@ -1021 +1350 @@
                                        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,
                           const float lambda)
 {
-    int start = 0, i, w, w2, g;
+    int start = 0, i, w, w2, g, sid_sf_boost;
     float M[128], S[128];
     float *L34 = s->scoefs, *R34 = s->scoefs + 128, *M34 = s->scoefs + 128*2, *S34 = s->scoefs + 128*3;
+    float mslambda = FFMIN(1.0f, lambda / 120.0f);
     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;
+            cpe->ms_mask[w*16+g] = 0;
             if (!cpe->ch[0].zeroes[w*16+g] && !cpe->ch[1].zeroes[w*16+g]) {
-                float dist1 = 0.0f, dist2 = 0.0f;
+                float Mmax = 0.0f, Smax = 0.0f;
+                
+                /* Must compute mid/side SF and book for the whole window group */
                 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);
                     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;
+                        M[i] = (sce0->pcoeffs[start+w2*128+i]
+                              + sce1->pcoeffs[start+w2*128+i]) * 0.5f;
                         S[i] =  M[i]
-                              - sce1->coeffs[start+w2*128+i];
+                              - sce1->pcoeffs[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]);
+                    }
+                }
+                
+                for (sid_sf_boost = 0; sid_sf_boost < 4; sid_sf_boost++) {
+                    float dist1 = 0.0f, dist2 = 0.0f;
+                    int B0 = 0, B1 = 0;
+                    int minidx;
+                    int mididx, sididx;
+                    int midcb, sidcb;
+                    
+                    minidx = FFMIN(sce0->sf_idx[w*16+g], sce1->sf_idx[w*16+g]);
+                    mididx = av_clip(minidx, min_sf_idx_mid, min_sf_idx_mid + SCALE_MAX_DIFF);
+                    sididx = av_clip(minidx - sid_sf_boost * 3, min_sf_idx_side, min_sf_idx_side + SCALE_MAX_DIFF);
+                    midcb = find_min_book(Mmax, mididx);
+                    sidcb = find_min_book(Smax, sididx);
+                    
+                    if ((mididx > minidx) || (sididx > minidx)) {
+                        /* scalefactor range violation, bad stuff, will decrease quality unacceptably */
+                        continue;
+                    }
+                    
+                    /* No CB can be zero */
+                    midcb = FFMAX(1,midcb);
+                    sidcb = FFMAX(1,sidcb);
+                    
+                    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);
+                        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;
+                            S[i] =  M[i]
+                                  - sce1->coeffs[start+w2*128+i];
+                        }
+                        abs_pow34_v(L34, sce0->coeffs+start+w2*128, sce0->ics.swb_sizes[g]);
+                        abs_pow34_v(R34, sce1->coeffs+start+w2*128, sce0->ics.swb_sizes[g]);
+                        abs_pow34_v(M34, M,                         sce0->ics.swb_sizes[g]);
+                        abs_pow34_v(S34, S,                         sce0->ics.swb_sizes[g]);
+                        dist1 += quantize_band_cost(s, sce0->coeffs + start + w2*128,
+                                                    L34,
+                                                    sce0->ics.swb_sizes[g],
+                                                    sce0->sf_idx[(w+w2)*16+g],
+                                                    sce0->band_type[w*16+g],
+                                                    lambda / band0->threshold, INFINITY, &b1, 0);
+                        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*16+g],
+                                                    lambda / band1->threshold, INFINITY, &b2, 0);
+                        dist2 += quantize_band_cost(s, M,
+                                                    M34,
+                                                    sce0->ics.swb_sizes[g],
+                                                    mididx,
+                                                    midcb,
+                                                    lambda / minthr, INFINITY, &b3, 0);
+                        dist2 += quantize_band_cost(s, S,
+                                                    S34,
+                                                    sce1->ics.swb_sizes[g],
+                                                    sididx,
+                                                    sidcb,
+                                                    lambda * mslambda / (minthr * bmax), INFINITY, &b4, 0);
+                        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]) {
+                        sce0->sf_idx[w*16+g] = mididx;
+                        sce1->sf_idx[w*16+g] = sididx;
+                        sce0->band_type[w*16+g] = midcb;
+                        sce1->band_type[w*16+g] = sidcb;
+                        break;
+                    } else if (B1 > B0) {
+                        /* More boost won't fix this */
+                        break;
                     }
-                    abs_pow34_v(L34, sce0->coeffs+start+w2*128, sce0->ics.swb_sizes[g]);
-                    abs_pow34_v(R34, sce1->coeffs+start+w2*128, sce0->ics.swb_sizes[g]);
-                    abs_pow34_v(M34, M,                         sce0->ics.swb_sizes[g]);
-                    abs_pow34_v(S34, S,                         sce0->ics.swb_sizes[g]);
-                    dist1 += quantize_band_cost(s, sce0->coeffs + start + w2*128,
-                                                L34,
-                                                sce0->ics.swb_sizes[g],
-                                                sce0->sf_idx[(w+w2)*16+g],
-                                                sce0->band_type[(w+w2)*16+g],
-                                                lambda / band0->threshold, INFINITY, NULL);
-                    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);
-                    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);
-                    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);
                 }
-                cpe->ms_mask[w*16+g] = dist2 < dist1;
             }
             start += sce0->ics.swb_sizes[g];
         }

libavcodec/aacenc.c

@@ -46 +46 @@
 #include "psymodel.h"
 
 #define AAC_MAX_CHANNELS 6
+#define CLIP_AVOIDANCE_FACTOR 0.95f
 
 #define ERROR_IF(cond, ...) \
     if (cond) { \
@@ -53 +54 @@
         return AVERROR(EINVAL); \
     }
 
+#define WARN_IF(cond, ...) \
+    if (cond) { \
+        av_log(avctx, AV_LOG_WARNING, __VA_ARGS__); \
+    }
+
 float ff_aac_pow34sf_tab[428];
 
 static const uint8_t swb_size_1024_96[] = {
@@ -102 +108 @@
     swb_size_1024_96, swb_size_1024_96, swb_size_1024_64,
     swb_size_1024_48, swb_size_1024_48, swb_size_1024_32,
     swb_size_1024_24, swb_size_1024_24, swb_size_1024_16,
-    swb_size_1024_16, swb_size_1024_16, swb_size_1024_8
+    swb_size_1024_16, swb_size_1024_16, swb_size_1024_8, 
+    swb_size_1024_8
 };
 
 static const uint8_t swb_size_128_96[] = {
@@ -131 +138 @@
     swb_size_128_96, swb_size_128_96, swb_size_128_96,
     swb_size_128_48, swb_size_128_48, swb_size_128_48,
     swb_size_128_24, swb_size_128_24, swb_size_128_16,
-    swb_size_128_16, swb_size_128_16, swb_size_128_8
+    swb_size_128_16, swb_size_128_16, swb_size_128_8, 
+    swb_size_128_8
 };
 
 /** default channel configurations */
@@ -260 +268 @@
         for (i = 0; i < 1024; i += 128)
             s->mdct128.mdct_calc(&s->mdct128, sce->coeffs + i, output + i*2);
     memcpy(audio, audio + 1024, sizeof(audio[0]) * 1024);
+    memcpy(sce->pcoeffs, sce->coeffs, sizeof(sce->pcoeffs));
 }
 
 /**
@@ -311 +320 @@
         start = 0;
         maxsfb = 0;
         cpe->ch[ch].pulse.num_pulse = 0;
-        for (w = 0; w < ics->num_windows*16; w += 16) {
-            for (g = 0; g < ics->num_swb; g++) {
-                //apply M/S
-                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];
+        for (w = 0; w < ics->num_windows; w += ics->group_len[w]) {
+            for (w2 = 0; w2 < ics->group_len[w]; w2++) {
+                start = (w+w2) * 128;
+                for (g = 0; g < ics->num_swb; g++) {
+                    //apply M/S
+                    if (cpe->common_window && !ch && cpe->ms_mask[w*16 + g]) {
+                        for (i = 0; i < ics->swb_sizes[g]; i++) {
+                            cpe->ch[0].coeffs[start+i] = (cpe->ch[0].pcoeffs[start+i] + cpe->ch[1].pcoeffs[start+i]) * 0.5f;
+                            cpe->ch[1].coeffs[start+i] = cpe->ch[0].coeffs[start+i] - cpe->ch[1].pcoeffs[start+i];
+                        }
                     }
+                    start += ics->swb_sizes[g];
                 }
-                start += ics->swb_sizes[g];
+                for (cmaxsfb = ics->num_swb; cmaxsfb > 0 && cpe->ch[ch].zeroes[w*16+cmaxsfb-1]; cmaxsfb--)
+                    ;
+                maxsfb = FFMAX(maxsfb, cmaxsfb);
             }
-            for (cmaxsfb = ics->num_swb; cmaxsfb > 0 && cpe->ch[ch].zeroes[w+cmaxsfb-1]; cmaxsfb--)
-                ;
-            maxsfb = FFMAX(maxsfb, cmaxsfb);
         }
         ics->max_sfb = maxsfb;
 
@@ -430 +442 @@
                                                    sce->ics.swb_sizes[i],
                                                    sce->sf_idx[w*16 + i],
                                                    sce->band_type[w*16 + i],
-                                                   s->lambda);
+                                                   s->lambda, sce->ics.window_clipping[w]);
+            start += sce->ics.swb_sizes[i];
+        }
+    }
+}
+
+/**
+ * Downscale spectral coefficients for near-clipping windows to avoid artifacts
+ */
+static void avoid_clipping(AACEncContext *s, SingleChannelElement *sce)
+{
+    int start, i, j, w, w2;
+
+    for (w = 0; w < sce->ics.num_windows; w += sce->ics.group_len[w]) {
+        start = 0;
+        for (i = 0; i < sce->ics.max_sfb; i++) {
+            if (sce->ics.window_clipping[w]) {
+                for (w2 = w; w2 < w + sce->ics.group_len[w]; w2++) {
+                    float *swb_coeffs = sce->coeffs + start + w2*128;
+                    for (j = 0; j < sce->ics.swb_sizes[i]; j++) 
+                        swb_coeffs[j] *= CLIP_AVOIDANCE_FACTOR;
+                }
+            }
             start += sce->ics.swb_sizes[i];
         }
     }
@@ -507 +541 @@
     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, ms_mode = 0;
     int chan_el_counter[4];
     FFPsyWindowInfo windows[AAC_MAX_CHANNELS];
 
@@ -565 +599 @@
             ics->num_swb            = tag == TYPE_LFE ? ics->num_swb : s->psy.num_bands[ics->num_windows == 8];
             for (w = 0; w < ics->num_windows; w++)
                 ics->group_len[w] = wi[ch].grouping[w];
+            for (w = 0; w < ics->num_windows; w++)
+                ics->window_clipping[w] = wi[ch].clipping[w];
 
             apply_window_and_mdct(s, &cpe->ch[ch], overlap);
+            avoid_clipping(s, &cpe->ch[ch]);
         }
         start_ch += chans;
     }
     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 +632 @@
             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.bitres.bits = avctx->frame_bits / s->channels;
             s->psy.model->analyze(&s->psy, start_ch, coeffs, wi);
+            if (s->psy.bitres.alloc > 0) {
+                /* Lambda unused here on purpose, we need to take psy's unscaled allocation */
+                target_bits += s->psy.bitres.alloc;
+                s->psy.bitres.alloc /= chans;
+            }
+            s->cur_type = tag;
             for (ch = 0; ch < chans; ch++) {
                 s->cur_channel = start_ch + ch;
                 s->coder->search_for_quantizers(avctx, s, &cpe->ch[ch], s->lambda);
@@ -626 +675 @@
                 if (cpe->common_window) {
                     put_ics_info(s, &cpe->ch[0].ics);
                     encode_ms_info(&s->pb, cpe);
+                    if (cpe->ms_mode) ms_mode = 1;
                 }
             }
             for (ch = 0; ch < chans; ch++) {
@@ -635 +685 @@
             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;
+        //fprintf(stderr, "l:%f\t%d\t%d\t%d\t%d\n", s->lambda, too_few_bits, frame_bits, target_bits, too_many_bits);
+        
+        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 ratio = ((float)target_bits) / frame_bits;
+            
+            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), but slowly
+                */
+                ratio = sqrtf(sqrtf(ratio));
+                ratio = av_clipf(ratio, 0.9f, 1.1f);
+            } else {
+                /* Not so fast though */
+                ratio = sqrtf(ratio);
+            }
+            s->lambda = FFMIN(s->lambda * ratio, 65536.f);
+            
+            // Keep iterating if we must reduce and lambda is in the sky
+            // NOT if we did M/S coding, it's undoable
+            if (s->lambda < 300.f || ratio > 0.9f)
+                break;
+            else {
+                if (ms_mode) {
+                    for (i = 0; i < s->chan_map[0]; i++) {
+                        // Must restore coeffs
+                        chans = tag == TYPE_CPE ? 2 : 1;
+                        cpe = &s->cpe[i];
+                        for (ch = 0; ch < chans; ch++)
+                            memcpy(cpe->ch[ch].coeffs, cpe->ch[ch].pcoeffs, sizeof(cpe->ch[ch].coeffs));
+                    }
+                }
+                its++;
+            }
+        } 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 (!frame)
         s->last_frame++;
 
@@ -733 +825 @@
 
     s->channels = avctx->channels;
 
-    ERROR_IF(i == 16,
+    ERROR_IF(i == 16 
+                || i >= (sizeof(swb_size_1024) / sizeof(*swb_size_1024))
+                || i >= (sizeof(swb_size_128) / sizeof(*swb_size_128)),
              "Unsupported sample rate %d\n", avctx->sample_rate);
     ERROR_IF(s->channels > AAC_MAX_CHANNELS,
              "Unsupported number of channels: %d\n", s->channels);
     ERROR_IF(avctx->profile != FF_PROFILE_UNKNOWN && avctx->profile != FF_PROFILE_AAC_LOW,
              "Unsupported profile %d\n", avctx->profile);
-    ERROR_IF(1024.0 * avctx->bit_rate / avctx->sample_rate > 6144 * s->channels,
-             "Too many bits per frame requested\n");
+    WARN_IF(1024.0 * avctx->bit_rate / avctx->sample_rate > 6144 * s->channels,
+             "Too many bits per frame requested, clamping to max\n");
+    
+    avctx->bit_rate = (int)FFMIN(
+        6144 * s->channels / 1024.0 * avctx->sample_rate, 
+        avctx->bit_rate);
 
     s->samplerate_index = i;
 
@@ -787 +885 @@
 
 #define AACENC_FLAGS AV_OPT_FLAG_ENCODING_PARAM | AV_OPT_FLAG_AUDIO_PARAM
 static const AVOption aacenc_options[] = {
-    {"stereo_mode", "Stereo coding method", offsetof(AACEncContext, options.stereo_mode), AV_OPT_TYPE_INT, {.i64 = 0}, -1, 1, AACENC_FLAGS, "stereo_mode"},
+    {"stereo_mode", "Stereo coding method", offsetof(AACEncContext, options.stereo_mode), AV_OPT_TYPE_INT, {.i64 = -1}, -1, 1, AACENC_FLAGS, "stereo_mode"},
         {"auto",     "Selected by the Encoder", 0, AV_OPT_TYPE_CONST, {.i64 = -1 }, INT_MIN, INT_MAX, AACENC_FLAGS, "stereo_mode"},
         {"ms_off",   "Disable Mid/Side coding", 0, AV_OPT_TYPE_CONST, {.i64 =  0 }, INT_MIN, INT_MAX, AACENC_FLAGS, "stereo_mode"},
         {"ms_force", "Force Mid/Side for the whole frame if possible", 0, AV_OPT_TYPE_CONST, {.i64 =  1 }, INT_MIN, INT_MAX, AACENC_FLAGS, "stereo_mode"},

libavcodec/aacenc.h

@@ -52 +52 @@
     void (*encode_window_bands_info)(struct AACEncContext *s, SingleChannelElement *sce,
                                      int win, int group_len, const float lambda);
     void (*quantize_and_encode_band)(struct AACEncContext *s, PutBitContext *pb, const float *in, int size,
-                                     int scale_idx, int cb, const float lambda);
+                                     int scale_idx, int cb, const float lambda, int rtz);
     void (*search_for_ms)(struct AACEncContext *s, ChannelElement *cpe, const float lambda);
 } AACCoefficientsEncoder;
 
@@ -78 +78 @@
     FFPsyContext psy;
     struct FFPsyPreprocessContext* psypp;
     AACCoefficientsEncoder *coder;
-    int cur_channel;
+    int cur_channel;                             ///< current channel for coder context
     int last_frame;
     float lambda;
+    enum RawDataBlockType cur_type;              ///< channel group type cur_channel belongs to
+    
     AudioFrameQueue afq;
     DECLARE_ALIGNED(16, int,   qcoefs)[96];      ///< quantized coefficients
     DECLARE_ALIGNED(32, float, scoefs)[1024];    ///< scaled coefficients

libavcodec/aacpsy.c

@@ -87 +87 @@
 };
 
 #define PSY_3GPP_BITS_TO_PE(bits) ((bits) * 1.18f)
+#define PSY_3GPP_PE_TO_BITS(bits) ((bits) / 1.18f)
 
 /* LAME psy model constants */
 #define PSY_LAME_FIR_LEN 21         ///< LAME psy model FIR order
@@ -95 +96 @@
 #define AAC_NUM_BLOCKS_SHORT 8      ///< number of blocks in a short sequence
 #define PSY_LAME_NUM_SUBBLOCKS 3    ///< Number of sub-blocks in each short block
 
+#define CLIPPING_THRESHOLD 0.98f
+
 /**
  * @}
  */
@@ -157 +160 @@
     } pe;
     AacPsyCoeffs psy_coef[2][64];
     AacPsyChannel *ch;
+    float global_quality; ///< normalized global quality taken from avctx
 }AacPsyContext;
 
 /**
@@ -255 +259 @@
 /**
  * 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 +302 @@
     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;
+    pctx->global_quality = (ctx->avctx->global_quality ? ctx->avctx->global_quality : 120) * 0.01f;
 
+    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 +401 @@
                                                  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];
@@ -501 +511 @@
     ctx->pe.max = FFMAX(pe, ctx->pe.max);
     ctx->pe.min = FFMIN(pe, ctx->pe.min);
 
-    return FFMIN(ctx->frame_bits * bit_factor, ctx->frame_bits + size - bits);
+    /* NOTE: allocate a minimum of 1/8th average frame bits, to avoid
+     *   reservoir starvation from producing zero-bit frames
+     */
+    return FFMIN(
+        ctx->frame_bits * bit_factor, 
+        FFMAX(ctx->frame_bits + size - bits, ctx->frame_bits / 8));
 }
 
 static float calc_pe_3gpp(AacPsyBand *band)
@@ -566 +581 @@
     return thr;
 }
 
+static av_always_inline float sqrf(float f)
+{
+    return f*f;
+}
+
 #ifndef calc_thr_3gpp
 static void calc_thr_3gpp(const FFPsyWindowInfo *wi, const int num_bands, AacPsyChannel *pch,
-                          const uint8_t *band_sizes, const float *coefs)
+                          const uint8_t *band_sizes, const float *coefs, float q)
 {
     int i, w, g;
-    int start = 0;
+    int start = 0, wstart = 0, end = 1024 / wi->num_windows;
+    float iq2 = sqrf(1.0f / FFMAX(0.125f, q));
     for (w = 0; w < wi->num_windows*16; w += 16) {
+        wstart = start;
         for (g = 0; g < num_bands; g++) {
             AacPsyBand *band = &pch->band[w+g];
 
             float form_factor = 0.0f;
-            float Temp;
+            float Temp, absthr;
             band->energy = 0.0f;
             for (i = 0; i < band_sizes[g]; i++) {
                 band->energy += coefs[start+i] * coefs[start+i];
                 form_factor  += sqrtf(fabs(coefs[start+i]));
             }
+            // Absolute threshold of hearing approx taking the low point 
+            // to be 16-bit quantization noise, scaled up to improve vbr
+            absthr = sqrf(sqrf(FFMAX(1.0f, (abs(start - wstart - end/16) - end/64) / end * 6.5f)) * end) * 3.0517578125e-5;
+            absthr *= band_sizes[g];
             Temp = band->energy > 0 ? sqrtf((float)band_sizes[g] / band->energy) : 0;
-            band->thr      = band->energy * 0.001258925f;
+            band->thr      = FFMAX(absthr * iq2, band->energy * 0.001258925f);
             band->nz_lines = form_factor * sqrtf(Temp);
 
             start += band_sizes[g];
@@ -628 +654 @@
     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);
+    calc_thr_3gpp(wi, num_bands, pch, band_sizes, coefs, pctx->global_quality);
 
     //modify thresholds and energies - spread, threshold in quiet, pre-echo control
     for (w = 0; w < wi->num_windows*16; w += 16) {
@@ -671 +697 @@
 
     /* 5.6.1.3.2 "Calculation of the desired perceptual entropy" */
     ctx->ch[channel].entropy = pe;
-    desired_bits = calc_bit_demand(pctx, pe, ctx->bitres.bits, ctx->bitres.size, wi->num_windows == 8);
-    desired_pe = PSY_3GPP_BITS_TO_PE(desired_bits);
-    /* NOTE: PE correction is kept simple. During initial testing it had very
-     *       little effect on the final bitrate. Probably a good idea to come
-     *       back and do more testing later.
-     */
-    if (ctx->bitres.bits > 0)
-        desired_pe *= av_clipf(pctx->pe.previous / PSY_3GPP_BITS_TO_PE(ctx->bitres.bits),
-                               0.85f, 1.15f);
+    if (ctx->avctx->flags & CODEC_FLAG_QSCALE) {
+        /* (2.5 * 120) achieves almost transparent rate, and we want to give
+         * ample room downwards, so we make that equivalent to QSCALE=2.4
+         */
+        desired_pe = pe * (ctx->avctx->global_quality ? ctx->avctx->global_quality : 120) / (2 * 2.5f * 120.0f);
+        desired_bits = FFMIN(2560, PSY_3GPP_PE_TO_BITS(desired_pe));
+        desired_pe = PSY_3GPP_BITS_TO_PE(desired_bits); // reflect clipping
+        
+        /* PE slope smoothing */
+        if (ctx->bitres.bits > 0) {
+            desired_bits = FFMIN(2560, PSY_3GPP_PE_TO_BITS(desired_pe));
+            desired_pe = PSY_3GPP_BITS_TO_PE(desired_bits); // reflect clipping
+        }
+        
+        pctx->pe.max = FFMAX(pe, pctx->pe.max);
+        pctx->pe.min = FFMIN(pe, pctx->pe.min);
+    } else {
+        desired_bits = calc_bit_demand(pctx, pe, ctx->bitres.bits, ctx->bitres.size, wi->num_windows == 8);
+        desired_pe = PSY_3GPP_BITS_TO_PE(desired_bits);
+        
+        /* NOTE: PE correction is kept simple. During initial testing it had very
+         *       little effect on the final bitrate. Probably a good idea to come
+         *       back and do more testing later.
+         */
+        if (ctx->bitres.bits > 0)
+            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" */
         for (w = 0; w < wi->num_windows*16; w += 16) {
@@ -717 +763 @@
             }
             desired_pe_no_ah = FFMAX(desired_pe - (pe - pe_no_ah), 0.0f);
             if (active_lines > 0.0f)
-                reduction += calc_reduction_3gpp(a, desired_pe_no_ah, pe_no_ah, active_lines);
+                reduction = calc_reduction_3gpp(a, desired_pe_no_ah, pe_no_ah, active_lines);
 
             pe = 0.0f;
             for (w = 0; w < wi->num_windows*16; w += 16) {
@@ -827 +873 @@
     int grouping     = 0;
     int uselongblock = 1;
     int attacks[AAC_NUM_BLOCKS_SHORT + 1] = { 0 };
+    uint8_t clippings[AAC_NUM_BLOCKS_SHORT];
     int i;
     FFPsyWindowInfo wi = { { 0 } };
 
@@ -841 +888 @@
 
         /* LAME comment: apply high pass filter of fs/4 */
         psy_hp_filter(firbuf, hpfsmpl, psy_fir_coeffs);
-
+        
         /* Calculate the energies of each sub-shortblock */
         for (i = 0; i < PSY_LAME_NUM_SUBBLOCKS; i++) {
             energy_subshort[i] = pch->prev_energy_subshort[i + ((AAC_NUM_BLOCKS_SHORT - 1) * PSY_LAME_NUM_SUBBLOCKS)];
@@ -916 +963 @@
 
     lame_apply_block_type(pch, &wi, uselongblock);
 
+    /* Calculate input sample maximums and evaluate clipping risk */
+    if (audio) {
+        for (i = 0; i < AAC_NUM_BLOCKS_SHORT; i++) {
+            const float *wbuf = audio + i * AAC_BLOCK_SIZE_SHORT;
+            float max = 0;
+            int j;
+            for (j = 0; j < AAC_BLOCK_SIZE_SHORT; j++)
+                max = FFMAX(max, fabsf(wbuf[j]));
+            clippings[i] = (max > CLIPPING_THRESHOLD) ? 1 : 0;
+        }
+    } else {
+        for (i = 0; i < 8; i++)
+            clippings[i] = 0;
+    }
+
     wi.window_type[1] = prev_type;
     if (wi.window_type[0] != EIGHT_SHORT_SEQUENCE) {
+        int clipping = 0;
+        
         wi.num_windows  = 1;
         wi.grouping[0]  = 1;
         if (wi.window_type[0] == LONG_START_SEQUENCE)
             wi.window_shape = 0;
         else
             wi.window_shape = 1;
+
+        for (i = 0; i < 8 && !clipping; i++)
+            clipping = clippings[i];
+        wi.clipping[0] = clipping;
     } else {
-        int lastgrp = 0;
+        int lastgrp = 0, anyclipped = 0;
 
         wi.num_windows = 8;
         wi.window_shape = 0;
@@ -934 +1002 @@
                 lastgrp = i;
             wi.grouping[lastgrp]++;
         }
+        
+        for (i = 0; i < 8; i += wi.grouping[i]) {
+            int w, clipping = 0;
+            for (w = 0; w < wi.grouping[i] && !clipping; w++)
+                clipping = clippings[i+w];
+            wi.clipping[i] = clipping;
+            if (clipping) anyclipped = 1;
+        }
     }
 
     /* Determine grouping, based on the location of the first attack, and save for
@@ -948 +1024 @@
             break;
         }
     }
+    
     pch->next_grouping = window_grouping[grouping];
 
     pch->prev_attack = attacks[8];

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);
@@ -101 +101 @@
     ctx        = av_mallocz(sizeof(FFPsyPreprocessContext));
     ctx->avctx = avctx;
 
-    if (avctx->cutoff > 0)
-        cutoff_coeff = 2.0 * avctx->cutoff / avctx->sample_rate;
-
-    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,
-                                             cutoff_coeff, 0.0, 0.0);
-    if (ctx->fcoeffs) {
-        ctx->fstate = av_mallocz_array(sizeof(ctx->fstate[0]), avctx->channels);
-        for (i = 0; i < avctx->channels; i++)
-            ctx->fstate[i] = ff_iir_filter_init_state(FILT_ORDER);
+    /* AAC has its own LP method */
+    if (avctx->codec_id != AV_CODEC_ID_AAC) {
+        if (avctx->cutoff > 0)
+            cutoff_coeff = 2.0 * avctx->cutoff / 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,
+                                                 cutoff_coeff, 0.0, 0.0);
+        if (ctx->fcoeffs) {
+            ctx->fstate = av_mallocz(sizeof(ctx->fstate[0]) * avctx->channels);
+            for (i = 0; i < avctx->channels; i++)
+                ctx->fstate[i] = ff_iir_filter_init_state(FILT_ORDER);
+        }
     }
 
     ff_iir_filter_init(&ctx->fiir);

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
@@ -67 +80 @@
     int window_shape;                 ///< window shape (sine/KBD/whatever)
     int num_windows;                  ///< number of windows in a frame
     int grouping[8];                  ///< window grouping (for e.g. AAC)
+    uint8_t clipping[8];              ///< set if a window group's samples are near clipping (for e.g. AAC)
     int *window_sizes;                ///< sequence of window sizes inside one frame (for eg. WMA)
 } FFPsyWindowInfo;
 
@@ -88 +102 @@
     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