Skip to content

Instantly share code, notes, and snippets.

@theWatchmen
Created April 13, 2024 13:40
Show Gist options
  • Save theWatchmen/746f35c349748525b412cfd9466608ce to your computer and use it in GitHub Desktop.
Save theWatchmen/746f35c349748525b412cfd9466608ce to your computer and use it in GitHub Desktop.
Resampling impementation in Python (2nd attempt)
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()
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment