swresample: support AV_SAMPLE_FMT_DSD input

Convert DSD to PCM with the same 96-tap symmetric lowpass filter the
libavcodec DSD decoders use. The implementation is copied from
libavcodec/dsd.c, minus the LSBF table variants, since AV_SAMPLE_FMT_DSD
is defined as MSBF.

The libavcodec copy will be removed in future commits as decoders are
transitioned to SWR.

The conversion is stateful (per-channel filter history), so the audio
conversion functions gain a per-channel state argument.

DSD silence is the 0x69 bit pattern instead of the 0x80 used for PCM.

Only a direct DSD to float conversion function is provided. For format
pairs without a direct conversion function (e.g. DSD to s16 with equal
rates), swr_init() now falls back from the full_convert shortcut to
the generic two-stage path instead of failing, and reports unsupported
conversions with a proper error message instead of ENOMEM.

Signed-off-by: Kacper Michajłow <kasper93@gmail.com>

libswresample/Makefile

@@ -7,6 +7,7 @@ HEADERS = swresample.h                       \
           version_major.h                    \
 
 OBJS = audioconvert.o                        \
+       dsd2pcm.o                             \
        dither.o                              \
        options.o                             \
        rematrix.o                            \
@@ -19,7 +20,7 @@ OBJS = audioconvert.o                        \
 OBJS-$(CONFIG_LIBSOXR) += soxr_resample.o
 
 # Objects duplicated from other libraries for shared builds
-SHLIBOBJS              += log2_tab.o
+SHLIBOBJS              += log2_tab.o reverse.o
 
 # Windows resource file
 SHLIBOBJS-$(HAVE_GNU_WINDRES) += swresampleres.o

libswresample/audioconvert.c

@@ -36,7 +36,7 @@
 
 //FIXME rounding ?
 #define CONV_FUNC(ofmt, otype, ifmt, expr)\
