secretrabbit/package.d

   1 /*
   2 ** Copyright (C) 2002-2011 Erik de Castro Lopo <erikd@mega-nerd.com>
   3 **
   4 ** This program is free software; you can redistribute it and/or modify
   5 ** it under the terms of the GNU General Public License as published by
   6 ** the Free Software Foundation; either version 2 of the License, or
   7 ** (at your option) any later version.
   8 **
   9 ** This program is distributed in the hope that it will be useful,
  10 ** but WITHOUT ANY WARRANTY; without even the implied warranty of
  11 ** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  12 ** GNU General Public License for more details.
  13 **
  14 ** You should have received a copy of the GNU General Public License
  15 ** along with this program; if not, write to the Free Software
  16 ** Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307, USA.
  17 */
  18 /*
  19 ** This code is part of Secret Rabbit Code aka libsamplerate. A commercial
  20 ** use license for this code is available, please see:
  21 **    http://www.mega-nerd.com/SRC/procedure.html
  22 */
  23 module iv.secretrabbit /*is aliced*/;
  24 import iv.alice;
  25
  26 //version = lib_secret_rabbit_allow_hq_filter;
  27 //version = lib_secret_rabbit_do_additional_checks;
  28
  29 version(lib_secret_rabbit_allow_hq_filter) enum rabbitHasHQ = true; else enum rabbitHasHQ = false;
  30
  31
  32 // ////////////////////////////////////////////////////////////////////////// //
  33 public struct SecretRabbit {
  34 public:
  35   enum Error {
  36     OK = 0,
  37     Memory,
  38     BadData,
  39     BadRatio,
  40     NotInitialized,
  41     BadConverter,
  42     BadChannelCount,
  43     DataAfterEOS,
  44     DataOverlaps,
  45   }
  46
  47   enum Type {
  48     Best,
  49     Medium,
  50     Fast,
  51   }
  52
  53   // Data is used to pass data to `process()`
  54   static struct Data {
  55     float[] dataIn, dataOut; // samples; `.length%channel` must be 0!
  56     // alas, we can't change `.length` in @nogc, so...
  57     usize inputUsed, outputUsed; // in array items, not frames!
  58     bool endOfInput; // set if this is the last input chunk
  59     double srcRatio; // destinationSampleRate/sourceSampleRate
  60   }
  61
  62 private:
  63   // private filter data
  64   double lastRatio, lastPosition;
  65
  66   Error lastError;
  67   int channels;
  68   bool wasEOI;
  69   Type interpType;
  70
  71   // private data for interpolators
  72   usize inCount, inUsed;
  73   usize outCount, outUsed;
  74
  75   // sinc data
  76   int indexInc;
  77
  78   double srcRatio, inputIndex;
  79
  80   immutable(float)[] coeffs;
  81
  82   int bCurrent, bEnd, bRealEnd, bLen;
  83
  84   // sure hope noone does more than 128 channels at once
  85   double[128] leftCalc;
  86   union {
  87     double[128] rightCalc;
  88     float[128] lastValue;
  89   }
  90
  91   float* buffer; // malloced
  92
  93   // other interpolators (linear and zoh)
  94   bool resetFlt;
  95   //float[128] lastValue;
  96   //alias lastValue = leftCalc;
  97
  98   // processing functions
  99   Error function (ref SecretRabbit filter, ref Data data) nothrow @trusted @nogc processFn; // process function
 100   void function (ref SecretRabbit filter) nothrow @trusted @nogc resetFn; // state reset
 101   void function (ref SecretRabbit filter) nothrow @trusted @nogc deinitFn; // free additional memory, etc
 102
 103 nothrow @trusted @nogc:
 104 public:
 105   this (Type flt, int chans) { setup(flt, chans); }
 106   ~this () { deinit(); }
 107
 108   @disable this (this); // no copies
 109
 110   @property Error error () pure const { return lastError; }
 111   @property Type type () pure const { return interpType; }
 112   @property bool wasComplete () pure const { return wasEOI; }
 113
 114   void setup (Type flt, int chans) {
 115     deinit();
 116     if (chans < 1 || chans > leftCalc.length) { lastError = Error.BadChannelCount; return; }
 117     interpType = flt;
 118     channels = chans;
 119     if ((lastError = setupSinc(flt)) != 0) return;
 120     reset();
 121   }
 122
 123   void deinit () {
 124     if (buffer !is null) {
 125       import core.stdc.stdlib : free;
 126       free(buffer);
 127       buffer = null;
 128     }
 129     processFn = null;
 130     lastError = Error.OK;
 131     channels = 0;
 132     wasEOI = false;
 133     coeffs = null;
 134   }
 135
 136   // `data.inputUsed` and `data.outputUsed` will be set
 137   Error process (ref Data data) {
 138     import std.math : isNaN;
 139
 140     // set the input and output counts to zero
 141     data.inputUsed = data.outputUsed = 0;
 142
 143     if (processFn is null) return (lastError = Error.NotInitialized);
 144
 145     // check for valid Data first
 146     if (data.dataIn.length > 0 && wasEOI) return (lastError = Error.DataAfterEOS);
 147     if (data.endOfInput) wasEOI = true;
 148
 149     // and that dataIn and dataOut are valid
 150     if (data.dataOut.length == 0 || data.dataIn.length%channels != 0 || data.dataOut.length%channels != 0) return (lastError = Error.BadData);
 151
 152     // check srcRatio is in range
 153     if (isBadSrcRatio(data.srcRatio)) return (lastError = Error.BadRatio);
 154
 155     if (data.dataIn.ptr < data.dataOut.ptr) {
 156       if (data.dataIn.ptr+data.dataIn.length > data.dataOut.ptr) return (lastError = Error.DataOverlaps);
 157     } else if (data.dataOut.ptr+data.dataOut.length > data.dataIn.ptr) {
 158       return (lastError = Error.DataOverlaps);
 159     }
 160
 161     // special case for when lastRatio has not been set
 162     if (isNaN(lastRatio) || lastRatio < (1.0/SRC_MAX_RATIO)) lastRatio = data.srcRatio;
 163
 164     // now process
 165     lastError = processFn(this, data);
 166
 167     data.inputUsed = inUsed/*/filter.channels*/;
 168     data.outputUsed = outUsed/*/filter.channels*/;
 169
 170     return lastError;
 171   }
 172
 173   // use this to immediately change resampling ratio instead of nicely sliding to it
 174   Error setRatio (double newRatio) {
 175     if (processFn is null) return (lastError = Error.NotInitialized);
 176     if (isBadSrcRatio(newRatio)) return (lastError = Error.BadRatio);
 177     lastRatio = newRatio;
 178     return (lastError = Error.OK);
 179   }
 180
 181   // if you want to start new conversion from the scratch, reset resampler
 182   Error reset () {
 183     import core.stdc.string : memset;
 184     lastPosition = 0.0;
 185     lastRatio = 0.0; // really?
 186     wasEOI = false;
 187     resetFlt = true;
 188     // sinc reset
 189     bCurrent = bEnd = 0;
 190     bRealEnd = -1;
 191     srcRatio = inputIndex = 0.0;
 192     if (buffer !is null) {
 193       if (bLen > 0) buffer[0..bLen] = 0;
 194       // set this for a sanity check
 195       memset(buffer+bLen, 0xAA, channels*buffer[0].sizeof);
 196     }
 197     leftCalc[] = 0;
 198     rightCalc[] = 0;
 199     //lastValue[] = 0;
 200     return (lastError = Error.OK);
 201   }
 202
 203 private:
 204   Error setupSinc (Type flt) nothrow @trusted @nogc {
 205     import core.stdc.math : lrint;
 206     import core.stdc.stdlib : malloc;
 207
 208     switch (channels) {
 209       case 1: processFn = &sincMonoProcessor; break;
 210       case 2: processFn = &sincStereoProcessor; break;
 211       default: processFn = &sincMultiChanProcessor; break;
 212     }
 213
 214     switch (flt) with (Type) {
 215       case Fast:
 216         coeffs = coeffsFT.coeffs;
 217         indexInc = coeffsFT.increment;
 218         break;
 219       case Medium:
 220         coeffs = coeffsMD.coeffs;
 221         indexInc = coeffsMD.increment;
 222         break;
 223       case Best:
 224         version(lib_secret_rabbit_allow_hq_filter) {
 225           coeffs = coeffsHQ.coeffs;
 226           indexInc = coeffsHQ.increment;
 227         } else {
 228           // use "medium" filter if we were compiled without "hq" one
 229           coeffs = coeffsMD.coeffs;
 230           indexInc = coeffsMD.increment;
 231         }
 232         break;
 233       default:
 234         return Error.BadConverter;
 235     }
 236
 237     /*
 238      * FIXME : This needs to be looked at more closely to see if there is
 239      * a better way. Need to look at prepareData() at the same time.
 240      */
 241
 242     bLen = lrint(2.5*(cast(int)coeffs.length-1)/(indexInc*1.0)*SRC_MAX_RATIO);
 243     bLen = max(bLen, 4096);
 244     bLen *= channels;
 245
 246     buffer = cast(float*)malloc(buffer[0].sizeof*(bLen+channels));
 247     if (buffer is null) return Error.Memory;
 248     //buffer = buf[0..bLen+channels];
 249
 250     version(lib_secret_rabbit_do_additional_checks) {
 251       import core.stdc.stdlib : free;
 252       int count = (cast(int)coeffs.length-1);
 253       int bits = void;
 254       for (bits = 0; (1<<bits) < count; ++bits) count |= (1<<bits);
 255       if (bits+SHIFT_BITS-1 >= int.sizeof*8) {
 256         free(buffer);
 257         assert(0, "SecretRabbit: corrupted filter data!");
 258       }
 259     }
 260
 261     return Error.OK;
 262   }
 263
 264 public:
 265 static:
 266   string errorStr (Error err) pure { static if (__VERSION__ > 2067) pragma(inline, true); return (err >= Error.min && err <= Error.max ? rabbitErrorStrings[err] : "rabbit-wtf"); }
 267   string name (Type flt) pure { static if (__VERSION__ > 2067) pragma(inline, true); return (flt >= Type.min && flt <= Type.max ? rabbitFilterNames[flt] : "invalid interpolator type"); }
 268   string description (Type flt) pure { static if (__VERSION__ > 2067) pragma(inline, true); return (flt >= Type.min && flt <= Type.max ? rabbitFilterDescs[flt] : "invalid interpolator type"); }
 269   bool isValidRatio (double ratio) pure { static if (__VERSION__ > 2067) pragma(inline, true); return !isBadSrcRatio(ratio); }
 270
 271   // will not resize output
 272   void short2float (in short[] input, float[] output) {
 273     if (output.length < input.length) assert(0, "invalid length");
 274     foreach (immutable idx, short v; input) output.ptr[idx] = cast(float)(v/(1.0*0x8000));
 275   }
 276
 277   // will not resize output
 278   void float2short (in float[] input, short[] output) {
 279     import core.stdc.math : lrint;
 280     double scaledValue = void;
 281     if (output.length < input.length) assert(0, "invalid length");
 282     foreach (immutable idx, float v; input) {
 283       scaledValue = v*(8.0*0x10000000);
 284       if (scaledValue >= 1.0*0x7FFFFFFF) output.ptr[idx] = 32767;
 285       else if (scaledValue <= -8.0*0x10000000) output.ptr[idx] = -32768;
 286       else output.ptr[idx] = cast(short)(lrint(scaledValue)>>16);
 287     }
 288   }
 289 }
 290
 291
 292 // ////////////////////////////////////////////////////////////////////////// //
 293 enum SRC_MAX_RATIO = 256;
 294 enum SRC_MIN_RATIO_DIFF = (1.0e-20);
 295
 296
 297 // ////////////////////////////////////////////////////////////////////////// //
 298 immutable string[SecretRabbit.Error.max+1] rabbitErrorStrings = [
 299   "no error",
 300   "out of memory",
 301   "invalid data passed",
 302   "ratio outside [1/256, 256] range",
 303   "trying to use uninitialized filter",
 304   "bad converter number",
 305   "channel count must be >= 1 and <= 128",
 306   "process() called without reset after end of input",
 307   "input and output data arrays overlap",
 308 ];
 309
 310 immutable string[SecretRabbit.Type.max+1] rabbitFilterNames = [
 311   "Best Sinc Interpolator",
 312   "Medium Sinc Interpolator",
 313   "Fastest Sinc Interpolator",
 314 ];
 315
 316 immutable string[SecretRabbit.Type.max+1] rabbitFilterDescs = [
 317   "Band limited sinc interpolation, fastest, 97dB SNR, 80% BW",
 318   "Band limited sinc interpolation, medium quality, 121dB SNR, 90% BW",
 319   "Band limited sinc interpolation, best quality, 145dB SNR, 96% BW",
 320 ];
 321
 322
 323 // ////////////////////////////////////////////////////////////////////////// //
 324 // sinc interpolator blobs
 325 static struct CoeffData {
 326   int increment;
 327   float[] coeffs;
 328 }
 329
 330 version(lib_secret_rabbit_allow_hq_filter) immutable CoeffData coeffsHQ;
 331 immutable CoeffData coeffsMD, coeffsFT;
 332
 333 shared static this () {
 334   version(lib_secret_rabbit_allow_hq_filter) immutable coeffs_hq_data = import("secretrabbit/coeffs_hq.bin");
 335   immutable coeffs_md_data = import("secretrabbit/coeffs_md.bin");
 336   immutable coeffs_ft_data = import("secretrabbit/coeffs_ft.bin");
 337
 338   version(lib_secret_rabbit_allow_hq_filter) {
 339     coeffsHQ.increment = *(cast(immutable(int)*)coeffs_hq_data);
 340     coeffsHQ.coeffs = (cast(immutable(float)*)(coeffs_hq_data.ptr+4))[0..(coeffs_hq_data.length-4)/4];
 341   }
 342
 343   coeffsMD.increment = *(cast(immutable(int)*)coeffs_md_data);
 344   coeffsMD.coeffs = (cast(immutable(float)*)(coeffs_md_data.ptr+4))[0..(coeffs_md_data.length-4)/4];
 345
 346   coeffsFT.increment = *(cast(immutable(int)*)coeffs_ft_data);
 347   coeffsFT.coeffs = (cast(immutable(float)*)(coeffs_ft_data.ptr+4))[0..(coeffs_ft_data.length-4)/4];
 348 }
 349
 350
 351 // ////////////////////////////////////////////////////////////////////////// //
 352 enum SHIFT_BITS = 12;
 353 enum FP_ONE = cast(double)(1<<SHIFT_BITS);
 354 enum INV_FP_ONE = (1.0/FP_ONE);
 355
 356 // quick sanity check
 357 static assert(SHIFT_BITS < int.sizeof*8-1, "internal error: SHIFT_BITS too large");
 358
 359 // ////////////////////////////////////////////////////////////////////////// //
 360 bool isBadSrcRatio() (double ratio) { static if (__VERSION__ > 2067) pragma(inline, true); import std.math : isNaN; return (isNaN(ratio) || ratio < (1.0/SRC_MAX_RATIO) || ratio > (1.0*SRC_MAX_RATIO)); }
 361
 362 double fmodOne() (double x) {
 363   static if (__VERSION__ > 2067) pragma(inline, true);
 364   import core.stdc.math : lrint;
 365   double res = x-lrint(x);
 366   return (res < 0.0 ? res+1.0 : res);
 367 }
 368
 369 T min(T) (T v0, T v1) { static if (__VERSION__ > 2067) pragma(inline, true); return (v0 < v1 ? v0 : v1); }
 370 T max(T) (T v0, T v1) { static if (__VERSION__ > 2067) pragma(inline, true); return (v0 < v1 ? v1 : v0); }
 371 T fabs(T) (T v) { static if (__VERSION__ > 2067) pragma(inline, true); return (v < 0 ? -v : v); }
 372
 373 int double2fp() (double x) { static if (__VERSION__ > 2067) pragma(inline, true); import core.stdc.math : lrint; return (lrint((x)*FP_ONE)); } /* double2fp */
 374 int int2fp() (int x) { static if (__VERSION__ > 2067) pragma(inline, true); return (x<<SHIFT_BITS); } /* int2fp */
 375 int fp2int() (int x) { static if (__VERSION__ > 2067) pragma(inline, true); return (x>>SHIFT_BITS); } /* fp2int */
 376 int fpfrac() (int x) { static if (__VERSION__ > 2067) pragma(inline, true); return (x&((1<<SHIFT_BITS)-1)); } /* fpfrac */
 377 double fp2double() (int x) { static if (__VERSION__ > 2067) pragma(inline, true); return fpfrac(x)*INV_FP_ONE; } /* fp2double */
 378
 379
 380 // ////////////////////////////////////////////////////////////////////////// //
 381 // Beware all ye who dare pass this point. There be dragons here.
 382 double sincCalcOutputMono (ref SecretRabbit filter, int increment, int startFilterIndex) nothrow @trusted @nogc {
 383   double fraction, left, right, icoeff;
 384   int filterIndex, maxFilterIndex;
 385   int dataIndex, coeffCount, indx;
 386
 387   // convert input parameters into fixed point
 388   maxFilterIndex = int2fp((cast(int)filter.coeffs.length-1));
 389
 390   // first apply the left half of the filter
 391   filterIndex = startFilterIndex;
 392   coeffCount = (maxFilterIndex-filterIndex)/increment;
 393   filterIndex = filterIndex+coeffCount*increment;
 394   dataIndex = filter.bCurrent-coeffCount;
 395
 396   left = 0.0;
 397   do {
 398     fraction = fp2double(filterIndex);
 399     indx = fp2int(filterIndex);
 400     icoeff = filter.coeffs.ptr[indx]+fraction*(filter.coeffs.ptr[indx+1]-filter.coeffs.ptr[indx]);
 401     left += icoeff*filter.buffer[dataIndex];
 402     filterIndex -= increment;
 403     dataIndex = dataIndex+1;
 404   } while (filterIndex >= 0);
 405
 406   // now apply the right half of the filter
 407   filterIndex = increment-startFilterIndex;
 408   coeffCount = (maxFilterIndex-filterIndex)/increment;
 409   filterIndex = filterIndex+coeffCount*increment;
 410   dataIndex = filter.bCurrent+1+coeffCount;
 411
 412   right = 0.0;
 413   do {
 414     fraction = fp2double(filterIndex);
 415     indx = fp2int(filterIndex);
 416     icoeff = filter.coeffs.ptr[indx]+fraction*(filter.coeffs.ptr[indx+1]-filter.coeffs.ptr[indx]);
 417     right += icoeff*filter.buffer[dataIndex];
 418     filterIndex -= increment;
 419     dataIndex = dataIndex-1;
 420   } while (filterIndex > 0);
 421
 422   return left+right;
 423 }
 424
 425
 426 SecretRabbit.Error sincMonoProcessor (ref SecretRabbit filter, ref SecretRabbit.Data data) nothrow @trusted @nogc {
 427   import core.stdc.math : lrint;
 428
 429   double inputIndex, srcRatio, count, floatIncrement, terminate, rem;
 430   int increment, startFilterIndex;
 431   int halfFilterChanLen, samplesInHand;
 432
 433   filter.inCount = data.dataIn.length/*mul chans*/;
 434   filter.outCount = data.dataOut.length/*mul chans*/;
 435   filter.inUsed = filter.outUsed = 0;
 436
 437   srcRatio = filter.lastRatio;
 438
 439   // check the sample rate ratio wrt the buffer len
 440   count = ((cast(int)filter.coeffs.length-1)+2.0)/filter.indexInc;
 441   if (min(filter.lastRatio, data.srcRatio) < 1.0) count /= min(filter.lastRatio, data.srcRatio);
 442
 443   // maximum coefficientson either side of center point
 444   halfFilterChanLen = filter.channels*(lrint(count)+1);
 445
 446   inputIndex = filter.lastPosition;
 447   floatIncrement = filter.indexInc;
 448
 449   rem = fmodOne (inputIndex);
 450   filter.bCurrent = (filter.bCurrent+filter.channels*lrint(inputIndex-rem))%filter.bLen;
 451   inputIndex = rem;
 452
 453   terminate = 1.0/srcRatio+1.0e-20;
 454
 455   // main processing loop
 456   while (filter.outUsed < filter.outCount) {
 457     // need to reload buffer?
 458     samplesInHand = (filter.bEnd-filter.bCurrent+filter.bLen)%filter.bLen;
 459
 460     if (samplesInHand <= halfFilterChanLen) {
 461       prepareData(filter, data, halfFilterChanLen);
 462       samplesInHand = (filter.bEnd-filter.bCurrent+filter.bLen)%filter.bLen;
 463       if (samplesInHand <= halfFilterChanLen) break;
 464     }
 465
 466     // this is the termination condition
 467     if (filter.bRealEnd >= 0 && filter.bCurrent+inputIndex+terminate >= filter.bRealEnd) break;
 468
 469     if (filter.outCount > 0 && fabs(filter.lastRatio-data.srcRatio) > 1.0e-10) {
 470       srcRatio = filter.lastRatio+filter.outUsed*(data.srcRatio-filter.lastRatio)/filter.outCount;
 471     }
 472
 473     floatIncrement = filter.indexInc*1.0;
 474     if (srcRatio < 1.0) floatIncrement = filter.indexInc*srcRatio;
 475
 476     increment = double2fp(floatIncrement);
 477
 478     startFilterIndex = double2fp(inputIndex*floatIncrement);
 479
 480     data.dataOut.ptr[filter.outUsed] = cast(float)((floatIncrement/filter.indexInc)*sincCalcOutputMono(filter, increment, startFilterIndex));
 481     ++filter.outUsed;
 482
 483     // figure out the next index
 484     inputIndex += 1.0/srcRatio;
 485     rem = fmodOne(inputIndex);
 486
 487     filter.bCurrent = (filter.bCurrent+filter.channels*lrint(inputIndex-rem))%filter.bLen;
 488     inputIndex = rem;
 489   }
 490
 491   filter.lastPosition = inputIndex;
 492
 493   // save current ratio rather then target ratio
 494   filter.lastRatio = srcRatio;
 495
 496   return SecretRabbit.Error.OK;
 497 }
 498
 499
 500 // ////////////////////////////////////////////////////////////////////////// //
 501 void sincCalcOutputStereo (ref SecretRabbit filter, int increment, int startFilterIndex, double scale, float* output) nothrow @trusted @nogc {
 502   double fraction, icoeff;
 503   double[2] left, right;
 504   int filterIndex, maxFilterIndex;
 505   int dataIndex, coeffCount, indx;
 506
 507   // convert input parameters into fixed point
 508   maxFilterIndex = int2fp((cast(int)filter.coeffs.length-1));
 509
 510   // first apply the left half of the filter
 511   filterIndex = startFilterIndex;
 512   coeffCount = (maxFilterIndex-filterIndex)/increment;
 513   filterIndex = filterIndex+coeffCount*increment;
 514   dataIndex = filter.bCurrent-filter.channels*coeffCount;
 515
 516   left.ptr[0] = left.ptr[1] = 0.0;
 517   do {
 518     fraction = fp2double(filterIndex);
 519     indx = fp2int(filterIndex);
 520     icoeff = filter.coeffs.ptr[indx]+fraction*(filter.coeffs.ptr[indx+1]-filter.coeffs.ptr[indx]);
 521     left.ptr[0] += icoeff*filter.buffer[dataIndex];
 522     left.ptr[1] += icoeff*filter.buffer[dataIndex+1];
 523     filterIndex -= increment;
 524     dataIndex = dataIndex+2;
 525   } while (filterIndex >= 0);
 526
 527   // now apply the right half of the filter
 528   filterIndex = increment-startFilterIndex;
 529   coeffCount = (maxFilterIndex-filterIndex)/increment;
 530   filterIndex = filterIndex+coeffCount*increment;
 531   dataIndex = filter.bCurrent+filter.channels*(1+coeffCount);
 532
 533   right.ptr[0] = right.ptr[1] = 0.0;
 534   do {
 535     fraction = fp2double (filterIndex);
 536     indx = fp2int (filterIndex);
 537     icoeff = filter.coeffs.ptr[indx]+fraction*(filter.coeffs.ptr[indx+1]-filter.coeffs.ptr[indx]);
 538     right.ptr[0] += icoeff*filter.buffer[dataIndex];
 539     right.ptr[1] += icoeff*filter.buffer[dataIndex+1];
 540     filterIndex -= increment;
 541     dataIndex = dataIndex-2;
 542   } while (filterIndex > 0);
 543
 544   output[0] = scale*(left.ptr[0]+right.ptr[0]);
 545   output[1] = scale*(left.ptr[1]+right.ptr[1]);
 546 }
 547
 548
 549 SecretRabbit.Error sincStereoProcessor (ref SecretRabbit filter, ref SecretRabbit.Data data) nothrow @trusted @nogc {
 550   import core.stdc.math : lrint;
 551
 552   double inputIndex, srcRatio, count, floatIncrement, terminate, rem;
 553   int increment, startFilterIndex;
 554   int halfFilterChanLen, samplesInHand;
 555
 556   filter.inCount = data.dataIn.length/*mul chans*/;
 557   filter.outCount = data.dataOut.length/*mul chans*/;
 558   filter.inUsed = filter.outUsed = 0;
 559
 560   srcRatio = filter.lastRatio;
 561
 562   // check the sample rate ratio wrt the buffer len
 563   count = ((cast(int)filter.coeffs.length-1)+2.0)/filter.indexInc;
 564   if (min(filter.lastRatio, data.srcRatio) < 1.0) count /= min(filter.lastRatio, data.srcRatio);
 565
 566   // maximum coefficientson either side of center point
 567   halfFilterChanLen = filter.channels*(lrint(count)+1);
 568
 569   inputIndex = filter.lastPosition;
 570   floatIncrement = filter.indexInc;
 571
 572   rem = fmodOne(inputIndex);
 573   filter.bCurrent = (filter.bCurrent+filter.channels*lrint(inputIndex-rem))%filter.bLen;
 574   inputIndex = rem;
 575
 576   terminate = 1.0/srcRatio+1e-20;
 577
 578   // main processing loop
 579   while (filter.outUsed < filter.outCount) {
 580     // need to reload buffer?
 581     samplesInHand = (filter.bEnd-filter.bCurrent+filter.bLen)%filter.bLen;
 582
 583     if (samplesInHand <= halfFilterChanLen) {
 584       prepareData(filter, data, halfFilterChanLen);
 585       samplesInHand = (filter.bEnd-filter.bCurrent+filter.bLen)%filter.bLen;
 586       if (samplesInHand <= halfFilterChanLen) break;
 587     }
 588
 589     // this is the termination condition
 590     if (filter.bRealEnd >= 0 && filter.bCurrent+inputIndex+terminate >= filter.bRealEnd) break;
 591
 592     if (filter.outCount > 0 && fabs(filter.lastRatio-data.srcRatio) > 1.0e-10) {
 593       srcRatio = filter.lastRatio+filter.outUsed*(data.srcRatio-filter.lastRatio)/filter.outCount;
 594     }
 595
 596     floatIncrement = filter.indexInc*1.0;
 597     if (srcRatio < 1.0) floatIncrement = filter.indexInc*srcRatio;
 598
 599     increment = double2fp(floatIncrement);
 600
 601     startFilterIndex = double2fp(inputIndex*floatIncrement);
 602
 603     sincCalcOutputStereo(filter, increment, startFilterIndex, floatIncrement/filter.indexInc, data.dataOut.ptr+filter.outUsed);
 604     filter.outUsed += 2;
 605
 606     // figure out the next index
 607     inputIndex += 1.0/srcRatio;
 608     rem = fmodOne(inputIndex);
 609
 610     filter.bCurrent = (filter.bCurrent+filter.channels*lrint(inputIndex-rem))%filter.bLen;
 611     inputIndex = rem;
 612   }
 613
 614   filter.lastPosition = inputIndex;
 615
 616   // save current ratio rather then target ratio
 617   filter.lastRatio = srcRatio;
 618
 619   return SecretRabbit.Error.OK;
 620 }
 621
 622
 623 // ////////////////////////////////////////////////////////////////////////// //
 624 void sincCalcOutputMultiChan (ref SecretRabbit filter, int increment, int startFilterIndex, int channels, double scale, float* output) nothrow @trusted @nogc {
 625   double fraction, icoeff;
 626   double* left, right;
 627   int filterIndex, maxFilterIndex;
 628   int dataIndex, coeffCount, indx, ch;
 629
 630   left = filter.leftCalc.ptr;
 631   right = filter.rightCalc.ptr;
 632
 633   // convert input parameters into fixed point
 634   maxFilterIndex = int2fp((cast(int)filter.coeffs.length-1));
 635
 636   // first apply the left half of the filter
 637   filterIndex = startFilterIndex;
 638   coeffCount = (maxFilterIndex-filterIndex)/increment;
 639   filterIndex = filterIndex+coeffCount*increment;
 640   dataIndex = filter.bCurrent-channels*coeffCount;
 641
 642   //memset(left, 0, left[0].sizeof*channels);
 643   left[0..channels] = 0;
 644   do {
 645     fraction = fp2double(filterIndex);
 646     indx = fp2int(filterIndex);
 647     icoeff = filter.coeffs.ptr[indx]+fraction*(filter.coeffs.ptr[indx+1]-filter.coeffs.ptr[indx]);
 648     /*
 649     **  Duff's Device.
 650     **  See : http://en.wikipedia.org/wiki/Duff's_device
 651     */
 652     ch = channels;
 653     do {
 654       switch (ch%8) {
 655         default:
 656           --ch;
 657           left[ch] += icoeff*filter.buffer[dataIndex+ch];
 658           goto case;
 659         case 7:
 660           --ch;
 661           left[ch] += icoeff*filter.buffer[dataIndex+ch];
 662           goto case;
 663         case 6:
 664           --ch;
 665           left[ch] += icoeff*filter.buffer[dataIndex+ch];
 666           goto case;
 667         case 5:
 668           --ch;
 669           left[ch] += icoeff*filter.buffer[dataIndex+ch];
 670           goto case;
 671         case 4:
 672           --ch;
 673           left[ch] += icoeff*filter.buffer[dataIndex+ch];
 674           goto case;
 675         case 3:
 676           --ch;
 677           left[ch] += icoeff*filter.buffer[dataIndex+ch];
 678           goto case;
 679         case 2:
 680           --ch;
 681           left[ch] += icoeff*filter.buffer[dataIndex+ch];
 682           goto case;
 683         case 1:
 684           --ch;
 685           left[ch] += icoeff*filter.buffer[dataIndex+ch];
 686           break;
 687       }
 688     } while (ch > 0);
 689
 690     filterIndex -= increment;
 691     dataIndex = dataIndex+channels;
 692   } while (filterIndex >= 0);
 693
 694   // now apply the right half of the filter
 695   filterIndex = increment-startFilterIndex;
 696   coeffCount = (maxFilterIndex-filterIndex)/increment;
 697   filterIndex = filterIndex+coeffCount*increment;
 698   dataIndex = filter.bCurrent+channels*(1+coeffCount);
 699
 700   //memset(right, 0, right[0].sizeof*channels);
 701   right[0..channels] = 0;
 702   do {
 703     fraction = fp2double (filterIndex);
 704     indx = fp2int (filterIndex);
 705     icoeff = filter.coeffs.ptr[indx]+fraction*(filter.coeffs.ptr[indx+1]-filter.coeffs.ptr[indx]);
 706     ch = channels;
 707     do {
 708       switch (ch%8) {
 709         default:
 710           --ch;
 711           right[ch] += icoeff*filter.buffer[dataIndex+ch];
 712           goto case;
 713         case 7:
 714           --ch;
 715           right[ch] += icoeff*filter.buffer[dataIndex+ch];
 716           goto case;
 717         case 6:
 718           --ch;
 719           right[ch] += icoeff*filter.buffer[dataIndex+ch];
 720           goto case;
 721         case 5:
 722           --ch;
 723           right[ch] += icoeff*filter.buffer[dataIndex+ch];
 724           goto case;
 725         case 4:
 726           --ch;
 727           right[ch] += icoeff*filter.buffer[dataIndex+ch];
 728           goto case;
 729         case 3:
 730           --ch;
 731           right[ch] += icoeff*filter.buffer[dataIndex+ch];
 732           goto case;
 733         case 2:
 734           --ch;
 735           right[ch] += icoeff*filter.buffer[dataIndex+ch];
 736           goto case;
 737         case 1:
 738           --ch;
 739           right[ch] += icoeff*filter.buffer[dataIndex+ch];
 740           break;
 741       }
 742     } while (ch > 0);
 743
 744     filterIndex -= increment;
 745     dataIndex = dataIndex-channels;
 746   } while (filterIndex > 0);
 747
 748   ch = channels;
 749   do {
 750     switch (ch%8) {
 751       default:
 752         --ch;
 753         output[ch] = scale*(left[ch]+right[ch]);
 754         goto case;
 755       case 7:
 756         --ch;
 757         output[ch] = scale*(left[ch]+right[ch]);
 758         goto case;
 759       case 6:
 760         --ch;
 761         output[ch] = scale*(left[ch]+right[ch]);
 762         goto case;
 763       case 5:
 764         --ch;
 765         output[ch] = scale*(left[ch]+right[ch]);
 766         goto case;
 767       case 4:
 768         --ch;
 769         output[ch] = scale*(left[ch]+right[ch]);
 770         goto case;
 771       case 3:
 772         --ch;
 773         output[ch] = scale*(left[ch]+right[ch]);
 774         goto case;
 775       case 2:
 776         --ch;
 777         output[ch] = scale*(left[ch]+right[ch]);
 778         goto case;
 779       case 1:
 780         --ch;
 781         output[ch] = scale*(left[ch]+right[ch]);
 782         break;
 783     }
 784   } while (ch > 0);
 785 }
 786
 787
 788 SecretRabbit.Error sincMultiChanProcessor (ref SecretRabbit filter, ref SecretRabbit.Data data) nothrow @trusted @nogc {
 789   import core.stdc.math : lrint;
 790
 791   double inputIndex, srcRatio, count, floatIncrement, terminate, rem;
 792   int increment, startFilterIndex;
 793   int halfFilterChanLen, samplesInHand;
 794
 795   filter.inCount = data.dataIn.length/*mul chans*/;
 796   filter.outCount = data.dataOut.length/*mul chans*/;
 797   filter.inUsed = filter.outUsed = 0;
 798
 799   srcRatio = filter.lastRatio;
 800
 801   // check the sample rate ratio wrt the buffer len
 802   count = ((cast(int)filter.coeffs.length-1)+2.0)/filter.indexInc;
 803   if (min(filter.lastRatio, data.srcRatio) < 1.0) count /= min(filter.lastRatio, data.srcRatio);
 804
 805   // maximum coefficientson either side of center point
 806   halfFilterChanLen = filter.channels*(lrint(count)+1);
 807
 808   inputIndex = filter.lastPosition;
 809   floatIncrement = filter.indexInc;
 810
 811   rem = fmodOne (inputIndex);
 812   filter.bCurrent = (filter.bCurrent+filter.channels*lrint(inputIndex-rem))%filter.bLen;
 813   inputIndex = rem;
 814
 815   terminate = 1.0/srcRatio+1e-20;
 816
 817   // main processing loop
 818   while (filter.outUsed < filter.outCount) {
 819     // need to reload buffer?
 820     samplesInHand = (filter.bEnd-filter.bCurrent+filter.bLen)%filter.bLen;
 821
 822     if (samplesInHand <= halfFilterChanLen) {
 823       prepareData(filter, data, halfFilterChanLen);
 824       samplesInHand = (filter.bEnd-filter.bCurrent+filter.bLen)%filter.bLen;
 825       if (samplesInHand <= halfFilterChanLen) break;
 826     }
 827
 828     // this is the termination condition
 829     if (filter.bRealEnd >= 0 && filter.bCurrent+inputIndex+terminate >= filter.bRealEnd) break;
 830
 831     if (filter.outCount > 0 && fabs(filter.lastRatio-data.srcRatio) > 1.0e-10) {
 832       srcRatio = filter.lastRatio+filter.outUsed*(data.srcRatio-filter.lastRatio)/filter.outCount;
 833     }
 834
 835     floatIncrement = filter.indexInc*1.0;
 836     if (srcRatio < 1.0) floatIncrement = filter.indexInc*srcRatio;
 837
 838     increment = double2fp(floatIncrement);
 839
 840     startFilterIndex = double2fp(inputIndex*floatIncrement);
 841
 842     sincCalcOutputMultiChan(filter, increment, startFilterIndex, filter.channels, floatIncrement/filter.indexInc, data.dataOut.ptr+filter.outUsed);
 843     filter.outUsed += filter.channels;
 844
 845     /* Figure out the next index. */
 846     inputIndex += 1.0/srcRatio;
 847     rem = fmodOne(inputIndex);
 848
 849     filter.bCurrent = (filter.bCurrent+filter.channels*lrint(inputIndex-rem))%filter.bLen;
 850     inputIndex = rem;
 851   }
 852
 853   filter.lastPosition = inputIndex;
 854
 855   // save current ratio rather then target ratio
 856   filter.lastRatio = srcRatio;
 857
 858   return SecretRabbit.Error.OK;
 859 }
 860
 861
 862 // ////////////////////////////////////////////////////////////////////////// //
 863 void prepareData (ref SecretRabbit filter, ref SecretRabbit.Data data, int halfFilterChanLen) nothrow @trusted @nogc {
 864   import core.stdc.string : memcpy, memmove, memset;
 865
 866   int len = 0;
 867   if (filter.bRealEnd >= 0) return; // should be terminating: just return
 868   if (filter.bCurrent == 0) {
 869     // initial state. Set up zeros at the start of the buffer and then load new data after that
 870     len = filter.bLen-2*halfFilterChanLen;
 871     filter.bCurrent = filter.bEnd = halfFilterChanLen;
 872   } else if (filter.bEnd+halfFilterChanLen+filter.channels < filter.bLen) {
 873     // load data at current end position
 874     len = max(filter.bLen-filter.bCurrent-halfFilterChanLen, 0);
 875   } else {
 876     // move data at end of buffer back to the start of the buffer
 877     len = filter.bEnd-filter.bCurrent;
 878     memmove(filter.buffer, filter.buffer+filter.bCurrent-halfFilterChanLen, (halfFilterChanLen+len)*filter.buffer[0].sizeof);
 879     filter.bCurrent = halfFilterChanLen;
 880     filter.bEnd = filter.bCurrent+len;
 881     // now load data at current end of buffer
 882     len = max(filter.bLen-filter.bCurrent-halfFilterChanLen, 0);
 883   }
 884   len = min(filter.inCount-filter.inUsed, len);
 885   len -= (len%filter.channels);
 886   if (len < 0 || filter.bEnd+len > filter.bLen) assert(0, "SecretRabbit internal error: bad length in prepareData()");
 887   memcpy(filter.buffer+filter.bEnd, data.dataIn.ptr+filter.inUsed, len*filter.buffer[0].sizeof);
 888   filter.bEnd += len;
 889   filter.inUsed += len;
 890   if (filter.inUsed == filter.inCount && filter.bEnd-filter.bCurrent < 2*halfFilterChanLen && data.endOfInput) {
 891     // handle the case where all data in the current buffer has been consumed and this is the last buffer
 892     if (filter.bLen-filter.bEnd < halfFilterChanLen+5) {
 893       // if necessary, move data down to the start of the buffer
 894       len = filter.bEnd-filter.bCurrent;
 895       memmove(filter.buffer, filter.buffer+filter.bCurrent-halfFilterChanLen, (halfFilterChanLen+len)*filter.buffer[0].sizeof);
 896       filter.bCurrent = halfFilterChanLen;
 897       filter.bEnd = filter.bCurrent+len;
 898     }
 899     filter.bRealEnd = filter.bEnd;
 900     len = halfFilterChanLen+5;
 901     if (len < 0 || filter.bEnd+len > filter.bLen) len = filter.bLen-filter.bEnd;
 902     memset(filter.buffer+filter.bEnd, 0, len*filter.buffer[0].sizeof);
 903     filter.bEnd += len;
 904   }
 905 }