# Help with resampling algorithm implementation

I am trying to implement the resampling algorithm as described in this paper: https://ccrma.stanford.edu/~jos/resample/resample.pdf

I managed to construct the $$h$$ (windowed sinc) table, and I have implemented the loop. I am having a hard time to compute the $$\eta$$ interpolation factor and the $$l$$ index index for the $$h$$ table. The "fixed point" implementation of the paper is a bit confusing to me.

Looking at figure 7 in the paper, I think the interpolation factor needs to be computed between the $$h$$ entries instead of between the original samples. However I don't understand how to compute the starting point ($$h(0)$$). As in the paper, the table only contains the half wing.

I have just started working with audio code, so I might be misunderstanding something fundamental about the paper.

Here's the main code:

Fs = 44100
Fs_prime = 96000

Nz = 13
L = 128
M = (Nz * 2 + 1) * L

def resample(original, windowed_sinc):
ratio = Fs_prime / Fs
t_increment = 1 / Fs_prime
inteval_orig = 1 / Fs
Lo = len(original)
l_half = M // 2

y = np.zeros(int(ratio * Lo + 0.5))

t = 0
l = 0
for n in range(y.shape[0]):
n0 = int(t / inteval_orig)
original_t = n0 * inteval_orig
eta0 = 1 - (t - original_t ) / inteval_orig
eta1 = 1 - eta0

v = 0
eta = eta0

for i in range(Nz):
ni = n0 - i
if ni < 0 or ni >= Lo:
continue

lL = l + (i * L)
print(lL)

if lL + 1 < 0 or lL + 1 >= l_half:
break

v += original[ni] * ( w[lL] + eta * (w[lL + 1] - w[lL]) )

eta = eta1
for i in range(Nz):
ni = n0 + 1 + i
if ni < 0 or ni >= Lo:
continue

lL = l + (1 + i * L)
print(lL)

if lL + 1 >= l_half:
break

v += original[ni] * ( w[lL] + eta * (w[lL + 1] - w[lL]) )

l = (l + n0) % l_half

y[n] = v

t += t_increment

return y


You have to first decide how many equally-spaced fractional-sample delays you want. You may have to linearly interpolate between adjacent fractional-sample delays. Call that number $$R$$. (For extremely high-quality audio resampling, I have found $$R=512$$ is good. But fewer discrete delays may be acceptable. Like $$R=128$$.)

Then you have to decide how many non-zero taps you use to compute your fractional-delayed sample. It should be an even number with an equal number of samples preceding the interpolated sample as the number of samples following the interpolated sample. Call this number $$N$$. I have found that $$N=16$$ to be good enough for waveform synthesis, but for high-quality audio, I have gone to $$N=32$$ (so the 32 original samples surrounding the interpolated sample on both sides are used).

Now the number of non-zero coefficients of your ideal brickwall LPF FIR is $$N \times R - 1$$ (the FIR order is $$NR-2$$) and the and the target cutoff frequency for the brickwall LPF is

$$f_\mathrm{cutoff}=\frac{f_\mathrm{Nyquist}}{R}$$

So, first design the impulse response for your FIR filter (sinc or Parks-McClellan or Least-squares). Append a zero preceding the $$NR-1$$ samples in the FIR. Now you have exactly $$NR$$ samples.

Now, you have to reshape (that's the MATLAB term) this $$1 \times NR$$ column matrix into an $$R \times N$$ matrix. $$R$$ columns each having $$N$$ samples. Those are the impulse responses for each of the $$R$$ fractional-sample delays.

• Thanks for your reply! And thanks also for dsp.stackexchange.com/a/37715/71569, I used it to implement the sinc window. I now understand the indexing of $h$ in the original paper, if you think of it in terms of a matrix. I guess I still don't get how to determine how to pick the right column. If I have $t_{new}$ that sits between $old_{t0}$ and $old_{t1}$, how do I compute $l$ to index into the matrix? Commented Mar 18 at 10:16
• Consider the "equally-spaced fractional-sample delays" as hypothetical samples. If $R=128$ it's like you upsampled by a factor of 128, but you haven't bothered to compute all 128 new samples for each input sample at the original sample rate. Now the fraction that $t_{new}$ sits between $old_{t0}$ and $old_{t1}$ is the same fraction of the 128 samples. You only need to compute two adjacent samples. Then linearly interpolate between those two. , Commented Mar 19 at 3:51