aboutsummaryrefslogtreecommitdiffstats
diff options
context:
space:
mode:
Diffstat
-rw-r--r--src/dsd2pcm/dsd2pcm.c188
-rw-r--r--src/dsd2pcm/dsd2pcm.h64
-rw-r--r--src/dsd2pcm/dsd2pcm.hpp41
-rw-r--r--src/dsd2pcm/info.txt38
-rw-r--r--src/dsd2pcm/main.cpp120
-rw-r--r--src/dsd2pcm/noiseshape.c83
-rw-r--r--src/dsd2pcm/noiseshape.h57
-rw-r--r--src/dsd2pcm/noiseshape.hpp46
8 files changed, 637 insertions, 0 deletions
diff --git a/src/dsd2pcm/dsd2pcm.c b/src/dsd2pcm/dsd2pcm.c
new file mode 100644
index 000000000..31f9a106d
--- /dev/null
+++ b/src/dsd2pcm/dsd2pcm.c
@@ -0,0 +1,188 @@
+#include <stdlib.h>
+#include <string.h>
+
+#include "dsd2pcm.h"
+
+#define HTAPS 48 /* number of FIR constants */
+#define FIFOSIZE 16 /* must be a power of two */
+#define FIFOMASK (FIFOSIZE-1) /* bit mask for FIFO offsets */
+#define CTABLES ((HTAPS+7)/8) /* number of "8 MACs" lookup tables */
+
+#if FIFOSIZE*8 < HTAPS*2
+#error "FIFOSIZE too small"
+#endif
+
+/*
+ * 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[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 float ctables[CTABLES][256];
+static unsigned char bitreverse[256];
+static int precalculated = 0;
+
+static void precalc()
+{
+ int t, e, m, k;
+ double acc;
+ if (precalculated) return;
+ for (t=0, e=0; t<256; ++t) {
+ bitreverse[t] = e;
+ for (m=128; m && !((e^=m)&m); m>>=1)
+ ;
+ }
+ for (t=0; t<CTABLES; ++t) {
+ k = HTAPS - t*8;
+ if (k>8) k=8;
+ for (e=0; e<256; ++e) {
+ acc = 0.0;
+ for (m=0; m<k; ++m) {
+ acc += (((e >> (7-m)) & 1)*2-1) * htaps[t*8+m];
+ }
+ ctables[CTABLES-1-t][e] = (float)acc;
+ }
+ }
+ precalculated = 1;
+}
+
+struct dsd2pcm_ctx_s
+{
+ unsigned char fifo[FIFOSIZE];
+ unsigned fifopos;
+};
+
+extern dsd2pcm_ctx* dsd2pcm_init()
+{
+ dsd2pcm_ctx* ptr;
+ if (!precalculated) precalc();
+ ptr = (dsd2pcm_ctx*) malloc(sizeof(dsd2pcm_ctx));
+ if (ptr) dsd2pcm_reset(ptr);
+ return ptr;
+}
+
+extern void dsd2pcm_destroy(dsd2pcm_ctx* ptr)
+{
+ free(ptr);
+}
+
+extern dsd2pcm_ctx* dsd2pcm_clone(dsd2pcm_ctx* ptr)
+{
+ dsd2pcm_ctx* p2;
+ p2 = (dsd2pcm_ctx*) malloc(sizeof(dsd2pcm_ctx));
+ if (p2) {
+ memcpy(p2,ptr,sizeof(dsd2pcm_ctx));
+ }
+ return p2;
+}
+
+extern void dsd2pcm_reset(dsd2pcm_ctx* ptr)
+{
+ int i;
+ for (i=0; i<FIFOSIZE; ++i)
+ ptr->fifo[i] = 0x69; /* my favorite silence pattern */
+ ptr->fifopos = 0;
+ /* 0x69 = 01101001
+ * This pattern "on repeat" makes a low energy 352.8 kHz tone
+ * and a high energy 1.0584 MHz tone which should be filtered
+ * out completely by any playback system --> silence
+ */
+}
+
+extern void dsd2pcm_translate(
+ dsd2pcm_ctx* ptr,
+ size_t samples,
+ const unsigned char *src, ptrdiff_t src_stride,
+ int lsbf,
+ float *dst, ptrdiff_t dst_stride)
+{
+ unsigned ffp;
+ unsigned i;
+ unsigned bite1, bite2;
+ unsigned char* p;
+ double acc;
+ ffp = ptr->fifopos;
+ lsbf = lsbf ? 1 : 0;
+ while (samples-- > 0) {
+ bite1 = *src & 0xFFu;
+ if (lsbf) bite1 = bitreverse[bite1];
+ ptr->fifo[ffp] = bite1; src += src_stride;
+ p = ptr->fifo + ((ffp-CTABLES) & FIFOMASK);
+ *p = bitreverse[*p & 0xFF];
+ acc = 0;
+ for (i=0; i<CTABLES; ++i) {
+ bite1 = ptr->fifo[(ffp -i) & FIFOMASK] & 0xFF;
+ bite2 = ptr->fifo[(ffp-(CTABLES*2-1)+i) & FIFOMASK] & 0xFF;
+ acc += ctables[i][bite1] + ctables[i][bite2];
+ }
+ *dst = (float)acc; dst += dst_stride;
+ ffp = (ffp + 1) & FIFOMASK;
+ }
+ ptr->fifopos = ffp;
+}
+
diff --git a/src/dsd2pcm/dsd2pcm.h b/src/dsd2pcm/dsd2pcm.h
new file mode 100644
index 000000000..26f6b6f10
--- /dev/null
+++ b/src/dsd2pcm/dsd2pcm.h
@@ -0,0 +1,64 @@
+#ifndef DSD2PCM_H_INCLUDED
+#define DSD2PCM_H_INCLUDED
+
+#include <stddef.h>
+#include <string.h>
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+struct dsd2pcm_ctx_s;
+
+typedef struct dsd2pcm_ctx_s dsd2pcm_ctx;
+
+/**
+ * initializes a "dsd2pcm engine" for one channel
+ * (precomputes tables and allocates memory)
+ *
+ * This is the only function that is not thread-safe in terms of the
+ * POSIX thread-safety definition because it modifies global state
+ * (lookup tables are computed during the first call)
+ */
+extern dsd2pcm_ctx* dsd2pcm_init();
+
+/**
+ * deinitializes a "dsd2pcm engine"
+ * (releases memory, don't forget!)
+ */
+extern void dsd2pcm_destroy(dsd2pcm_ctx *ctx);
+
+/**
+ * clones the context and returns a pointer to the
+ * newly allocated copy
+ */
+extern dsd2pcm_ctx* dsd2pcm_clone(dsd2pcm_ctx *ctx);
+
+/**
+ * resets the internal state for a fresh new stream
+ */
+extern void dsd2pcm_reset(dsd2pcm_ctx *ctx);
+
+/**
+ * "translates" a stream of octets to a stream of floats
+ * (8:1 decimation)
+ * @param ctx -- pointer to abstract context (buffers)
+ * @param samples -- number of octets/samples to "translate"
+ * @param src -- pointer to first octet (input)
+ * @param src_stride -- src pointer increment
+ * @param lsbitfirst -- bitorder, 0=msb first, 1=lsbfirst
+ * @param dst -- pointer to first float (output)
+ * @param dst_stride -- dst pointer increment
+ */
+extern void dsd2pcm_translate(dsd2pcm_ctx *ctx,
+ size_t samples,
+ const unsigned char *src, ptrdiff_t src_stride,
+ int lsbitfirst,
+ float *dst, ptrdiff_t dst_stride);
+
+#ifdef __cplusplus
+} /* extern "C" */
+#endif
+
+#endif /* include guard DSD2PCM_H_INCLUDED */
+
diff --git a/src/dsd2pcm/dsd2pcm.hpp b/src/dsd2pcm/dsd2pcm.hpp
new file mode 100644
index 000000000..b1b2ae1c5
--- /dev/null
+++ b/src/dsd2pcm/dsd2pcm.hpp
@@ -0,0 +1,41 @@
+#ifndef DSD2PCM_HXX_INCLUDED
+#define DSD2PCM_HXX_INCLUDED
+
+#include <algorithm>
+#include <stdexcept>
+#include "dsd2pcm.h"
+
+/**
+ * C++ PImpl Wrapper for the dsd2pcm C library
+ */
+
+class dxd
+{
+ dsd2pcm_ctx *handle;
+public:
+ dxd() : handle(dsd2pcm_init())
+ { if (!handle) throw std::runtime_error("wtf?!"); }
+
+ dxd(dxd const& x) : handle(dsd2pcm_clone(x.handle))
+ { if (!handle) throw std::runtime_error("wtf?!"); }
+
+ ~dxd() { dsd2pcm_destroy(handle); }
+
+ friend void swap(dxd & a, dxd & b)
+ { std::swap(a.handle,b.handle); }
+
+ dxd& operator=(dxd x)
+ { swap(*this,x); return *this; }
+
+ void translate(size_t samples,
+ const unsigned char *src, ptrdiff_t src_stride,
+ bool lsbitfirst,
+ float *dst, ptrdiff_t dst_stride)
+ {
+ dsd2pcm_translate(handle,samples,src,src_stride,
+ lsbitfirst,dst,dst_stride);
+ }
+};
+
+#endif // DSD2PCM_HXX_INCLUDED
+
diff --git a/src/dsd2pcm/info.txt b/src/dsd2pcm/info.txt
new file mode 100644
index 000000000..15ff29245
--- /dev/null
+++ b/src/dsd2pcm/info.txt
@@ -0,0 +1,38 @@
+You downloaded the source code for "dsd2pcm" which is a simple little
+"filter" program, that takes a DSD data stream on stdin and converts
+it to a PCM stream (352.8 kHz, either 16 or 24 bits) and writes it to
+stdout. The code is split into three modules:
+
+ (1) dsd2pcm
+
+ This is where the 8:1 decimation magic happens. It's an
+ implementation of a symmetric 96-taps FIR lowpass filter
+ optimized for DSD inputs. If you feed this converter with
+ DSD64 you get a PCM stream at 352.8 kHz and floating point
+ samples. This module is independent and can be reused.
+
+ (2) noiseshape
+
+ A module for applying generic noise shaping filters. It's
+ used for the 16-bit output mode in "main" to preserve the
+ dynamic range. This module is independent and can be reused.
+
+ (3) main.cpp (file contains the main function and handles I/O)
+
+The first two modules are pure C for maximum portability. In addition,
+there are C++ wrapper headers for convenient use of these modules in
+C++. The main application is a C++ application and makes use of the
+C++ headers to access the functionality of the first two modules.
+
+
+Under Linux this program is easily compiled by typing
+
+ g++ *.c *.cpp -O3 -o dsd2pcm
+
+provided you have GCC installed. That's why I didn't bother writing
+any makefiles. :-p
+
+
+Cheers!
+SG
+
diff --git a/src/dsd2pcm/main.cpp b/src/dsd2pcm/main.cpp
new file mode 100644
index 000000000..0b58888a8
--- /dev/null
+++ b/src/dsd2pcm/main.cpp
@@ -0,0 +1,120 @@
+#include <iostream>
+#include <vector>
+#include <cstring>
+
+#include "dsd2pcm.hpp"
+#include "noiseshape.hpp"
+
+namespace {
+
+const float my_ns_coeffs[] = {
+// b1 b2 a1 a2
+ -1.62666423, 0.79410094, 0.61367127, 0.23311013, // section 1
+ -1.44870017, 0.54196219, 0.03373857, 0.70316556 // section 2
+};
+
+const int my_ns_soscount = sizeof(my_ns_coeffs)/(sizeof(my_ns_coeffs[0])*4);
+
+inline long myround(float x)
+{
+ return static_cast<long>(x + (x>=0 ? 0.5f : -0.5f));
+}
+
+template<typename T>
+struct id { typedef T type; };
+
+template<typename T>
+inline T clip(
+ typename id<T>::type min,
+ T v,
+ typename id<T>::type max)
+{
+ if (v<min) return min;
+ if (v>max) return max;
+ return v;
+}
+
+inline void write_intel16(unsigned char * ptr, unsigned word)
+{
+ ptr[0] = word & 0xFF;
+ ptr[1] = (word >> 8) & 0xFF;
+}
+
+inline void write_intel24(unsigned char * ptr, unsigned long word)
+{
+ ptr[0] = word & 0xFF;
+ ptr[1] = (word >> 8) & 0xFF;
+ ptr[2] = (word >> 16) & 0xFF;
+}
+
+} // anonymous namespace
+
+using std::vector;
+using std::cin;
+using std::cout;
+using std::cerr;
+
+int main(int argc, char *argv[])
+{
+ const int block = 16384;
+ int channels = -1;
+ int lsbitfirst = -1;
+ int bits = -1;
+ if (argc==4) {
+ if ('1'<=argv[1][0] && argv[1][0]<='9') channels = 1 + (argv[1][0]-'1');
+ if (argv[2][0]=='m' || argv[2][0]=='M') lsbitfirst=0;
+ if (argv[2][0]=='l' || argv[2][0]=='L') lsbitfirst=1;
+ if (!strcmp(argv[3],"16")) bits = 16;
+ if (!strcmp(argv[3],"24")) bits = 24;
+ }
+ if (channels<1 || lsbitfirst<0 || bits<0) {
+ cerr << "\n"
+ "DSD2PCM filter (raw DSD64 --> 352 kHz raw PCM)\n"
+ "(c) 2009 Sebastian Gesemann\n\n"
+ "(filter as in \"reads data from stdin and writes to stdout\")\n\n"
+ "Syntax: dsd2pcm <channels> <bitorder> <bitdepth>\n"
+ "channels = 1,2,3,...,9 (number of channels in DSD stream)\n"
+ "bitorder = L (lsb first), M (msb first) (DSD stream option)\n"
+ "bitdepth = 16 or 24 (intel byte order, output option)\n\n"
+ "Note: At 16 bits/sample a noise shaper kicks in that can preserve\n"
+ "a dynamic range of 135 dB below 30 kHz.\n\n";
+ return 1;
+ }
+ int bytespersample = bits/8;
+ vector<dxd> dxds (channels);
+ vector<noise_shaper> ns;
+ if (bits==16) {
+ ns.resize(channels, noise_shaper(my_ns_soscount, my_ns_coeffs) );
+ }
+ vector<unsigned char> dsd_data (block * channels);
+ vector<float> float_data (block);
+ vector<unsigned char> pcm_data (block * channels * bytespersample);
+ char * const dsd_in = reinterpret_cast<char*>(&dsd_data[0]);
+ char * const pcm_out = reinterpret_cast<char*>(&pcm_data[0]);
+ while (cin.read(dsd_in,block * channels)) {
+ for (int c=0; c<channels; ++c) {
+ dxds[c].translate(block,&dsd_data[0]+c,channels,
+ lsbitfirst,
+ &float_data[0],1);
+ unsigned char * out = &pcm_data[0] + c*bytespersample;
+ if (bits==16) {
+ for (int s=0; s<block; ++s) {
+ float r = float_data[s]*32768 + ns[c].get();
+ long smp = clip(-32768,myround(r),32767);
+ ns[c].update( clip(-1,smp-r,1) );
+ write_intel16(out,smp);
+ out += channels*bytespersample;
+ }
+ } else {
+ for (int s=0; s<block; ++s) {
+ float r = float_data[s]*8388608;
+ long smp = clip(-8388608,myround(r),8388607);
+ write_intel24(out,smp);
+ out += channels*bytespersample;
+ }
+ }
+ }
+ cout.write(pcm_out,block*channels*bytespersample);
+ }
+}
+
diff --git a/src/dsd2pcm/noiseshape.c b/src/dsd2pcm/noiseshape.c
new file mode 100644
index 000000000..ecd2f251d
--- /dev/null
+++ b/src/dsd2pcm/noiseshape.c
@@ -0,0 +1,83 @@
+#include <stdlib.h>
+#include <string.h>
+
+#include "noiseshape.h"
+
+extern int noise_shape_init(
+ noise_shape_ctx *ctx,
+ int sos_count,
+ const float *coeffs)
+{
+ int i;
+ ctx->sos_count = sos_count;
+ ctx->bbaa = coeffs;
+ ctx->t1 = (float*) malloc(sizeof(float)*sos_count);
+ if (!ctx->t1) goto escape1;
+ ctx->t2 = (float*) malloc(sizeof(float)*sos_count);
+ if (!ctx->t2) goto escape2;
+ for (i=0; i<sos_count; ++i) {
+ ctx->t1[i] = 0.f;
+ ctx->t2[i] = 0.f;
+ }
+ return 0;
+escape2:
+ free(ctx->t1);
+escape1:
+ return -1;
+}
+
+extern void noise_shape_destroy(
+ noise_shape_ctx *ctx)
+{
+ free(ctx->t1);
+ free(ctx->t2);
+}
+
+extern int noise_shape_clone(
+ const noise_shape_ctx *from,
+ noise_shape_ctx *to)
+{
+ to->sos_count = from->sos_count;
+ to->bbaa = from->bbaa;
+ to->t1 = (float*) malloc(sizeof(float)*to->sos_count);
+ if (!to->t1) goto error1;
+ to->t2 = (float*) malloc(sizeof(float)*to->sos_count);
+ if (!to->t2) goto error2;
+ memcpy(to->t1,from->t1,sizeof(float)*to->sos_count);
+ memcpy(to->t2,from->t2,sizeof(float)*to->sos_count);
+ return 0;
+error2:
+ free(to->t1);
+error1:
+ return -1;
+}
+
+extern float noise_shape_get(noise_shape_ctx *ctx)
+{
+ int i;
+ float acc;
+ const float *c;
+ acc = 0.0;
+ c = ctx->bbaa;
+ for (i=0; i<ctx->sos_count; ++i) {
+ float t1i = ctx->t1[i];
+ float t2i = ctx->t2[i];
+ ctx->t2[i] = acc -= t1i * c[2] + t2i * c[3];
+ acc += t1i * c[0] + t2i * c[1];
+ c += 4;
+ }
+ return acc;
+}
+
+extern void noise_shape_update(noise_shape_ctx *ctx, float qerror)
+{
+ float *p;
+ int i;
+ for (i=0; i<ctx->sos_count; ++i) {
+ ctx->t2[i] += qerror;
+ }
+ p = ctx->t1;
+ ctx->t1 = ctx->t2;
+ ctx->t2 = p;
+}
+
diff --git a/src/dsd2pcm/noiseshape.h b/src/dsd2pcm/noiseshape.h
new file mode 100644
index 000000000..bb6b5bf5d
--- /dev/null
+++ b/src/dsd2pcm/noiseshape.h
@@ -0,0 +1,57 @@
+#ifndef NOISE_SHAPE_H_INCLUDED
+#define NOISE_SHAPE_H_INCLUDED
+
+#ifdef __cpluspluc
+extern "C" {
+#endif
+
+typedef struct noise_shape_ctx_s {
+ int sos_count; /* number of second order sections */
+ const float *bbaa; /* filter coefficients, owned by user */
+ float *t1, *t2; /* filter state, owned by ns library */
+} noise_shape_ctx;
+
+/**
+ * initializes a noise_shaper context
+ * returns an error code or 0
+ */
+extern int noise_shape_init(
+ noise_shape_ctx *ctx,
+ int sos_count,
+ const float *coeffs);
+
+/**
+ * destroys a noise_shaper context
+ */
+extern void noise_shape_destroy(
+ noise_shape_ctx *ctx);
+
+/**
+ * initializes a noise_shaper context so that its state
+ * is a copy of a given context
+ * returns an error code or 0
+ */
+extern int noise_shape_clone(
+ const noise_shape_ctx *from, noise_shape_ctx *to);
+
+/**
+ * computes the next "noise shaping sample". Note: This call
+ * alters the internal state. xxx_get and xxx_update must be
+ * called in an alternating manner.
+ */
+extern float noise_shape_get(
+ noise_shape_ctx *ctx);
+
+/**
+ * updates the noise shaper's state with the
+ * last quantization error
+ */
+extern void noise_shape_update(
+ noise_shape_ctx *ctx, float qerror);
+
+#ifdef __cpluspluc
+} /* extern "C" */
+#endif
+
+#endif /* NOISE_SHAPE_H_INCLUDED */
+
diff --git a/src/dsd2pcm/noiseshape.hpp b/src/dsd2pcm/noiseshape.hpp
new file mode 100644
index 000000000..726272f91
--- /dev/null
+++ b/src/dsd2pcm/noiseshape.hpp
@@ -0,0 +1,46 @@
+#ifndef NOISE_SHAPE_HXX_INCLUDED
+#define NOISE_SHAPE_HXX_INCLUDED
+
+#include <stdexcept>
+#include "noiseshape.h"
+
+/**
+ * C++ wrapper for the noiseshape C library
+ */
+
+class noise_shaper
+{
+ noise_shape_ctx ctx;
+public:
+ noise_shaper(int sos_count, const float *bbaa)
+ {
+ if (noise_shape_init(&ctx,sos_count,bbaa))
+ throw std::runtime_error("noise shaper initialization failed");
+ }
+
+ noise_shaper(noise_shaper const& x)
+ {
+ if (noise_shape_clone(&x.ctx,&ctx))
+ throw std::runtime_error("noise shaper initialization failed");
+ }
+
+ ~noise_shaper()
+ { noise_shape_destroy(&ctx); }
+
+ noise_shaper& operator=(noise_shaper const& x)
+ {
+ if (this != &x) {
+ noise_shape_destroy(&ctx);
+ if (noise_shape_clone(&x.ctx,&ctx))
+ throw std::runtime_error("noise shaper initialization failed");
+ }
+ return *this;
+ }
+
+ float get() { return noise_shape_get(&ctx); }
+
+ void update(float error) { noise_shape_update(&ctx,error); }
+};
+
+#endif /* NOISE_SHAPE_HXX_INCLUDED */
+
generated by cgit v1.2.3 (git 2.25.1) at 2024-11-15 01:56:47 +0000