Creating and using windowed sinc filters (low-pass and high pass) from scratch

windowed_sinc_filter.py

import math
import numpy as np
import matplotlib.pyplot as plt

# Generate a lowpass filter using a windowed sinc function
def sinc(x):
    if x == 0:
        return 1
    return math.sin(math.pi * x) / (math.pi * x)

Fs = 2000 # Hz
Ts = 1 / Fs

Fc = (Fs * 0.1) # Hz, not normalised
Bt = (Fs / 40) # Hz, transition bandwidth, not normalised

# Calculate filter length based on the transition bandwidth (b = 4/N)
N = int(3.1 / (Bt / Fs)) # Note: This is independent of the actual sampling frequency, since Bt is used normalised
N = N if (N % 2 == 0) else N + 1

print(f"Generating a {N} tap filter")

# Compute the filter
lpf = []
for i in range(N):
    # Note: This filter is somewhat independent of the sampling frequency, since Fc is normalised
    #       This means the filter will be valid for any sampling frequency where the cutoff is *proportionally* the same
    lpf.append(sinc(2 * (Fc / Fs) * (i - (N - 1) // 2)))

# Compute the window function
window = []
for i in range(N):
    # Regular hamming window
    window.append(0.54 - (1 - 0.54) * math.cos((2 * math.pi * i) / (N - 1)))

# Apply the window to the sinc filter
for i, w in enumerate(window):
    lpf[i] *= w

# Normalise the filter coefficients to get a final gain of 1
coefficients_total = sum(lpf)
for i in range(N):
    lpf[i] /= coefficients_total

# Convert to a highpass filter with the same cutoff, by spectral inversion
hpf = []
for v in lpf:
    hpf.append(-v)
hpf[(N - 1) // 2] += 1

# Generate a signal with low and high frequency components
signal = []
signal_n = N * 2
for i in range(signal_n):
    components = 0
    components += math.sin(2 * math.pi * 40 * i * Ts)
    components += math.sin(2 * math.pi * 500 * i * Ts)
    signal.append(components)

# Apply the filter on an entire signal, end to end
output_lp = np.convolve(signal, lpf, mode="same")
output_hp = np.convolve(signal, hpf, mode="same")

# Apply the filter on a "streaming" buffer, sample by sample
output_lp2 = []
for i in range(N):
    # We need at least one full filter lengths worth of samples to use, so consider our "buffer index" to be N-1
    # Since we're applying this across a bunch of samples, add an offset of i each time
    buffer_index = N-1 + i

    # Convolving by hand, start off with a result of 0
    sample = 0
    for n in range(len(lpf)):
        sample += lpf[n] * signal[buffer_index - n]
    output_lp2.append(sample)

fig, ax = plt.subplots()
ax.plot(signal)
ax.plot(output_lp2)
# ax.plot(output_hp)
# ax.plot(window)
# ax.plot(lpf)
ax.grid()
plt.show()