[FFmpeg-devel] [RFC] AAC Encoder

Kostya kostya.shishkov
Fri Aug 15 18:59:52 CEST 2008


On Thu, Aug 14, 2008 at 03:38:17PM +0200, Michael Niedermayer wrote:
> viterbi for determining band_types ...
> look this isnt hard, its not even slow in this paricular case,let me explain
[explanation skipped]

Hmm, I have not understood it at the beginning but then I found out it's
strikingly similar to the one-pass almost-optimal LZ matching scheme.
So here it is with other comments taken care of too.

P.S. I was surprised to find out that I won't be near computer next week
so I'll try to make encoder fit for committing ASAP.

[...]
> -- 
> Michael     GnuPG fingerprint: 9FF2128B147EF6730BADF133611EC787040B0FAB
> 
> Republics decline into democracies and democracies degenerate into
> despotisms. -- Aristotle
-------------- next part --------------
--- /home/kst/cvs-get/ffmpeg/libavcodec/aacenc.c	2008-08-14 19:01:30.000000000 +0300
+++ aacenc.c	2008-08-15 19:50:22.000000000 +0300
@@ -27,8 +27,7 @@
 /***********************************
  *              TODOs:
  * psy model selection with some option
- * change greedy codebook search into something more optimal, like Viterbi algorithm
- * determine run lengths along with codebook
+ * add sane pulse detection
  ***********************************/
 
 #include "avcodec.h"
@@ -119,6 +118,34 @@
     swb_size_128_16, swb_size_128_16, swb_size_128_8
 };
 
+#define CB_UNSIGNED 0x01    ///< coefficients are coded as absolute values
+#define CB_PAIRS    0x02    ///< coefficients are grouped into pairs before coding (quads by default)
+#define CB_ESCAPE   0x04    ///< codebook allows escapes
+
+/** spectral coefficients codebook information */
+static const struct {
+    int16_t maxval;         ///< maximum possible value
+     int8_t cb_num;         ///< codebook number
+    uint8_t flags;          ///< codebook features
+} aac_cb_info[] = {
+    {    0, -1, CB_UNSIGNED }, // zero codebook
+    {    1,  0, 0 },
+    {    1,  1, 0 },
+    {    2,  2, CB_UNSIGNED },
+    {    2,  3, CB_UNSIGNED },
+    {    4,  4, CB_PAIRS },
+    {    4,  5, CB_PAIRS },
+    {    7,  6, CB_PAIRS | CB_UNSIGNED },
+    {    7,  7, CB_PAIRS | CB_UNSIGNED },
+    {   12,  8, CB_PAIRS | CB_UNSIGNED },
+    {   12,  9, CB_PAIRS | CB_UNSIGNED },
+    { 8191, 10, CB_PAIRS | CB_UNSIGNED | CB_ESCAPE },
+    {   -1, -1, 0 }, // reserved
+    {   -1, -1, 0 }, // perceptual noise substitution
+    {   -1, -1, 0 }, // intensity out-of-phase
+    {   -1, -1, 0 }, // intensity in-phase
+};
+
 /** default channel configurations */
 static const uint8_t aac_chan_configs[6][5] = {
  {1, ID_SCE},                         // 1 channel  - single channel element
@@ -130,6 +157,39 @@
 };
 
 /**
+ * structure used in optimal codebook search
+ */
+typedef struct BandCodingPath {
+    int prev_idx; ///< pointer to the previous path point
+    int codebook; ///< codebook for coding band run
+    int bits;     ///< number of bit needed to code given number of bands
+} BandCodingPath;
+
+/**
+ * AAC encoder context
+ */
+typedef struct {
+    PutBitContext pb;
+    MDCTContext mdct1024;                        ///< long (1024 samples) frame transform context
+    MDCTContext mdct128;                         ///< short (128 samples) frame transform context
+    DSPContext  dsp;
+    DECLARE_ALIGNED_16(FFTSample, output[2048]); ///< temporary buffer for MDCT input coefficients
+    int16_t* samples;                            ///< saved preprocessed input
+
+    int samplerate_index;                        ///< MPEG-4 samplerate index
+    const uint8_t *swb_sizes1024;                ///< scalefactor band sizes for long frame
+    int swb_num1024;                             ///< number of scalefactor bands for long frame
+    const uint8_t *swb_sizes128;                 ///< scalefactor band sizes for short frame
+    int swb_num128;                              ///< number of scalefactor bands for short frame
+
+    ChannelElement *cpe;                         ///< channel elements
+    AACPsyContext psy;                           ///< psychoacoustic model context
+    int last_frame;
+    BandCodingPath path[64];                     ///< auxiliary data needed for optimal band info coding
+    int band_bits[64][12];                       ///< bits needed to encode each band with each codebook
+} AACEncContext;
+
+/**
  * Make AAC audio config object.
  * @see 1.6.2.1 "Syntax - AudioSpecificConfig"
  */
@@ -176,6 +236,11 @@
     dsputil_init(&s->dsp, avctx);
     ff_mdct_init(&s->mdct1024, 11, 0);
     ff_mdct_init(&s->mdct128,   8, 0);
+    // window init
+    ff_kbd_window_init(ff_aac_kbd_long_1024, 4.0, 1024);
+    ff_kbd_window_init(ff_aac_kbd_short_128, 6.0, 128);
+    ff_sine_window_init(ff_sine_1024, 1024);
+    ff_sine_window_init(ff_sine_128, 128);
 
     s->samples = av_malloc(2 * 1024 * avctx->channels * sizeof(s->samples[0]));
     s->cpe = av_mallocz(sizeof(ChannelElement) * aac_chan_configs[avctx->channels-1][0]);
@@ -189,6 +254,55 @@
     return 0;
 }
 
