For years, whenever I wanted to resample an audio signal, for example to play it faster and at a higher pitch, I used simple linear interpolation. From a theoretical point of view, linear interpolation is pretty awful, but in practice it works surprisingly well with many signals.

Ideal resampling requires the use of an ideal low-pass filter, whose impulse response is a sinc function; since its impulse response has infinite length, in practice we use a finite length windowed sinc function.

I won’t go in to any further detail on the theory, as it’s all been explained very well in numerous papers and textbooks. This particular implementation is based on Stanford professor Julius O. Smith’s excellent Digital Audio Resampling Home Page.

Numerous simplifications and optimizations are possible, especially if you take advantage of fast platform-specific libraries such the vDSP functions in Apple’s Accelerate framework.

// =============================================================================================================
// Bandlimited Audio Interpolation Demo
//
// The code is a demonstration of bandlimited audio interpolation, and follows as closely as possible
// the method and notation as in [1], Julius O. Smith's excellent tutorial on the subject. While this code
// is reasonably well-commented, it will be very hard to understand how it works without refering to [1],
// unless you happen to already be intimately familiar with how to implement bandlimited interpolation.
//
// For simplicity, everything is in a single source file, including generation of a simple test signal.
// There are no dependencies on libraries other than standard ones included with any modern C++ compiler.
//
// The code could easily be repackaged into an Interpolator class, or perhaps turned into a command line
// tool that takes input and output wav files as arguments. If that's what you need, please note that since
// I'm releasing the code under the MIT License, you're free to do it yourself if you like. There's also plenty
// of scope for various optimizations and extensions.
//
// References:
// [1] Julius O. Smith, "Digital Audio Resampling Home Page"
// https://ccrma.stanford.edu/~jos/resample/
// For printing, the PDF is more convenient:
// https://ccrma.stanford.edu/~jos/resample/resample.pdf
//
//
// Released under the MIT License
//
// The MIT License (MIT)
//
// Copyright (c) 2015 Gerald T Beauregard
//
// Permission is hereby granted, free of charge, to any person obtaining a copy
// of this software and associated documentation files (the "Software"), to deal
// in the Software without restriction, including without limitation the rights
// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
// copies of the Software, and to permit persons to whom the Software is
// furnished to do so, subject to the following conditions:
//
// The above copyright notice and this permission notice shall be included in
// all copies or substantial portions of the Software.
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
// THE SOFTWARE.
// =============================================================================================================
#include <cstdio>
#include <cmath>
#include <cstdint>
#include <cfloat>
#include <vector>
#include <cassert>
using namespace std;
// Test signal generation function. The sample number i is a float
// so that we can generate in-between samples to compare with those
// obtained from the interpolator.
float genFunc(float i)
{
const float SIN_FREQ = 1000;
const float SAMPLE_RATE = 44100;
return sin(2*M_PI*i*SIN_FREQ/SAMPLE_RATE);
}
// -- Resampler variables and constants --
std::vector<float> h; // Filter coefficient table
std::vector<float> hb; // Differences in successive table values (h-bar in [1])
// Time register elements (as per Fig 6)
//
// Binary point -->|<------------ P [0,1)----------->|
//
// +----------------+----------------+----------------+
// | Input | Filter Table | Interpolation |
// | Signal | Initial | Factor eta |
// | Index n | Index l | |
// +---------------------------------+----------------+
//
// |<----- nn ----->|<----- nl ----->|<----- ne ----->|
//
const int MAX_TIME_BITS = 64; // Maximum total number of bits for time register
const int nn = 32; // Number of bits for input signal sample index
const int nl = 9; // Log2 of number of filter table elements per zero crossing
const int ne = 16; // Number of bits for filter table interpolation factor (Greek letter eta). [1] uses 8; more bits is better!
const int nP = nl + ne; // Total number of bits in P
const int L = 1 << nl; // Number of filter table elements per zero-crossing
const int Nz = 13; // Number of zero crossings per "wing" of windowed sinc
const int Nh = L*Nz; // Maximum index in tables h and hb
const uint64_t nMask = ((1L<<nn)-1) << nP; // Bit mask to extract input signal sample index n
const uint64_t eMask = ((1L<<ne)-1); // Bit mask to extract filter table interpolation factor eta
// Forward declaration of resampling functions
void initInterpolator();
void interpolate(vector<float> &x, vector<float> &y, float rho);
// ==============================================================================================
// main
// ==============================================================================================
int main(int argc, const char * argv[])
{
// Generate the input signal
const int N = 1024;
std::vector<float> x(N);
for (int i = 0; i < N; i++)
{
x[i] = genFunc(i);
}
// Output signal vector
std::vector<float> y;
// Initialize the resampling filter
initInterpolator();
// Do the interpolation, with resampling factor rho = output sample rate / input sample rate.
// rho > 1 means more output samples than input samples, rho < 1 means fewer output samples than input.
float rho = 1.2;
interpolate(x, y, rho);
// Compare interpolated output with expected output.
printf("%10s\t%10s\t%10s\t%10s\n", "index", "expected", "y[i]", "error");
for (int i = 0; i < y.size(); i++)
{
float expected = genFunc(i/rho);
float error = y[i] - expected;
printf("%10d\t%10.5f\t%10.5f\t%10.5f\n", i, expected, y[i], error);
}
return 0;
}
// ==============================================================================================
// initInterpolator
// ==============================================================================================
void initInterpolator()
{
assert(nn + nl + ne <= MAX_TIME_BITS);
// Resize the tables.
h.resize(Nh+1);
hb.resize(Nh+1);
// Compute the "right wing" of the filter table.
// We use letter l for the index for consistency with notation in [1],
// despite the possibility for confusion with the digit 1. Sorry.
for (int l = 0; l <= Nh; l++)
{
// Sinc function
float t = (float)l/L;
float sinc = (t == 0) ? 1.0 : sin(M_PI*t)/(M_PI*t);
// Blackman window. A Kaiser window is used in [1], but Blackman is
// so much simpler to calculate, and is good enough for my purposes.
float w = 0.42 + 0.5*cos(M_PI*l/Nh) + 0.08*cos(2*M_PI*l/Nh);
// Windowed sinc
h[l] = w * sinc;
}
// Get difference between successive elements of h. In [1] this
// is notated as h with a bar on top; we abbreviate "h-bar" as hb.
// Note that we scale these to compensate for fixed-point multiplication by eta below.
float scale = 1.0 / (1 << ne);
for (int l = 0; l < Nh; l++)
{
hb[l] = (h[l+1] - h[l]) * scale;
}
hb[Nh] = 0.0;
}
// ==============================================================================================
// interpolate
//
// @param x Buffer of input samples
// @param y Buffer for output samples
// @param rho Resampling factor Fsout / Fsin (Greek letter rho, looks like "p")
// ==============================================================================================
void interpolate(vector<float> &x, vector<float> &y, float rho)
{
int lenX = (int)x.size();
uint64_t maxT = (uint64_t)(lenX-1) << nP;
// Initialize fixed point input time register
uint64_t t = 0;
// Input time register increment per output sample
uint64_t dt = (1L << nP)/rho;
// Size of increment through the filter table h. When rho >= 1, we step through
// the filter table with a step-size of L. When rho < 1, the step-size through
// the filter table is reduced to rho L insead of L; this lowers the filter
// cutoff to avoid aliasing.
uint64_t dP = (L << ne);
if (rho < 1)
dP *= rho;
// For each output sample
y.resize(0);
while (t <= maxT)
//for (int yPos = 0; yPos < yLen; yPos++)
{
uint64_t P; // Position in h table, as fixed point number l.eta
uint64_t l; // Integer position in h table.
uint64_t eta; // Fractional position in h table (looks a bit like lower-case n)
uint64_t n; // Integer source sample position
// Extract integer source sample index from time register
n = (t & nMask) >> nP;
float v = 0;
// Get initial value for P as a fixed-point number in range [0,1) by
// masking out the integer sample number. When rho < 1, the initial
// P is replaced by P' = rho P.
P = t & ~nMask;
if (rho < 1)
P = rho*P;
// Compute sum for "left wing" as per equation (3)
// h end
// v <- Sum x(n-i)[h(l + iL) + eta hb(l + iL)
// i=0
int i = 0;
while (n >= i && (l = P >> ne) <= Nh)
{
eta = P & eMask;
v += x[n-i] * (h[l] + eta*hb[l]);
P += dP;
i++;
}
// Compute P <- 1 - P, as per equation (4).
// When rho < 1, 1-P becomes rho (1-P)
P = t & ~nMask;
P = ((1L << nP) - P);
if (rho < 1)
P = rho*P;
// Equation (5). Add in the "right wing".
// h end
// y(t) <- v + Sum x(n+1+i)[h(l + iL) + eta hb(l + iL)
// i=0
i = 0;
while (((n+1+i) < lenX) && (l = P >> ne) <= Nh)
{
eta = P & eMask;
v += x[n+1+i] * (h[l] + eta*hb[l]);
P += dP;
i++;
}
// When rho < 1, the sinc filter function should have been rescaled to maintain
// unity gain in the passband. Since the sinc is precaculated and not scaled, we
// compensate by scaling here.
if (rho < 1)
v *= rho;
// We have a new output sample!
y.push_back(v);
// Increment time register
t += dt;
}
}