plectrum

Plectrum: instrument tuner for Android
Log | Files | Refs | README | LICENSE

Resampler.java (17695B)


      1 /*
      2 *      _______                       _____   _____ _____  
      3 *     |__   __|                     |  __ \ / ____|  __ \ 
      4 *        | | __ _ _ __ ___  ___  ___| |  | | (___ | |__) |
      5 *        | |/ _` | '__/ __|/ _ \/ __| |  | |\___ \|  ___/ 
      6 *        | | (_| | |  \__ \ (_) \__ \ |__| |____) | |     
      7 *        |_|\__,_|_|  |___/\___/|___/_____/|_____/|_|     
      8 *                                                         
      9 * -------------------------------------------------------------
     10 *
     11 * TarsosDSP is developed by Joren Six at IPEM, University Ghent
     12 *  
     13 * -------------------------------------------------------------
     14 *
     15 *  Info: http://0110.be/tag/TarsosDSP
     16 *  Github: https://github.com/JorenSix/TarsosDSP
     17 *  Releases: http://0110.be/releases/TarsosDSP/
     18 *  
     19 *  TarsosDSP includes modified source code by various authors,
     20 *  for credits and info, see README.
     21 * 
     22 */
     23 
     24 /******************************************************************************
     25  *
     26  * libresample4j
     27  * Copyright (c) 2009 Laszlo Systems, Inc. All Rights Reserved.
     28  *
     29  * libresample4j is a Java port of Dominic Mazzoni's libresample 0.1.3,
     30  * which is in turn based on Julius Smith's Resample 1.7 library.
     31  *      http://www-ccrma.stanford.edu/~jos/resample/
     32  *
     33  * License: LGPL -- see the file LICENSE.txt for more information
     34  *
     35  *****************************************************************************/
     36 package be.tarsos.dsp.resample;
     37 
     38 import java.nio.FloatBuffer;
     39 
     40 public class Resampler {
     41 
     42     public static class Result {
     43         public final int inputSamplesConsumed;
     44         public final int outputSamplesGenerated;
     45 
     46         public Result(int inputSamplesConsumed, int outputSamplesGenerated) {
     47             this.inputSamplesConsumed = inputSamplesConsumed;
     48             this.outputSamplesGenerated = outputSamplesGenerated;
     49         }
     50     }
     51 
     52     // number of values per 1/delta in impulse response
     53     protected static final int Npc = 4096;
     54 
     55     private final float[] Imp;
     56     private final float[] ImpD;
     57     private final float LpScl;
     58     private final int Nmult;
     59     private final int Nwing;
     60     private final double minFactor;
     61     private final double maxFactor;
     62     private final int XSize;
     63     private final float[] X;
     64     private int Xp; // Current "now"-sample pointer for input
     65     private int Xread; // Position to put new samples
     66     private final int Xoff;
     67     private final float[] Y;
     68     private int Yp;
     69     private double Time;
     70 
     71     /**
     72      * Clone an existing resampling session. Faster than creating one from scratch.
     73      *
     74      * @param other
     75      */
     76     public Resampler(Resampler other) {
     77         this.Imp = other.Imp.clone();
     78         this.ImpD = other.ImpD.clone();
     79         this.LpScl = other.LpScl;
     80         this.Nmult = other.Nmult;
     81         this.Nwing = other.Nwing;
     82         this.minFactor = other.minFactor;
     83         this.maxFactor = other.maxFactor;
     84         this.XSize = other.XSize;
     85         this.X = other.X.clone();
     86         this.Xp = other.Xp;
     87         this.Xread = other.Xread;
     88         this.Xoff = other.Xoff;
     89         this.Y = other.Y.clone();
     90         this.Yp = other.Yp;
     91         this.Time = other.Time;
     92     }
     93 
     94     /**
     95      * Create a new resampling session.
     96      *
     97      * @param highQuality true for better quality, slower processing time
     98      * @param minFactor   lower bound on resampling factor for this session
     99      * @param maxFactor   upper bound on resampling factor for this session
    100      * @throws IllegalArgumentException if minFactor or maxFactor is not
    101      *                                  positive, or if maxFactor is less than minFactor
    102      */
    103     public Resampler(boolean highQuality, double minFactor, double maxFactor) {
    104         if (minFactor <= 0.0 || maxFactor <= 0.0) {
    105             throw new IllegalArgumentException("minFactor and maxFactor must be positive");
    106         }
    107         if (maxFactor < minFactor) {
    108             throw new IllegalArgumentException("minFactor must be <= maxFactor");
    109         }
    110 
    111         this.minFactor = minFactor;
    112         this.maxFactor = maxFactor;
    113         this.Nmult = highQuality ? 35 : 11;
    114         this.LpScl = 1.0f;
    115         this.Nwing = Npc * (this.Nmult - 1) / 2; // # of filter coeffs in right wing
    116 
    117         double Rolloff = 0.90;
    118         double Beta = 6;
    119 
    120         double[] Imp64 = new double[this.Nwing];
    121 
    122         FilterKit.lrsLpFilter(Imp64, this.Nwing, 0.5 * Rolloff, Beta, Npc);
    123         this.Imp = new float[this.Nwing];
    124         this.ImpD = new float[this.Nwing];
    125 
    126         for (int i = 0; i < this.Nwing; i++) {
    127             this.Imp[i] = (float) Imp64[i];
    128         }
    129 
    130         // Storing deltas in ImpD makes linear interpolation
    131         // of the filter coefficients faster
    132         for (int i = 0; i < this.Nwing - 1; i++) {
    133             this.ImpD[i] = this.Imp[i + 1] - this.Imp[i];
    134         }
    135 
    136         // Last coeff. not interpolated
    137         this.ImpD[this.Nwing - 1] = -this.Imp[this.Nwing - 1];
    138 
    139         // Calc reach of LP filter wing (plus some creeping room)
    140         int Xoff_min = (int) (((this.Nmult + 1) / 2.0) * Math.max(1.0, 1.0 / minFactor) + 10);
    141         int Xoff_max = (int) (((this.Nmult + 1) / 2.0) * Math.max(1.0, 1.0 / maxFactor) + 10);
    142         this.Xoff = Math.max(Xoff_min, Xoff_max);
    143 
    144         // Make the inBuffer size at least 4096, but larger if necessary
    145         // in order to store the minimum reach of the LP filter and then some.
    146         // Then allocate the buffer an extra Xoff larger so that
    147         // we can zero-pad up to Xoff zeros at the end when we reach the
    148         // end of the input samples.
    149         this.XSize = Math.max(2 * this.Xoff + 10, 4096);
    150         this.X = new float[this.XSize + this.Xoff];
    151         this.Xp = this.Xoff;
    152         this.Xread = this.Xoff;
    153 
    154         // Make the outBuffer long enough to hold the entire processed
    155         // output of one inBuffer
    156         int YSize = (int) (((double) this.XSize) * maxFactor + 2.0);
    157         this.Y = new float[YSize];
    158         this.Yp = 0;
    159 
    160         this.Time = (double) this.Xoff; // Current-time pointer for converter
    161     }
    162 
    163     public int getFilterWidth() {
    164         return this.Xoff;
    165     }
    166 
    167     /**
    168      * Process a batch of samples. There is no guarantee that the input buffer will be drained.
    169      *
    170      * @param factor    factor at which to resample this batch
    171      * @param buffers   sample buffer for producing input and consuming output
    172      * @param lastBatch true if this is known to be the last batch of samples
    173      * @return true iff resampling is complete (ie. no input samples consumed and no output samples produced)
    174      */
    175     public boolean process(double factor, SampleBuffers buffers, boolean lastBatch) {
    176         if (factor < this.minFactor || factor > this.maxFactor) {
    177             throw new IllegalArgumentException("factor " + factor + " is not between minFactor=" + minFactor
    178                     + " and maxFactor=" + maxFactor);
    179         }
    180 
    181         int outBufferLen = buffers.getOutputBufferLength();
    182         int inBufferLen = buffers.getInputBufferLength();
    183 
    184         float[] Imp = this.Imp;
    185         float[] ImpD = this.ImpD;
    186         float LpScl = this.LpScl;
    187         int Nwing = this.Nwing;
    188         boolean interpFilt = false; // TRUE means interpolate filter coeffs
    189 
    190         int inBufferUsed = 0;
    191         int outSampleCount = 0;
    192 
    193         // Start by copying any samples still in the Y buffer to the output
    194         // buffer
    195         if ((this.Yp != 0) && (outBufferLen - outSampleCount) > 0) {
    196             int len = Math.min(outBufferLen - outSampleCount, this.Yp);
    197 
    198             buffers.consumeOutput(this.Y, 0, len);
    199             //for (int i = 0; i < len; i++) {
    200             //    outBuffer[outBufferOffset + outSampleCount + i] = this.Y[i];
    201             //}
    202 
    203             outSampleCount += len;
    204             for (int i = 0; i < this.Yp - len; i++) {
    205                 this.Y[i] = this.Y[i + len];
    206             }
    207             this.Yp -= len;
    208         }
    209 
    210         // If there are still output samples left, return now - we need
    211         // the full output buffer available to us...
    212         if (this.Yp != 0) {
    213             return inBufferUsed == 0 && outSampleCount == 0;
    214         }
    215 
    216         // Account for increased filter gain when using factors less than 1
    217         if (factor < 1) {
    218             LpScl = (float) (LpScl * factor);
    219         }
    220 
    221         while (true) {
    222 
    223             // This is the maximum number of samples we can process
    224             // per loop iteration
    225 
    226             /*
    227              * #ifdef DEBUG
    228              * printf("XSize: %d Xoff: %d Xread: %d Xp: %d lastFlag: %d\n",
    229              * this.XSize, this.Xoff, this.Xread, this.Xp, lastFlag); #endif
    230              */
    231 
    232             // Copy as many samples as we can from the input buffer into X
    233             int len = this.XSize - this.Xread;
    234 
    235             if (len >= inBufferLen - inBufferUsed) {
    236                 len = inBufferLen - inBufferUsed;
    237             }
    238 
    239             buffers.produceInput(this.X, this.Xread, len);
    240             //for (int i = 0; i < len; i++) {
    241             //    this.X[this.Xread + i] = inBuffer[inBufferOffset + inBufferUsed + i];
    242             //}
    243 
    244             inBufferUsed += len;
    245             this.Xread += len;
    246 
    247             int Nx;
    248             if (lastBatch && (inBufferUsed == inBufferLen)) {
    249                 // If these are the last samples, zero-pad the
    250                 // end of the input buffer and make sure we process
    251                 // all the way to the end
    252                 Nx = this.Xread - this.Xoff;
    253                 for (int i = 0; i < this.Xoff; i++) {
    254                     this.X[this.Xread + i] = 0;
    255                 }
    256             } else {
    257                 Nx = this.Xread - 2 * this.Xoff;
    258             }
    259 
    260             /*
    261              * #ifdef DEBUG fprintf(stderr, "new len=%d Nx=%d\n", len, Nx);
    262              * #endif
    263              */
    264 
    265             if (Nx <= 0) {
    266                 break;
    267             }
    268 
    269             // Resample stuff in input buffer
    270             int Nout;
    271             if (factor >= 1) { // SrcUp() is faster if we can use it */
    272                 Nout = lrsSrcUp(this.X, this.Y, factor, /* &this.Time, */Nx, Nwing, LpScl, Imp, ImpD, interpFilt);
    273             } else {
    274                 Nout = lrsSrcUD(this.X, this.Y, factor, /* &this.Time, */Nx, Nwing, LpScl, Imp, ImpD, interpFilt);
    275             }
    276 
    277             /*
    278              * #ifdef DEBUG
    279              * printf("Nout: %d\n", Nout);
    280              * #endif
    281              */
    282 
    283             this.Time -= Nx; // Move converter Nx samples back in time
    284             this.Xp += Nx; // Advance by number of samples processed
    285 
    286             // Calc time accumulation in Time
    287             int Ncreep = (int) (this.Time) - this.Xoff;
    288             if (Ncreep != 0) {
    289                 this.Time -= Ncreep; // Remove time accumulation
    290                 this.Xp += Ncreep; // and add it to read pointer
    291             }
    292 
    293             // Copy part of input signal that must be re-used
    294             int Nreuse = this.Xread - (this.Xp - this.Xoff);
    295 
    296             for (int i = 0; i < Nreuse; i++) {
    297                 this.X[i] = this.X[i + (this.Xp - this.Xoff)];
    298             }
    299 
    300             /*
    301             #ifdef DEBUG
    302             printf("New Xread=%d\n", Nreuse);
    303             #endif */
    304 
    305             this.Xread = Nreuse; // Pos in input buff to read new data into
    306             this.Xp = this.Xoff;
    307 
    308             this.Yp = Nout;
    309 
    310             // Copy as many samples as possible to the output buffer
    311             if (this.Yp != 0 && (outBufferLen - outSampleCount) > 0) {
    312                 len = Math.min(outBufferLen - outSampleCount, this.Yp);
    313 
    314                 buffers.consumeOutput(this.Y, 0, len);
    315                 //for (int i = 0; i < len; i++) {
    316                 //    outBuffer[outBufferOffset + outSampleCount + i] = this.Y[i];
    317                 //}
    318 
    319                 outSampleCount += len;
    320                 for (int i = 0; i < this.Yp - len; i++) {
    321                     this.Y[i] = this.Y[i + len];
    322                 }
    323                 this.Yp -= len;
    324             }
    325 
    326             // If there are still output samples left, return now,
    327             //   since we need the full output buffer available
    328             if (this.Yp != 0) {
    329                 break;
    330             }
    331         }
    332 
    333         return inBufferUsed == 0 && outSampleCount == 0;
    334     }
    335 
    336     /**
    337      * Process a batch of samples. Convenience method for when the input and output are both floats.
    338      *
    339      * @param factor       factor at which to resample this batch
    340      * @param inputBuffer  contains input samples in the range -1.0 to 1.0
    341      * @param outputBuffer output samples will be deposited here
    342      * @param lastBatch    true if this is known to be the last batch of samples
    343      * @return true iff resampling is complete (ie. no input samples consumed and no output samples produced)
    344      */
    345     public boolean process(double factor, final FloatBuffer inputBuffer, boolean lastBatch, final FloatBuffer outputBuffer) {
    346         SampleBuffers sampleBuffers = new SampleBuffers() {
    347             public int getInputBufferLength() {
    348                 return inputBuffer.remaining();
    349             }
    350 
    351             public int getOutputBufferLength() {
    352                 return outputBuffer.remaining();
    353             }
    354 
    355             public void produceInput(float[] array, int offset, int length) {
    356                 inputBuffer.get(array, offset, length);
    357             }
    358 
    359             public void consumeOutput(float[] array, int offset, int length) {
    360                 outputBuffer.put(array, offset, length);
    361             }
    362         };
    363         return process(factor, sampleBuffers, lastBatch);
    364     }
    365 
    366     /**
    367      * Process a batch of samples. Alternative interface if you prefer to work with arrays.
    368      *
    369      * @param factor         resampling rate for this batch
    370      * @param inBuffer       array containing input samples in the range -1.0 to 1.0
    371      * @param inBufferOffset offset into inBuffer at which to start processing
    372      * @param inBufferLen    number of valid elements in the inputBuffer
    373      * @param lastBatch      pass true if this is the last batch of samples
    374      * @param outBuffer      array to hold the resampled data
    375      * @param outBufferOffset Offset in the output buffer.
    376      * @param outBufferLen    Output buffer length.
    377      * @return the number of samples consumed and generated
    378      */
    379     public Result process(double factor, float[] inBuffer, int inBufferOffset, int inBufferLen, boolean lastBatch, float[] outBuffer, int outBufferOffset, int outBufferLen) {
    380         FloatBuffer inputBuffer = FloatBuffer.wrap(inBuffer, inBufferOffset, inBufferLen);
    381         FloatBuffer outputBuffer = FloatBuffer.wrap(outBuffer, outBufferOffset, outBufferLen);
    382 
    383         process(factor, inputBuffer, lastBatch, outputBuffer);
    384 
    385         return new Result(inputBuffer.position() - inBufferOffset, outputBuffer.position() - outBufferOffset);
    386     }
    387 
    388 
    389 
    390     /*
    391      * Sampling rate up-conversion only subroutine; Slightly faster than
    392      * down-conversion;
    393      */
    394     private int lrsSrcUp(float X[], float Y[], double factor, int Nx, int Nwing, float LpScl, float Imp[],
    395                          float ImpD[], boolean Interp) {
    396 
    397         float[] Xp_array = X;
    398         int Xp_index;
    399 
    400         float[] Yp_array = Y;
    401         int Yp_index = 0;
    402 
    403         float v;
    404 
    405         double CurrentTime = this.Time;
    406         double dt; // Step through input signal
    407         double endTime; // When Time reaches EndTime, return to user
    408 
    409         dt = 1.0 / factor; // Output sampling period
    410 
    411         endTime = CurrentTime + Nx;
    412         while (CurrentTime < endTime) {
    413             double LeftPhase = CurrentTime - Math.floor(CurrentTime);
    414             double RightPhase = 1.0 - LeftPhase;
    415 
    416             Xp_index = (int) CurrentTime; // Ptr to current input sample
    417             // Perform left-wing inner product
    418             v = FilterKit.lrsFilterUp(Imp, ImpD, Nwing, Interp, Xp_array, Xp_index++, LeftPhase, -1);
    419             // Perform right-wing inner product
    420             v += FilterKit.lrsFilterUp(Imp, ImpD, Nwing, Interp, Xp_array, Xp_index, RightPhase, 1);
    421 
    422             v *= LpScl; // Normalize for unity filter gain
    423 
    424             Yp_array[Yp_index++] = v; // Deposit output
    425             CurrentTime += dt; // Move to next sample by time increment
    426         }
    427 
    428         this.Time = CurrentTime;
    429         return Yp_index; // Return the number of output samples
    430     }
    431 
    432     private int lrsSrcUD(float X[], float Y[], double factor, int Nx, int Nwing, float LpScl, float Imp[],
    433                          float ImpD[], boolean Interp) {
    434 
    435         float[] Xp_array = X;
    436         int Xp_index;
    437 
    438         float[] Yp_array = Y;
    439         int Yp_index = 0;
    440 
    441         float v;
    442 
    443         double CurrentTime = this.Time;
    444         double dh; // Step through filter impulse response
    445         double dt; // Step through input signal
    446         double endTime; // When Time reaches EndTime, return to user
    447 
    448         dt = 1.0 / factor; // Output sampling period
    449 
    450         dh = Math.min(Npc, factor * Npc); // Filter sampling period
    451 
    452         endTime = CurrentTime + Nx;
    453         while (CurrentTime < endTime) {
    454             double LeftPhase = CurrentTime - Math.floor(CurrentTime);
    455             double RightPhase = 1.0 - LeftPhase;
    456 
    457             Xp_index = (int) CurrentTime; // Ptr to current input sample
    458             // Perform left-wing inner product
    459             v = FilterKit.lrsFilterUD(Imp, ImpD, Nwing, Interp, Xp_array, Xp_index++, LeftPhase, -1, dh);
    460             // Perform right-wing inner product
    461             v += FilterKit.lrsFilterUD(Imp, ImpD, Nwing, Interp, Xp_array, Xp_index, RightPhase, 1, dh);
    462 
    463             v *= LpScl; // Normalize for unity filter gain
    464 
    465             Yp_array[Yp_index++] = v; // Deposit output
    466 
    467             CurrentTime += dt; // Move to next sample by time increment
    468         }
    469 
    470         this.Time = CurrentTime;
    471         return Yp_index; // Return the number of output samples
    472     }
    473 
    474 }