+static void apply_window_and_mdct(AVCodecContext *avctx, AACEncContext *s, ChannelElement *cpe, short *audio, int channel)
+{
+    int i, j, k;
+    const float * lwindow = cpe->ch[channel].ics.use_kb_window[0] ? ff_aac_kbd_long_1024 : ff_sine_1024;
+    const float * swindow = cpe->ch[channel].ics.use_kb_window[0] ? ff_aac_kbd_short_128 : ff_sine_128;
+    const float * pwindow = cpe->ch[channel].ics.use_kb_window[1] ? ff_aac_kbd_short_128 : ff_sine_128;
+
+    if (cpe->ch[channel].ics.window_sequence[0] != EIGHT_SHORT_SEQUENCE) {
+        memcpy(s->output, cpe->ch[channel].saved, sizeof(float)*1024);
+        if(cpe->ch[channel].ics.window_sequence[0] == LONG_STOP_SEQUENCE){
+            memset(s->output, 0, sizeof(s->output[0]) * 448);
+            for(i = 448; i < 576; i++)
+                s->output[i] = cpe->ch[channel].saved[i] * pwindow[i - 448];
+            for(i = 576; i < 704; i++)
+                s->output[i] = cpe->ch[channel].saved[i];
+        }
+        if(cpe->ch[channel].ics.window_sequence[0] != LONG_START_SEQUENCE){
+            j = channel;
+            for (i = 0; i < 1024; i++, j += avctx->channels){
+                s->output[i+1024]         = audio[j] * lwindow[1024 - i - 1];
+                cpe->ch[channel].saved[i] = audio[j] * lwindow[i];
+            }
+        }else{
+            j = channel;
+            for(i = 0; i < 448; i++, j += avctx->channels)
+                s->output[i+1024]         = audio[j];
+            for(i = 448; i < 576; i++, j += avctx->channels)
+                s->output[i+1024]         = audio[j] * swindow[576 - i - 1];
+            memset(s->output+1024+576, 0, sizeof(s->output[0]) * 448);
+            j = channel;
+            for(i = 0; i < 1024; i++, j += avctx->channels)
+                cpe->ch[channel].saved[i] = audio[j];
+        }
+        ff_mdct_calc(&s->mdct1024, cpe->ch[channel].coeffs, s->output);
+    }else{
+        j = channel;
+        for (k = 0; k < 1024; k += 128) {
+            for(i = 448 + k; i < 448 + k + 256; i++)
+                s->output[i - 448 - k] = (i < 1024) ? cpe->ch[channel].saved[i] : audio[channel + (i-1024)*avctx->channels] / 512.0;
+            s->dsp.vector_fmul        (s->output,     k ?  swindow : pwindow, 128);
+            s->dsp.vector_fmul_reverse(s->output+128, s->output+128, swindow, 128);
+            ff_mdct_calc(&s->mdct128, cpe->ch[channel].coeffs + k, s->output);
+        }
+        j = channel;
+        for(i = 0; i < 1024; i++, j += avctx->channels)
+            cpe->ch[channel].saved[i] = audio[j];
+    }
+}
+
 /**
  * Encode ics_info element.
  * @see Table 4.6 (syntax of ics_info)
@@ -196,7 +310,7 @@
 static void put_ics_info(AVCodecContext *avctx, IndividualChannelStream *info)
 {
     AACEncContext *s = avctx->priv_data;
-    int i;
+    int wg;
 
     put_bits(&s->pb, 1, 0);                // ics_reserved bit
     put_bits(&s->pb, 2, info->window_sequence[0]);
@@ -206,12 +320,425 @@
         put_bits(&s->pb, 1, 0);            // no prediction
     }else{
         put_bits(&s->pb, 4, info->max_sfb);
-        for(i = 1; i < info->num_windows; i++)
-            put_bits(&s->pb, 1, info->group_len[i]);
+        for(wg = 0; wg < info->num_window_groups; wg++){
+            if(wg)
+                put_bits(&s->pb, 1, 0);
+            if(info->group_len[wg] > 1)
+                put_sbits(&s->pb, info->group_len[wg] - 1, 0xFF);
+        }
+    }
+}
+
+/**
+ * Encode MS data.
+ * @see 4.6.8.1 "Joint Coding - M/S Stereo"
+ */
+static void encode_ms_info(PutBitContext *pb, ChannelElement *cpe)
+{
+    int i, w, wg;
+
+    put_bits(pb, 2, cpe->ms.present);
+    if(cpe->ms.present == 1){
+        w = 0;
+        for(wg = 0; wg < cpe->ch[0].ics.num_window_groups; wg++){
+            for(i = 0; i < cpe->ch[0].ics.max_sfb; i++)
+                put_bits(pb, 1, cpe->ms.mask[w][i]);
+            w += cpe->ch[0].ics.group_len[wg];
+        }
     }
 }
 
 /**
+ * Return number of bits needed to write codebook run length value.
+ *
+ * @param run     run length
+ * @param bits    number of bits used to code value (5 for long frames, 3 for short frames)
+ */
+static av_always_inline int calculate_run_bits(int run, const int bits)
+{
+    int esc = (1 << bits) - 1;
+    return (1 + (run >= esc)) * bits;
+}
+
+/**
+ * Calculate the number of bits needed to code given band with given codebook.
+ *
+ * @param s       encoder context
+ * @param cpe     channel element
+ * @param channel channel number inside channel pair
+ * @param win     window group start number
+ * @param start   scalefactor band position in spectral coefficients
+ * @param size    scalefactor band size
+ * @param cb      codebook number
+ */
+static int calculate_band_bits(AACEncContext *s, ChannelElement *cpe, int channel, int win, int group_len, int start, int size, int cb)
+{
+    int i, j, w;
+    int score = 0, dim, idx, start2;
+    int range;
+
+    if(!cb) return 0;
+    cb--;
+    dim = (aac_cb_info[cb].flags & CB_PAIRS) ? 2 : 4;
+    if(aac_cb_info[cb].flags & CB_UNSIGNED)
+        range = aac_cb_info[cb].maxval + 1;
+    else
+        range = aac_cb_info[cb].maxval*2 + 1;
+
+    start2 = start;
+    if(aac_cb_info[cb].flags & CB_ESCAPE){
+        int coef_abs[2];
+        for(w = win; w < win + group_len; w++){
+            for(i = start2; i < start2 + size; i += dim){
+                idx = 0;
+                for(j = 0; j < dim; j++)
+                    coef_abs[j] = FFABS(cpe->ch[channel].icoefs[i+j]);
+                for(j = 0; j < dim; j++)
+                    idx = idx*17 + FFMIN(coef_abs[j], 16);
+                score += ff_aac_spectral_bits[cb][idx];
+                for(j = 0; j < dim; j++)
+                    if(cpe->ch[channel].icoefs[i+j])
+                        score++;
+                for(j = 0; j < dim; j++)
+                    if(coef_abs[j] > 15)
+                        score += av_log2(coef_abs[j]) * 2 - 4 + 1;
+            }
+            start2 += 128;
+       }
+    }else if(aac_cb_info[cb].flags & CB_UNSIGNED){
+        for(w = win; w < win + group_len; w++){
+            for(i = start2; i < start2 + size; i += dim){
+                idx = 0;
+                for(j = 0; j < dim; j++)
+                    idx = idx * range + FFABS(cpe->ch[channel].icoefs[i+j]);
+                score += ff_aac_spectral_bits[cb][idx];
+                for(j = 0; j < dim; j++)
+                     if(cpe->ch[channel].icoefs[i+j])
+                         score++;
+            }
+            start2 += 128;
+        }
+    }else{
+        for(w = win; w < win + group_len; w++){
+            for(i = start2; i < start2 + size; i += dim){
+                idx = 0;
+                for(j = 0; j < dim; j++)
+                    idx = idx * range + cpe->ch[channel].icoefs[i+j] + aac_cb_info[cb].maxval;
+                score += ff_aac_spectral_bits[cb][idx];
+            }
+            start2 += 128;
+        }
+    }
+    return score;
+}
+
+/**
+ * Encode band info for single window group bands.
+ */
+static void encode_window_bands_info(AACEncContext *s, ChannelElement *cpe, int channel, int win, int group_len){
+    int maxval;
+    int w, swb, cb, ccb, start, start2, size;
+    int i, j, k;
+    const int max_sfb = cpe->ch[channel].ics.max_sfb;
+    const int run_bits = cpe->ch[channel].ics.num_windows == 1 ? 5 : 3;
+    const int run_esc = (1 << run_bits) - 1;
+    int bits, idx, count;
+    int stack[64], stack_len;
+
+    start = win*128;
+    for(swb = 0; swb < max_sfb; swb++){
+        maxval = 0;
+        start2 = start;
+        size = cpe->ch[channel].ics.swb_sizes[swb];
+        if(cpe->ch[channel].zeroes[win][swb])
+            maxval = 0;
+        else{
+            for(w = win; w < win + group_len; w++){
+                for(i = start2; i < start2 + size; i++){
+                    maxval = FFMAX(maxval, FFABS(cpe->ch[channel].icoefs[i]));
+                }
+                start2 += 128;
+            }
+        }
+        for(cb = 0; cb < 12; cb++){
+            if(aac_cb_info[cb].maxval < maxval)
+                s->band_bits[swb][cb] = INT_MAX;
+            else
+                s->band_bits[swb][cb] = calculate_band_bits(s, cpe, channel, win, group_len, start, size, cb);
+        }
+        start += cpe->ch[channel].ics.swb_sizes[swb];
+    }
+    s->path[0].bits = 0;
+    for(i = 1; i <= max_sfb; i++)
+        s->path[i].bits = INT_MAX;
+    for(i = 0; i < max_sfb; i++){
+        for(j = 1; j <= max_sfb - i; j++){
+            bits = INT_MAX;
+            ccb = 0;
+            for(cb = 0; cb < 12; cb++){
+                int sum = 0;
+                for(k = 0; k < j; k++){
+                    if(s->band_bits[i + k][cb] == INT_MAX){
+                        sum = INT_MAX;
+                        break;
+                    }
+                    sum += s->band_bits[i + k][cb];
+                }
+                if(sum < bits){
+                    bits = sum;
+                    ccb  = cb;
+                }
+            }
+            assert(bits != INT_MAX);
+            bits += s->path[i].bits + calculate_run_bits(j, run_bits);
+            if(bits < s->path[i+j].bits){
+                s->path[i+j].bits     = bits;
+                s->path[i+j].codebook = ccb;
+                s->path[i+j].prev_idx = i;
+            }
+        }
+    }
+
+    //convert resulting path from backward-linked list
+    stack_len = 0;
+    idx = max_sfb;
+    while(idx > 0){
+        stack[stack_len++] = idx;
+        idx = s->path[idx].prev_idx;
+    }
+
+    //perform actual band info encoding
+    start = 0;
+    for(i = stack_len - 1; i >= 0; i--){
+        put_bits(&s->pb, 4, s->path[stack[i]].codebook);
+        count = stack[i] - s->path[stack[i]].prev_idx;
+        for(j = 0; j < count; j++){
+            cpe->ch[channel].band_type[win][start] =  s->path[stack[i]].codebook;
+            cpe->ch[channel].zeroes[win][start]    = !s->path[stack[i]].codebook;
+            start++;
+        }
+        while(count >= run_esc){
+            put_bits(&s->pb, run_bits, run_esc);
+            count -= run_esc;
+        }
+        put_bits(&s->pb, run_bits, count);
+    }
+}
+
+/**
+ * Encode one scalefactor band with selected codebook.
+ */
+static void encode_band_coeffs(AACEncContext *s, ChannelElement *cpe, int channel, int start, int size, int cb)
+{
+    const uint8_t  *bits  = ff_aac_spectral_bits [aac_cb_info[cb].cb_num];
+    const uint16_t *codes = ff_aac_spectral_codes[aac_cb_info[cb].cb_num];
+    const int dim = (aac_cb_info[cb].flags & CB_PAIRS) ? 2 : 4;
+    int i, j, idx, range;
+
+    if(!bits) return;
+
+    if(aac_cb_info[cb].flags & CB_UNSIGNED)
+        range = aac_cb_info[cb].maxval + 1;
+    else
+        range = aac_cb_info[cb].maxval*2 + 1;
+
+    if(aac_cb_info[cb].flags & CB_ESCAPE){
+        int coef_abs[2];
+        for(i = start; i < start + size; i += dim){
+            idx = 0;
+            for(j = 0; j < dim; j++)
+                coef_abs[j] = FFABS(cpe->ch[channel].icoefs[i+j]);
+            for(j = 0; j < dim; j++)
+                idx = idx*17 + FFMIN(coef_abs[j], 16);
+            put_bits(&s->pb, bits[idx], codes[idx]);
+            //output signs
+            for(j = 0; j < dim; j++)
+                if(cpe->ch[channel].icoefs[i+j])
+                    put_bits(&s->pb, 1, cpe->ch[channel].icoefs[i+j] < 0);
+            //output escape values
+            for(j = 0; j < dim; j++)
+                if(coef_abs[j] > 15){
+                    int len = av_log2(coef_abs[j]);
+
+                    put_bits(&s->pb, len - 4 + 1, (1 << (len - 4 + 1)) - 2);
+                    put_bits(&s->pb, len, coef_abs[j] & ((1 << len) - 1));
+                }
+        }
+    }else if(aac_cb_info[cb].flags & CB_UNSIGNED){
+        for(i = start; i < start + size; i += dim){
+            idx = 0;
+            for(j = 0; j < dim; j++)
+                idx = idx * range + FFABS(cpe->ch[channel].icoefs[i+j]);
+            put_bits(&s->pb, bits[idx], codes[idx]);
+            //output signs
+            for(j = 0; j < dim; j++)
+                if(cpe->ch[channel].icoefs[i+j])
+                    put_bits(&s->pb, 1, cpe->ch[channel].icoefs[i+j] < 0);
+        }
+    }else{
+        for(i = start; i < start + size; i += dim){
+            idx = 0;
+            for(j = 0; j < dim; j++)
+                idx = idx * range + cpe->ch[channel].icoefs[i+j] + aac_cb_info[cb].maxval;
+            put_bits(&s->pb, bits[idx], codes[idx]);
+        }
+    }
+}
+
+/**
+ * Encode scalefactor band coding type.
+ */
+static void encode_band_info(AVCodecContext *avctx, AACEncContext *s, ChannelElement *cpe, int channel)
+{
+    int w, wg;
+
+    w = 0;
+    for(wg = 0; wg < cpe->ch[channel].ics.num_window_groups; wg++){
+        encode_window_bands_info(s, cpe, channel, w, cpe->ch[channel].ics.group_len[wg]);
+        w += cpe->ch[channel].ics.group_len[wg];
+    }
+}
+
+/**
+ * Encode scalefactors.
+ */
+static void encode_scale_factors(AVCodecContext *avctx, AACEncContext *s, ChannelElement *cpe, int channel, int global_gain)
+{
+    int off = global_gain, diff;
+    int i, w, wg;
+
+    w = 0;
+    for(wg = 0; wg < cpe->ch[channel].ics.num_window_groups; wg++){
+        for(i = 0; i < cpe->ch[channel].ics.max_sfb; i++){
+            if(!cpe->ch[channel].zeroes[w][i]){
+                if(cpe->ch[channel].sf_idx[w][i] == 256) cpe->ch[channel].sf_idx[w][i] = off;
+                diff = cpe->ch[channel].sf_idx[w][i] - off + SCALE_DIFF_ZERO;
+                if(diff < 0 || diff > 120) av_log(avctx, AV_LOG_ERROR, "Scalefactor difference is too big to be coded\n");
+                off = cpe->ch[channel].sf_idx[w][i];
+                put_bits(&s->pb, ff_aac_scalefactor_bits[diff], ff_aac_scalefactor_code[diff]);
+            }
+        }
+        w += cpe->ch[channel].ics.group_len[wg];
+    }
+}
+
+/**
+ * Encode pulse data.
+ */
+static void encode_pulses(AVCodecContext *avctx, AACEncContext *s, Pulse *pulse, int channel)
+{
+    int i;
+
+    put_bits(&s->pb, 1, !!pulse->num_pulse);
+    if(!pulse->num_pulse) return;
+
+    put_bits(&s->pb, 2, pulse->num_pulse - 1);
+    put_bits(&s->pb, 6, pulse->start);
+    for(i = 0; i < pulse->num_pulse; i++){
+        put_bits(&s->pb, 5, pulse->offset[i]);
+        put_bits(&s->pb, 4, pulse->amp[i]);
+    }
+}
+
+/**
+ * Encode temporal noise shaping data.
+ */
+static void encode_tns_data(AVCodecContext *avctx, AACEncContext *s, ChannelElement *cpe, int channel)
+{
+    int i, w;
+
+    put_bits(&s->pb, 1, cpe->ch[channel].tns.present);
+    if(!cpe->ch[channel].tns.present) return;
+    if(cpe->ch[channel].ics.window_sequence[0] == EIGHT_SHORT_SEQUENCE){
+        for(w = 0; w < cpe->ch[channel].ics.num_windows; w++){
+            put_bits(&s->pb, 1, cpe->ch[channel].tns.n_filt[w]);
+            if(!cpe->ch[channel].tns.n_filt[w]) continue;
+            put_bits(&s->pb, 1, cpe->ch[channel].tns.coef_res[w] - 3);
+            put_bits(&s->pb, 4, cpe->ch[channel].tns.length[w][0]);
+            put_bits(&s->pb, 3, cpe->ch[channel].tns.order[w][0]);
+            if(cpe->ch[channel].tns.order[w][0]){
+                put_bits(&s->pb, 1, cpe->ch[channel].tns.direction[w][0]);
+                put_bits(&s->pb, 1, cpe->ch[channel].tns.coef_compress[w][0]);
+                for(i = 0; i < cpe->ch[channel].tns.order[w][0]; i++)
+                     put_bits(&s->pb, cpe->ch[channel].tns.coef_len[w][0], cpe->ch[channel].tns.coef[w][0][i]);
+            }
+        }
+    }else{
+        put_bits(&s->pb, 1, cpe->ch[channel].tns.n_filt[0]);
+        if(!cpe->ch[channel].tns.n_filt[0]) return;
+        put_bits(&s->pb, 1, cpe->ch[channel].tns.coef_res[0] - 3);
+        for(w = 0; w < cpe->ch[channel].tns.n_filt[0]; w++){
+            put_bits(&s->pb, 6, cpe->ch[channel].tns.length[0][w]);
+            put_bits(&s->pb, 5, cpe->ch[channel].tns.order[0][w]);
+            if(cpe->ch[channel].tns.order[0][w]){
+                put_bits(&s->pb, 1, cpe->ch[channel].tns.direction[0][w]);
+                put_bits(&s->pb, 1, cpe->ch[channel].tns.coef_compress[0][w]);
+                for(i = 0; i < cpe->ch[channel].tns.order[0][w]; i++)
+                     put_bits(&s->pb, cpe->ch[channel].tns.coef_len[0][w], cpe->ch[channel].tns.coef[0][w][i]);
+            }
+        }
+    }
+}
+
+/**
+ * Encode spectral coefficients processed by psychoacoustic model.
+ */
+static void encode_spectral_coeffs(AVCodecContext *avctx, AACEncContext *s, ChannelElement *cpe, int channel)
+{
+    int start, i, w, w2, wg;
+
+    w = 0;
+    for(wg = 0; wg < cpe->ch[channel].ics.num_window_groups; wg++){
+        start = 0;
+        for(i = 0; i < cpe->ch[channel].ics.max_sfb; i++){
+            if(cpe->ch[channel].zeroes[w][i]){
+                start += cpe->ch[channel].ics.swb_sizes[i];
+                continue;
+            }
+            for(w2 = w; w2 < w + cpe->ch[channel].ics.group_len[wg]; w2++){
+                encode_band_coeffs(s, cpe, channel, start + w2*128, cpe->ch[channel].ics.swb_sizes[i], cpe->ch[channel].band_type[w][i]);
+            }
+            start += cpe->ch[channel].ics.swb_sizes[i];
+        }
+        w += cpe->ch[channel].ics.group_len[wg];
+    }
+}
+
+/**
+ * Encode one channel of audio data.
+ */
+static int encode_individual_channel(AVCodecContext *avctx, ChannelElement *cpe, int channel)
+{
+    AACEncContext *s = avctx->priv_data;
+    int g, w, wg;
+    int global_gain;
+
+    //determine global gain as standard recommends - the first scalefactor value
+    global_gain = 0;
+    w = 0;
+    for(wg = 0; wg < cpe->ch[channel].ics.num_window_groups; wg++){
+        for(g = 0; g < cpe->ch[channel].ics.max_sfb; g++){
+            if(!cpe->ch[channel].zeroes[w][g]){
+                global_gain = cpe->ch[channel].sf_idx[w][g];
+                break;
+            }
+        }
+        if(global_gain) break;
+        w += cpe->ch[channel].ics.group_len[wg];
+    }
+
+    put_bits(&s->pb, 8, global_gain);
+    if(!cpe->common_window) put_ics_info(avctx, &cpe->ch[channel].ics);
+    encode_band_info(avctx, s, cpe, channel);
+    encode_scale_factors(avctx, s, cpe, channel, global_gain);
+    encode_pulses(avctx, s, &cpe->ch[channel].pulse, channel);
+    encode_tns_data(avctx, s, cpe, channel);
+    put_bits(&s->pb, 1, 0); //ssr
+    encode_spectral_coeffs(avctx, s, cpe, channel);
+    return 0;
+}
+
+/**
  * Write some auxiliary information about the created AAC file.
  */
 static void put_bitstream_info(AVCodecContext *avctx, AACEncContext *s, const char *name)
