aaccoder: move the Intensity Stereo implementation out

This commit moves the intensity stereo implementation
out from aaccoder and into a separate file. This was
possible using the previous commits.

This commit also drastically improves the IS implementation
by making it phase invariant e.g. it will always choose the
best possible phase regardless of whether M/S coding is on
or most of the coefficients have identical phases.
This also increases the quality and reduces any distortions
introduced by enablind intensity stereo.

Users are encouraged to test it out using the -aac_is 1
parameter as it has always been.

Signed-off-by: Rostislav Pehlivanov <atomnuker@gmail.com>

5 files changed, 192 insertions, 102 deletions

diff --git a/libavcodec/Makefile b/libavcodec/Makefile
index d595fe178e..3d6913a368 100644
--- a/libavcodec/Makefile
+++ b/libavcodec/Makefile
@@ -132,6 +132,7 @@ OBJS-$(CONFIG_AAC_FIXED_DECODER)       += aacdec_fixed.o aactab.o aacsbr_fixed.o
                                           sbrdsp_fixed.o aacpsdsp_fixed.o
 OBJS-$(CONFIG_AAC_ENCODER)             += aacenc.o aaccoder.o aacenctab.o    \
                                           aacpsy.o aactab.o      \
+                                          aacenc_is.o \
                                           psymodel.o mpeg4audio.o kbdwin.o
 OBJS-$(CONFIG_AASC_DECODER)            += aasc.o msrledec.o
 OBJS-$(CONFIG_AC3_DECODER)             += ac3dec_float.o ac3dec_data.o ac3.o kbdwin.o
diff --git a/libavcodec/aaccoder.c b/libavcodec/aaccoder.c
index e325788560..5df30a61ef 100644
--- a/libavcodec/aaccoder.c
+++ b/libavcodec/aaccoder.c
@@ -44,6 +44,8 @@
 #include "aacenc_quantization.h"
 #include "aac_tablegen_decl.h"
 
+#include "aacenc_is.h"
+
 /** Frequency in Hz for lower limit of noise substitution **/
 #define NOISE_LOW_LIMIT 4500
 
@@ -889,104 +891,6 @@ static void search_for_pns(AACEncContext *s, AVCodecContext *avctx, SingleChanne
     }
 }
 
