```
import numpy as np
from scipy import signal
import matplotlib.pyplot as plt
# Create low-pass filter (8th order Butterworth, f_cutoff = 20 kHz)
fs = 192000
N = 8
f_cutoff = 20000
butter_sos = signal.butter(N=N, Wn=f_cutoff, btype='low', output='sos', fs=fs)
# 2 kHz square wave band-limited to its first 5 terms
f_signal = 2000
duration = 1.0
omega = 2 * np.pi * f_signal
t = np.linspace(0.0, duration, int(duration * fs), endpoint=False)
n_pts = t.shape[0]
x = np.empty([5, t.size])
for idx, k in enumerate([1, 3, 5, 7, 9]):
x[idx, :] = (4 / np.pi) * np.sin(k * omega * t) / k
x_sum = np.sum(x, axis=0)
# Band-limited square wave filtered by the low-pass filter
x_butter = signal.sosfilt(butter_sos, x_sum)
# Plot original and filtered square waves
plot_pts = np.arange(int(3 * fs / f_signal)) + int(3 * fs / f_signal)
fig2, ax = plt.subplots(figsize=(8, 5))
ax.plot(t[plot_pts], x_sum[plot_pts], label='$x_{orig}$')
ax.plot(t[plot_pts], x_butter[plot_pts], label='$x_{butter}$')
ax.legend(loc='upper right')
plt.tight_layout()
plt.show()
```