@@ -231,6 +758,80 @@
     put_bits(&s->pb, 12 - padbits, 0);
 }
 
+static int aac_encode_frame(AVCodecContext *avctx,
+                            uint8_t *frame, int buf_size, void *data)
+{
+    AACEncContext *s = avctx->priv_data;
+    int16_t *samples = s->samples, *samples2, *la;
+    ChannelElement *cpe;
+    int i, j, chans, tag, start_ch;
+    const uint8_t *chan_map = aac_chan_configs[avctx->channels-1];
+    int chan_el_counter[4];
+
+    if(s->last_frame)
+        return 0;
+    if(data){
+        if((s->psy.flags & PSY_MODEL_NO_PREPROC) == PSY_MODEL_NO_PREPROC){
+            memcpy(s->samples + 1024 * avctx->channels, data, 1024 * avctx->channels * sizeof(s->samples[0]));
+        }else{
+            start_ch = 0;
+            samples2 = s->samples + 1024 * avctx->channels;
+            for(i = 0; i < chan_map[0]; i++){
+                tag = chan_map[i+1];
+                chans = tag == ID_CPE ? 2 : 1;
+                ff_aac_psy_preprocess(&s->psy, (uint16_t*)data + start_ch, samples2 + start_ch, i, tag);
+                start_ch += chans;
+            }
+        }
+    }
+    if(!avctx->frame_number){
+        memcpy(s->samples, s->samples + 1024 * avctx->channels, 1024 * avctx->channels * sizeof(s->samples[0]));
+        return 0;
+    }
+
+    init_put_bits(&s->pb, frame, buf_size*8);
+    if(avctx->frame_number==1 && !(avctx->flags & CODEC_FLAG_BITEXACT)){
+        put_bitstream_info(avctx, s, LIBAVCODEC_IDENT);
+    }
+    start_ch = 0;
+    memset(chan_el_counter, 0, sizeof(chan_el_counter));
+    for(i = 0; i < chan_map[0]; i++){
+        tag = chan_map[i+1];
+        chans = tag == ID_CPE ? 2 : 1;
+        cpe = &s->cpe[i];
+        samples2 = samples + start_ch;
+        la = samples2 + 1024 * avctx->channels + start_ch;
+        if(!data) la = NULL;
+        ff_aac_psy_suggest_window(&s->psy, samples2, la, i, tag, cpe);
+        for(j = 0; j < chans; j++){
+            apply_window_and_mdct(avctx, s, cpe, samples2, j);
+        }
+        ff_aac_psy_analyze(&s->psy, i, tag, cpe);
+        put_bits(&s->pb, 3, tag);
+        put_bits(&s->pb, 4, chan_el_counter[tag]++);
+        if(chans == 2){
+            put_bits(&s->pb, 1, cpe->common_window);
+            if(cpe->common_window){
+                put_ics_info(avctx, &cpe->ch[0].ics);
+                encode_ms_info(&s->pb, cpe);
+            }
+        }
+        for(j = 0; j < chans; j++){
+            encode_individual_channel(avctx, cpe, j);
+        }
+        start_ch += chans;
+    }
+
+    put_bits(&s->pb, 3, ID_END);
+    flush_put_bits(&s->pb);
+    avctx->frame_bits = put_bits_count(&s->pb);
+
+    if(!data)
+        s->last_frame = 1;
+    memcpy(s->samples, s->samples + 1024 * avctx->channels, 1024 * avctx->channels * sizeof(s->samples[0]));
+    return put_bits_count(&s->pb)>>3;
+}
+
 static av_cold int aac_encode_end(AVCodecContext *avctx)
 {
     AACEncContext *s = avctx->priv_data;
-------------- next part --------------
/*
 * AAC encoder psychoacoustic model
 * Copyright (C) 2008 Konstantin Shishkov
 *
 * This file is part of FFmpeg.
 *
 * FFmpeg is free software; you can redistribute it and/or
 * modify it under the terms of the GNU Lesser General Public
 * License as published by the Free Software Foundation; either
 * version 2.1 of the License, or (at your option) any later version.
 *
 * FFmpeg is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * Lesser General Public License for more details.
 *
 * You should have received a copy of the GNU Lesser General Public
 * License along with FFmpeg; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
 */

#ifndef FFMPEG_AACPSY_H
#define FFMPEG_AACPSY_H

#include "avcodec.h"
#include "aac.h"
#include "lowpass.h"

enum AACPsyModelType{
    AAC_PSY_NULL,              ///< do nothing with frequencies
    AAC_PSY_NULL8,             ///< do nothing with frequencies but work with short windows
    AAC_PSY_3GPP,              ///< model following recommendations from 3GPP TS 26.403

    AAC_NB_PSY_MODELS          ///< total number of psychoacoustic models, since it's not a part of the ABI new models can be added freely
};

enum AACPsyModelMode{
    PSY_MODE_CBR,              ///< follow bitrate as closely as possible
    PSY_MODE_ABR,              ///< try to achieve bitrate but actual bitrate may differ significantly
    PSY_MODE_QUALITY,          ///< try to achieve set quality instead of bitrate
};

#define PSY_MODEL_MODE_MASK  0x0000000F ///< bit fields for storing mode (CBR, ABR, VBR)
#define PSY_MODEL_NO_SWITCH  0x00000020 ///< disable window switching
#define PSY_MODEL_NO_ST_ATT  0x00000040 ///< disable stereo attenuation
#define PSY_MODEL_NO_LOWPASS 0x00000080 ///< disable low-pass filtering

#define PSY_MODEL_NO_PREPROC (PSY_MODEL_NO_ST_ATT | PSY_MODEL_NO_LOWPASS)

#define PSY_MODEL_MODE(a)  ((a) & PSY_MODEL_MODE_MASK)

/**
 * context used by psychoacoustic model
 */
typedef struct AACPsyContext {
    AVCodecContext *avctx;            ///< encoder context

    int flags;                        ///< model flags
    const uint8_t *bands1024;         ///< scalefactor band sizes for long (1024 samples) frame
    int num_bands1024;                ///< number of scalefactor bands for long frame
    const uint8_t *bands128;          ///< scalefactor band sizes for short (128 samples) frame
    int num_bands128;                 ///< number of scalefactor bands for short frame

    const struct AACPsyModel *model;  ///< pointer to the psychoacoustic model implementation
    void* model_priv_data;            ///< psychoacoustic model implementation private data

    float stereo_att;                 ///< stereo attenuation factor
    LPFilterCoeffs lp_coeffs;         ///< lowpass filter coefficients
    LPFilterState *lp_state;          ///< lowpass filter state
}AACPsyContext;

typedef struct AACPsyModel {
    const char *name;
    int   (*init)   (AACPsyContext *apc, int elements);
    void  (*window) (AACPsyContext *apc, int16_t *audio, int16_t *la, int tag, int type, ChannelElement *cpe);
    void  (*process)(AACPsyContext *apc, int tag, int type, ChannelElement *cpe);
    void  (*end)    (AACPsyContext *apc);
}AACPsyModel;

/**
 * Initialize psychoacoustic model.
 *
 * @param ctx           model context
 * @param avctx         codec context
 * @param model         model implementation that will be used
 * @param elements      number of channel elements (single channel or channel pair) to handle by  model
 * @param flags         model flags, may be ignored by model if unsupported
 * @param bands1024     scalefactor band lengths for long (1024 samples) frame
 * @param num_bands1024 number of scalefactor bands for long frame
 * @param bands128      scalefactor band lengths for short (128 samples) frame
 * @param num_bands128  number of scalefactor bands for short frame
 *
 * @return zero if successful, a negative value if not
 */
int ff_aac_psy_init(AACPsyContext *ctx, AVCodecContext *avctx,
                    enum AACPsyModelType model, int elements, int flags,
                    const uint8_t *bands1024, int num_bands1024,
                    const uint8_t *bands128,  int num_bands128);

/**
 * Preprocess audio frame in order to compress it better.
 *
 * @param ctx   model context
 * @param audio samples to preprocess
 * @param dest  place to put filtered samples
 * @param tag   number of channel element to analyze
 * @param type  channel element type (e.g. ID_SCE or ID_CPE)
 */
void ff_aac_psy_preprocess(AACPsyContext *ctx, int16_t *audio, int16_t *dest, int tag, int type);

/**
 * Set window sequence and related parameters for channel element.
 *
 * @param ctx   model context
 * @param audio samples for the current frame
 * @param la    lookahead samples (NULL when unavailable)
 * @param tag   number of channel element to analyze
 * @param type  channel element type (e.g. ID_SCE or ID_CPE)
 * @param cpe   pointer to the current channel element
 */
void ff_aac_psy_suggest_window(AACPsyContext *ctx, int16_t *audio, int16_t *la, int tag, int type, ChannelElement *cpe);

/**
 * Perform psychoacoustic analysis and output coefficients in integer form
 * along with scalefactors, M/S flags, etc.
 *
 * @param ctx   model context
 * @param tag   number of channel element to analyze
 * @param type  channel element type (e.g. ID_SCE or ID_CPE)
 * @param cpe   pointer to the current channel element
 */
void ff_aac_psy_analyze(AACPsyContext *ctx, int tag, int type, ChannelElement *cpe);

/**
 * Cleanup model context at the end.
 *
 * @param ctx model context
 */
void ff_aac_psy_end(AACPsyContext *ctx);

#endif /* FFMPEG_AACPSY_H */

-------------- next part --------------
/*
 * AAC encoder psychoacoustic model
 * Copyright (C) 2008 Konstantin Shishkov
 *
 * This file is part of FFmpeg.
 *
 * FFmpeg is free software; you can redistribute it and/or
 * modify it under the terms of the GNU Lesser General Public
 * License as published by the Free Software Foundation; either
 * version 2.1 of the License, or (at your option) any later version.
 *
 * FFmpeg is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * Lesser General Public License for more details.
 *
 * You should have received a copy of the GNU Lesser General Public
 * License along with FFmpeg; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
 */

/**
 * @file aacpsy.c
 * AAC encoder psychoacoustic model
 */

#include "avcodec.h"
#include "aacpsy.h"
#include "aactab.h"

/***********************************
 *              TODOs:
 * General:
 * better audio preprocessing (add DC highpass filter?)
 * more psy models
 *
 * 3GPP-based psy model:
 * thresholds linearization after their modifications for attaining given bitrate
 * try other bitrate controlling mechanism (maybe use ratecontrol.c?)
 * control quality for quality-based output
 **********************************/

/**
 * 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)
{
    return av_clip((int)(pow(fabsf(coef) * Q, 0.75) + 0.4054), 0, 8191);
}

/**
 * Convert coefficients to integers.
 * @return sum of coefficients
 */
static inline int quantize_coeffs(float *in, int *out, int size, int scale_idx)
{
    int i, sign, sum = 0;
    const float Q = ff_aac_pow2sf_tab[200 - scale_idx + SCALE_ONE_POS - SCALE_DIV_512];
    for(i = 0; i < size; i++){
        sign = in[i] > 0.0;
        out[i] = quant(in[i], Q);
        sum += out[i];
        if(sign) out[i] = -out[i];
    }
    return sum;
}

static inline float calc_distortion(float *c, int size, int scale_idx)
{
    int i;
    int q;
    float coef, unquant, sum = 0.0f;
    const float Q  = ff_aac_pow2sf_tab[200 - scale_idx + SCALE_ONE_POS - SCALE_DIV_512];
    const float IQ = ff_aac_pow2sf_tab[200 + scale_idx - SCALE_ONE_POS + SCALE_DIV_512];
    for(i = 0; i < size; i++){
        coef = fabs(c[i]);
        q = quant(c[i], Q);
        unquant = (q * cbrt(q)) * IQ;
        sum += (coef - unquant) * (coef - unquant);
    }
    return sum;
}

/**
 * Produce integer coefficients from scalefactors provided by the model.
 */
static void psy_create_output(AACPsyContext *apc, ChannelElement *cpe, int chans)
{
    int i, w, w2, wg, g, ch;
    int start, sum, maxsfb, cmaxsfb;

    for(ch = 0; ch < chans; ch++){
        start = 0;
        maxsfb = 0;
        cpe->ch[ch].pulse.num_pulse = 0;
        for(w = 0; w < cpe->ch[ch].ics.num_windows; w++){
            for(g = 0; g < cpe->ch[ch].ics.num_swb; g++){
                sum = 0;
                //apply M/S
                if(!ch && cpe->ms.mask[w][g]){
                    for(i = 0; i < cpe->ch[ch].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];
                    }
                }
                if(!cpe->ch[ch].zeroes[w][g])
                    sum = quantize_coeffs(cpe->ch[ch].coeffs + start, cpe->ch[ch].icoefs + start, cpe->ch[ch].ics.swb_sizes[g], cpe->ch[ch].sf_idx[w][g]);
                else
                    memset(cpe->ch[ch].icoefs + start, 0, cpe->ch[ch].ics.swb_sizes[g] * sizeof(cpe->ch[0].icoefs[0]));
                cpe->ch[ch].zeroes[w][g] = !sum;
                start += cpe->ch[ch].ics.swb_sizes[g];
            }
            for(cmaxsfb = cpe->ch[ch].ics.num_swb; cmaxsfb > 0 && cpe->ch[ch].zeroes[w][cmaxsfb-1]; cmaxsfb--);
            maxsfb = FFMAX(maxsfb, cmaxsfb);
        }
        cpe->ch[ch].ics.max_sfb = maxsfb;

        //adjust zero bands for window groups
        w = 0;
        for(wg = 0; wg < cpe->ch[ch].ics.num_window_groups; wg++){
            for(g = 0; g < cpe->ch[ch].ics.max_sfb; g++){
                i = 1;
                for(w2 = 0; w2 < cpe->ch[ch].ics.group_len[wg]; w2++){
                    if(!cpe->ch[ch].zeroes[w + w2][g]){
                        i = 0;
                        break;
                    }
                }
                cpe->ch[ch].zeroes[w][g] = i;
            }
            w += cpe->ch[ch].ics.group_len[wg];
        }
    }

    if(chans > 1 && cpe->common_window){
        int msc = 0;
        cpe->ch[0].ics.max_sfb = FFMAX(cpe->ch[0].ics.max_sfb, cpe->ch[1].ics.max_sfb);
        cpe->ch[1].ics.max_sfb = cpe->ch[0].ics.max_sfb;
        for(w = 0; w < cpe->ch[0].ics.num_windows; w++)
            for(i = 0; i < cpe->ch[0].ics.max_sfb; i++)
                if(cpe->ms.mask[w][i]) msc++;
        if(msc == 0 || cpe->ch[0].ics.max_sfb == 0) cpe->ms.present = 0;
        else cpe->ms.present = msc < cpe->ch[0].ics.max_sfb ? 1 : 2;
    }
}

