avcodec/dcaenc: Initial implementation of ADPCM encoding for DCA encoder

libavcodec/Makefile

@@ -244,7 +244,8 @@ OBJS-$(CONFIG_CYUV_DECODER)            += cyuv.o
 OBJS-$(CONFIG_DCA_DECODER)             += dcadec.o dca.o dcadata.o dcahuff.o \
                                           dca_core.o dca_exss.o dca_xll.o dca_lbr.o \
                                           dcadsp.o dcadct.o synth_filter.o
-OBJS-$(CONFIG_DCA_ENCODER)             += dcaenc.o dca.o dcadata.o dcahuff.o
+OBJS-$(CONFIG_DCA_ENCODER)             += dcaenc.o dca.o dcadata.o dcahuff.o \
+                                          dcaadpcm.o
 OBJS-$(CONFIG_DDS_DECODER)             += dds.o
 OBJS-$(CONFIG_DIRAC_DECODER)           += diracdec.o dirac.o diracdsp.o diractab.o \
                                           dirac_arith.o dirac_dwt.o dirac_vlc.o

libavcodec/dca_core.c

@@ -18,6 +18,7 @@
  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
  */
 
+#include "dcaadpcm.h"
 #include "dcadec.h"
 #include "dcadata.h"
 #include "dcahuff.h"
@@ -670,46 +671,21 @@ static inline int extract_audio(DCACoreDecoder *s, int32_t *audio, int abits, in
     return 0;
 }
 
-static inline void dequantize(int32_t *output, const int32_t *input,
-                              int32_t step_size, int32_t scale, int residual)
-{
-    // Account for quantizer step size
-    int64_t step_scale = (int64_t)step_size * scale;
-    int n, shift = 0;
-
-    // Limit scale factor resolution to 22 bits
-    if (step_scale > (1 << 23)) {
-        shift = av_log2(step_scale >> 23) + 1;
-        step_scale >>= shift;
-    }
-
-    // Scale the samples
-    if (residual) {
-        for (n = 0; n < DCA_SUBBAND_SAMPLES; n++)
-            output[n] += clip23(norm__(input[n] * step_scale, 22 - shift));
-    } else {
-        for (n = 0; n < DCA_SUBBAND_SAMPLES; n++)
-            output[n]  = clip23(norm__(input[n] * step_scale, 22 - shift));
-    }
-}
-
 static inline void inverse_adpcm(int32_t **subband_samples,
                                  const int16_t *vq_index,
                                  const int8_t *prediction_mode,
                                  int sb_start, int sb_end,
                                  int ofs, int len)
 {
-    int i, j, k;
+    int i, j;
 
     for (i = sb_start; i < sb_end; i++) {
         if (prediction_mode[i]) {
-            const int16_t *coeff = ff_dca_adpcm_vb[vq_index[i]];
+            const int pred_id = vq_index[i];
             int32_t *ptr = subband_samples[i] + ofs;
             for (j = 0; j < len; j++) {
-                int64_t err = 0;
-                for (k = 0; k < DCA_ADPCM_COEFFS; k++)
-                    err += (int64_t)ptr[j - k - 1] * coeff[k];
-                ptr[j] = clip23(ptr[j] + clip23(norm13(err)));
+                int32_t x = ff_dcaadpcm_predict(pred_id, ptr + j - DCA_ADPCM_COEFFS);
+                ptr[j] = clip23(ptr[j] + x);
             }
         }
     }
@@ -817,8 +793,8 @@ static int parse_subframe_audio(DCACoreDecoder *s, int sf, enum HeaderType heade
                     scale = clip23(adj * scale >> 22);
                 }
 
-                dequantize(s->subband_samples[ch][band] + ofs,
-                           audio, step_size, scale, 0);
+                ff_dca_core_dequantize(s->subband_samples[ch][band] + ofs,
+                           audio, step_size, scale, 0, DCA_SUBBAND_SAMPLES);
             }
         }
 
@@ -1146,8 +1122,8 @@ static int parse_xbr_subframe(DCACoreDecoder *s, int xbr_base_ch, int xbr_nchann
                 else
                     scale = xbr_scale_factors[ch][band][1];
 
-                dequantize(s->subband_samples[ch][band] + ofs,
-                           audio, step_size, scale, 1);
+                ff_dca_core_dequantize(s->subband_samples[ch][band] + ofs,
+                           audio, step_size, scale, 1, DCA_SUBBAND_SAMPLES);
             }
         }
 
@@ -1326,8 +1302,8 @@ static int parse_x96_subframe_audio(DCACoreDecoder *s, int sf, int xch_base, int
                 // Get the scale factor
                 scale = s->scale_factors[ch][band >> 1][band & 1];
 
-                dequantize(s->x96_subband_samples[ch][band] + ofs,
-                           audio, step_size, scale, 0);
+                ff_dca_core_dequantize(s->x96_subband_samples[ch][band] + ofs,
+                           audio, step_size, scale, 0, DCA_SUBBAND_SAMPLES);
             }
         }
 

