Bandlimited Audio Resampling

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"
//		For printing, the PDF is more convenient:
//	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.
// =============================================================================================================

#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
	// 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.
	// 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
	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;
		// 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;
		// 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!
		// Increment time register
		t += dt;
Posted in Uncategorized | Leave a comment

Jordan Rudess is a fan!

Jordan Rudess, Dream Theater‘s extraordinary keyboard player, is a big fan of AudioStretch. He’s made it his “Featured App of the Week“, and did a nice YouTube demo. Thanks Jordan!

Posted in Uncategorized | Leave a comment

AudioStretch now supports Audiobus!

I’m happy to report that AudioStretch now supports Audiobus. Audiobus is the de-facto standard on iOS for sending audio from one app to another. It’s in version 3.6, which was released to the iOS App Store just minutes ago. (Note that Audiobus support is not available in the free edition, AudioStretch Lite).

As most readers of this blog know already, AudioStretch is my audio time-stretching app for iOS. It’s unique among apps of this sort in that you can scrub the audio by dragging the waveform, without changing the pitch even at zero speed, effectively allowing you to ‘freeze’ the sound. This makes it incredibly useful for music transcription as you can basically drag through the waveform, note-by-note, as you work out complicated passages (say, Charlie Parker or Joe Satriani solos).

While I originally designed AudioStretch mainly for the music transcription use case, many of my customers find the scrubbing feature makes it great as a sound design and manipulation tool. Folks in this sound designer crowd have been asking for Audiobus for months, on forums, reviews, and via support email.

Their wish is now granted :-)

Here’s an example created by one of my beta testers, Peter Schaufuss, scrubbing a bit of recorded speech in AudioStretch, feeding that into effects in IK Multimedia’s AmpliTube, and recording the resulting other-worldly sounds using Kymatica’s AudioShare.

Posted in Uncategorized | Leave a comment

AudioStretch ad in Guitar Techniques Magazine

I’m really excited to have a quarter-page ad running in the next three issue of UK-based Guitar Techniques magazine, starting with issue 237 (December 2014) which should appear on newsstands any day now.

The ad concept and design are by Toronto-based Aestheticize Media. The aim was to focus on the most unique feature of AudioStretch, namely the ability to ‘scrub’ the audio without changing the pitch, and do that in a way that excites and arouses curiosity. Aestheticize nailed it!


It’s a huge improvement over the crowded little 1/8 page ad I ran in the Classifieds section of Guitar Techniques last year. Since it was designed by an engineer (me!), it’s not surprising that it looks like something pulled from an engineering specification: basically a screenshot and list of features. Informative, but hardly exciting!  What was I thinking?




Posted in Uncategorized | 1 Comment

AudioStretch: “One Of The Essential Music Apps”

Nice video review of my app, AudioStretch for iOS.

Posted in Uncategorized | Leave a comment

No iCloud Drive on iOS? No problem

For quite a while, AudioStretch has been able to open audio files that are not in an iOS device’s Music library via the “Open In…” mechanism. For example, I can open the Dropbox app, choose an mp3, then click on the “Action” icon, choose “Open In…” then choose AudioStretch from the list of compatible apps that appears. After a second or two, the file opens up in AudioStretch.

I find this really handy for keeping a small selection of test files that I don’t want to have cluttering up my Music library, e.g. brief excerpts of various types of music, sine tones, etc.

Now that iOS 10.10 (aka Yosemite) and iOS 8 are out, I thought I’d try iCloud Drive instead. So I dragged a few audio files into iCloud Drive on my Mac, which is really easy, as iCloud Drive is nicely integrated with the Finder (as is Dropbox).

I figured I’d then open up Apple’s iCloud Drive app on my iPhone, choose one of the files, and use “Open In…” just like I do with Dropbox. Problem: Apple, rather incredibly, doesn’t have a iCloud Drive app!

With iOS 8, it’s now possible for an app to access files you put in iCloud Drive, but only if the app has been set up to do so. Specifically, it needs to use the UIDocumentMenuViewController class. I’m experimenting with it right now in my AudioStretch app. Works well, but before I can ship this feature, I have to figure out where to add yet another button to the UI, test it thoroughly, submit it to the app store, etc. which means it’ll be easily several weeks before it’s deployed.

In the meantime, there’s a stop gap measure: the Cloud Grab app. Ridiculously simple app: it just shows the documents you have in iCloud, and allows you to “Grab” and “Use” them (via Open In) in other apps. For $0.99 it was a no-brainer for me to buy it.

Seems a little crazy that Apple didn’t provide something like this out of the box, but hey, the developer of Cloud Grab was smart enough to spot this gap, so good on him! That’s entrepreneurship: spotting some little niche and filling it fast.


Posted in Uncategorized | Leave a comment

Predictable PIN numbers

An interesting analysis of PIN numbers in this article in Data Genetics shows astonishing lack of randomness in how people choose PIN numbers.   For example, just 20 4-digit numbers account for a staggering 26% of 4-digit PIN numbers, with “1234” alone accounting for over 10%!

Posted in Uncategorized | Leave a comment