static void psy_null_window(AACPsyContext *apc, int16_t *audio, int16_t *la, int tag, int type, ChannelElement *cpe)
{
    int ch;
    int chans = type == ID_CPE ? 2 : 1;

    for(ch = 0; ch < chans; ch++){
        cpe->ch[ch].ics.window_sequence[0] = ONLY_LONG_SEQUENCE;
        cpe->ch[ch].ics.use_kb_window[0] = 1;
        cpe->ch[ch].ics.num_windows = 1;
        cpe->ch[ch].ics.swb_sizes = apc->bands1024;
        cpe->ch[ch].ics.num_swb = apc->num_bands1024;
        cpe->ch[ch].ics.num_window_groups = 1;
        cpe->ch[ch].ics.group_len[0] = 1;
    }
    cpe->common_window = cpe->ch[0].ics.use_kb_window[0] == cpe->ch[1].ics.use_kb_window[0];
}

static void psy_null_process(AACPsyContext *apc, int tag, int type, ChannelElement *cpe)
{
    int start;
    int ch, g, i;
    int minscale;
    int chans = type == ID_CPE ? 2 : 1;

    for(ch = 0; ch < chans; ch++){
        start = 0;
        for(g = 0; g < apc->num_bands1024; g++){
            float energy = 0.0f, ffac = 0.0f, thr, dist;

            for(i = 0; i < apc->bands1024[g]; i++){
                energy += cpe->ch[ch].coeffs[start+i]*cpe->ch[ch].coeffs[start+i];
                ffac += sqrt(FFABS(cpe->ch[ch].coeffs[start+i]));
            }
            thr = energy * 0.001258925f;
            cpe->ch[ch].sf_idx[ch][g] = 136;
            cpe->ch[ch].zeroes[ch][g] = (energy == 0.0);
            if(cpe->ch[ch].zeroes[ch][g]) continue;
            minscale = (int)(2.66667 * (log2(6.75*thr) - log2(ffac)));
            cpe->ch[ch].sf_idx[ch][g] = SCALE_ONE_POS - minscale;
            while(cpe->ch[ch].sf_idx[ch][g] > 3){
                dist = calc_distortion(cpe->ch[ch].coeffs + start, apc->bands1024[g], cpe->ch[ch].sf_idx[ch][g]);
                if(dist < thr) break;
                cpe->ch[ch].sf_idx[ch][g] -= 3;
            }
        }
    }
    for(ch = 0; ch < chans; ch++){
        minscale = 255;
        for(g = 0; g < apc->num_bands1024; g++)
            if(!cpe->ch[ch].zeroes[0][g])
                minscale = FFMIN(minscale, cpe->ch[ch].sf_idx[0][g]);
        for(g = 0; g < apc->num_bands1024; g++)
            if(!cpe->ch[ch].zeroes[0][g])
                cpe->ch[ch].sf_idx[0][g] = FFMIN(minscale + SCALE_MAX_DIFF, cpe->ch[ch].sf_idx[0][g]);
    }
    psy_create_output(apc, cpe, chans);
}