libavcodec/dca_core.h

@@ -33,6 +33,7 @@
 #include "dca_exss.h"
 #include "dcadsp.h"
 #include "dcadct.h"
+#include "dcamath.h"
 #include "dcahuff.h"
 #include "fft.h"
 #include "synth_filter.h"
@@ -43,7 +44,6 @@
 #define DCA_SUBFRAMES           16
 #define DCA_SUBBAND_SAMPLES     8
 #define DCA_PCMBLOCK_SAMPLES    32
-#define DCA_ADPCM_COEFFS        4
 #define DCA_LFE_HISTORY         8
 #define DCA_ABITS_MAX           26
 
@@ -195,6 +195,29 @@ static inline int ff_dca_core_map_spkr(DCACoreDecoder *core, int spkr)
     return -1;
 }
 
+static inline void ff_dca_core_dequantize(int32_t *output, const int32_t *input,
+                                          int32_t step_size, int32_t scale, int residual, int len)
+{
+    // Account for quantizer step size
+    int64_t step_scale = (int64_t)step_size * scale;
+    int n, shift = 0;
+
+    // Limit scale factor resolution to 22 bits
+    if (step_scale > (1 << 23)) {
+        shift = av_log2(step_scale >> 23) + 1;
+        step_scale >>= shift;
+    }
+
+    // Scale the samples
+    if (residual) {
+        for (n = 0; n < len; n++)
+            output[n] += clip23(norm__(input[n] * step_scale, 22 - shift));
+    } else {
+        for (n = 0; n < len; n++)
+            output[n]  = clip23(norm__(input[n] * step_scale, 22 - shift));
+    }
+}
+
 int ff_dca_core_parse(DCACoreDecoder *s, uint8_t *data, int size);
 int ff_dca_core_parse_exss(DCACoreDecoder *s, uint8_t *data, DCAExssAsset *asset);
 int ff_dca_core_filter_fixed(DCACoreDecoder *s, int x96_synth);

libavcodec/dcaadpcm.c