-static void CONV_FUNC_NAME(ofmt, ifmt)(uint8_t *po, const uint8_t *pi, int is, int os, uint8_t *end)\
+static void CONV_FUNC_NAME(ofmt, ifmt)(DSDContext *st, uint8_t *po, const uint8_t *pi, int is, int os, uint8_t *end)\
 {\
     uint8_t *end2 = end - 3*os;\
     while(po < end2){\
@@ -87,6 +87,12 @@ CONV_FUNC(AV_SAMPLE_FMT_S32, int32_t, AV_SAMPLE_FMT_DBL, av_clipl_int32(llrint(*
 CONV_FUNC(AV_SAMPLE_FMT_S64, int64_t, AV_SAMPLE_FMT_DBL, llrint(*(const double*)pi * (UINT64_C(1)<<63)))
 CONV_FUNC(AV_SAMPLE_FMT_FLT, float  , AV_SAMPLE_FMT_DBL, *(const double*)pi)
 CONV_FUNC(AV_SAMPLE_FMT_DBL, double , AV_SAMPLE_FMT_DBL, *(const double*)pi)
+CONV_FUNC(AV_SAMPLE_FMT_DSD, uint8_t, AV_SAMPLE_FMT_DSD, *(const uint8_t*)pi)
+
+static void CONV_FUNC_NAME(AV_SAMPLE_FMT_FLT, AV_SAMPLE_FMT_DSD)(DSDContext *st, uint8_t *po, const uint8_t *pi, int is, int os, uint8_t *end)
+{
+    swri_dsd2pcm_translate(st, (end - po) / os, pi, is, (float *)po, os / sizeof(float));
+}
 
 #define FMT_PAIR_FUNC(out, in) [(out) + AV_SAMPLE_FMT_NB*(in)] = CONV_FUNC_NAME(out, in)
 
@@ -127,6 +133,8 @@ static conv_func_type * const fmt_pair_to_conv_functions[AV_SAMPLE_FMT_NB*AV_SAM
     FMT_PAIR_FUNC(AV_SAMPLE_FMT_FLT, AV_SAMPLE_FMT_S64),
     FMT_PAIR_FUNC(AV_SAMPLE_FMT_DBL, AV_SAMPLE_FMT_S64),
     FMT_PAIR_FUNC(AV_SAMPLE_FMT_S64, AV_SAMPLE_FMT_S64),
+    FMT_PAIR_FUNC(AV_SAMPLE_FMT_FLT, AV_SAMPLE_FMT_DSD),
+    FMT_PAIR_FUNC(AV_SAMPLE_FMT_DSD, AV_SAMPLE_FMT_DSD),
 };
 
 static void cpy1(uint8_t **dst, const uint8_t **src, int len){
@@ -166,6 +174,12 @@ AudioConvert *swri_audio_convert_alloc(enum AVSampleFormat out_fmt,
     ctx->ch_map   = ch_map;
     if (in_fmt == AV_SAMPLE_FMT_U8 || in_fmt == AV_SAMPLE_FMT_U8P)
         memset(ctx->silence, 0x80, sizeof(ctx->silence));
+    if (in_fmt == AV_SAMPLE_FMT_DSD) {
+        swri_dsd2pcm_init();
+        memset(ctx->silence, 0x69, sizeof(ctx->silence));
+        for (int ch = 0; ch < FF_ARRAY_ELEMS(ctx->dsd_state); ch++)
+            memset(ctx->dsd_state[ch].buf, 0x69, sizeof(ctx->dsd_state[ch].buf));
+    }
 
     if(out_fmt == in_fmt && !ch_map) {
         switch(av_get_bytes_per_sample(in_fmt)){
@@ -245,7 +259,7 @@ int swri_audio_convert(AudioConvert *ctx, AudioData *out, AudioData *in, int len
         if(!po)
             continue;
         end = po + os * len;
-        ctx->conv_f(po+off*os, pi+off*is, is, os, end);
+        ctx->conv_f(&ctx->dsd_state[ch], po+off*os, pi+off*is, is, os, end);
     }
     return 0;
 }

libswresample/audioconvert.h

@@ -30,9 +30,10 @@
 
 
 #include "swresample_internal.h"
+#include "dsd2pcm.h"
 
 
-typedef void (conv_func_type)(uint8_t *po, const uint8_t *pi, int is, int os, uint8_t *end);
+typedef void (conv_func_type)(DSDContext *st, uint8_t *po, const uint8_t *pi, int is, int os, uint8_t *end);
 typedef void (simd_func_type)(uint8_t **dst, const uint8_t **src, int len);
 
 typedef struct AudioConvert {
@@ -43,6 +44,7 @@ typedef struct AudioConvert {
     simd_func_type *simd_f;
     const int *ch_map;
     uint8_t silence[8]; ///< silence input sample
+    DSDContext dsd_state[SWR_CH_MAX]; ///< per-channel state for DSD input
 }AudioConvert;
 
 /**

libswresample/dsd2pcm.c

@@ -0,0 +1,130 @@
+/*
+ * DSD to PCM conversion
+ * based on BSD licensed dsd2pcm by Sebastian Gesemann
+ * Copyright (c) 2009, 2011 Sebastian Gesemann. All rights reserved.
+ * Copyright (c) 2014 Peter Ross
+ *
+ * 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 <string.h>
+#include "libavutil/attributes.h"
+#include "libavutil/reverse.h"
+#include "libavutil/thread.h"
+#include "dsd2pcm.h"
+
+#define DSD_CTABLES ((DSD_HTAPS + 7) / 8) /** number of "8 MACs" lookup tables */
+
+/*
+ * Properties of this 96-tap lowpass filter when applied on a signal
+ * with sampling rate of 44100*64 Hz:
+ *
+ * () has a delay of 17 microseconds.
+ *
+ * () flat response up to 48 kHz
+ *
+ * () if you downsample afterwards by a factor of 8, the
+ *    spectrum below 70 kHz is practically alias-free.
+ *
+ * () stopband rejection is about 160 dB
+ *
+ * The coefficient tables ("ctables") take only 6 Kibi Bytes and
+ * should fit into a modern processor's fast cache.
+ */
+
+/**
+ * The 2nd half (48 coeffs) of a 96-tap symmetric lowpass filter
+ */
+static const double htaps[DSD_HTAPS] = {
+     0.09950731974056658,    0.09562845727714668,    0.08819647126516944,
+     0.07782552527068175,    0.06534876523171299,    0.05172629311427257,
+     0.0379429484910187,     0.02490921351762261,    0.0133774746265897,
+     0.003883043418804416,  -0.003284703416210726,  -0.008080250212687497,
+    -0.01067241812471033,   -0.01139427235000863,   -0.0106813877974587,
+    -0.009007905078766049,  -0.006828859761015335,  -0.004535184322001496,
+    -0.002425035959059578,  -0.0006922187080790708,  0.0005700762133516592,
+     0.001353838005269448,   0.001713709169690937,   0.001742046839472948,
+     0.001545601648013235,   0.001226696225277855,   0.0008704322683580222,
+     0.0005381636200535649,  0.000266446345425276,   7.002968738383528e-05,
+    -5.279407053811266e-05, -0.0001140625650874684, -0.0001304796361231895,
+    -0.0001189970287491285, -9.396247155265073e-05, -6.577634378272832e-05,
+    -4.07492895872535e-05,  -2.17407957554587e-05,  -9.163058931391722e-06,
+    -2.017460145032201e-06,  1.249721855219005e-06,  2.166655190537392e-06,
+     1.930520892991082e-06,  1.319400334374195e-06,  7.410039764949091e-07,
+     3.423230509967409e-07,  1.244182214744588e-07,  3.130441005359396e-08
+};
+
+static double ctables[DSD_CTABLES][256];
+
+static av_cold void dsd2pcm_ctables_tableinit(void)
+{
+    int t, e, m, sign;
+    double acc[DSD_CTABLES];
+    for (e = 0; e < 256; ++e) {
+        memset(acc, 0, sizeof(acc));
+        for (m = 0; m < 8; ++m) {
+            sign = (((e >> (7 - m)) & 1) * 2 - 1);
+            for (t = 0; t < DSD_CTABLES; ++t)
+                acc[t] += sign * htaps[t * 8 + m];
+        }
+        for (t = 0; t < DSD_CTABLES; ++t)
+            ctables[DSD_CTABLES - 1 - t][e] = acc[t];
+    }
+}
+
+av_cold void swri_dsd2pcm_init(void)
+{
+    static AVOnce init_static_once = AV_ONCE_INIT;
+    ff_thread_once(&init_static_once, dsd2pcm_ctables_tableinit);
+}
+
+void swri_dsd2pcm_translate(DSDContext *s, size_t samples,
+                            const uint8_t *src, ptrdiff_t src_stride,
+                            float *dst, ptrdiff_t dst_stride)
+{
+    uint8_t buf[DSD_FIFOSIZE];
+    unsigned pos, i;
+    uint8_t* p;
+    double sum;
+
+    pos = s->pos;
+
+    memcpy(buf, s->buf, sizeof(buf));
+
+    while (samples-- > 0) {
+        buf[pos] = *src;
+        src += src_stride;
+
+        p = buf + ((pos - DSD_CTABLES) & DSD_FIFOMASK);
+        *p = ff_reverse[*p];
+
+        sum = 0.0;
+        for (i = 0; i < DSD_CTABLES; i++) {
+            uint8_t a = buf[(pos                   - i) & DSD_FIFOMASK];
+            uint8_t b = buf[(pos - (DSD_CTABLES*2 - 1) + i) & DSD_FIFOMASK];
+            sum += ctables[i][a] + ctables[i][b];
+        }
+
+        *dst = (float)sum;
+        dst += dst_stride;
+
+        pos = (pos + 1) & DSD_FIFOMASK;
+    }
+
+    s->pos = pos;
+    memcpy(s->buf, buf, sizeof(buf));
+}

libswresample/dsd2pcm.h

@@ -0,0 +1,56 @@
+/*
+ * DSD to PCM conversion
+ * based on BSD licensed dsd2pcm by Sebastian Gesemann
+ * Copyright (c) 2009, 2011 Sebastian Gesemann. All rights reserved.
+ * Copyright (c) 2014 Peter Ross
+ *
+ * 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 SWRESAMPLE_DSD2PCM_H
+#define SWRESAMPLE_DSD2PCM_H
+
+#include <stddef.h>
+#include <stdint.h>
+
+#define DSD_HTAPS   48               /** number of FIR constants */
+#define DSD_FIFOSIZE 16              /** must be a power of two */
+#define DSD_FIFOMASK (DSD_FIFOSIZE - 1)  /** bit mask for FIFO offsets */
+
+#if DSD_FIFOSIZE * 8 < DSD_HTAPS * 2
+#error "DSD_FIFOSIZE too small"
+#endif
+
+/**
+ * Per-channel buffer
+ */
+typedef struct DSDContext {
+    uint8_t buf[DSD_FIFOSIZE];
+    unsigned pos;
+} DSDContext;
+
+void swri_dsd2pcm_init(void);
+
+/**
+ * Convert one channel of MSB-first DSD data (one byte = 8 samples) to
+ * float PCM at 1/8th of the DSD bit rate. Strides are in elements.
+ */
+void swri_dsd2pcm_translate(DSDContext *s, size_t samples,
+                            const uint8_t *src, ptrdiff_t src_stride,
+                            float *dst, ptrdiff_t dst_stride);
+
+#endif /* SWRESAMPLE_DSD2PCM_H */

libswresample/reverse.c

@@ -0,0 +1 @@
+#include "libavutil/reverse.c"

libswresample/swresample.c

@@ -177,6 +177,14 @@ av_cold int swr_init(struct SwrContext *s){
         return AVERROR(EINVAL);
     }
 
+    if (s->out_sample_fmt == AV_SAMPLE_FMT_DSD &&
+        !(s->in_sample_fmt == AV_SAMPLE_FMT_DSD &&
+          s->in_sample_rate == s->out_sample_rate &&
+          !(s->flags & SWR_FLAG_RESAMPLE))) {
+        av_log(s, AV_LOG_ERROR, "Conversion to DSD is not supported\n");
+        return AVERROR(EINVAL);
+    }
+
     s->out.ch_count  = s-> user_out_chlayout.nb_channels;
     s-> in.ch_count  = s->  user_in_chlayout.nb_channels;
 
@@ -225,8 +233,12 @@ av_cold int swr_init(struct SwrContext *s){
                  s->rematrix_custom;
 
     if(s->int_sample_fmt == AV_SAMPLE_FMT_NONE){
+        // DSD to PCM conversion is done in floating point
+        if(   s->in_sample_fmt == AV_SAMPLE_FMT_DSD
+           && s->out_sample_fmt != AV_SAMPLE_FMT_DSD) {
+            s->int_sample_fmt= AV_SAMPLE_FMT_FLTP;
         // 16bit or less to 16bit or less with the same sample rate
-        if(   av_get_bytes_per_sample(s-> in_sample_fmt) <= 2
+        } else if(   av_get_bytes_per_sample(s-> in_sample_fmt) <= 2
            && av_get_bytes_per_sample(s->out_sample_fmt) <= 2
            && s->out_sample_rate==s->in_sample_rate) {
             s->int_sample_fmt= AV_SAMPLE_FMT_S16P;
@@ -346,7 +358,10 @@ av_assert0(s->out.ch_count);
     if(!s->resample && !s->rematrix && !s->channel_map && !s->dither.method){
         s->full_convert = swri_audio_convert_alloc(s->out_sample_fmt,
                                                    s-> in_sample_fmt, s-> in.ch_count, NULL, 0);
-        return 0;
+        // fall through to the generic path for conversions that have no
+        // direct implementation (e.g. DSD input to non-float output)
+        if (s->full_convert)
+            return 0;
     }
 
     s->in_convert = swri_audio_convert_alloc(s->int_sample_fmt,
@@ -355,7 +370,10 @@ av_assert0(s->out.ch_count);
                                              s->int_sample_fmt, s->out.ch_count, NULL, 0);
 
     if (!s->in_convert || !s->out_convert) {
-        ret = AVERROR(ENOMEM);
+        av_log(s, AV_LOG_ERROR, "Cannot convert %s sample format to %s sample format\n",
+               av_get_sample_fmt_name(!s->in_convert ? s->in_sample_fmt : s->int_sample_fmt),
+               av_get_sample_fmt_name(!s->in_convert ? s->int_sample_fmt : s->out_sample_fmt));
+        ret = AVERROR(EINVAL);
         goto fail;
     }
 
@@ -865,10 +883,14 @@ int swr_inject_silence(struct SwrContext *s, int count){
     if((ret=swri_realloc_audio(&s->silence, count))<0)
         return ret;
 
-    if(s->silence.planar) for(i=0; i<s->silence.ch_count; i++) {
-        memset(s->silence.ch[i], s->silence.bps==1 ? 0x80 : 0, count*s->silence.bps);
-    } else
-        memset(s->silence.ch[0], s->silence.bps==1 ? 0x80 : 0, count*s->silence.bps*s->silence.ch_count);
+    {
+        int fill = s->silence.fmt == AV_SAMPLE_FMT_DSD ? 0x69 :
+                   s->silence.bps == 1                 ? 0x80 : 0;
+        if(s->silence.planar) for(i=0; i<s->silence.ch_count; i++) {
+            memset(s->silence.ch[i], fill, count*s->silence.bps);
+        } else
+            memset(s->silence.ch[0], fill, count*s->silence.bps*s->silence.ch_count);
+    }
 
     reversefill_audiodata(&s->silence, tmp_arg);
     av_log(s, AV_LOG_VERBOSE, "adding %d audio samples of silence\n", count);

libswresample/version.h

@@ -30,7 +30,7 @@
 
 #include "version_major.h"
 
-#define LIBSWRESAMPLE_VERSION_MINOR   2
+#define LIBSWRESAMPLE_VERSION_MINOR   3
 #define LIBSWRESAMPLE_VERSION_MICRO 100
 
 #define LIBSWRESAMPLE_VERSION_INT  AV_VERSION_INT(LIBSWRESAMPLE_VERSION_MAJOR, \

tests/ref/fate/source

@@ -16,6 +16,7 @@ libavformat/log2_tab.c
 libavformat/rangecoder_dec.c
 libavformat/riscv/cpu_common.c
 libswresample/log2_tab.c
+libswresample/reverse.c
 libswscale/aarch64/ops_entries.c
 libswscale/log2_tab.c
 libswscale/riscv/cpu_common.c