Bank-switched Farrow resampler
A Farrow resampler variant that requires fewer multiplications per output sample, at the expense of more coefficients.
A conventional Farrow structure uses piecewise polynomial approximation on a continuous-time impulse response. In this variant, the impulse response is broken down further into sub-segments, allowing a reduction of approximation order and thus reduced computational complexity (fewer multiplications per output sample). The price is a higher number of coefficients and additional bank-switching logic.
To understand the code, it may be a good idea to start with the conventional Farrow resampler, which can be found here.
Comparing both implementations, the modifications needed to add bank-switching are relatively small.
Further discussion of the coefficient bank-switching concept can be found in the related blog entry.
It is mostly useful for impulse responses that require high (i.e. 5th) order approximation. Reducing the order of cubic (or even linear) approximation is usually not practical, since an excessive number of sub-segments is required to maintain accuracy, and the number of coefficients explodes.
The program contains two alternative Farrow resampler implementations that approximate the same impulse response with a length of 14 input samples:
- conventional:
14 segments, 3rd order polynomials.
4 multiplications / tap and output sample - bank-switched:
42 segments (3 / tap), 2nd order polynomial.
3 multiplications / tap and output sample
Code tested on OctaveForge and Matlab
% **************************************************************
% bank-switched Farrow resampler
% M. Nentwig, 2011
% Note: Uses cyclic signals (wraps around)
% **************************************************************
close all; clear all;
% inData contains input to the resampling process (instead of function arguments)
inData = struct();
% **************************************************************
% example coefficients.
% Each column [c0; c1; c2; ...] describes a polynomial for one tap coefficent in fractional time ft [0, 1]:
% tapCoeff = c0 + c1 * ft + c2 * ft ^ 2 + ...
% Each column corresponds to one tap.
% the matrix size may be changed arbitrarily.
%
% The example filter is based on a 6th order Chebyshev Laplace-domain prototype.
% **************************************************************
if false
% for comparison, this is a conventional design (no bank switching)
inData.cMatrix =[ -8.57738278e-3 7.82989032e-1 7.19303539e+000 6.90955718e+000 -2.62377450e+000 -6.85327127e-1 1.44681608e+000 -8.79147907e-1 7.82633997e-2 1.91318985e-1 -1.88573400e-1 6.91790782e-2 3.07723786e-3 -6.74800912e-3
2.32448021e-1 2.52624309e+000 7.67543936e+000 -8.83951796e+000 -5.49838636e+000 6.07298348e+000 -2.16053205e+000 -7.59142947e-1 1.41269409e+000 -8.17735712e-1 1.98119464e-1 9.15904145e-2 -9.18092030e-2 2.74136108e-2
-1.14183319e+000 6.86126458e+000 -6.86015957e+000 -6.35135894e+000 1.10745051e+001 -3.34847578e+000 -2.22405694e+000 3.14374725e+000 -1.68249886e+000 2.54083065e-1 3.22275037e-1 -3.04794927e-1 1.29393976e-1 -3.32026332e-2
1.67363115e+000 -2.93090391e+000 -1.13549165e+000 5.65274939e+000 -3.60291782e+000 -6.20715544e-1 2.06619782e+000 -1.42159644e+000 3.75075865e-1 1.88433333e-1 -2.64135123e-1 1.47117661e-1 -4.71871047e-2 1.24921920e-2] / 231.46 * 20;
inData.nBanks = 1;
else
% same example filter as above, but now the matrix contains three alternative coefficient banks for each tap.
% The order was reduced from cubic to quadratic.
% column 1: first bank, tap 1
% column 2: second bank, tap 1
% column 3: third bank, tap 1
% column 4: first bank, tap 2
% and so on
inData.cMatrix =[ 2.87810386e-4 4.70096244e-3 7.93412570e-2 4.39824536e-1 1.31192924e+000 2.67892232e+000 4.16465421e+000 5.16499621e+000 5.15592605e+000 3.99000369e+000 2.00785470e+000 -7.42377060e-2 -1.52569354e+000 -1.94402804e+000 -1.40915797e+000 -3.86484652e-1 5.44712939e-1 9.77559688e-1 8.32191447e-1 3.22691788e-1 -2.13133045e-1 -5.08501962e-1 -4.82928807e-1 -2.36313854e-1 4.76034568e-2 2.16891966e-1 2.20894063e-1 1.08361553e-1 -2.63421832e-2 -1.06276015e-1 -1.07491548e-1 -5.53793711e-2 4.86314061e-3 3.94357182e-2 4.06217506e-2 2.17199064e-2 1.60318761e-3 -8.40370106e-3 -8.10525279e-3 -3.62112499e-3 -4.13413072e-4 2.33101911e-4
-3.26760325e-3 -6.46028234e-3 1.46793247e-1 5.90235537e-1 1.18931309e+000 1.57853546e+000 1.40402774e+000 5.76506323e-1 -6.33522788e-1 -1.74564700e+000 -2.24153717e+000 -1.91309453e+000 -9.55568978e-1 1.58239169e-1 9.36193787e-1 1.10969783e+000 7.33284446e-1 1.06542194e-1 -4.15412084e-1 -6.06616434e-1 -4.54898908e-1 -1.20841199e-1 1.82941623e-1 3.12543429e-1 2.49935829e-1 8.05376898e-2 -7.83213666e-2 -1.47769751e-1 -1.18735248e-1 -3.70656555e-2 3.72608374e-2 6.71425397e-2 5.17812605e-2 1.55564930e-2 -1.40896327e-2 -2.35058137e-2 -1.59635057e-2 -3.44701792e-3 4.14108065e-3 4.56234829e-3 1.59503132e-3 -3.17301882e-4
5.64310141e-3 7.74786707e-2 2.11791763e-1 2.84703201e-1 1.85158633e-1 -8.41118142e-2 -3.98497442e-1 -5.86821615e-1 -5.40397941e-1 -2.47558080e-1 1.50864737e-1 4.59312895e-1 5.41539400e-1 3.84673917e-1 9.39576331e-2 -1.74932542e-1 -3.01635463e-1 -2.56239225e-1 -9.87146864e-2 6.82216764e-2 1.59795852e-1 1.48668245e-1 6.62563431e-2 -2.71234898e-2 -8.07045577e-2 -7.76841351e-2 -3.55333136e-2 1.23206602e-2 3.88535040e-2 3.64199073e-2 1.54608563e-2 -6.59814558e-3 -1.72735099e-2 -1.46307777e-2 -5.04363288e-3 3.31049461e-3 6.01267607e-3 3.83904192e-3 3.92549958e-4 -1.36315264e-3 -9.76017430e-4 7.46699178e-5] / 133.64 * 20;
inData.nBanks = 3;
end
% **************************************************************
% Create example signal
% **************************************************************
nIn = 50; % used for test signal generation only
if false
% complex signal
inData.signal = cos(2*pi*(0:(nIn-1)) / nIn) + cos(2*2*pi*(0:(nIn-1)) / nIn) + 1.5;
else
% impulse response
inData.signal = zeros(1, nIn); inData.signal(1) = 1; %
% inData.cMatrix = 0 * inData.cMatrix; inData.cMatrix(1, 1) = 1; % enable to show constant c in first tap, first bank
% inData.cMatrix = 0 * inData.cMatrix; inData.cMatrix(2, 1) = 1; % enable to show linear c in first tap, first bank
% inData.cMatrix = 0 * inData.cMatrix; inData.cMatrix(3, 1) = 1; % enable to show quadratic c in first tap, first bank
% inData.cMatrix = 0 * inData.cMatrix; inData.cMatrix(1, 2) = 1; % enable to show constant c in first tap, second bank
end
% **************************************************************
% Resample to the following number of output samples
% must be integer, otherwise arbitrary
% **************************************************************
inData.nSamplesOut = floor(nIn * 3 * 6.28);
% **************************************************************
% Set up Farrow resampling
% **************************************************************
nSamplesIn = size(inData.signal, 2);
nSamplesOut = inData.nSamplesOut;
order = size(inData.cMatrix, 1) - 1; % polynomial order
% number of input samples that contribute to one output sample (FIR size)
nTaps = size(inData.cMatrix, 2);
assert(mod(nTaps, inData.nBanks) == 0);
nTaps = nTaps / inData.nBanks; % only one out of nBanks coefficients contributes at any time
% pointer to the position in the input stream for each output sample (row vector, real numbers), starting at 0
inputIndex = (0:nSamplesOut-1) / nSamplesOut * nSamplesIn;
% split into integer part (0..nSamplesIn - 1) ...
inputIndexIntegerPart = floor(inputIndex);
% ... and fractional part [0, 1[
inputIndexFractionalPart = inputIndex - inputIndexIntegerPart;
% bank switching
% the fractional part is again split into an integer and fractional part
inputIndexFractionalPart = inputIndexFractionalPart * inData.nBanks;
inputIndexFractionalPart_int = floor(inputIndexFractionalPart); % coefficient bank index
inputIndexFractionalPart_frac = inputIndexFractionalPart - inputIndexFractionalPart_int; % fractional time 0..1 within each bank
% **************************************************************
% Calculate output stream
% First constant term (conventional FIR), then linear, quadratic, cubic, ...
% **************************************************************
outStream = zeros(1, inData.nSamplesOut);
for ixOrder = 0 : order
% note: fractional time is now defined for each sub-segment [0, 1[
x = inputIndexFractionalPart_frac .^ ixOrder;
% **************************************************************
% Add the contribution of each tap one-by-one
% **************************************************************
for ixTap = 0 : nTaps - 1
% coefficient bank switching: There are inData.nBanks alternative coefficients for each tap
c = inData.cMatrix(ixOrder+1, ixTap * inData.nBanks + inputIndexFractionalPart_int + 1);
% index of input sample that contributes to output via the current tap
% higher tap index => longer delay => older input sample => smaller data index
dataIx = inputIndexIntegerPart - ixTap;
% wrap around
dataIx = mod(dataIx, nSamplesIn);
% array indexing starts at 1
dataIx = dataIx + 1;
delayed = inData.signal(dataIx);
% for each individual output sample (index in row vector),
% - evaluate f = c(order, tapindex) * fracPart .^ order
% - scale the delayed input sample with f
% - accumulate the contribution of all taps
% this implementation performs the operation for all output samples in parallel (row vector)
outStream = outStream + c .* delayed .* x;
end % for ixTap
end % for ixOrder
% **************************************************************
% plot
% **************************************************************
xIn = linspace(0, 1, nSamplesIn + 1); xIn = xIn(1:end-1);
xOut = linspace(0, 1, nSamplesOut + 1); xOut = xOut(1:end-1);
figure(); grid on; hold on;
stem(xIn, inData.signal, 'k+-');
plot(xOut, outStream, 'b+-');
legend('input', 'output');
title('bank-switched Farrow resampling. Signals are cyclic.');