+/*
+ * DCA ADPCM engine
+ * Copyright (C) 2017 Daniil Cherednik
+ *
+ * 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
+ */
+
+
+#include "dcaadpcm.h"
+#include "dcaenc.h"
+#include "dca_core.h"
+#include "mathops.h"
+
+typedef int32_t premultiplied_coeffs[10];
+
+//assume we have DCA_ADPCM_COEFFS values before x
+static inline int64_t calc_corr(const int32_t *x, int len, int j, int k)
+{
+    int n;
+    int64_t s = 0;
+    for (n = 0; n < len; n++)
+        s += MUL64(x[n-j], x[n-k]);
+    return s;
+}
+
+static inline int64_t apply_filter(const int16_t a[DCA_ADPCM_COEFFS], const int64_t corr[15], const int32_t aa[10])
+{
+    int64_t err = 0;
+    int64_t tmp = 0;
+
+    err = corr[0];
+
+    tmp += MUL64(a[0], corr[1]);
+    tmp += MUL64(a[1], corr[2]);
+    tmp += MUL64(a[2], corr[3]);
+    tmp += MUL64(a[3], corr[4]);
+
+    tmp = norm__(tmp, 13);
+    tmp += tmp;
+
+    err -= tmp;
+    tmp = 0;
+
+    tmp += MUL64(corr[5], aa[0]);
+    tmp += MUL64(corr[6], aa[1]);
+    tmp += MUL64(corr[7], aa[2]);
+    tmp += MUL64(corr[8], aa[3]);
+
+    tmp += MUL64(corr[9], aa[4]);
+    tmp += MUL64(corr[10], aa[5]);
+    tmp += MUL64(corr[11], aa[6]);
+
+    tmp += MUL64(corr[12], aa[7]);
+    tmp += MUL64(corr[13], aa[8]);
+
+    tmp += MUL64(corr[14], aa[9]);
+
+    tmp = norm__(tmp, 26);
+
+    err += tmp;
+
+    return llabs(err);
+}
+
+static int64_t find_best_filter(const DCAADPCMEncContext *s, const int32_t *in, int len)
+{
+    const premultiplied_coeffs *precalc_data = s->private_data;
+    int i, j, k = 0;
+    int vq;
+    int64_t err;
+    int64_t min_err = 1ll << 62;
+    int64_t corr[15];
+
+    for (i = 0; i <= DCA_ADPCM_COEFFS; i++)
+        for (j = i; j <= DCA_ADPCM_COEFFS; j++)
+            corr[k++] = calc_corr(in+4, len, i, j);
+
+    for (i = 0; i < DCA_ADPCM_VQCODEBOOK_SZ; i++) {
+        err = apply_filter(ff_dca_adpcm_vb[i], corr, *precalc_data);
+        if (err < min_err) {
+            min_err = err;
+            vq = i;
+        }
+        precalc_data++;
+    }
+
+    return vq;
+}
+
+static inline int64_t calc_prediction_gain(int pred_vq, const int32_t *in, int32_t *out, int len)
+{
+    int i;
+    int32_t error;
+
+    int64_t signal_energy = 0;
+    int64_t error_energy = 0;
+
+    for (i = 0; i < len; i++) {
+        error = in[DCA_ADPCM_COEFFS + i] - ff_dcaadpcm_predict(pred_vq, in + i);
+        out[i] = error;
+        signal_energy += MUL64(in[DCA_ADPCM_COEFFS + i], in[DCA_ADPCM_COEFFS + i]);
+        error_energy += MUL64(error, error);
+    }
+
+    if (!error_energy)
+        return -1;
+
+    return signal_energy / error_energy;
+}
+
+int ff_dcaadpcm_subband_analysis(const DCAADPCMEncContext *s, const int32_t *in, int len, int *diff)
+{
+    int pred_vq, i;
+    int32_t input_buffer[16 + DCA_ADPCM_COEFFS];
+    int32_t input_buffer2[16 + DCA_ADPCM_COEFFS];
+
+    int32_t max = 0;
+    int shift_bits;
+    uint64_t pg = 0;
+
+    for (i = 0; i < len + DCA_ADPCM_COEFFS; i++)
+        max |= FFABS(in[i]);
+
+    // normalize input to simplify apply_filter
+    shift_bits = av_log2(max) - 11;
+
+    for (i = 0; i < len + DCA_ADPCM_COEFFS; i++) {
+        input_buffer[i] = norm__(in[i], 7);
+        input_buffer2[i] = norm__(in[i], shift_bits);
+    }
+
+    pred_vq = find_best_filter(s, input_buffer2, len);
+
+    if (pred_vq < 0)
+        return -1;
+
+    pg = calc_prediction_gain(pred_vq, input_buffer, diff, len);
+
+    // Greater than 10db (10*log(10)) prediction gain to use ADPCM.
+    // TODO: Tune it.
+    if (pg < 10)
+        return -1;
+
+    for (i = 0; i < len; i++)
+        diff[i] <<= 7;
+
+    return pred_vq;
+}
+
+static void precalc(premultiplied_coeffs *data)
+{
+    int i, j, k;
+
+    for (i = 0; i < DCA_ADPCM_VQCODEBOOK_SZ; i++) {
+        int id = 0;
+        int32_t t = 0;
+        for (j = 0; j < DCA_ADPCM_COEFFS; j++) {
+            for (k = j; k < DCA_ADPCM_COEFFS; k++) {
+                t = (int32_t)ff_dca_adpcm_vb[i][j] * (int32_t)ff_dca_adpcm_vb[i][k];
+                if (j != k)
+                    t *= 2;
+                (*data)[id++] = t;
+             }
+        }
+        data++;
+    }
+}
+
+int ff_dcaadpcm_do_real(int pred_vq_index,
+                        softfloat quant, int32_t scale_factor, int32_t step_size,
+                        const int32_t *prev_hist, const int32_t *in, int32_t *next_hist, int32_t *out,
+                        int len, int32_t peak)
+{
+    int i;
+    int64_t delta;
+    int32_t dequant_delta;
+    int32_t work_bufer[16 + DCA_ADPCM_COEFFS];
+
+    memcpy(work_bufer, prev_hist, sizeof(int32_t) * DCA_ADPCM_COEFFS);
+
+    for (i = 0; i < len; i++) {
+        work_bufer[DCA_ADPCM_COEFFS + i] = ff_dcaadpcm_predict(pred_vq_index, &work_bufer[i]);
+
+        delta = (int64_t)in[i] - ((int64_t)work_bufer[DCA_ADPCM_COEFFS + i] << 7);
+
+        out[i] = quantize_value(av_clip64(delta, -peak, peak), quant);
+
+        ff_dca_core_dequantize(&dequant_delta, &out[i], step_size, scale_factor, 0, 1);
+
+        work_bufer[DCA_ADPCM_COEFFS+i] += dequant_delta;
+    }
+
+    memcpy(next_hist, &work_bufer[len], sizeof(int32_t) * DCA_ADPCM_COEFFS);
+
+    return 0;
+}
+
+av_cold int ff_dcaadpcm_init(DCAADPCMEncContext *s)
+{
+    if (!s)
+        return -1;
+
+    s->private_data = av_malloc(sizeof(premultiplied_coeffs) * DCA_ADPCM_VQCODEBOOK_SZ);
+    precalc(s->private_data);
+    return 0;
+}
+
+av_cold void ff_dcaadpcm_free(DCAADPCMEncContext *s)
+{
+    if (!s)
+        return;
+
+    av_freep(&s->private_data);
+}

libavcodec/dcaadpcm.h

+/*
+ * DCA ADPCM engine
+ * Copyright (C) 2017 Daniil Cherednik
+ *
+ * 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 AVCODEC_DCAADPCM_H
+#define AVCODEC_DCAADPCM_H
+
+#include "dcamath.h"
+#include "dcadata.h"
+#include "dcaenc.h"
+
+typedef struct DCAADPCMEncContext {
+    void *private_data;
+} DCAADPCMEncContext;
+
+static inline int64_t ff_dcaadpcm_predict(int pred_vq_index, const int32_t *input)
+{
+    int i;
+    const int16_t *coeff = ff_dca_adpcm_vb[pred_vq_index];
+    int64_t pred = 0;
+    for (i = 0; i < DCA_ADPCM_COEFFS; i++)
+        pred += (int64_t)input[DCA_ADPCM_COEFFS - 1 - i] * coeff[i];
+
+    return clip23(norm13(pred));
+}
+
+int ff_dcaadpcm_subband_analysis(const DCAADPCMEncContext *s, const int32_t *input, int len, int *diff);
+
+int ff_dcaadpcm_do_real(int pred_vq_index,
+                        softfloat quant, int32_t scale_factor, int32_t step_size,
+                        const int32_t *prev_hist, const int32_t *in, int32_t *next_hist, int32_t *out,
+                        int len, int32_t peak);
+
+av_cold int ff_dcaadpcm_init(DCAADPCMEncContext *s);
+av_cold void ff_dcaadpcm_free(DCAADPCMEncContext *s);
+
+#endif /* AVCODEC_DCAADPCM_H */

libavcodec/dcadata.c

@@ -61,7 +61,7 @@ const uint8_t ff_dca_quant_index_group_size[DCA_CODE_BOOKS] = {
 /* ADPCM data */
 
 /* 16 bits signed fractional Q13 binary codes */
-const int16_t ff_dca_adpcm_vb[4096][4] = {
+const int16_t ff_dca_adpcm_vb[DCA_ADPCM_VQCODEBOOK_SZ][DCA_ADPCM_COEFFS] = {
     {   9928,  -2618,  -1093, -1263 },
     {  11077,  -2876,  -1747,  -308 },
     {  10503,  -1082,  -1426, -1167 },

libavcodec/dcadata.h

@@ -25,6 +25,9 @@
 
 #include "dcahuff.h"
 
+#define DCA_ADPCM_COEFFS        4
+#define DCA_ADPCM_VQCODEBOOK_SZ 4096
+
 extern const uint32_t ff_dca_bit_rates[32];
 
 extern const uint8_t ff_dca_channels[16];
@@ -36,7 +39,7 @@ extern const uint8_t ff_dca_dmix_primary_nch[8];
 extern const uint8_t ff_dca_quant_index_sel_nbits[DCA_CODE_BOOKS];
 extern const uint8_t ff_dca_quant_index_group_size[DCA_CODE_BOOKS];
 
-extern const int16_t ff_dca_adpcm_vb[4096][4];
+extern const int16_t ff_dca_adpcm_vb[DCA_ADPCM_VQCODEBOOK_SZ][DCA_ADPCM_COEFFS];
 
 extern const uint32_t ff_dca_scale_factor_quant6[64];
 extern const uint32_t ff_dca_scale_factor_quant7[128];

libavcodec/dcaenc.h

@@ -24,6 +24,8 @@
 
 #include <stdint.h>
 
+#include "dcamath.h"
+
 typedef struct {
     int32_t m;
     int32_t e;
@@ -144,4 +146,13 @@ static const int8_t channel_reorder_nolfe[16][9] = {
     { 3,  2,  4,  0,  1,  5,  7,  6, -1 },
 };
 
+static inline int32_t quantize_value(int32_t value, softfloat quant)
+{
+    int32_t offset = 1 << (quant.e - 1);
+
+    value = mul32(value, quant.m) + offset;
+    value = value >> quant.e;
+    return value;
+}
+
 #endif /* AVCODEC_DCAENC_H */

libavcodec/dcamath.h

@@ -49,6 +49,7 @@ static inline int32_t mul17(int32_t a, int32_t b) { return mul__(a, b, 17); }
 static inline int32_t mul22(int32_t a, int32_t b) { return mul__(a, b, 22); }
 static inline int32_t mul23(int32_t a, int32_t b) { return mul__(a, b, 23); }
 static inline int32_t mul31(int32_t a, int32_t b) { return mul__(a, b, 31); }
+static inline int32_t mul32(int32_t a, int32_t b) { return mul__(a, b, 32); }
 
 static inline int32_t clip23(int32_t a) { return av_clip_intp2(a, 23); }