-static void search_for_is(AACEncContext *s, AVCodecContext *avctx, ChannelElement *cpe)
-{
-    float IS[128];
-    float *L34  = s->scoefs + 128*0, *R34  = s->scoefs + 128*1;
-    float *I34  = s->scoefs + 128*2;
-    SingleChannelElement *sce0 = &cpe->ch[0];
-    SingleChannelElement *sce1 = &cpe->ch[1];
-    int start = 0, count = 0, i, w, w2, g;
-    const float freq_mult = avctx->sample_rate/(1024.0f/sce0->ics.num_windows)/2.0f;
-    const float lambda = s->lambda;
-
-    for (w = 0; w < 128; w++)
-        if (sce1->band_type[w] >= INTENSITY_BT2)
-            sce1->band_type[w] = 0;
-
-    if (!cpe->common_window)
-        return;
-    for (w = 0; w < sce0->ics.num_windows; w += sce0->ics.group_len[w]) {
-        start = 0;
-        for (g = 0;  g < sce0->ics.num_swb; g++) {
-            if (start*freq_mult > INT_STEREO_LOW_LIMIT*(lambda/170.0f) &&
-                cpe->ch[0].band_type[w*16+g] != NOISE_BT && !cpe->ch[0].zeroes[w*16+g] &&
-                cpe->ch[1].band_type[w*16+g] != NOISE_BT && !cpe->ch[1].zeroes[w*16+g]) {
-                int phase = 0;
-                float ener0 = 0.0f, ener1 = 0.0f, ener01 = 0.0f;
-                float dist1 = 0.0f, dist2 = 0.0f;
-                for (w2 = 0; w2 < sce0->ics.group_len[w]; w2++) {
-                    for (i = 0; i < sce0->ics.swb_sizes[g]; i++) {
-                        float coef0 = sce0->pcoeffs[start+(w+w2)*128+i];
-                        float coef1 = sce1->pcoeffs[start+(w+w2)*128+i];
-                        phase += coef0*coef1 >= 0.0f ? 1 : -1;
-                        ener0 += coef0*coef0;
-                        ener1 += coef1*coef1;
-                        ener01 += (coef0 + coef1)*(coef0 + coef1);
-                    }
-                }
-                if (!phase) { /* Too much phase difference between channels */
-                    start += sce0->ics.swb_sizes[g];
-                    continue;
-                }
-                phase = av_clip(phase, -1, 1);
-                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];
-                    int is_band_type, is_sf_idx = FFMAX(1, sce0->sf_idx[(w+w2)*16+g]-4);
-                    float e01_34 = phase*pow(sqrt(ener1/ener0), 3.0/4.0);
-                    float maxval, dist_spec_err = 0.0f;
-                    float minthr = FFMIN(band0->threshold, band1->threshold);
-                    for (i = 0; i < sce0->ics.swb_sizes[g]; i++)
-                        IS[i] = (sce0->pcoeffs[start+(w+w2)*128+i] + phase*sce1->pcoeffs[start+(w+w2)*128+i]) * sqrt(ener0/ener01);
-                    abs_pow34_v(L34, sce0->coeffs+start+(w+w2)*128, sce0->ics.swb_sizes[g]);
-                    abs_pow34_v(R34, sce1->coeffs+start+(w+w2)*128, sce0->ics.swb_sizes[g]);
-                    abs_pow34_v(I34, IS,                            sce0->ics.swb_sizes[g]);
-                    maxval = find_max_val(1, sce0->ics.swb_sizes[g], I34);
-                    is_band_type = find_min_book(maxval, is_sf_idx);
-                    dist1 += quantize_band_cost(s, sce0->coeffs + start + (w+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, 0);
-                    dist1 += quantize_band_cost(s, sce1->coeffs + start + (w+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, 0);
-                    dist2 += quantize_band_cost(s, IS,
-                                                I34,
-                                                sce0->ics.swb_sizes[g],
-                                                is_sf_idx,
-                                                is_band_type,
-                                                lambda / minthr, INFINITY, NULL, 0);
-                    for (i = 0; i < sce0->ics.swb_sizes[g]; i++) {
-                        dist_spec_err += (L34[i] - I34[i])*(L34[i] - I34[i]);
-                        dist_spec_err += (R34[i] - I34[i]*e01_34)*(R34[i] - I34[i]*e01_34);
-                    }
-                    dist_spec_err *= lambda / minthr;
-                    dist2 += dist_spec_err;
-                }
-                if (dist2 <= dist1) {
-                    cpe->is_mask[w*16+g] = 1;
-                    cpe->ms_mask[w*16+g] = 0;
-                    cpe->ch[0].is_ener[w*16+g] = sqrt(ener0/ener01);
-                    cpe->ch[1].is_ener[w*16+g] = ener0/ener1;
-                    if (phase)
-                        cpe->ch[1].band_type[w*16+g] = INTENSITY_BT;
-                    else
-                        cpe->ch[1].band_type[w*16+g] = INTENSITY_BT2;
-                    count++;
-                }
-            }
-            start += sce0->ics.swb_sizes[g];
-        }
-    }
-    cpe->is_mode = !!count;
-}
-
 static void search_for_ms(AACEncContext *s, ChannelElement *cpe)
 {
     int start = 0, i, w, w2, g;
@@ -1000,7 +904,7 @@ static void search_for_ms(AACEncContext *s, ChannelElement *cpe)
     for (w = 0; w < sce0->ics.num_windows; w += sce0->ics.group_len[w]) {
         start = 0;
         for (g = 0;  g < sce0->ics.num_swb; g++) {
-            if (!cpe->ch[0].zeroes[w*16+g] && !cpe->ch[1].zeroes[w*16+g] && !cpe->is_mask[w*16+g]) {
+            if (!cpe->ch[0].zeroes[w*16+g] && !cpe->ch[1].zeroes[w*16+g]) {
                 float dist1 = 0.0f, dist2 = 0.0f;
                 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];
diff --git a/libavcodec/aacenc.c b/libavcodec/aacenc.c
index 924064823f..4ade340bf7 100644
--- a/libavcodec/aacenc.c
+++ b/libavcodec/aacenc.c
@@ -459,7 +459,9 @@ static int aac_encode_frame(AVCodecContext *avctx, AVPacket *avpkt,
     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, ms_mode = 0, is_mode = 0;
+    SingleChannelElement *sce;
+    int i, ch, w, g, chans, tag, start_ch, ret;
+    int ms_mode = 0, is_mode = 0, tns_mode = 0, pred_mode = 0;
     int chan_el_counter[4];
     FFPsyWindowInfo windows[AAC_MAX_CHANNELS];
 
@@ -608,7 +610,7 @@ static int aac_encode_frame(AVCodecContext *avctx, AVPacket *avpkt,
                     s->coder->search_for_ms(s, cpe);
                 }
             }
-            if (chans > 1 && s->options.intensity_stereo && s->coder->search_for_is) {
+            if (s->options.intensity_stereo && s->coder->search_for_is) {
                 s->coder->search_for_is(s, avctx, cpe);
                 if (cpe->is_mode) is_mode = 1;
             }
@@ -636,7 +638,7 @@ static int aac_encode_frame(AVCodecContext *avctx, AVPacket *avpkt,
             s->psy.bitres.bits = frame_bits / s->channels;
             break;
         }
-        if (is_mode || ms_mode) {
+        if (is_mode || ms_mode || tns_mode || pred_mode) {
             for (i = 0; i < s->chan_map[0]; i++) {
                 // Must restore coeffs
                 chans = tag == TYPE_CPE ? 2 : 1;
diff --git a/libavcodec/aacenc_is.c b/libavcodec/aacenc_is.c
new file mode 100644
index 0000000000..4ae10542c4
--- /dev/null
+++ b/libavcodec/aacenc_is.c
@@ -0,0 +1,136 @@
+/*
+ * AAC encoder intensity stereo
+ * Copyright (C) 2015 Rostislav Pehlivanov
+ *
+ * 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
+ * AAC encoder Intensity Stereo
+ * @author Rostislav Pehlivanov ( atomnuker gmail com )
+ */
+
+#include "aacenc.h"
+#include "aacenc_utils.h"
+#include "aacenc_is.h"
+#include "aacenc_quantization.h"
+
+struct is_error calc_is_encoding_err(AACEncContext *s, ChannelElement *cpe,
+                                     int start, int w, int g, float ener0,
+                                     float ener1, float ener01, int phase)
+{
+    int i, w2;
+    float *L34 = &s->scoefs[256*0], *R34 = &s->scoefs[256*1];
+    float *IS  = &s->scoefs[256*2], *I34 = &s->scoefs[256*3];
+    float dist1 = 0.0f, dist2 = 0.0f;
+    struct is_error is_error = {0};
+    SingleChannelElement *sce0 = &cpe->ch[0];
+    SingleChannelElement *sce1 = &cpe->ch[1];
+
+    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];
+        int is_band_type, is_sf_idx = FFMAX(1, sce0->sf_idx[(w+w2)*16+g]-4);
+        float e01_34 = phase*pow(sqrt(ener1/ener0), 3.0/4.0);
+        float maxval, dist_spec_err = 0.0f;
+        float minthr = FFMIN(band0->threshold, band1->threshold);
+        for (i = 0; i < sce0->ics.swb_sizes[g]; i++) {
+            IS[i] = (sce0->pcoeffs[start+(w+w2)*128+i]+
+                     phase*sce1->pcoeffs[start+(w+w2)*128+i])*
+                     sqrt(ener0/ener01);
+        }
+        abs_pow34_v(L34, &sce0->coeffs[start+(w+w2)*128], sce0->ics.swb_sizes[g]);
+        abs_pow34_v(R34, &sce1->coeffs[start+(w+w2)*128], sce0->ics.swb_sizes[g]);
+        abs_pow34_v(I34, IS,                            sce0->ics.swb_sizes[g]);
+        maxval = find_max_val(1, sce0->ics.swb_sizes[g], I34);
+        is_band_type = find_min_book(maxval, is_sf_idx);
+        dist1 += quantize_band_cost(s, &sce0->coeffs[start + (w+w2)*128], L34,
+                                    sce0->ics.swb_sizes[g],
+                                    sce0->sf_idx[(w+w2)*16+g],
+                                    sce0->band_type[(w+w2)*16+g],
+                                    s->lambda / band0->threshold, INFINITY, NULL, 0);
+        dist1 += quantize_band_cost(s, &sce1->coeffs[start + (w+w2)*128], R34,
+                                    sce1->ics.swb_sizes[g],
+                                    sce1->sf_idx[(w+w2)*16+g],
+                                    sce1->band_type[(w+w2)*16+g],
+                                    s->lambda / band1->threshold, INFINITY, NULL, 0);
+        dist2 += quantize_band_cost(s, IS, I34, sce0->ics.swb_sizes[g],
+                                    is_sf_idx, is_band_type,
+                                    s->lambda / minthr, INFINITY, NULL, 0);
+        for (i = 0; i < sce0->ics.swb_sizes[g]; i++) {
+            dist_spec_err += (L34[i] - I34[i])*(L34[i] - I34[i]);
+            dist_spec_err += (R34[i] - I34[i]*e01_34)*(R34[i] - I34[i]*e01_34);
+        }
+        dist_spec_err *= s->lambda / minthr;
+        dist2 += dist_spec_err;
+    }
+
+    is_error.pass = dist2 <= dist1;
+    is_error.phase = phase;
+    is_error.error = fabsf(dist1 - dist2);
+    is_error.dist1 = dist1;
+    is_error.dist2 = dist2;
+
+    return is_error;
+}
+
+void search_for_is(AACEncContext *s, AVCodecContext *avctx, ChannelElement *cpe)
+{
+    SingleChannelElement *sce0 = &cpe->ch[0];
+    SingleChannelElement *sce1 = &cpe->ch[1];
+    int start = 0, count = 0, w, w2, g, i;
+    const float freq_mult = avctx->sample_rate/(1024.0f/sce0->ics.num_windows)/2.0f;
+
+    if (!cpe->common_window)
+        return;
+
+    for (w = 0; w < sce0->ics.num_windows; w += sce0->ics.group_len[w]) {
+        start = 0;
+        for (g = 0;  g < sce0->ics.num_swb; g++) {
+            if (start*freq_mult > INT_STEREO_LOW_LIMIT*(s->lambda/170.0f) &&
+                cpe->ch[0].band_type[w*16+g] != NOISE_BT && !cpe->ch[0].zeroes[w*16+g] &&
+                cpe->ch[1].band_type[w*16+g] != NOISE_BT && !cpe->ch[1].zeroes[w*16+g]) {
+                float ener0 = 0.0f, ener1 = 0.0f, ener01 = 0.0f;
+                struct is_error ph_err1, ph_err2, *erf;
+                for (w2 = 0; w2 < sce0->ics.group_len[w]; w2++) {
+                    for (i = 0; i < sce0->ics.swb_sizes[g]; i++) {
+                        float coef0 = sce0->pcoeffs[start+(w+w2)*128+i];
+                        float coef1 = sce1->pcoeffs[start+(w+w2)*128+i];
+                        ener0 += coef0*coef0;
+                        ener1 += coef1*coef1;
+                        ener01 += (coef0 + coef1)*(coef0 + coef1);
+                    }
+                }
+                ph_err1 = calc_is_encoding_err(s, cpe, start, w, g,
+                                               ener0, ener1, ener01, -1);
+                ph_err2 = calc_is_encoding_err(s, cpe, start, w, g,
+                                               ener0, ener1, ener01, +1);
+                erf = ph_err1.error < ph_err2.error ? &ph_err1 : &ph_err2;
+                if (erf->pass) {
+                    cpe->is_mask[w*16+g] = 1;
+                    cpe->ch[0].is_ener[w*16+g] = sqrt(ener0/ener01);
+                    cpe->ch[1].is_ener[w*16+g] = ener0/ener1;
+                    cpe->ch[1].band_type[w*16+g] = erf->phase ? INTENSITY_BT : INTENSITY_BT2;
+                    count++;
+                }
+            }
+            start += sce0->ics.swb_sizes[g];
+        }
+    }
+    cpe->is_mode = !!count;
+}
diff --git a/libavcodec/aacenc_is.h b/libavcodec/aacenc_is.h
new file mode 100644
index 0000000000..abb9b32daf
--- /dev/null
+++ b/libavcodec/aacenc_is.h
@@ -0,0 +1,47 @@
+/*
+ * AAC encoder intensity stereo
+ * Copyright (C) 2015 Rostislav Pehlivanov
+ *
+ * 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
+ * AAC encoder Intensity Stereo
+ * @author Rostislav Pehlivanov ( atomnuker gmail com )
+ */
+
+#ifndef AVCODEC_AACENC_INTENSITY_STEREO_H
+#define AVCODEC_AACENC_INTENSITY_STEREO_H
+
+/** Frequency in Hz for lower limit of intensity stereo **/
+#define INT_STEREO_LOW_LIMIT 6100
+
+struct is_error {
+    int pass;    /* 1 if dist2 <= dist1  */
+    int phase;   /* -1 or +1             */
+    float error; /* fabs(dist1 - dist2)  */
+    float dist1; /* From original coeffs */
+    float dist2; /* From IS'd coeffs     */
+};
+
+struct is_error calc_is_encoding_err(AACEncContext *s, ChannelElement *cpe,
+                                     int start, int g, int w, float ener0,
+                                     float ener1, float ener01, int phase);
+void search_for_is(AACEncContext *s, AVCodecContext *avctx, ChannelElement *cpe);
+
+#endif /* AVCODEC_AACENC_INTENSITY_STEREO_H */