Resampling impementation in Python (2nd attempt)

resampling2.py

import numpy as np
import matplotlib.pyplot as plt

Fs = 44100
Fs_prime = 96000
A = 90
TEST_SAMPLES = 1000

Nz = 13
Nzd = Nz * 2 + 1
# Cutoff frequency = pi / L
Nn = 7
Neta = 6
Np = Nn + Neta
L = 1 << Nn
M = Nzd * L

K = 32

def beta():
    if A > 50:
        return 0.1102 * (A - 8.7)
    elif A > 21:
        return 0.5842 * (A - 21)**(2./5.) + 0.07886 * (A - 21)
    else:
        return 0.0

# From https://dsp.stackexchange.com/questions/37714/kaiser-window-approximation
def bessel_coeff():
    result = np.zeros(K)
    fact = 1

    for k in range(K):
        fact *= (k + 1)
        exp = 2 ** k
        result[k] = 1 / (fact * exp)**2

    return result

def bessel(x, bc):
    k = K - 1
    x2 = x**2
    y = x2 * bc[k]
    for k in range(K - 2, 0, -1):
        y = x2 * (bc[k] + y)

    return y

def window(bc):
    b = beta()

    b_beta = bessel(b, bc)

    m2 =  M//2
    y = np.zeros(M)

    for n in range(m2 + 1):
        bn = b * ( 1 - ( n / m2 )**2 )**0.5
        bb = bessel(bn, bc)
        y[n + m2 - 1] = (1 / b_beta) * bb
        y[m2 - n] = y[n + m2 - 1]

    return y

def resample_window_sinc(o, r, w):
    rho = Fs_prime / Fs
    input_increment = 1 / rho

    y = np.zeros(r.shape[0])

    o_padded = np.concatenate([np.zeros(Nz), o, np.zeros(Nz)])
    Lo = len(o_padded)

    for n in range(y.shape[0]):
        input_delay = Nz + n * input_increment
        n0 = int(input_delay)

        delay_start = int(input_delay - Nz)
        additional_delay = input_delay - delay_start
        integral_delay = int(additional_delay)
        fractional_delay = additional_delay - integral_delay
        l_offset = fractional_delay * L
        l0 = int(l_offset)

        eta0 = 1 - (l_offset - int(l_offset))
        eta1 = 1 - eta0

        v = 0
        nn = (l0) % L
        for z in range(-Nz, Nz + 1):
            ni = n0 + z
            if ni < 0 or ni >= Lo:
                continue
            v += o_padded[ni] * eta0 * w[nn][z + Nz]

        nn = (l0 + 1) % L
        for z in range(-Nz, Nz + 1):
            ni = n0 + 1 + z
            if ni < 0 or ni >= Lo:
                continue

            v += o_padded[ni] * eta1 * w[nn][z + Nz]

        y[n] = v

    return y

def resample_box(o, r):
    t_increment = 1 / Fs_prime
    inteval_orig = 1 / Fs

    y = np.zeros(r.shape[0])

    t = 0
    for n in range(y.shape[0]):
        n0 = int(t / inteval_orig)
        ot = n0 * inteval_orig
        eta = (t - ot) / inteval_orig

        if eta <= 0.5:
            y[n] = o[n0]
        else:
            new_n = min(len(o) - 1, n0 + 1)
            y[n] = o[new_n]

        t += t_increment

    return y

if __name__ == '__main__':
    bc = bessel_coeff()

    x = np.array(np.arange(-16, 16, 1/M))
    y = np.zeros(x.shape)

    for i in range(y.shape[0]):
        y[i] = bessel(x[i], bc)

    w = window(bc)

    scale = min(1, Fs_prime/Fs)
    omega = min(Fs,Fs_prime)
    fc = 0.5 / L
    sinc = [np.sin(2*np.pi * fc * i) / (2*np.pi * fc * i) if i != 0 else 1 for i in range(-M//2, M//2)]
    sinc = [w[i] * sinc[i] for i in range(M)]
    sinc = [0]*34 + sinc # TODO(marco): compute correct number of zeros
    sinc = np.array(sinc)

    sincr = np.zeros((L, Nzd))
    for l in range(L):
        for i in range(Nzd):
            sincr[l][i] = sinc[L*i + l]

    o_count = TEST_SAMPLES
    xo = np.array(np.arange(0, 1, 1/o_count))
    original = np.sin(400 * np.pi * 2 * xo)

    r_count = TEST_SAMPLES * Fs_prime/Fs
    xr = np.array(np.arange(0, 1, 1/r_count))
    r = resample_window_sinc(original, xr, sincr)

    fft_o = np.fft.fft(original)
    fft_r = np.fft.fft(r)

    fig, (ax1, ax2, ax3, ax4) = plt.subplots(4, 1, layout='constrained')
    ax1.set_title('Individual Filters')
    for i in range(sincr.shape[0]):
        ax1.plot(sincr[i])

    ax2.set_title('Original')
    ax2.plot(original[:TEST_SAMPLES // 4])

    ax3.set_title('Output')
    ax3.plot(r[:r.shape[0] // 4])

    ax4.set_title('FFT Output')
    ax4.plot(np.abs(fft_r[:fft_r.shape[0]//2]))
    ax4.plot(np.abs(fft_o[:fft_o.shape[0]//2]))

    plt.show()