static void psy_null8_window(AACPsyContext *apc, int16_t *audio, int16_t *la, int tag, int type, ChannelElement *cpe)
{
    int ch, i;
    int chans = type == ID_CPE ? 2 : 1;

    for(ch = 0; ch < chans; ch++){
        int prev_seq = cpe->ch[ch].ics.window_sequence[1];
        cpe->ch[ch].ics.use_kb_window[1] = cpe->ch[ch].ics.use_kb_window[0];
        cpe->ch[ch].ics.window_sequence[1] = cpe->ch[ch].ics.window_sequence[0];
        switch(cpe->ch[ch].ics.window_sequence[0]){
        case ONLY_LONG_SEQUENCE:   if(prev_seq == ONLY_LONG_SEQUENCE)cpe->ch[ch].ics.window_sequence[0] = LONG_START_SEQUENCE;   break;
        case LONG_START_SEQUENCE:  cpe->ch[ch].ics.window_sequence[0] = EIGHT_SHORT_SEQUENCE; break;
        case EIGHT_SHORT_SEQUENCE: if(prev_seq == EIGHT_SHORT_SEQUENCE)cpe->ch[ch].ics.window_sequence[0] = LONG_STOP_SEQUENCE;  break;
        case LONG_STOP_SEQUENCE:   cpe->ch[ch].ics.window_sequence[0] = ONLY_LONG_SEQUENCE;   break;
        }

        if(cpe->ch[ch].ics.window_sequence[0] != EIGHT_SHORT_SEQUENCE){
            cpe->ch[ch].ics.use_kb_window[0] = 1;
            cpe->ch[ch].ics.num_windows = 1;
            cpe->ch[ch].ics.swb_sizes = apc->bands1024;
            cpe->ch[ch].ics.num_swb = apc->num_bands1024;
            cpe->ch[ch].ics.num_window_groups = 1;
            cpe->ch[ch].ics.group_len[0] = 1;
        }else{
            cpe->ch[ch].ics.use_kb_window[0] = 1;
            cpe->ch[ch].ics.num_windows = 8;
            cpe->ch[ch].ics.swb_sizes = apc->bands128;
            cpe->ch[ch].ics.num_swb = apc->num_bands128;
            cpe->ch[ch].ics.num_window_groups = 4;
            for(i = 0; i < 4; i++)
                cpe->ch[ch].ics.group_len[i] = 2;
        }
    }
    cpe->common_window = cpe->ch[0].ics.use_kb_window[0] == cpe->ch[1].ics.use_kb_window[0];
}

static void psy_null8_process(AACPsyContext *apc, int tag, int type, ChannelElement *cpe)
{
    int start;
    int w, ch, g, i;
    int chans = type == ID_CPE ? 2 : 1;

    //detect M/S
    if(chans > 1 && cpe->common_window){
        start = 0;
        for(w = 0; w < cpe->ch[0].ics.num_windows; w++){
            for(g = 0; g < cpe->ch[0].ics.num_swb; g++){
                float diff = 0.0f;

                for(i = 0; i < cpe->ch[0].ics.swb_sizes[g]; i++)
                    diff += fabs(cpe->ch[0].coeffs[start+i] - cpe->ch[1].coeffs[start+i]);
                cpe->ms.mask[w][g] = diff == 0.0;
            }
        }
    }
    for(ch = 0; ch < chans; ch++){
        for(w = 0; w < cpe->ch[ch].ics.num_windows; w++){
            for(g = 0; g < cpe->ch[ch].ics.num_swb; g++){
                cpe->ch[ch].sf_idx[w][g] = SCALE_ONE_POS;
                cpe->ch[ch].zeroes[w][g] = 0;
            }
        }
    }
    psy_create_output(apc, cpe, chans);
}

/**
 * constants for 3GPP AAC psychoacoustic model
 * @{
 */
#define PSY_3GPP_C1 3.0f                    // log2(8.0)
#define PSY_3GPP_C2 1.32192809488736234787f // log2(2.5)
#define PSY_3GPP_C3 0.55935730170421255071f // 1 - C2/C1

#define PSY_3GPP_SPREAD_LOW  1.5f // spreading factor for ascending threshold spreading  (15 dB/Bark)
#define PSY_3GPP_SPREAD_HI   3.0f // spreading factor for descending threshold spreading (30 dB/Bark)

#define PSY_3GPP_RPEMIN      0.01f
#define PSY_3GPP_RPELEV      2.0f
/**
 * @}
 */

/**
 * information for single band used by 3GPP TS26.403-inspired psychoacoustic model
 */
typedef struct Psy3gppBand{
    float energy;    ///< band energy
    float ffac;      ///< form factor
    float thr;       ///< energy threshold
    float pe;        ///< perceptual entropy
    float a;         ///< constant part in perceptual entropy
    float b;         ///< variable part in perceptual entropy
    float nl;        ///< predicted number of lines left after quantization
    float min_snr;   ///< minimal SNR
    float thr_quiet; ///< threshold in quiet
}Psy3gppBand;

/**
 * single/pair channel context for psychoacoustic model
 */
typedef struct Psy3gppChannel{
    float       a[2];                       ///< parameter used for perceptual entropy - constant part
    float       b[2];                       ///< parameter used for perceptual entropy - variable part
    float       pe[2];                      ///< channel perceptual entropy
    float       thr[2];                     ///< channel thresholds sum
    Psy3gppBand band[2][128];               ///< bands information
    Psy3gppBand prev_band[2][128];          ///< bands information from the previous frame

    float       win_nrg[2];                 ///< sliding average of channel energy
    float       iir_state[2][2];            ///< hi-pass IIR filter state
    uint8_t     next_grouping[2];           ///< stored grouping scheme for the next frame (in case of 8 short window sequence)
    enum WindowSequence next_window_seq[2]; ///< window sequence to be used in the next frame
}Psy3gppChannel;

/**
 * 3GPP TS26.403-inspired psychoacoustic model specific data
 */
