import dsd2pcm_src.zip

[this is the code from dsd2pcm_src.zip, published on a forum by
Sebastian Gesemann.  Upon request, he has given permission to
redistribute and modify his code, without referring to a specific
license. - mk]

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 */
+