Help with resampling algorithm implementation

Question

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 &lt; 0 or ni &gt;= Lo:
            continue

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

        if lL + 1 &lt; 0 or lL + 1 &gt;= 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 &lt; 0 or ni &gt;= Lo:
            continue

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

        if lL + 1 &gt;= 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

score 2 · Answer 1 · answered Mar 17 '24 at 22:40

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 https://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? — Marco Castorina, Mar 18 '24 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. , — robert bristow-johnson, Mar 19 '24 at 03:51

Help with resampling algorithm implementation

1 Answers1