typedef struct Psy3gppContext{
    float       barks [1024]; ///< Bark value for each spectral line
    float       bark_l[64];   ///< Bark value for each spectral band in long frame
    float       bark_s[16];   ///< Bark value for each spectral band in short frame
    float       s_low_l[64];  ///< spreading factor for low-to-high threshold spreading in long frame
    float       s_low_s[16];  ///< spreading factor for low-to-high threshold spreading in short frame
    float       s_hi_l [64];  ///< spreading factor for high-to-low threshold spreading in long frame
    float       s_hi_s [16];  ///< spreading factor for high-to-low threshold spreading in short frame
    int         reservoir;    ///< bit reservoir fullness
    int         avg_bits;     ///< average frame size of bits for CBR
    float       ath_l[64];    ///< absolute threshold of hearing per bands in long frame
    float       ath_s[16];    ///< absolute threshold of hearing per bands in short frame
    Psy3gppChannel *ch;
}Psy3gppContext;

/**
 * Calculate Bark value for given line.
 */
static inline float calc_bark(float f)
{
    return 13.3f * atanf(0.00076f * f) + 3.5f * atanf((f / 7500.0f) * (f / 7500.0f));
}

#define ATH_ADD 4
/**
 * Calculate ATH value for given frequency.
 * Borrowed from Lame.
 */
static inline float ath(float f, float add)
{
    f /= 1000.0f;
    return   3.64 * pow(f, -0.8)
            - 6.8  * exp(-0.6  * (f - 3.4) * (f - 3.4))
            + 6.0  * exp(-0.15 * (f - 8.7) * (f - 8.7))
            + (0.6 + 0.04 * add) * 0.001 * f * f * f * f;
}

static av_cold int psy_3gpp_init(AACPsyContext *apc, int elements)
{
    Psy3gppContext *pctx;
    int i, g, start;
    float prev, minscale, minath;
    apc->model_priv_data = av_mallocz(sizeof(Psy3gppContext));
    pctx = (Psy3gppContext*) apc->model_priv_data;

    for(i = 0; i < 1024; i++)
        pctx->barks[i] = calc_bark(i * apc->avctx->sample_rate / 2048.0);
    i = 0;
    prev = 0.0;
    for(g = 0; g < apc->num_bands1024; g++){
        i += apc->bands1024[g];
        pctx->bark_l[g] = (pctx->barks[i - 1] + prev) / 2.0;
        prev = pctx->barks[i - 1];
    }
    for(g = 0; g < apc->num_bands1024 - 1; g++){
        pctx->s_low_l[g] = pow(10.0, -(pctx->bark_l[g+1] - pctx->bark_l[g]) * PSY_3GPP_SPREAD_LOW);
        pctx->s_hi_l [g] = pow(10.0, -(pctx->bark_l[g+1] - pctx->bark_l[g]) * PSY_3GPP_SPREAD_HI);
    }
    i = 0;
    prev = 0.0;
    for(g = 0; g < apc->num_bands128; g++){
        i += apc->bands128[g];
        pctx->bark_s[g] = (pctx->barks[i - 1] + prev) / 2.0;
        prev = pctx->barks[i - 1];
    }
    for(g = 0; g < apc->num_bands128 - 1; g++){
        pctx->s_low_s[g] = pow(10.0, -(pctx->bark_s[g+1] - pctx->bark_s[g]) * PSY_3GPP_SPREAD_LOW);
        pctx->s_hi_s [g] = pow(10.0, -(pctx->bark_s[g+1] - pctx->bark_s[g]) * PSY_3GPP_SPREAD_HI);
    }
    start = 0;
    minath = ath(3410, ATH_ADD);
    for(g = 0; g < apc->num_bands1024; g++){
        minscale = ath(apc->avctx->sample_rate * start / 1024.0, ATH_ADD);
        for(i = 1; i < apc->bands1024[g]; i++){
            minscale = fminf(minscale, ath(apc->avctx->sample_rate * (start + i) / 1024.0 / 2.0, ATH_ADD));
        }
        pctx->ath_l[g] = minscale - minath;
        start += apc->bands1024[g];
    }
    start = 0;
    for(g = 0; g < apc->num_bands128; g++){
        minscale = ath(apc->avctx->sample_rate * start / 1024.0, ATH_ADD);
        for(i = 1; i < apc->bands128[g]; i++){
            minscale = fminf(minscale, ath(apc->avctx->sample_rate * (start + i) / 1024.0 / 2.0, ATH_ADD));
        }
        pctx->ath_s[g] = minscale - minath;
        start += apc->bands128[g];
    }

    pctx->avg_bits = apc->avctx->bit_rate * 1024 / apc->avctx->sample_rate;
    pctx->ch = av_mallocz(sizeof(Psy3gppChannel) * elements);
    return 0;
}

/**
 * IIR filter used in block switching decision
 */
static float iir_filter(int in, float state[2])
{
    float ret;

    ret = 0.7548f * (in - state[0]) + 0.5095f * state[1];
    state[0] = in;
    state[1] = ret;
    return ret;
}

/**
 * window grouping information stored as bits (0 - new group, 1 - group continues)
 */
static const uint8_t window_grouping[9] = {
    0xB6, 0x6C, 0xD8, 0xB2, 0x66, 0xC6, 0x96, 0x36, 0x36
};

/**
 * Tell encoder which window types to use.
 * @see 3GPP TS26.403 5.4.1 "Blockswitching"
 */
static void psy_3gpp_window(AACPsyContext *apc, int16_t *audio, int16_t *la, int tag, int type, ChannelElement *cpe)
{
    int ch;
    int chans = type == ID_CPE ? 2 : 1;
    int i, j;
    int br = apc->avctx->bit_rate / apc->avctx->channels;
    int attack_ratio = (br <= 16000 + 8000*chans) ? 18 : 10;
    Psy3gppContext *pctx = (Psy3gppContext*) apc->model_priv_data;
    Psy3gppChannel *pch = &pctx->ch[tag];
    uint8_t grouping[2];
    enum WindowSequence win[2];

    if(la && !(apc->flags & PSY_MODEL_NO_SWITCH)){
        float s[8], v;
        for(ch = 0; ch < chans; ch++){
            enum WindowSequence last_window_sequence = cpe->ch[ch].ics.window_sequence[0];
            int switch_to_eight = 0;
            float sum = 0.0, sum2 = 0.0;
            int attack_n = 0;
            for(i = 0; i < 8; i++){
                for(j = 0; j < 128; j++){
                    v = iir_filter(audio[(i*128+j)*apc->avctx->channels+ch], pch->iir_state[ch]);
                    sum += v*v;
                }
                s[i] = sum;
                sum2 += sum;
            }
            for(i = 0; i < 8; i++){
                if(s[i] > pch->win_nrg[ch] * attack_ratio){
                    attack_n = i + 1;
                    switch_to_eight = 1;
                    break;
                }
            }
            pch->win_nrg[ch] = pch->win_nrg[ch]*7/8 + sum2/64;

            switch(last_window_sequence){
            case ONLY_LONG_SEQUENCE:
                win[ch] = switch_to_eight ? LONG_START_SEQUENCE : ONLY_LONG_SEQUENCE;
                grouping[ch] = 0;
                break;
            case LONG_START_SEQUENCE:
                win[ch] = EIGHT_SHORT_SEQUENCE;
                grouping[ch] = pch->next_grouping[ch];
                break;
            case LONG_STOP_SEQUENCE:
                win[ch] = ONLY_LONG_SEQUENCE;
                grouping[ch] = 0;
                break;
            case EIGHT_SHORT_SEQUENCE:
                win[ch] = switch_to_eight ? EIGHT_SHORT_SEQUENCE : LONG_STOP_SEQUENCE;
                grouping[ch] = switch_to_eight ? pch->next_grouping[ch] : 0;
                break;
            }
            pch->next_grouping[ch] = window_grouping[attack_n];
        }
    }else{
        for(ch = 0; ch < chans; ch++){
            win[ch] = (cpe->ch[ch].ics.window_sequence[0] == EIGHT_SHORT_SEQUENCE) ? EIGHT_SHORT_SEQUENCE : ONLY_LONG_SEQUENCE;
            grouping[ch] = (cpe->ch[ch].ics.window_sequence[0] == EIGHT_SHORT_SEQUENCE) ? window_grouping[0] : 0;
        }
    }

    for(ch = 0; ch < chans; ch++){
        cpe->ch[ch].ics.window_sequence[0] = win[ch];
        cpe->ch[ch].ics.use_kb_window[0] = 1;
        if(win[ch] != EIGHT_SHORT_SEQUENCE){
            cpe->ch[ch].ics.num_windows = 1;
            cpe->ch[ch].ics.swb_sizes = apc->bands1024;
            cpe->ch[ch].ics.num_swb = apc->num_bands1024;
            cpe->ch[ch].ics.num_window_groups = 1;
            cpe->ch[ch].ics.group_len[0] = 1;
        }else{
            cpe->ch[ch].ics.num_windows = 8;
            cpe->ch[ch].ics.swb_sizes = apc->bands128;
            cpe->ch[ch].ics.num_swb = apc->num_bands128;
            cpe->ch[ch].ics.num_window_groups = 0;
            cpe->ch[ch].ics.group_len[0] = 1;
            for(i = 0; i < 8; i++){
                if((grouping[ch] >> i) & 1){
                    cpe->ch[ch].ics.group_len[cpe->ch[ch].ics.num_window_groups - 1]++;
                }else{
                    cpe->ch[ch].ics.num_window_groups++;
                    cpe->ch[ch].ics.group_len[cpe->ch[ch].ics.num_window_groups - 1] = 1;
                }
            }
        }
    }
    cpe->common_window = chans > 1 && cpe->ch[0].ics.window_sequence[0] == cpe->ch[1].ics.window_sequence[0] && cpe->ch[0].ics.use_kb_window[0] == cpe->ch[1].ics.use_kb_window[0];
    if(cpe->common_window && cpe->ch[0].ics.window_sequence[0] == EIGHT_SHORT_SEQUENCE && grouping[0] != grouping[1])
        cpe->common_window = 0;
    if(PSY_MODEL_MODE(apc->flags) > PSY_MODE_QUALITY){
        av_log(apc->avctx, AV_LOG_ERROR, "Unknown mode %d, defaulting to CBR\n", PSY_MODEL_MODE(apc->flags));
    }
}

