The book "Digital Signal Processing - Fundamentals and applications" (Li Tan, Jean Jiang, 2nd Edition) contains an example in chapter 8, section 5 of a digital audio equaliser. Unfortunately, this example is written in Matlab; here the same concepts are implemented in a Python script, using scipy.signal. The script uses cascaded second-order section (SOS) filters; a SOS filter is made by taking a high-order IIR filter and splitting it into several 2nd-order IIR filters (AKA digital biquad filters) which are applied to an input signal sequentially. The motivation for using SOS filters is that high-order IIR filters can contain both very large and very small filter coefficients, and using such a filter can be numerically unstable when using floating-point arithmetic. Below is the Python script which implements a simple digital audio equaliser, followed by the resulting graphs:
import numpy as np
import matplotlib.pyplot as plt
from scipy import signal
def get_sos_filter_coeffs(f_c_list, bw_list, f_n, butter_order=1, verbose=True):
""" get_sos_filter_coeffs: Get a list of SOS (second-order section) coefficients
for a bank of band-pass Butterworth filters. """
# Create list containing the coefficients for each filter
filter_coefs = []
for f_c, bw in zip(f_c_list, bw_list):
# Half-power frequencies (f_c = sqrt(hpf_lo * hpf_hi))
hpf_lo = np.sqrt(f_c*f_c + bw*bw/4) - bw/2
hpf_hi = np.sqrt(f_c*f_c + bw*bw/4) + bw/2
# Generate butterworth filter
sos = signal.butter(
butter_order, [hpf_lo/f_n, hpf_hi/f_n], btype="bandpass",
output="sos"
)
filter_coefs.append(sos)
if verbose: print(hpf_lo, hpf_hi, sos, "-"*50, sep="\n")
return filter_coefs
def equalise_signal(x, filter_coefs, gains):
""" equalise_signal: apply individual gains to different frequency bands within
a signal. """
y = x.copy()
for sos, g in zip(filter_coefs, gains):
y += signal.sosfilt(sos, x) * g
return y
def plot_filter_coeffs(
filter_coefs, f_plot, f_s, filename="Filter responses",
xlims=[1e1, 1e5], ylims=[1e-5, 1e0], plot_func=plt.loglog
):
""" plot_filter_coeffs: plot the frequency response for a list of filter
coefficients, and save the image file to disk """
# Check that plot_func is valid
valid_plot_funcs = [plt.plot, plt.semilogx, plt.semilogy, plt.loglog]
if plot_func not in valid_plot_funcs:
raise ValueError("plot func must be in {}".format(
[("plt." + pf.__name__) for pf in valid_plot_funcs]
))
# Initialise plot
plt.figure(figsize=[8, 6])
colours = plt.get_cmap("hsv")(
np.linspace(0, 1, len(filter_coefs), endpoint=False)
)
# Plot filter responses
for i, sos in enumerate(filter_coefs):
_, h = signal.sosfreqz(sos, worN=f_plot, fs=f_s)
plot_func(f_plot, abs(h), c=colours[i])
# Format, save and close
plt.title("Filter responses")
plt.xlabel("Frequency (Hz)")
plt.ylabel("Gain")
plt.grid(which="major", ls="-")
plt.grid(which="minor", ls=":")
plt.tight_layout()
plt.xlim(*xlims)
plt.ylim(*ylims)
plt.savefig(filename)
plt.close()
def plot_equalised_signal(
x, y, f_s, filename="Signal spectra", xlims=[1e1, 1e5], ylims=[1e-4, 2e1],
plot_func_name="loglog"
):
""" plot_equalised_signal: plot the frequency spectrum of a signal before and
after equalisation, and save the image file to disk """
# Check that plot_func is valid
valid_plot_func_names = ["plot", "semilogx", "semilogy", "loglog"]
if plot_func_name not in valid_plot_func_names:
raise ValueError("plot func must be in " + str(valid_plot_func_names))
# Generate frequencies and frequency response
N = x.size
f = np.linspace(0, f_s, N)
Ax = 2 * np.abs(np.fft.fft(x)) / N
Ay = 2 * np.abs(np.fft.fft(y)) / N
fig, axes = plt.subplots(2, 1, sharex=True)
fig.set_size_inches(8, 6)
# Only plot frequencies less than or equal to the Nyquist frequency
plot_inds = (f <= (f_s / 2))
getattr(axes[0], plot_func_name)(f[plot_inds], Ax[plot_inds], "b")
axes[0].set_title("Original signal spectrum")
getattr(axes[1], plot_func_name)(f[plot_inds], Ay[plot_inds], "b")
axes[1].set_title("Equalised signal spectrum")
# Format, save and close
axes[1].set(xlabel="Frequency (Hz)", xlim=xlims)
for ax in axes:
ax.set(ylabel="Gain", ylim=ylims)
ax.grid(which="major", ls="-")
ax.grid(which="minor", ls=":")
fig.tight_layout()
plt.savefig(filename)
plt.close()
if __name__ == "__main__":
# Sampling and Nyquist frequencies
f_s = 44.1e3
f_n = f_s / 2
# Define the centre frequency and bandwidth for each filter
f_c_list = [100, 200, 400, 1000, 2500, 6000, 15000]
bw_list = [50, 100, 200, 500, 1250, 3000, 7500]
# Generate band-pass filters and plot the responses
filter_coefs = get_sos_filter_coeffs(
f_c_list, bw_list, f_n, butter_order=4
)
f_plot = np.logspace(1, np.log10(f_n), 500)
plot_filter_coeffs(filter_coefs, f_plot, f_s)
# Generate input signal to equalise
t = np.arange(20481) / f_s
x = np.zeros(t.shape)
periods = np.arange(7) * np.pi / 14
for f, p in zip(f_c_list, periods):
x += np.sin(2*np.pi*f*t + p)
# Equalise the signal and plot the spectrum
gains = [10, 10, 0, 0, 0, 10, 10]
y = equalise_signal(x, filter_coefs, gains)
plot_equalised_signal(x, y, f_s)
