240 lines
8.2 KiB
Python
240 lines
8.2 KiB
Python
import numpy as np
|
|
import matplotlib.pyplot as plt
|
|
import scipy.signal as sig
|
|
|
|
# optionally, use an alternate style for plotting.
|
|
# this reconfigures colors, fonts, padding, etc.
|
|
# i personally prefer the look of ggplot, but the contrast is generally worse.
|
|
plt.style.use('ggplot')
|
|
|
|
# create log-spaced points from 20 to 20480 Hz.
|
|
# we'll compute our y values using these later.
|
|
# increasing precision will result in a smoother plot.
|
|
precision = 4096
|
|
xs = np.arange(0, precision) / precision
|
|
xs = 20 * 1024**xs
|
|
|
|
# decide on a sample-rate and a center-frequency (the -3dB point).
|
|
fs = 44100
|
|
fc = 1000
|
|
|
|
|
|
def response_setup(ax, ymin=-24, ymax=24, major_ticks_every=6, minor_ticks_split=2):
|
|
"""configure axes for magnitude plotting within the range of human hearing."""
|
|
ax.set_xlim(20, 20000)
|
|
ax.set_ylim(ymin, ymax)
|
|
ax.grid(True, 'both') # enable major and minor gridlines on both axes
|
|
ax.set_xlabel('frequency (Hz)')
|
|
ax.set_ylabel('magnitude (dB)')
|
|
# manually set y tick positions.
|
|
# this is usually better than relying on matplotlib to choose good values.
|
|
from matplotlib import ticker
|
|
ax.set_yticks(tuple(range(ymin, ymax + 1, major_ticks_every)))
|
|
minor_ticker = ticker.AutoMinorLocator(minor_ticks_split)
|
|
ax.yaxis.set_minor_locator(minor_ticker)
|
|
|
|
|
|
def modified_bilinear(b, a, w0):
|
|
"""
|
|
the modified bilinear transform defers specifying the center frequency
|
|
of an analog filter to the conversion to a digital filter:
|
|
this s-plane to z-plane transformation.
|
|
it also compensates for the frequency-warping effect.
|
|
|
|
b and a are the normalized numerator and denominator coefficients (respectively)
|
|
of the polynomial in the s-plane, and w0 is the angular center frequency.
|
|
effectively, w0 = 2 * pi * fc / fs.
|
|
|
|
the modified bilinear transform is specified by RBJ [1] as:
|
|
s := (1 - z^-1) / (1 + z^-1) / tan(w0 / 2)
|
|
as opposed to the traditional bilinear transform:
|
|
s := (1 - z^-1) / (1 + z^-1)
|
|
|
|
[1]: http://musicdsp.org/files/Audio-EQ-Cookbook.txt
|
|
"""
|
|
|
|
# in my opinion, this is a more convenient way of creating digital filters.
|
|
# scipy doesn't implement the modified bilinear transform,
|
|
# so i'll just write the expanded form for second-order filters here.
|
|
# TODO: include the equation that expanded/simplified to this
|
|
|
|
# use padding to allow for first-order filters.
|
|
# untested
|
|
if len(b) == 2:
|
|
b = (b[0], b[1], 0)
|
|
if len(a) == 2:
|
|
a = (a[0], a[1], 0)
|
|
|
|
assert len(b) == 3 and len(a) == 3, "unsupported order of filter"
|
|
|
|
cw = np.cos(w0)
|
|
sw = np.sin(w0)
|
|
zb = np.array((
|
|
b[0]*(1 - cw) + b[2]*(1 + cw) - b[1]*sw,
|
|
2*(b[0]*(1 - cw) - b[2]*(1 + cw)),
|
|
b[0]*(1 - cw) + b[2]*(1 + cw) + b[1]*sw,
|
|
))
|
|
za = np.array((
|
|
a[0]*(1 - cw) + a[2]*(1 + cw) - a[1]*sw,
|
|
2*(a[0]*(1 - cw) - a[2]*(1 + cw)),
|
|
a[0]*(1 - cw) + a[2]*(1 + cw) + a[1]*sw,
|
|
))
|
|
return zb, za
|
|
|
|
|
|
def digital_magnitude(xs, cascade, fc, fs, decibels=True):
|
|
"""
|
|
compute the digital magnitude response
|
|
of an analog SOS filter (cascade)
|
|
at the points given at xs
|
|
and the center frequency fc
|
|
for the sampling rate of fs.
|
|
"""
|
|
|
|
to_magnitude_decibels = lambda x: np.real(np.log10(np.abs(x)**2) * 10)
|
|
|
|
# compute angular frequencies
|
|
ws = 2 * np.pi * xs / fs
|
|
w0 = 2 * np.pi * fc / fs
|
|
|
|
digital_cascade = [modified_bilinear(*ba, w0=w0) for ba in cascade]
|
|
|
|
# we don't need the first result freqz returns
|
|
# because it'll just be the worN we're passing.
|
|
responses = [sig.freqz(*ba, worN=ws)[1] for ba in digital_cascade]
|
|
|
|
if decibels:
|
|
mags = [to_magnitude_decibels(response) for response in responses]
|
|
mags = np.array(mags)
|
|
composite = np.sum(mags, axis=0)
|
|
else:
|
|
# untested
|
|
mags = [np.abs(response) for response in responses]
|
|
mags = np.array(mags)
|
|
composite = np.prod(mags, axis=0)
|
|
|
|
# we could also compute phase using (untested):
|
|
# -np.arctan2(np.imag(response), np.real(response))
|
|
|
|
return composite, mags
|
|
|
|
|
|
# compute the butterworth coefficients for later.
|
|
# we'll request them in SOS format: second-order sections.
|
|
# that means we'll get an array of second-order filters
|
|
# that combine to produce one higher-order filter.
|
|
# note that we specify 1 as the frequency; we'll give the actual frequency
|
|
# later, when we transform it to a digital filter.
|
|
butterworth_2 = sig.butter(N=2, Wn=1, analog=True, output='sos')
|
|
butterworth_6 = sig.butter(N=6, Wn=1, analog=True, output='sos')
|
|
# scipy returns the b and a coefficients in a flat array,
|
|
# but we need them separate.
|
|
butterworth_2 = butterworth_2.reshape(-1, 2, 3)
|
|
butterworth_6 = butterworth_6.reshape(-1, 2, 3)
|
|
butterworth_2_times3 = butterworth_2.repeat(3, axis=0)
|
|
|
|
|
|
# set up our first plot.
|
|
fig, ax = plt.subplots()
|
|
response_setup(ax, -30, 3)
|
|
ax.set_title("comparison of digital butterworth filters (fc=1kHz, fs=44.1kHz)")
|
|
|
|
# compute the magnitude of the filters at our x positions.
|
|
# we don't want the individual magnitudes here,
|
|
# so assign them to the dummy variable _.
|
|
ys1, _ = digital_magnitude(xs, butterworth_2, fc, fs)
|
|
ys2, _ = digital_magnitude(xs, butterworth_2_times3, fc, fs)
|
|
ys3, _ = digital_magnitude(xs, butterworth_6, fc, fs)
|
|
|
|
# our x axis is frequency, so it should be spaced logarithmically.
|
|
# our y axis is decibels, which is already a logarithmic unit.
|
|
# thus, semilogx is the plotting function to use here.
|
|
ax.semilogx(xs, ys1, label="second order")
|
|
ax.semilogx(xs, ys2, label="second order, three times")
|
|
ax.semilogx(xs, ys3, label="sixth order")
|
|
|
|
ax.legend() # display the legend, given by labels when plotting.
|
|
fig.show()
|
|
fig.savefig('butterworth_comparison.png')
|
|
|
|
|
|
# set up our second plot.
|
|
fig, ax = plt.subplots()
|
|
response_setup(ax, -12, 12)
|
|
ax.set_title("digital butterworth decomposition (fc=1kHz, fs=44.1kHz)")
|
|
|
|
cascade = butterworth_6
|
|
composite, mags = digital_magnitude(xs, cascade, fc, fs)
|
|
for i, mag in enumerate(mags):
|
|
# Q has an inverse relation to the 1st-degree coefficient
|
|
# of the analog filter's denominator, so we can compute it that way.
|
|
Q = 1 / cascade[i][1][1]
|
|
ax.semilogx(xs, mag, label="second order (Q = {:.6f})".format(Q))
|
|
ax.semilogx(xs, composite, label="composite")
|
|
|
|
ax.legend()
|
|
fig.show()
|
|
fig.savefig('butterworth_composite.png')
|
|
|
|
|
|
# as a bonus, here's how the modified bilinear transform is especially useful:
|
|
# we can specify all the basic second-order analog filter types
|
|
# with minimal code!
|
|
|
|
lowpass2 = lambda A, Q: ((1, 0, 0),
|
|
(1, 1/Q, 1))
|
|
highpass2 = lambda A, Q: ((0, 0, 1),
|
|
(1, 1/Q, 1))
|
|
|
|
allpass2 = lambda A, Q: ((1, 1/Q, 1),
|
|
(1, 1/Q, 1))
|
|
|
|
# peaking filters are also (perhaps for the better) known as bell filters.
|
|
peaking2 = lambda A, Q: ((1, A/Q, 1),
|
|
(1, 1/A/Q, 1))
|
|
# TODO: add classical "analog" peaking filter, where bandwidth and amplitude are interlinked
|
|
|
|
# there are two common bandpass variants:
|
|
# one where the bandwidth and amplitude are interlinked,
|
|
bandpass2a = lambda A, Q: ((0, -A/Q, 0),
|
|
(1, 1/A/Q, 1))
|
|
# and one where they aren't.
|
|
bandpass2b = lambda A, Q: ((0,-A*A/Q, 0),
|
|
(1, 1/Q, 1))
|
|
|
|
notch2 = lambda A, Q: ((1, 0, 1),
|
|
(1, 1/Q, 1))
|
|
|
|
lowshelf2 = lambda A, Q: (( A, np.sqrt(A)/Q, 1),
|
|
(1/A, 1/np.sqrt(A)/Q, 1))
|
|
|
|
highshelf2 = lambda A, Q: ((1, np.sqrt(A)/Q, A),
|
|
(1, 1/np.sqrt(A)/Q, 1/A))
|
|
|
|
|
|
# here's an example of how to plot these.
|
|
# we specify the gain in decibels and bandwidth in octaves,
|
|
# then convert these to A and Q.
|
|
fig, ax = plt.subplots()
|
|
response_setup(ax, -30, 30)
|
|
|
|
invsqrt2 = 1 / np.sqrt(2) # for bandwidth <-> Q calculation
|
|
# note that A is squared, so the decibel math here is a little different.
|
|
designer = lambda f, decibels, bandwidth: f(10**(decibels / 40.), invsqrt2 / bandwidth)
|
|
|
|
ax.set_title("example of interlinked bandpass amplitude/bandwidth")
|
|
cascade = [
|
|
designer(bandpass2a, 0, 2),
|
|
designer(bandpass2a, -18, 2),
|
|
designer(bandpass2a, 18, 2),
|
|
]
|
|
|
|
composite, mags = digital_magnitude(xs, cascade, fc, fs)
|
|
for mag in mags:
|
|
ax.semilogx(xs, mag)
|
|
#ax.semilogx(xs, composite)
|
|
|
|
#ax.legend()
|
|
fig.show()
|
|
fig.savefig('filter_example.png')
|