/**
 * Modify threshold by adding some value in loudness domain.
 * @see 3GPP TS26.403 5.6.1.1.1 "Addition of noise with equal loudness"
 */
static inline float modify_thr(float thr, float r){
    float t;
    t = pow(thr, 0.25) + r;
    return t*t*t*t;
}

/**
 * Calculate perceptual entropy and its corresponding values for one band.
 * @see 3GPP TS26.403 5.6.1.3 "Calculation of the reduction value"
 */
static void calc_pe(Psy3gppBand *band, int band_width)
{
    if(band->energy <= band->thr){
        band->a  = 0.0f;
        band->b  = 0.0f;
        band->nl = 0.0f;
        return;
    }
    band->nl = band->ffac / pow(band->energy/band_width, 0.25);
    if(band->energy >= band->thr * 8.0){
        band->a = band->nl * log2(band->energy);
        band->b = band->nl;
    }else{
        band->a = band->nl * (PSY_3GPP_C2 + PSY_3GPP_C3 * log2(band->energy));
        band->b = band->nl * PSY_3GPP_C3;
    }
    band->pe = band->a - band->b * log2(band->thr);
    band->min_snr = 1.0 / (pow(2.0, band->pe / band_width) - 1.5);
    if(band->min_snr < 1.26f)     band->min_snr = 1.26f;
    if(band->min_snr > 316.2277f) band->min_snr = 316.2277f;
}

/**
 * Determine scalefactors and prepare coefficients for encoding.
 * @see 3GPP TS26.403 5.4 "Psychoacoustic model"
 */
static void psy_3gpp_process(AACPsyContext *apc, int tag, int type, ChannelElement *cpe)
{
    int start;
    int ch, w, wg, g, g2, i;
    int prev_scale;
    Psy3gppContext *pctx = (Psy3gppContext*) apc->model_priv_data;
    float pe_target;
    int bits_avail;
    int chans = type == ID_CPE ? 2 : 1;
    Psy3gppChannel *pch = &pctx->ch[tag];

    //calculate energies, initial thresholds and related values - 5.4.2 "Threshold Calculation"
    memset(pch->band, 0, sizeof(pch->band));
    for(ch = 0; ch < chans; ch++){
        start = 0;
        for(w = 0; w < cpe->ch[ch].ics.num_windows; w++){
            for(g = 0; g < cpe->ch[ch].ics.num_swb; g++){
                g2 = w*16 + g;
                for(i = 0; i < cpe->ch[ch].ics.swb_sizes[g]; i++)
                    pch->band[ch][g2].energy +=  cpe->ch[ch].coeffs[start+i] *  cpe->ch[ch].coeffs[start+i];
                pch->band[ch][g2].energy /= 262144.0f;
                pch->band[ch][g2].thr = pch->band[ch][g2].energy * 0.001258925f;
                start += cpe->ch[ch].ics.swb_sizes[g];
                if(pch->band[ch][g2].energy != 0.0){
                    float ffac = 0.0;

                    for(i = 0; i < cpe->ch[ch].ics.swb_sizes[g]; i++)
                        ffac += sqrt(FFABS(cpe->ch[ch].coeffs[start+i]));
                    pch->band[ch][g2].ffac = ffac / sqrt(512.0);
                }
            }
        }
    }

    //modify thresholds - spread, threshold in quiet - 5.4.3 "Spreaded Energy Calculation"
    for(ch = 0; ch < chans; ch++){
        for(w = 0; w < cpe->ch[ch].ics.num_windows; w++){
            for(g = 1; g < cpe->ch[ch].ics.num_swb; g++){
                g2 = w*16 + g;
                if(cpe->ch[ch].ics.num_swb == apc->num_bands1024)
                    pch->band[ch][g2].thr = FFMAX(pch->band[ch][g2].thr, pch->band[ch][g2-1].thr * pctx->s_low_l[g-1]);
                else
                    pch->band[ch][g2].thr = FFMAX(pch->band[ch][g2].thr, pch->band[ch][g2-1].thr * pctx->s_low_s[g-1]);
            }
            for(g = cpe->ch[ch].ics.num_swb - 2; g >= 0; g--){
                g2 = w*16 + g;
                if(cpe->ch[ch].ics.num_swb == apc->num_bands1024)
                    pch->band[ch][g2].thr = FFMAX(pch->band[ch][g2].thr, pch->band[ch][g2+1].thr * pctx->s_hi_l[g+1]);
                else
                    pch->band[ch][g2].thr = FFMAX(pch->band[ch][g2].thr, pch->band[ch][g2+1].thr * pctx->s_hi_s[g+1]);
            }
            for(g = 0; g < cpe->ch[ch].ics.num_swb; g++){
                g2 = w*16 + g;
                if(cpe->ch[ch].ics.num_swb == apc->num_bands1024)
                    pch->band[ch][g2].thr_quiet = FFMAX(pch->band[ch][g2].thr, pctx->ath_l[g]);
                else
                    pch->band[ch][g2].thr_quiet = FFMAX(pch->band[ch][g2].thr, pctx->ath_s[g]);
                pch->band[ch][g2].thr_quiet = fmaxf(PSY_3GPP_RPEMIN*pch->band[ch][g2].thr_quiet, fminf(pch->band[ch][g2].thr_quiet, PSY_3GPP_RPELEV*pch->prev_band[ch][g2].thr_quiet));
                pch->band[ch][g2].thr = FFMAX(pch->band[ch][g2].thr, pch->band[ch][g2].thr_quiet * 0.25);
            }
        }
    }

    // M/S detection - 5.5.2 "Mid/Side Stereo"
    if(chans > 1 && cpe->common_window){
        start = 0;
        for(w = 0; w < cpe->ch[0].ics.num_windows; w++){
            for(g = 0; g < cpe->ch[0].ics.num_swb; g++){
                double en_m = 0.0, en_s = 0.0, ff_m = 0.0, ff_s = 0.0, l1;
                float m, s;

                g2 = w*16 + g;
                cpe->ms.mask[w][g] = 0;
                if(pch->band[0][g2].energy == 0.0 || pch->band[1][g2].energy == 0.0)
                    continue;
                for(i = 0; i < cpe->ch[0].ics.swb_sizes[g]; i++){
                    m = (cpe->ch[0].coeffs[start+i] + cpe->ch[1].coeffs[start+i]) / 2.0;
                    s = (cpe->ch[0].coeffs[start+i] - cpe->ch[1].coeffs[start+i]) / 2.0;
                    en_m += m*m;
                    en_s += s*s;
                    ff_m += sqrt(FFABS(m));
                    ff_s += sqrt(FFABS(s));
                }
                en_m /= 262144.0;
                en_s /= 262144.0;
                ff_m /= sqrt(512.0);
                ff_s /= sqrt(512.0);
                l1 = FFMIN(pch->band[0][g2].thr, pch->band[1][g2].thr);
                if(en_m == 0.0 || en_s == 0.0 || l1*l1 / (en_m * en_s) >= (pch->band[0][g2].thr * pch->band[1][g2].thr / (pch->band[0][g2].energy * pch->band[1][g2].energy))){
                    cpe->ms.mask[w][g] = 1;
                    pch->band[0][g2].energy = en_m;
                    pch->band[1][g2].energy = en_s;
                    pch->band[0][g2].ffac = ff_m;
                    pch->band[1][g2].ffac = ff_s;
                    pch->band[0][g2].thr = en_m * 0.001258925f;
                    pch->band[1][g2].thr = en_s * 0.001258925f;
                }
            }
        }
    }

    for(ch = 0; ch < chans; ch++){
        pch->a[ch] = pch->b[ch] = pch->pe[ch] = pch->thr[ch] = 0.0f;
        for(w = 0; w < cpe->ch[ch].ics.num_windows; w++){
            for(g = 0; g < cpe->ch[ch].ics.num_swb; g++){
                g2 = w*16 + g;
                if(pch->band[ch][g2].energy != 0.0)
                    calc_pe(&pch->band[ch][g2], cpe->ch[ch].ics.swb_sizes[g]);
                if(pch->band[ch][g2].thr < pch->band[ch][g2].energy){
                    pch->a[ch]   += pch->band[ch][g2].a;
                    pch->b[ch]   += pch->band[ch][g2].b;
                    pch->pe[ch]  += pch->band[ch][g2].pe;
                    pch->thr[ch] += pch->band[ch][g2].thr;
                }
            }
        }
    }

    switch(PSY_MODEL_MODE(apc->flags)){
    case PSY_MODE_CBR:
    case PSY_MODE_ABR:
        //bitrate reduction - 5.6.1 "Reduction of psychoacoustic requirements"
        if(PSY_MODEL_MODE(apc->flags) != PSY_MODE_ABR){
            pctx->reservoir += pctx->avg_bits - apc->avctx->frame_bits;
            bits_avail = pctx->avg_bits + pctx->reservoir;
            bits_avail = FFMIN(bits_avail, pctx->avg_bits * 1.5);
            pe_target = 1.18f * bits_avail / apc->avctx->channels * chans;
        }else{
            pe_target = pctx->avg_bits / apc->avctx->channels * chans;
        }
        for(i = 0; i < 2; i++){
            float t0, pe, r, a0 = 0.0f, pe0 = 0.0f, b0 = 0.0f;
            for(ch = 0; ch < chans; ch++){
                a0  += pch->a[ch];
                b0  += pch->b[ch];
                pe0 += pch->pe[ch];
            }
            if(pe0 == 0.0f) break;
            t0 = pow(2.0, (a0 - pe0)       / (4.0 * b0));
            r  = pow(2.0, (a0 - pe_target) / (4.0 * b0)) - t0;

            //add correction factor to thresholds and recalculate perceptual entropy
            for(ch = 0; ch < chans; ch++){
                pch->a[ch] = pch->b[ch] = pch->pe[ch] = pch->thr[ch] = 0.0;
                pe = 0.0f;
                for(w = 0; w < cpe->ch[ch].ics.num_windows; w++){
                    for(g = 0; g < cpe->ch[ch].ics.num_swb; g++){
                        g2 = w*16 + g;
                        pch->band[ch][g2].thr = modify_thr(pch->band[ch][g2].thr, r);
                        calc_pe(&pch->band[ch][g2], cpe->ch[ch].ics.swb_sizes[g]);
                        if(pch->band[ch][g2].thr < pch->band[ch][g2].energy){
                            pch->a[ch]   += pch->band[ch][g2].a;
                            pch->b[ch]   += pch->band[ch][g2].b;
                            pch->pe[ch]  += pch->band[ch][g2].pe;
                            pch->thr[ch] += pch->band[ch][g2].thr;
                        }
                    }
                }
            }
        }

        //determine scalefactors - 5.6.2 "Scalefactor determination"
        for(ch = 0; ch < chans; ch++){
            prev_scale = -1;
            for(w = 0; w < cpe->ch[ch].ics.num_windows; w++){
                for(g = 0; g < cpe->ch[ch].ics.num_swb; g++){
                    g2 = w*16 + g;
                    cpe->ch[ch].zeroes[w][g] = pch->band[ch][g2].thr >= pch->band[ch][g2].energy;
                    if(cpe->ch[ch].zeroes[w][g]) continue;
                    //spec gives constant for lg() but we scaled it for log2()
                    cpe->ch[ch].sf_idx[w][g] = (int)(2.66667 * (log2(6.75*pch->band[ch][g2].thr) - log2(pch->band[ch][g2].ffac)));
                    if(prev_scale != -1)
                        cpe->ch[ch].sf_idx[w][g] = av_clip(cpe->ch[ch].sf_idx[w][g], prev_scale - SCALE_MAX_DIFF, prev_scale + SCALE_MAX_DIFF);
                    prev_scale = cpe->ch[ch].sf_idx[w][g];
                }
            }
        }
        break;
    case PSY_MODE_QUALITY:
        for(ch = 0; ch < chans; ch++){
            start = 0;
            for(w = 0; w < cpe->ch[ch].ics.num_windows; w++){
                for(g = 0; g < cpe->ch[ch].ics.num_swb; g++){
                    g2 = w*16 + g;
                    if(pch->band[ch][g2].thr >= pch->band[ch][g2].energy){
                        cpe->ch[ch].sf_idx[w][g] = 0;
                        cpe->ch[ch].zeroes[w][g] = 1;
                    }else{
                        cpe->ch[ch].zeroes[w][g] = 0;
                        cpe->ch[ch].sf_idx[w][g] = (int)(2.66667 * (log2(6.75*pch->band[ch][g2].thr) - log2(pch->band[ch][g2].ffac)));
                        while(cpe->ch[ch].sf_idx[ch][g] > 3){
                            float dist = calc_distortion(cpe->ch[ch].coeffs + start, cpe->ch[ch].ics.swb_sizes[g], SCALE_ONE_POS + cpe->ch[ch].sf_idx[ch][g]);
                            if(dist < pch->band[ch][g2].thr) break;
                            cpe->ch[ch].sf_idx[ch][g] -= 3;
                        }
                    }
                    start += cpe->ch[ch].ics.swb_sizes[g];
                }
            }
        }
        break;
    }

    //limit scalefactors
    for(ch = 0; ch < chans; ch++){
        int min_scale = 256;
        for(w = 0; w < cpe->ch[ch].ics.num_windows; w++)
            for(g = 0; g < cpe->ch[ch].ics.num_swb; g++){
                if(cpe->ch[ch].zeroes[w][g]) continue;
                min_scale = FFMIN(min_scale, cpe->ch[ch].sf_idx[w][g]);
            }
        for(w = 0; w < cpe->ch[ch].ics.num_windows; w++)
            for(g = 0; g < cpe->ch[ch].ics.num_swb; g++){
                if(cpe->ch[ch].zeroes[w][g]) continue;
                cpe->ch[ch].sf_idx[w][g] = FFMIN(cpe->ch[ch].sf_idx[w][g], min_scale + SCALE_MAX_DIFF);
            }
        for(w = 0; w < cpe->ch[ch].ics.num_windows; w++)
            for(g = 0; g < cpe->ch[ch].ics.num_swb; g++){
                if(cpe->ch[ch].zeroes[w][g])
                    cpe->ch[ch].sf_idx[w][g] = 256;
                else
                    cpe->ch[ch].sf_idx[w][g] = av_clip(SCALE_ONE_POS + cpe->ch[ch].sf_idx[w][g], 0, SCALE_MAX_POS);
            }

        //adjust scalefactors for window groups
        w = 0;
        for(wg = 0; wg < cpe->ch[ch].ics.num_window_groups; wg++){
            int min_scale = 256;

            for(g = 0; g < cpe->ch[ch].ics.num_swb; g++){
                for(i = w; i < w + cpe->ch[ch].ics.group_len[wg]; i++){
                    if(cpe->ch[ch].zeroes[i][g]) continue;
                    min_scale = FFMIN(min_scale, cpe->ch[ch].sf_idx[i][g]);
                }
                for(i = w; i < w + cpe->ch[ch].ics.group_len[wg]; i++)
                    cpe->ch[ch].sf_idx[i][g] = min_scale;
            }
            w += cpe->ch[ch].ics.group_len[wg];
        }
    }

    memcpy(pch->prev_band, pch->band, sizeof(pch->band));
    psy_create_output(apc, cpe, chans);
}

