diff options
Diffstat (limited to 'src/dsd2pcm')
| -rw-r--r-- | src/dsd2pcm/dsd2pcm.c | 184 | ||||
| -rw-r--r-- | src/dsd2pcm/dsd2pcm.h | 64 | ||||
| -rw-r--r-- | src/dsd2pcm/dsd2pcm.hpp | 41 | ||||
| -rw-r--r-- | src/dsd2pcm/info.txt | 38 | ||||
| -rw-r--r-- | src/dsd2pcm/main.cpp | 120 | ||||
| -rw-r--r-- | src/dsd2pcm/noiseshape.c | 83 | ||||
| -rw-r--r-- | src/dsd2pcm/noiseshape.h | 57 | ||||
| -rw-r--r-- | src/dsd2pcm/noiseshape.hpp | 46 |
8 files changed, 633 insertions, 0 deletions
diff --git a/src/dsd2pcm/dsd2pcm.c b/src/dsd2pcm/dsd2pcm.c new file mode 100644 index 000000000..4c7640853 --- /dev/null +++ b/src/dsd2pcm/dsd2pcm.c @@ -0,0 +1,184 @@ +#include "util/bit_reverse.h" + +#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 int precalculated = 0; + +static void precalc(void) +{ + int t, e, m, k; + double acc; + if (precalculated) return; + 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(void) +{ + 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 = bit_reverse(bite1); + ptr->fifo[ffp] = bite1; src += src_stride; + p = ptr->fifo + ((ffp-CTABLES) & FIFOMASK); + *p = bit_reverse(*p); + 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..80e8ce0cc --- /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(void); + +/** + * 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..6075f0d88 --- /dev/null +++ b/src/dsd2pcm/noiseshape.h @@ -0,0 +1,57 @@ +#ifndef NOISE_SHAPE_H_INCLUDED +#define NOISE_SHAPE_H_INCLUDED + +#ifdef __cplusplus +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 __cplusplus +} /* 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 */ + |