src/dither.c

   1 /* Effect: dither/noise-shape   Copyright (c) 2008-9 robs@users.sourceforge.net
   2  *
   3  * This library is free software; you can redistribute it and/or modify it
   4  * under the terms of the GNU Lesser General Public License as published by
   5  * the Free Software Foundation; either version 2.1 of the License, or (at
   6  * your option) any later version.
   7  *
   8  * This library is distributed in the hope that it will be useful, but
   9  * WITHOUT ANY WARRANTY; without even the implied warranty of
  10  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser
  11  * General Public License for more details.
  12  *
  13  * You should have received a copy of the GNU Lesser General Public License
  14  * along with this library; if not, write to the Free Software Foundation,
  15  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
  16  */
  17
  18 #ifdef NDEBUG /* Enable assert always. */
  19 #undef NDEBUG /* Must undef above assert.h or other that might include it. */
  20 #endif
  21
  22 #include "sox_i.h"
  23 #include <assert.h>
  24
  25 #undef RANQD1
  26 #define RANQD1 ranqd1(p->ranqd1)
  27
  28 typedef enum { /* Collection of various filters from the net */
  29   Shape_none, Shape_lipshitz, Shape_f_weighted, Shape_modified_e_weighted,
  30   Shape_improved_e_weighted, Shape_gesemann, Shape_shibata, Shape_low_shibata, Shape_high_shibata
  31 } filter_name_t;
  32 static lsx_enum_item const filter_names[] = {
  33   LSX_ENUM_ITEM(Shape_,none)
  34   LSX_ENUM_ITEM(Shape_,lipshitz)
  35   {"f-weighted", Shape_f_weighted},
  36   {"modified-e-weighted", Shape_modified_e_weighted},
  37   {"improved-e-weighted", Shape_improved_e_weighted},
  38   LSX_ENUM_ITEM(Shape_,gesemann)
  39   LSX_ENUM_ITEM(Shape_,shibata)
  40   {"low-shibata", Shape_low_shibata},
  41   {"high-shibata", Shape_high_shibata},
  42   {0, 0}};
  43
  44 typedef struct {
  45   sox_rate_t  rate;
  46   enum {fir, iir} type;
  47   size_t len;
  48   int gain_cB; /* Chosen so clips are few if any, but not guaranteed none. */
  49   double const * coefs;
  50   filter_name_t name;
  51 } filter_t;
  52
  53 static double const lip44[] = {2.033, -2.165, 1.959, -1.590, .6149};
  54 static double const fwe44[] = {
  55   2.412, -3.370, 3.937, -4.174, 3.353, -2.205, 1.281, -.569, .0847};
  56 static double const mew44[] = {
  57   1.662, -1.263, .4827, -.2913, .1268, -.1124, .03252, -.01265, -.03524};
  58 static double const iew44[] = {
  59   2.847, -4.685, 6.214, -7.184, 6.639, -5.032, 3.263, -1.632, .4191};
  60 static double const ges44[] = {
  61   2.2061, -.4706, -.2534, -.6214, 1.0587, .0676, -.6054, -.2738};
  62 static double const ges48[] = {
  63   2.2374, -.7339, -.1251, -.6033, .903, .0116, -.5853, -.2571};
  64
  65 static double const shi48[] = {
  66   2.8720729351043701172,  -5.0413231849670410156,   6.2442994117736816406,
  67   -5.8483986854553222656, 3.7067542076110839844,  -1.0495119094848632812,
  68   -1.1830236911773681641,   2.1126792430877685547, -1.9094531536102294922,
  69   0.99913084506988525391, -0.17090806365013122559, -0.32615602016448974609,
  70   0.39127644896507263184, -0.26876461505889892578,  0.097676105797290802002,
  71   -0.023473845794796943665,
  72 };
  73 static double const shi44[] = {
  74   2.6773197650909423828,  -4.8308925628662109375,   6.570110321044921875,
  75   -7.4572014808654785156, 6.7263274192810058594,  -4.8481650352478027344,
  76   2.0412089824676513672,   0.7006359100341796875, -2.9537565708160400391,
  77   4.0800385475158691406,  -4.1845216751098632812,   3.3311812877655029297,
  78   -2.1179926395416259766,   0.879302978515625,      -0.031759146600961685181,
  79   -0.42382788658142089844, 0.47882103919982910156, -0.35490813851356506348,
  80   0.17496839165687561035, -0.060908168554306030273,
  81 };
  82 static double const shi38[] = {
  83   1.6335992813110351562,  -2.2615492343902587891,   2.4077029228210449219,
  84   -2.6341717243194580078, 2.1440362930297851562,  -1.8153258562088012695,
  85   1.0816224813461303711,  -0.70302653312683105469, 0.15991993248462677002,
  86   0.041549518704414367676, -0.29416576027870178223,  0.2518316805362701416,
  87   -0.27766478061676025391,  0.15785403549671173096, -0.10165894031524658203,
  88   0.016833892092108726501,
  89 };
  90 static double const shi32[] =
  91 { /* dmaker 32000: bestmax=4.99659 (inverted) */
  92 0.82118552923202515,
  93 -1.0063692331314087,
  94 0.62341964244842529,
  95 -1.0447187423706055,
  96 0.64532512426376343,
  97 -0.87615132331848145,
  98 0.52219754457473755,
  99 -0.67434263229370117,
 100 0.44954317808151245,
 101 -0.52557498216629028,
 102 0.34567299485206604,
 103 -0.39618203043937683,
 104 0.26791760325431824,
 105 -0.28936097025871277,
 106 0.1883765310049057,
 107 -0.19097308814525604,
 108 0.10431359708309174,
 109 -0.10633844882249832,
 110 0.046832218766212463,
 111 -0.039653312414884567,
 112 };
 113 static double const shi22[] =
 114 { /* dmaker 22050: bestmax=5.77762 (inverted) */
 115 0.056581053882837296,
 116 -0.56956905126571655,
 117 -0.40727734565734863,
 118 -0.33870288729667664,
 119 -0.29810553789138794,
 120 -0.19039161503314972,
 121 -0.16510021686553955,
 122 -0.13468159735202789,
 123 -0.096633769571781158,
 124 -0.081049129366874695,
 125 -0.064953058958053589,
 126 -0.054459091275930405,
 127 -0.043378707021474838,
 128 -0.03660014271736145,
 129 -0.026256965473294258,
 130 -0.018786206841468811,
 131 -0.013387725688517094,
 132 -0.0090983230620622635,
 133 -0.0026585909072309732,
 134 -0.00042083300650119781,
 135 };
 136 static double const shi16[] =
 137 { /* dmaker 16000: bestmax=5.97128 (inverted) */
 138 -0.37251132726669312,
 139 -0.81423574686050415,
 140 -0.55010956525802612,
 141 -0.47405767440795898,
 142 -0.32624706625938416,
 143 -0.3161766529083252,
 144 -0.2286367267370224,
 145 -0.22916607558727264,
 146 -0.19565616548061371,
 147 -0.18160104751586914,
 148 -0.15423151850700378,
 149 -0.14104481041431427,
 150 -0.11844276636838913,
 151 -0.097583092749118805,
 152 -0.076493598520755768,
 153 -0.068106919527053833,
 154 -0.041881654411554337,
 155 -0.036922425031661987,
 156 -0.019364040344953537,
 157 -0.014994367957115173,
 158 };
 159 static double const shi11[] =
 160 { /* dmaker 11025: bestmax=5.9406 (inverted) */
 161 -0.9264228343963623,
 162 -0.98695987462997437,
 163 -0.631156325340271,
 164 -0.51966935396194458,
 165 -0.39738872647285461,
 166 -0.35679301619529724,
 167 -0.29720726609230042,
 168 -0.26310476660728455,
 169 -0.21719355881214142,
 170 -0.18561814725399017,
 171 -0.15404847264289856,
 172 -0.12687471508979797,
 173 -0.10339745879173279,
 174 -0.083688631653785706,
 175 -0.05875682458281517,
 176 -0.046893671154975891,
 177 -0.027950936928391457,
 178 -0.020740609616041183,
 179 -0.009366452693939209,
 180 -0.0060260160826146603,
 181 };
 182 static double const shi08[] =
 183 { /* dmaker 8000: bestmax=5.56234 (inverted) */
 184 -1.202863335609436,
 185 -0.94103097915649414,
 186 -0.67878556251525879,
 187 -0.57650017738342285,
 188 -0.50004476308822632,
 189 -0.44349345564842224,
 190 -0.37833768129348755,
 191 -0.34028723835945129,
 192 -0.29413089156150818,
 193 -0.24994957447052002,
 194 -0.21715600788593292,
 195 -0.18792112171649933,
 196 -0.15268312394618988,
 197 -0.12135542929172516,
 198 -0.099610626697540283,
 199 -0.075273610651493073,
 200 -0.048787496984004974,
 201 -0.042586319148540497,
 202 -0.028991291299462318,
 203 -0.011869125068187714,
 204 };
 205 static double const shl48[] = {
 206   2.3925774097442626953,  -3.4350297451019287109,   3.1853709220886230469,
 207   -1.8117271661758422852, -0.20124770700931549072,  1.4759907722473144531,
 208   -1.7210904359817504883,   0.97746700048446655273, -0.13790138065814971924,
 209   -0.38185903429985046387,  0.27421241998672485352,  0.066584214568138122559,
 210   -0.35223302245140075684,  0.37672343850135803223, -0.23964276909828186035,
 211   0.068674825131893157959,
 212 };
 213 static double const shl44[] = {
 214   2.0833916664123535156,  -3.0418450832366943359,   3.2047898769378662109,
 215   -2.7571926116943359375, 1.4978630542755126953,  -0.3427594602108001709,
 216   -0.71733748912811279297,  1.0737057924270629883, -1.0225815773010253906,
 217   0.56649994850158691406, -0.20968692004680633545, -0.065378531813621520996,
 218   0.10322438180446624756, -0.067442022264003753662, -0.00495197344571352005,
 219 };
 220 static double const shh44[] = {
 221    3.0259189605712890625, -6.0268716812133789062,   9.195003509521484375,
 222    -11.824929237365722656, 12.767142295837402344, -11.917946815490722656,
 223    9.1739168167114257812,  -5.3712320327758789062, 1.1393624544143676758,
 224    2.4484779834747314453,  -4.9719839096069335938,   6.0392003059387207031,
 225    -5.9359521865844726562,  4.903278350830078125,   -3.5527443885803222656,
 226    2.1909697055816650391, -1.1672389507293701172,  0.4903914332389831543,
 227    -0.16519790887832641602,  0.023217858746647834778,
 228 };
 229
 230 static const filter_t filters[] = {
 231   {44100, fir,  5, 210, lip44, Shape_lipshitz},
 232   {46000, fir,  9, 276, fwe44, Shape_f_weighted},
 233   {46000, fir,  9, 160, mew44, Shape_modified_e_weighted},
 234   {46000, fir,  9, 321, iew44, Shape_improved_e_weighted},
 235   {48000, iir,  4, 220, ges48, Shape_gesemann},
 236   {44100, iir,  4, 230, ges44, Shape_gesemann},
 237   {48000, fir, 16, 301, shi48, Shape_shibata},
 238   {44100, fir, 20, 333, shi44, Shape_shibata},
 239   {37800, fir, 16, 240, shi38, Shape_shibata},
 240   {32000, fir, 20, 240/*TBD*/, shi32, Shape_shibata},
 241   {22050, fir, 20, 240/*TBD*/, shi22, Shape_shibata},
 242   {16000, fir, 20, 240/*TBD*/, shi16, Shape_shibata},
 243   {11025, fir, 20, 240/*TBD*/, shi11, Shape_shibata},
 244   { 8000, fir, 20, 240/*TBD*/, shi08, Shape_shibata},
 245   {48000, fir, 16, 250, shl48, Shape_low_shibata},
 246   {44100, fir, 15, 250, shl44, Shape_low_shibata},
 247   {44100, fir, 20, 383, shh44, Shape_high_shibata},
 248   {    0, fir,  0,   0,  NULL, Shape_none},
 249 };
 250
 251 #define MAX_N 20
 252
 253 typedef struct {
 254   filter_name_t filter_name;
 255   sox_bool      auto_detect, alt_tpdf;
 256   double        dummy;
 257
 258   double        previous_errors[MAX_N * 2];
 259   double        previous_outputs[MAX_N * 2];
 260   size_t        pos, prec;
 261   uint64_t      num_output;
 262   int32_t       history, ranqd1, r;
 263   double const  * coefs;
 264   sox_bool      dither_off;
 265   sox_effect_handler_flow flow;
 266 } priv_t;
 267
 268 #define CONVOLVE _ _ _ _
 269 #define NAME flow_iir_4
 270 #define IIR
 271 #define N 4
 272 #include "dither.h"
 273 #undef IIR
 274 #define CONVOLVE _ _ _ _ _
 275 #define NAME flow_fir_5
 276 #define N 5
 277 #include "dither.h"
 278 #define CONVOLVE _ _ _ _ _ _ _ _ _
 279 #define NAME flow_fir_9
 280 #define N 9
 281 #include "dither.h"
 282 #define CONVOLVE _ _ _ _ _ _ _ _ _ _ _ _ _ _ _
 283 #define NAME flow_fir_15
 284 #define N 15
 285 #include "dither.h"
 286 #define CONVOLVE _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _
 287 #define NAME flow_fir_16
 288 #define N 16
 289 #include "dither.h"
 290 #define CONVOLVE _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _
 291 #define NAME flow_fir_20
 292 #define N 20
 293 #include "dither.h"
 294
 295 static int flow_no_shape(sox_effect_t * effp, const sox_sample_t * ibuf,
 296     sox_sample_t * obuf, size_t * isamp, size_t * osamp)
 297 {
 298   priv_t * p = (priv_t *)effp->priv;
 299   size_t len = *isamp = *osamp = min(*isamp, *osamp);
 300
 301   while (len--) {
 302     if (p->auto_detect) {
 303       p->history = (p->history << 1) +
 304           !!(*ibuf & (((unsigned)-1) >> p->prec));
 305       if (p->history && p->dither_off) {
 306         p->dither_off = sox_false;
 307         lsx_debug("flow %" PRIuPTR ": on  @ %" PRIu64, effp->flow, p->num_output);
 308       } else if (!p->history && !p->dither_off) {
 309         p->dither_off = sox_true;
 310         lsx_debug("flow %" PRIuPTR ": off @ %" PRIu64, effp->flow, p->num_output);
 311       }
 312     }
 313
 314     if (!p->dither_off) {
 315       int32_t r = RANQD1 >> p->prec;
 316       double d = ((double)*ibuf++ + r + (p->alt_tpdf? -p->r : (RANQD1 >> p->prec))) / (1 << (32 - p->prec));
 317       int i = d < 0? d - .5 : d + .5;
 318       p->r = r;
 319       if (i <= (-1 << (p->prec-1)))
 320         ++effp->clips, *obuf = SOX_SAMPLE_MIN;
 321       else if (i > (int)SOX_INT_MAX(p->prec))
 322         ++effp->clips, *obuf = SOX_INT_MAX(p->prec) << (32 - p->prec);
 323       else *obuf = i << (32 - p->prec);
 324       ++obuf;
 325     }
 326     else
 327       *obuf++ = *ibuf++;
 328     ++p->num_output;
 329   }
 330   return SOX_SUCCESS;
 331 }
 332
 333 static int getopts(sox_effect_t * effp, int argc, char * * argv)
 334 {
 335   priv_t * p = (priv_t *)effp->priv;
 336   int c;
 337   lsx_getopt_t optstate;
 338   lsx_getopt_init(argc, argv, "+aSsf:p:", NULL, lsx_getopt_flag_none, 1, &optstate);
 339
 340   while ((c = lsx_getopt(&optstate)) != -1) switch (c) {
 341     case 'a': p->auto_detect = sox_true; break;
 342     case 'S': p->alt_tpdf = sox_true; break;
 343     case 's': p->filter_name = Shape_shibata; break;
 344     case 'f':
 345       p->filter_name = lsx_enum_option(c, optstate.arg, filter_names);
 346       if (p->filter_name == INT_MAX)
 347         return SOX_EOF;
 348       break;
 349     GETOPT_NUMERIC(optstate, 'p', prec, 1, 24)
 350     default: lsx_fail("invalid option `-%c'", optstate.opt); return lsx_usage(effp);
 351   }
 352   argc -= optstate.ind, argv += optstate.ind;
 353   return argc? lsx_usage(effp) : SOX_SUCCESS;
 354 }
 355
 356 static int start(sox_effect_t * effp)
 357 {
 358   priv_t * p = (priv_t *)effp->priv;
 359   double mult = 1; /* Amount the noise shaping multiplies up the TPDF (+/-1) */
 360
 361   if (p->prec == 0)
 362     p->prec = effp->out_signal.precision;
 363
 364   if (effp->in_signal.precision <= p->prec || p->prec > 24)
 365     return SOX_EFF_NULL;   /* Dithering not needed at this resolution */
 366
 367   if (p->prec == 1) {
 368     /* The general dither routines don't work in this case, so notify
 369        user and leave it at that for now.
 370        TODO: Some special-case treatment of 1-bit noise shaping will be
 371          needed for meaningful DSD write support. */
 372     lsx_warn("Dithering/noise-shaping to 1 bit is currently not supported.");
 373     return SOX_EFF_NULL;
 374   }
 375
 376   effp->out_signal.precision = p->prec;
 377
 378   p->flow = flow_no_shape;
 379   if (p->filter_name) {
 380     filter_t const * f;
 381
 382     for (f = filters; f->len && (f->name != p->filter_name || fabs(effp->in_signal.rate - f->rate) / f->rate > .05); ++f); /* 5% leeway on frequency */
 383     if (!f->len) {
 384       p->alt_tpdf |= effp->in_signal.rate >= 22050;
 385       if (!effp->flow)
 386         lsx_warn("no `%s' filter is available for rate %g; using %s TPDF",
 387             lsx_find_enum_value(p->filter_name, filter_names)->text,
 388             effp->in_signal.rate, p->alt_tpdf? "sloped" : "plain");
 389     }
 390     else {
 391       assert(f->len <= MAX_N);
 392       if (f->type == fir) switch(f->len) {
 393         case  5: p->flow = flow_fir_5 ; break;
 394         case  9: p->flow = flow_fir_9 ; break;
 395         case 15: p->flow = flow_fir_15; break;
 396         case 16: p->flow = flow_fir_16; break;
 397         case 20: p->flow = flow_fir_20; break;
 398         default: assert(sox_false);
 399       } else switch(f->len) {
 400         case  4: p->flow = flow_iir_4 ; break;
 401         default: assert(sox_false);
 402       }
 403       p->coefs = f->coefs;
 404       mult = dB_to_linear(f->gain_cB * 0.1);
 405     }
 406   }
 407   p->ranqd1 = ranqd1(sox_globals.ranqd1) + effp->flow;
 408   if (effp->in_signal.mult) /* (Takes account of ostart mult (sox.c). */
 409     *effp->in_signal.mult *= (SOX_SAMPLE_MAX - (1 << (31 - p->prec)) *
 410         (2 * mult + 1)) / (SOX_SAMPLE_MAX - (1 << (31 - p->prec)));
 411   return SOX_SUCCESS;
 412 }
 413
 414 static int flow(sox_effect_t * effp, const sox_sample_t * ibuf,
 415     sox_sample_t * obuf, size_t * isamp, size_t * osamp)
 416 {
 417   priv_t * p = (priv_t *)effp->priv;
 418   return p->flow(effp, ibuf, obuf, isamp, osamp);
 419 }
 420
 421 sox_effect_handler_t const * lsx_dither_effect_fn(void)
 422 {
 423   static sox_effect_handler_t handler = {
 424     "dither", "[-S|-s|-f filter] [-a] [-p precision]"
 425     "\n  (none)   Use TPDF"
 426     "\n  -S       Use sloped TPDF (without noise shaping)"
 427     "\n  -s       Shape noise (with shibata filter)"
 428     "\n  -f name  Set shaping filter to one of: lipshitz, f-weighted,"
 429     "\n           modified-e-weighted, improved-e-weighted, gesemann,"
 430     "\n           shibata, low-shibata, high-shibata."
 431     "\n  -a       Automatically turn on & off dithering as needed (use with caution!)"
 432     "\n  -p bits  Override the target sample precision",
 433     SOX_EFF_PREC, getopts, start, flow, 0, 0, 0, sizeof(priv_t)
 434   };
 435   return &handler;
 436 }