static av_cold void psy_3gpp_end(AACPsyContext *apc)
{
    Psy3gppContext *pctx = (Psy3gppContext*) apc->model_priv_data;
    av_freep(&pctx->ch);
    av_freep(&apc->model_priv_data);
}

static const AACPsyModel psy_models[AAC_NB_PSY_MODELS] =
{
    {
       "Null model",
        NULL,
        psy_null_window,
        psy_null_process,
        NULL,
    },
    {
       "Null model - short windows",
        NULL,
        psy_null8_window,
        psy_null8_process,
        NULL,
    },
    {
       "3GPP TS 26.403-inspired model",
        psy_3gpp_init,
        psy_3gpp_window,
        psy_3gpp_process,
        psy_3gpp_end,
    },
};

int av_cold ff_aac_psy_init(AACPsyContext *ctx, AVCodecContext *avctx,
                            enum AACPsyModelType model, int elements, int flags,
                            const uint8_t *bands1024, int num_bands1024,
                            const uint8_t *bands128,  int num_bands128)
{
    int i;

    if(model >= AAC_NB_PSY_MODELS || !psy_models[model].window || !psy_models[model].process){
         av_log(avctx, AV_LOG_ERROR, "Invalid psy model\n");
         return -1;
    }

#ifndef CONFIG_HARDCODED_TABLES
   for (i = 0; i < 316; i++)
        ff_aac_pow2sf_tab[i] = pow(2, (i - 200)/4.);
#endif /* CONFIG_HARDCODED_TABLES */

    ctx->avctx = avctx;
    ctx->flags = flags;
    ctx->bands1024 = bands1024;
    ctx->num_bands1024 = num_bands1024;
    ctx->bands128 = bands128;
    ctx->num_bands128 = num_bands128;
    ctx->model = &psy_models[model];

    if(ctx->flags & PSY_MODEL_NO_ST_ATT || PSY_MODEL_MODE(ctx->flags) == PSY_MODE_QUALITY){
        ctx->flags |= PSY_MODEL_NO_ST_ATT;
        ctx->stereo_att = 0.5f;
    }else{
        ctx->stereo_att = av_clipf(avctx->bit_rate / elements / 192000.0, 0.0f, 0.5f);
    }
    if(ctx->flags & PSY_MODEL_NO_LOWPASS || PSY_MODEL_MODE(ctx->flags) == PSY_MODE_QUALITY){
        ctx->flags |= PSY_MODEL_NO_LOWPASS;
    }else{
        int cutoff;
        cutoff = avctx->bit_rate / elements / 8;
        if(ff_lowpass_filter_init_coeffs(&ctx->lp_coeffs, avctx->sample_rate/2, cutoff) < 0){
            ctx->flags |= PSY_MODEL_NO_LOWPASS;
        }else{
            ctx->lp_state = av_mallocz(sizeof(LPFilterState) * elements * 2);
        }
    }
    if(ctx->model->init)
        return ctx->model->init(ctx, elements);
    return 0;
}

void ff_aac_psy_suggest_window(AACPsyContext *ctx, int16_t *audio, int16_t *la, int tag, int type, ChannelElement *cpe)
{
    ctx->model->window(ctx, audio, la, tag, type, cpe);
}

void ff_aac_psy_analyze(AACPsyContext *ctx, int tag, int type, ChannelElement *cpe)
{
    ctx->model->process(ctx, tag, type, cpe);
}

void av_cold ff_aac_psy_end(AACPsyContext *ctx)
{
    av_freep(&ctx->lp_state);
    if(ctx->model->end)
        return ctx->model->end(ctx);
}

void ff_aac_psy_preprocess(AACPsyContext *ctx, int16_t *audio, int16_t *dest, int tag, int type)
{
    int chans = type == ID_CPE ? 2 : 1;
    const int chstride = ctx->avctx->channels;
    int i, ch;
    float t[2];

    if(chans == 1){
        for(ch = 0; ch < chans; ch++){
            for(i = 0; i < 1024; i++){
                dest[i * chstride + ch] = audio[i * chstride + ch];
            }
        }
    }else{
        for(i = 0; i < 1024; i++){
            if(ctx->flags & PSY_MODEL_NO_ST_ATT){
                for(ch = 0; ch < 2; ch++)
                    t[ch] = audio[i * chstride + ch];
            }else{
                t[0] = audio[i * chstride + 0] * (0.5 + ctx->stereo_att) + audio[i * chstride + 1] * (0.5 - ctx->stereo_att);
                t[1] = audio[i * chstride + 0] * (0.5 - ctx->stereo_att) + audio[i * chstride + 1] * (0.5 + ctx->stereo_att);
            }
            if(!(ctx->flags & PSY_MODEL_NO_LOWPASS)){
                LPFilterState *is = (LPFilterState*)ctx->lp_state + tag*2;
                for(ch = 0; ch < 2; ch++)
                    t[ch] = ff_lowpass_filter(&ctx->lp_coeffs, is + ch, t[ch]);
            }
            for(ch = 0; ch < 2; ch++)
                dest[i * chstride + ch] = av_clip_int16(t[ch]);
        }
    }
}




More information about the ffmpeg-devel mailing list