Farrow resampler implementation
Farrow resampler implementation in C, including optional bank switching.
The functionality is mainly identical to the matlab code snippet. A discussion of the optional bank switching variant (which may require fewer multiplications per output sample for a given signal quality) can be found in my blog entry.
While the matlab version is vectorized (it calculates all output samples at the same time), the code below operates sample-by-sample and is therefore closer to a hardware implementation.
The signal is not treated as cyclic. Instead, processing stops after the last input sample.
The code is meant to be readable, and should be optimized as needed.
The two example filters are identical to those in the matlab version, with the exception of the constant gain factor that was omitted.
Run from the Unix command line as
a.exe < inputdata.txt > outputdata.txt
On Windows, MinGW can be used, for example.
/* ****************************************************
* Farrow resampler (with optional bank switching)
* M. Nentwig, 2011
* Input stream is taken from stdin
* Output stream goes to stdout
* Main target was readability, is not optimized for efficiency.
* ****************************************************/
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <assert.h>
#if 1
/* **************************************************************
* 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 example filters is based on a 6th order Chebyshev Laplace-domain prototype.
* This version uses three sub-segments per tap (NBANKS = 3)
* **************************************************************/
const double 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};
#define NTAPS (14)
#define NBANKS (3)
#define ORDER (2)
#else
/* Alternative example
* Similar impulse response as above
* A conventional Farrow design (no subdivisions => NBANKS = 1), order 3
*/
const double 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};
#define NTAPS (14)
#define NBANKS (1)
#define ORDER (3)
#endif
/* Set here the ratio between output and input sample rate.
* It could be changed even during runtime! */
#define INCR (1.0 / 6.28 / 3)
/* delay line storage */
double delayLine[NTAPS];
/* Coefficient lookup "table" */
static double c(int ixTap, int ixBank, int ixPower){
return cMatrix[ixPower * (NTAPS * NBANKS) + ixTap * NBANKS + ixBank];
}
int main (void){
/* clear delay line */
int ix;
for (ix = 0; ix < NTAPS; ++ix)
delayLine[ix] = 0;
/* Index of last input sample that was read
* First valid sample starts at 0 */
int sampleIndexInput = -1;
/* Position of next output sample in input stream */
double sampleIndexOutput = 0.0;
/* loop forever. Terminate, once out of input data. */
while (1){
/* Split output sample location into integer and fractional part:
* Integer part corresponds to sample index in input stream
* fractional part [0, 1[ spans one tap (covering NBANKS segments) */
int sio_int = floor(sampleIndexOutput);
double sio_fracInTap = sampleIndexOutput - (double)sio_int;
assert((sio_fracInTap >= 0) && (sio_fracInTap < 1));
/* Further split the fractional part into
* - bank index
* - fractional part within the bank */
int sio_intBank = floor(sio_fracInTap * (double) NBANKS);
double sio_fracInBank = sio_fracInTap * (double) NBANKS - (double)sio_intBank;
assert((sio_fracInBank >= 0) && (sio_fracInBank < 1));
/* ****************************************************
* load new samples into the delay line, until the last
* processed input sample (sampleIndexInput) catches
* up with the integer part of the output stream position (sio_int)
* ***************************************************/
while (sampleIndexInput < sio_int){
/* Advance the delay line one step */
++sampleIndexInput;
for (ix = NTAPS-1; ix > 0; --ix)
delayLine[ix] = delayLine[ix-1];
/* Read one input sample */
int flag = scanf("%lf", &delayLine[0]);
if (flag != 1)
goto done; /* there's nothing wrong with "goto" as "break" for multiple loops ... */
} /* while delay line behind output data */
/* ****************************************************
* Calculate one output sample:
* "out" sums up the contribution of each tap
* ***************************************************/
double out = 0;
int ixTap;
for (ixTap = 0; ixTap < NTAPS; ++ixTap){
/* ****************************************************
* Contribution of tap "ixTap" to output:
* ***************************************************/
/* Evaluate polynomial in sio_fracInBank:
* c(ixTap, sio_intBank, 0) is the constant coefficient
* c(ixTap, sio_intBank, 1) is the linear coefficient etc
*/
double hornerSum = c(ixTap, sio_intBank, ORDER);
int ixPower;
for (ixPower = ORDER-1; ixPower >= 0; --ixPower){
hornerSum *= sio_fracInBank;
hornerSum += c(ixTap, sio_intBank, ixPower);
} /* for ixPower */
/* ****************************************************
* Weigh the delay line sample of this tap with the
* polynomial result and add to output
* ***************************************************/
out += hornerSum * delayLine[ixTap];
} /* for ixTap */
/* ****************************************************
* Generate output sample and advance the position of
* the next output sample in the input data stream
* ***************************************************/
printf("%1.12le\n", out);
sampleIndexOutput += INCR;
} /* loop until out of input data */
done: /* out of input data => break loops and continue here */
return 0;
} /* main */
