"importing Jupyter notebook from FilterUtils.ipynb\n"
]
}
],
"source": [
"%matplotlib inline\n",
"import matplotlib\n",
"import matplotlib.pyplot as plt\n",
"import numpy as np\n",
"from IPython.display import Audio\n",
"from scipy import signal\n",
"\n",
"import import_ipynb\n",
"from FilterUtils import *"
]
},
{
"cell_type": "code",
- "execution_count": 2,
+ "execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"plt.rcParams['figure.figsize'] = 14, 4 \n",
"matplotlib.rcParams.update({'font.size': 14})"
]
},
{
"cell_type": "code",
- "execution_count": 3,
+ "execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"DEFAULT_SF = 16000"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Introduction\n",
"\n",
"In this notebook we will explore a complete set of \"recipes\" to design second-order digital IIR filters. The transfer function of a generic second-order section (also known as a **biquad**) has the canonical form\n",
"The routines defined in the rest of this notebook will allow you to compute the values of the five biquad parameters in order to implement a variety of different filter prototypes according to the desired specifications. We will also explore how to cascade second-order sections to implement higher-order filters with improved characteristics. "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Common practices\n",
"\n",
"Although ultimately we will design digital filters, in audio applications it is important to become familiar with the main ideas behind _analog_ filter design and analysis; indeed, audio recording and production techniques have been developed and fine-tuned well before the advent of DSP and, even in today's world of DAWs, the language and many of the practices in current use still reflect the analog conventions of yore.\n",
"\n",
"In particular:\n",
" * filter specifications are almost always expressed in terms of real-world frequencies in Hz rather than as normalized frequencies over $[-\\pi, \\pi]$; this of course implies that the underlying sampling frequency is known\n",
" * plots of the magnitude response will usually be shown on a log-log graph, generally using a decibel scale for the amplitde and a decade scale for frequencies; this mirrors the way in which human audio perception is approximately logarithmic both in frequency and in scale.\n",
" \n",
"The companion notebook ``FilterUtils.ipynb`` implements a set of plotting routines that take these conventions into account. For example, this is how the magnitude response of the same leaky integrator looks like in the typical representations used in communication systems vs. audio equalization:"
"The next section provides a set of functions to compute the five biquad filter coefficients for a desired response (lowpass, bandpass, etc) and an associated set of specifications (cutoff, attenuation, etc). Rather than tackling the search for the coefficients as an abstract optimization problem, each recipe starts from a well-known _analog_ second-order filter with the desired characteristic, and then converts it to a discrete-time filter using a mapping called the _bilinear transform_ .\n",
"\n",
"The reason for this approach is that the classic analog filter prototypes (Butterworth, Chebyshev, etc.) (and the topologies used for their implementation) are extremely well understood and have proven extremely reliable over more than a century of research and practical experimentation. "
"Historically, the development of electronic filters began with the design of **passive** analog filters, that is, filters using only resistors, capacitors and inductors; indeed, an RLC circuit, that is, a network containing a resistor, a capacitor and an inductor, can implement the prototypical analog second-order section. Since these filters have no active elements that can provide signal amplification, the power of their output is at most equal to (but, in practice, smaller than) the power of the input. This implicitly guarantees the stability of these systems although, at least in theory, strong resonances can appear in the frequency response.\n",
"\n",
"Analog filters work by exploiting the frequency-dependent reactance of capacitors and inductors. The input-output characteristic for linear circuits using these electronic components is described by linear differential equations, which implies the existence of some form of _feedback_ in the circuits themselves. As a consequence, when we convert the analog prototypes to digital realizations, we invariably end up with IIR filters. \n",
"\n",
"Passive filters operating in the frequency range of audio applications require the use of bulky inductors and therefore **active** filters are usually preferred in the analog domain. In the digital domain, on the other hand, we are in fact free to use arbitrary gain factors (it's just multiplications!) and so the resulting transfer functions can approximate either type of design. "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## The bilinear transform\n",
"\n",
"The cookbook recipes below are obtained by mapping second-order analog filters prototypes to equivalent digital realization via the _bilinear transform_ . We will not go into full details but, as a quick reference, here are the main ideas behind the method. \n",
"\n",
"An analog filter is described by a transfer function $H(s)$, with $s\\in \\mathbb{C}$, which is the Laplace transform of the filter's continuous-time impulse response. The key facts about $H(s)$ are:\n",
" * filter stability requires that all the poles of $H(s)$ lie in the left half of the complex plane (i.e. their real part must be negative)\n",
" * the filter's frequency response is given by $H(j\\Omega)$, that is, by the values of $H(s)$ along the imaginary axis\n",
"The bilinear transform maps the complex $z$-plane (discrete time) to the complex $s$-plane (continuous time) as\n",
"\n",
"$$\n",
" s \\leftarrow c \\frac{1 - z^{-1}}{1 + z^{-1}} = \\Phi_{c}(z)\n",
"$$\n",
"\n",
"where $c$ is a real-valued constant. Given a stable analog filter $H(s)$, the transfer function of its discrete-time version is $H_d(z) = H(\\Phi_{c}(z))$ and it is relatively easy to verify that\n",
" * the inside of the unit circle on the $z$-plane is mapped to the left half of the $s$-plane, which preserves stability\n",
" * the unit circle on the $z$-plane is mapped to the imaginary axis of the $s$-plane\n",
" \n",
"The last property allows us to determine the frequency response of the digital filter as $H_d(e^{j\\omega}) = H(\\Phi_{c}(e^{j\\omega})) = H(j\\,c\\tan(\\omega/2))$, or: \n",
"we can see that $\\omega=0$ is mapped to $\\Omega=0$, $\\omega=\\pi/2$ is mapped to $\\Omega=c$, and $\\omega=\\pi$ is mapped to $\\Omega=\\infty$, which reveals the high nonlinearity of the frequency mapping. We usually need to precisely control a least one notable point $H_d(e^{j\\omega_0})$ in the frequency response of the discrete-time filter; for example, in a resonator, we need to place the magnitude peak at a specific frequency $\\omega_0$. To achieve this, we design the analog filter so that $H(1j) = H_d(e^{j\\omega_0})$ and then we set $c = 1/\\tan(\\omega_0/2)$ in the bilinear operator; this adjustment, called **pre-warping** , is used in the recipes below. "
]
},
{
"cell_type": "markdown",
- "metadata": {},
+ "metadata": {
+ "tags": []
+ },
"source": [
- "To illustrate the principle, here is a simple transfer function $H(s)$ that provide a Gaussian response along the imaginary axis; we can parametrize the width of the bell and its center position is at $\\Omega=1$ by default:"
+ "To illustrate the principle, here is a simple transfer function $H(s)$ that provide a triangular response centered at $\\Omega=1$; the default width is $1/2$ and the width can be optionally scaled.\n",
+ "\n",
+ "along the imaginary axis; we can parametrize the width of the bell and its center position is at by default:"
"Note that the nonlinear mapping between frequency axes has two consequences on the the discrete-time frequency response:\n",
- " * at low frequencies the response becomes more narrow; this can be compensated for by scaling the analog prototype\n",
- " * at high frequencies, the response is no longer exactly symmetric; this is much harder to compensate for because it would require a different analog design and it is therefore an accepted tradeoff."
+ " * at low and high digital frequencies the response becomes more narrow; this can be compensated for by scaling the analog prototype\n",
+ " * as we move to higher frequencies, the response is less and less symmetric; this is much harder to compensate for because it would require a different analog design and it is therefore an accepted tradeoff.\n",
+ " \n",
+ "The following example tries to keep the width of the response uniform"
+ "**Exercise**: how was the scaling factor derived? Can you improve on it?"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
- "## What about FIRs\n",
+ "Using the bilinear transform with pre-warping, we can move the equivalent discrete-time frequency response over the $[0, \\pi]$ interval"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "# What about FIRs?\n",
+ "\n",
+ "FIR filters are a great tool in digital signal processing; as opposed to IIR (which can be seen as a digital adaptation of electronic filters) FIRs offer:\n",
+ " * unconditional stability \n",
+ " * the possibility of a linear phase response\n",
+ " * a great design algorithm (Parks-McClellan) even for arbitrary responses\n",
+ " \n",
+ "The price for stability and linear phase is a much higher computational cost: for the same specifications, an FIR filter will require up to a hundred times more operations per sample with respect to an IIR implementation. Linear phase, however, is not terribly relevant in audio applications because of the limited phase sensitivity in the human auditory system. On the other hand, especially in real-time applications, the primary goal in audio processing is to minimize the overall processing delay; since linear phase FIRs have a symmetric impulse response, and since a well-performing filter will have a very long impulse response, the associated delay often makes FIRs difficult to use. Even if we give up linear phase and implement an asymmetric, minimum-phase FIR, the computational cost may be too high.\n",
+ "\n",
+ "\n",
+ "There are countless audiophile blogs debating the merits and demerits of FIRs in audio applications. Some of the purported negatives that are often quoted include\n",
+ " * FIRs sound \"cold\"\n",
+ " * linear phase FIRs cause pre-echos (because of their symmetric impulse response) \n",
+ " * minimum-phase FIRs exhibit excessive ringing in the impulse response\n",
"\n",
- "why not FIR? delay. pre-echo\n",
- "demonstrate pre-echo"
+ "It must be said that these artefacts, if they can be noticed at all, are anyway extremely subtle and unlikely to compromise overall sound quality in a significant way. The major obstacle to the use of FIRs remains their inherent processing delay."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# The cookbook\n",
"\n",
"In the following, we define a set of functions that return the five biquad coefficients for the most common types of audio filtering applications. Many of the formulas have been adapted from Robert Bristow-Johnson's famous [cookbook](https://webaudio.github.io/Audio-EQ-Cookbook/audio-eq-cookbook.html). \n",
"\n",
"Each function returns ``b`` and ``a``, two arrays of three floats each containing the coefficients of the transfer function\n",
"A second-order lowpass filter section will have a passband with approximately unit gain (0 dB) and a monotonically decreasing stopband. It is defined by two parameters:\n",
"\n",
" 1. the \"quality factor\" $Q$, which determines the shape of the magnitude response; by default $Q = \\sqrt{1/2}$, which yields a Butterworth characteristic (i.e. a monotonically decreasing response). \n",
" 1. the _corner frequency_ $f_c$ (also called the _cutoff_ frequency); the magnitude response will be equal to the quality factor $Q$ at $f_c$ and will decrease monotonically afterwards. For $Q = \\sqrt{1/2}$, the attenuation at $f_c$ is equal to $20\\log_{10}(\\sqrt{1/2}) \\approx -3$ dB, which yields a Butterworth (maximally flat) characteristic.\n"
"When $Q = 1/\\sqrt{2}$, as we said, the lowpass section corresponds to a Butterworth filter, that is, a filter with a maximally flat passband and a monotonically decreasing stopband. For higher $Q$ values the magnitude response exhibits a peak around $f_c$ which, in the time domain, corresponds to a damped oscillatory impulse response as shown in the following examples; for lower $Q$ values, the roll-off of the magnitude response will be less steep.\n",
"\n",
"While these $Q$ values are clearly not a good choice for a single-stage lowpass, values other than $1/\\sqrt{2}$ become useful when cascading multiple sections, as we will see later."
"A second-order bandpass filter section will have approximately unit gain (0 dB) in the passband and will decrease monotonically to zero in the stopband. It is defined by two parameters:\n",
"\n",
" 1. the center frequency $f_c$, where the gain is unitary\n",
" 1. the bandwidth $b = (f_+ - f_-)$, where $f_- < f_c < f_+$ are the first frequencies, left and right of $f_c$ where the attenuation reaches $-3$ dB. For the reasons explained above, note that the passband is almost but not exactly symmetric around $f_c$, with the asymmetry more pronounced towards the high end of the spectrum."
" frequency_response(b, a, dB=-50, half=True, axis=ax)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Notch\n",
"\n",
"A notch filter is the complementary filter to a resonator; its attenuation reaches $-\\infty$ at $f_c$ and its bandwidth is usually kept very small in order to selectively remove only a given frequency; this is achieved by placing a pair of complex-conjugate zeros _on_ the unit circle and by placing two poles very close to the zeros."
" plt.gcf().get_axes()[0].axvline(x=(2 * np.pi * CENTER / DEFAULT_SF), linewidth=0.5, color='r')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Shelves\n",
"\n",
"Shelving filters are used to amplify either the low or the high end of a signal's spectrum. A high shelf, for instanance, provides an arbitrary gain for high frequencies and has approximately unit gain in the low end of the spectrum. Shelving filters, high or low, are defined by the following parameters:\n",
"\n",
" 1. the desired _shelf gain_ in dB\n",
" 1. the midpoint frequency $f_c$, which corresponds to the frequency in the transition band where the gain reaches half its value.\n",
" 1. the \"quality factor\" $Q$, which determines the steepnes off the transition band; as for lowpass filters, the default value $Q = 1/\\sqrt{2}$ yields the steepest transition band while avoiding resonances.\n",
" \n",
"A common use case for shelving filters is in consumer audio appliances, where the standard \"Bass\" and \"Treble\" tone knobs control the gain of two complementary shelves with fixed midpoint frequency. "
]
},
{
"cell_type": "code",
"execution_count": 23,
"metadata": {},
"outputs": [],
"source": [
"def LSH(fc, gain, sf, Q=(1/np.sqrt(2))):\n",
" \"\"\"Biquad low shelf\"\"\"\n",
" w = 2 * np.pi * fc / sf\n",
" A = 10 ** (gain / 40) \n",
" alpha = np.sin(w) / (2 * Q)\n",
" c = np.cos(w)\n",
" b = np.array([A * ((A + 1) - (A - 1) * c + 2 * np.sqrt(A) * alpha),\n",
" 2 * A * ((A - 1) - (A + 1) * c),\n",
" A * ((A + 1) - (A - 1) * c - 2 * np.sqrt(A) * alpha)])\n",
" a = np.array([(A + 1) + (A - 1) * c + 2 * np.sqrt(A) * alpha,\n",
" -2 * ((A - 1) + (A + 1) * c),\n",
" (A + 1) + (A - 1) * c - 2 * np.sqrt(A) * alpha])\n",
"A peaking equalizer filter is the fundamental ingrediend in multiband parametric equalization. Each filter provides an arbitrary boost or attenuation for a given frequency band centered around a peak freqency and flattens to unit gain elsewhere. The filter is defined by the following parameters:\n",
"\n",
" 1. the desired gain in dB (which can be negative)\n",
" 1. the peak frequency $f_c$, where the desired gain is attained\n",
" 1. the bandwidth of the filter, defined as the interval around $f_c$ where the gain is greater (or smaller, for attenuators) than half the desired gain in dB; for instance, if the desired gain is 40dB, all frequencies within the filter's bandwidth will be boosted by at least 20dB. Note that the bandwdidth is not exactly symmetrical around $f_c$ "
]
},
{
"cell_type": "code",
"execution_count": 29,
"metadata": {},
"outputs": [],
"source": [
"def PEQ(fc, bw, gain, sf):\n",
" \"\"\"Biquad bandpass filter \"\"\"\n",
" w = 2 * np.pi * fc / sf\n",
" A = 10 ** (gain / 40) \n",
" alpha = np.tan(np. pi * bw / sf)\n",
" c = np.cos(w)\n",
" b = np.array([1 + alpha * A, -2 * c, 1 - alpha * A])\n",
" a = np.array([1 + alpha / A, -2 * c, 1 - alpha / A])\n",
" y = signal.lfilter(b, a, np.r_[1, np.zeros(200)])\n",
" plt.plot(y)\n",
" b, a = PEQ(CENTER, BW, -GAIN_DB, DEFAULT_SF)\n",
" y = signal.lfilter(b, a, y)\n",
" plt.plot(y)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Cascades of biquads\n",
"\n",
"The performance of a single biquad filter may not be adequate for a given application: a second-order lowpass filter, for instance, may not provide a sufficent amount of rejection in the stopband because of its rather slow roll-off characteristic; or we may want to design an equalizer with multiple peaks and dips. In all cases, we usually want to implement the final design as a cascade of biquad sections, because of their inherent numerical robustness."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Factorization of higher-order filters\n",
"\n",
"The first solution if a filter does not meet the requires specifications is to design a higher-order filter, possibly using different filter \"recipes\"; in the case of bandpass filters, for instance, we could try a Chebyshev or elliptic design. The resulting high-order transfer function can be then factored into a cascade of second-order sections (or, in the case of an odd-order filter, a cascade of second order-sections followed by a first-order filter):\n",
"The biquad elements returned by the factorization are not related to the \"cookbook\" prototypes of the previous section and therefore this method is simply an implementation strategy that focuses on second-order structures; the design algorithm, in other words, will be dependent on the particular type of filter. Nevertheless, both the design and the factorization are usually available in numerical packages such as Scipy, for instance and, in the following example, we illustrate the difference between a 6th-order elliptic lowpass and a single second-order butterworth. First we will use the full high-order realization and then we will show how a cascade of three second-order sections implements the same characteristic.\n",
"\n",
"Note that, when cascading transfer functions, the equivalent higher-order filter coefficients can be obtained simply by polynomial multiplication."
" # this returns an array of second-order filter coefficients. Each row corresponds to a section, \n",
" # with the first three columns providing the numerator coefficients and the last three providing the denominator\n",
" soe = signal.ellip(6, 1, 40, CUTOFF, fs=DEFAULT_SF, output='sos')\n",
" cb, ca = [1], [1]\n",
" for n in range(0, 3):\n",
" b, a = soe[n][0:3], soe[n][3:6]\n",
" analog_response(b, a, DEFAULT_SF, dB=-60, axis=ax, color=f'C{n}:')\n",
" cb = np.polymul(b, cb)\n",
" ca = np.polymul(a, ca)\n",
" analog_response(cb, ca, DEFAULT_SF, dB=-60, axis=ax, color='C3')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Cascading lowpass and highpass biquads\n",
"\n",
"A cascade of $N$ identical sections with transfer function $H(z)$ will yield the overall transfer function $H_c(z) = H^N(z)$ and thus the stopband attenuation in decibels will increase $N$-fold. For instance, the following example shows the cumulative magnitude responses obtained by cascading up to five identical second-order Butterworth lowpass sections:"
"As shown by the previous plot, a cascade of identical maximally flat lowpass sections yields a steeper roll-off and preserves the monotonicity of the response. However, since the passband of each filter is not perfectly flat, the $-3~\\mathrm{dB}$ cutoff frequency of the cascade becomes smaller with each added section and the effective bandwidth of the filter is reduced. In the previous example, the original $-3~\\mathrm{dB}$ cutoff frequency was $f_c = 1000~\\mathrm{Hz}$ but the magnitude response of the cascade at $f_c$ is $-15~\\mathrm{dB}$ whereas the actual $-3~\\mathrm{dB}$ point has shifted close to $600~\\mathrm{Hz}$.\n",
"\n",
"If our goal is to obtain a cascade with a maximally flat (Butterworth) response with a given $f_c$, an obvious approach is simply to factorize the transfer function of a high-order Butterworth as explained in the previous section. There is however a clever and simpler design strategy that is based on the geometric arrangement of the poles of an analog Butterworth filter of order $N$:\n",
" * the $N$ complex-conjugate poles are equally spaced along a circular contour centered on the origin of the $s$-plane\n",
" * the angle between poles is equal to $\\pi/N$\n",
" \n",
"With this, the pole angles in the upper $s$-plane are given by \n",
"Now, a generic second-order analog filter will have a single pair of complex-conjugate poles at $p_{1,2} = \\rho e^{\\pm \\theta}$ on the $s$-plane and, by cascading identical sections, we will only manage to increase the poles' multiplicity but we will not be able to change their position. In order to achieve a Butterworth pole configuration we will thus need to adjust the pole angle for each section; this is a simple task because it turns out that a second-order filter's quality factor $Q$ is related to the pole angle as \n",
"$$\n",
" 1/Q = 2\\cos \\theta\n",
"$$ \n",
"which means that we can choose the suitable pole angle for each section simply by setting $Q_n = 1/(2\\cos \\theta_n)$. We can now design $N/2$ discrete-time biquads with the same $Q_n$ values to obtain the desired result.\n",
"\n",
"Below is the example for a cascade of five lowpass sections (i.e. a 10th-order filter) compared to a single biquad, both with cutoff $f_c = 1000~\\mathrm{Hz}$; notice how the $-3~\\mathrm{dB}$ point has not moved in spite of the much steeper rolloff."
"The resulting digital filter has its poles arranged on a circular contour centered in $z=1$ if the cutoff frequency is less than $\\pi/2$ and centered on $z=-1$ otherwise. "
"Finally, the following plot shows the individual magnitude responses of the five sections. You can observe that the required $Q_n$ values lead to some biquad sections with a clear peak at the cutoff frequency, although the overall response is monotonic:"
" analog_response(cb / 10, ca, DEFAULT_SF, dB=-50, points=10001)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Parametric equalization\n",
"\n",
"Peaking equalizers with distinct bandwidths can be cascaded to obtain an arbitrary equalization curve for the entire range of input frequencies; indeed, this is the technique behind so-called _parametric equalizers_ where a bank of logarithmically spaced peaking eq's with independent gain controls allow the user to easily define a global equalization response."
"NB: these plotting functions are available in the companion notebook ``FilterUtils.ipynb``"
]
},
{
"cell_type": "markdown",
- "metadata": {},
+ "metadata": {
+ "slideshow": {
+ "slide_type": "slide"
+ }
+ },
"source": [
"## The frequency axis\n",
" \n",
" * audio DSP works within a very specific frequency range: 16 Hz to 20 kHz\n",
" * filter specifications use frequency values in Hz (instead of normalized values over $[-\\pi, \\pi]$)\n",
" * the frequency axis is logarithmic, because musical intervals are determined by frequency _ratios_ , not differences\n",
+ " * using $\\log_{10}$ splits the frequency axis into **decades**\n",
" \n",
- "Using real-world frequencies in filter specs implies that the underlying sampling frequency must always be specified"
+ "Using real-world frequencies in digital filter design requires knowing the sampling frequency used in the final application: **filter coefficients depend on sampling frequency**"
"Frequency mapping is $\\Omega \\leftarrow c\\tan(\\omega/2)$\n",
"\n",
"Suppose something important needs to happen at $\\omega_0$ in the digital filter:\n",
"\n",
" * we want $H_d(e^{j\\omega_0}) = v_0$\n",
- " * we design analog prototype so that $H(1j) = v_0$\n",
- " * we choose $c = 1/\\tan(\\omega_0/2)$ "
+ " * we design analog prototype so that $H(j) = v_0$\n",
+ " * we choose $c = 1/\\tan(\\omega_0/2)$ \n",
+ " \n",
+ "\n",
+ "**Caution with bandpass** : we have $H_d(e^{j\\omega_0}) = H(j) = v_0$ but the shape of the response may be very distorted! "
]
},
{
"cell_type": "markdown",
- "metadata": {},
+ "metadata": {
+ "slideshow": {
+ "slide_type": "slide"
+ }
+ },
"source": [
- "To illustrate the principle, here is a simple transfer function $H(s)$ that provide a Gaussian response along the imaginary axis; we can parametrize the width of the bell and its center position is at $\\Omega=1$ by default:"
+ "Example: simple continuous-time frequency response with triangular shape"
- "Note that the nonlinear mapping between frequency axes has two consequences on the the discrete-time frequency response:\n",
- " * at low frequencies the response becomes more narrow; this can be compensated for by scaling the analog prototype\n",
- " * at high frequencies, the response is no longer exactly symmetric; this is much harder to compensate for because it would require a different analog design and it is therefore an accepted tradeoff."
+ "We can try to compensate a bit by changing the width of the analog prototype"
+ "**Exercise**: try to formulate a full derivation for the scaling factor"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "slideshow": {
+ "slide_type": "slide"
+ }
+ },
"source": [
"# What about FIRs ?\n",
"\n",
- "why not FIR? delay. pre-echo\n",
- "demonstrate pre-echo"
+ "FIRs are great:\n",
+ " * you can have linear phase\n",
+ " * always stable\n",
+ " * great design algorithm (Parks-McClellan) even for arbitrary responses\n",
+ " \n",
+ "However:\n",
+ " * we don't care so much for linear phase in audio\n",
+ " * they are computationally much more expensive\n",
+ " * the long impulse response causes significant processing delay: problem in real-time applications"
]
},
{
"cell_type": "markdown",
- "metadata": {},
+ "metadata": {
+ "slideshow": {
+ "slide_type": "slide"
+ }
+ },
"source": [
- "# The cookbook\n",
+ "## Audiophiles' myths about FIRs\n",
+ " \n",
+ " * they sound \"cold\"\n",
+ " * you can hear pre-echos for linear phase FIRs\n",
+ " * excessive ringing in the impulse response for minimum-phase FIRs\n",
+ " * ...\n",
+ " \n",
+ "Honestly, I don't know. Never heard anything disturbing. Major problem remains the delay."
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "slideshow": {
+ "slide_type": "slide"
+ }
+ },
+ "source": [
+ "# The cookbook (finally!)\n",
"\n",
"The recipes have been adapted from Robert Bristow-Johnson's famous [cookbook](https://webaudio.github.io/Audio-EQ-Cookbook/audio-eq-cookbook.html). \n",
"\n",
"Each function returns ``b`` and ``a``, two arrays of three floats each containing the coefficients of the transfer function\n",
"Approximately unit gain (0 dB) in the passband. Specs:\n",
"\n",
" 1. center frequency $f_c$\n",
- " 1. bandwidth $b = (f_+ - f_-)$, where attenuation reaches $-3$ dB. Remember that the passband is almost but not exactly symmetric around $f_c$ (because of bilinear transform)"
+ " 1. bandwidth $b = (f_+ - f_-)$, frequencies where attenuation reaches $-3$ dB. Passband is almost but not exactly symmetric around $f_c$ (because of bilinear transform)"
"plt.gcf().get_axes()[0].axvline(x=(2 * np.pi * CENTER / DEFAULT_SF), linewidth=0.5, color='r');"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Shelves\n",
"\n",
- "Shelving filters are used to amplify either the low or the high end of a signal's spectrum. A high shelf, for instanance, provides an arbitrary gain for high frequencies and has approximately unit gain in the low end of the spectrum. Shelving filters, high or low, are defined by the following parameters:\n",
+ "Shelving filters are used to amplify either the low or the high end of a signal's spectrum; common in consumer audio appliances, behind the standard \"Bass\" and \"Treble\" tone knobs. Specs:\n",
"\n",
- " 1. the desired _shelf gain_ in dB\n",
- " 1. the midpoint frequency $f_c$, which corresponds to the frequency in the transition band where the gain reaches half its value.\n",
- " 1. the \"quality factor\" $Q$, which determines the steepnes off the transition band; as for lowpass filters, the default value $Q = 1/\\sqrt{2}$ yields the steepest transition band while avoiding resonances.\n",
- " \n",
- "A common use case for shelving filters is in consumer audio appliances, where the standard \"Bass\" and \"Treble\" tone knobs control the gain of two complementary shelves with fixed midpoint frequency. "
+ " 1. _shelf gain_ in dB\n",
+ " 1. the midpoint frequency $f_c$, where gain reaches half its value\n",
+ " 1. \"quality factor\" $Q$, as in lowpass filters"
]
},
{
"cell_type": "code",
"execution_count": 23,
"metadata": {},
"outputs": [],
"source": [
"def LSH(fc, gain, sf, Q=(1/np.sqrt(2))):\n",
" \"\"\"Biquad low shelf\"\"\"\n",
" w = 2 * np.pi * fc / sf\n",
" A = 10 ** (gain / 40) \n",
" alpha = np.sin(w) / (2 * Q)\n",
" c = np.cos(w)\n",
" b = np.array([A * ((A + 1) - (A - 1) * c + 2 * np.sqrt(A) * alpha),\n",
" 2 * A * ((A - 1) - (A + 1) * c),\n",
" A * ((A + 1) - (A - 1) * c - 2 * np.sqrt(A) * alpha)])\n",
" a = np.array([(A + 1) + (A - 1) * c + 2 * np.sqrt(A) * alpha,\n",
" -2 * ((A - 1) + (A + 1) * c),\n",
" (A + 1) + (A - 1) * c - 2 * np.sqrt(A) * alpha])\n",
- "A peaking equalizer filter is the fundamental ingrediend in multiband parametric equalization. Each filter provides an arbitrary boost or attenuation for a given frequency band centered around a peak freqency and flattens to unit gain elsewhere. The filter is defined by the following parameters:\n",
+ "Fundamental ingrediend in multiband equalization. Provides an arbitrary boost or attenuation for a given frequency band centered around a peak freqency and flattens to unit gain elsewhere. Specs:\n",
"\n",
- " 1. the desired gain in dB (which can be negative)\n",
- " 1. the peak frequency $f_c$, where the desired gain is attained\n",
- " 1. the bandwidth of the filter, defined as the interval around $f_c$ where the gain is greater (or smaller, for attenuators) than half the desired gain in dB; for instance, if the desired gain is 40dB, all frequencies within the filter's bandwidth will be boosted by at least 20dB. Note that the bandwdidth is not exactly symmetrical around $f_c$ "
+ " 1. desired gain in dB (can be negative)\n",
+ " 1. peak frequency $f_c$, where gain is attained\n",
+ " 1. bandwidth: interval around $f_c$ where gain is more than half of final gain (in dB); again, bandwdidth is not perfectly symmetrical around $f_c$ "
]
},
{
"cell_type": "code",
"execution_count": 29,
- "metadata": {},
+ "metadata": {
+ "slideshow": {
+ "slide_type": "slide"
+ }
+ },
"outputs": [],
"source": [
"def PEQ(fc, bw, gain, sf):\n",
" \"\"\"Biquad bandpass filter \"\"\"\n",
" w = 2 * np.pi * fc / sf\n",
" A = 10 ** (gain / 40) \n",
" alpha = np.tan(np. pi * bw / sf)\n",
" c = np.cos(w)\n",
" b = np.array([1 + alpha * A, -2 * c, 1 - alpha * A])\n",
" a = np.array([1 + alpha / A, -2 * c, 1 - alpha / A])\n",
- " b, a = PEQ(CENTER, BW, GAIN_DB, DEFAULT_SF)\n",
- " y = signal.lfilter(b, a, np.r_[1, np.zeros(200)])\n",
- " plt.plot(y)\n",
- " b, a = PEQ(CENTER, BW, -GAIN_DB, DEFAULT_SF)\n",
- " y = signal.lfilter(b, a, y)\n",
- " plt.plot(y)"
+ "CENTER, BW, GAIN_DB = 800, 400, 40\n",
+ "b, a = PEQ(CENTER, BW, GAIN_DB, DEFAULT_SF)\n",
+ "y = signal.lfilter(b, a, np.r_[1, np.zeros(200)])\n",
+ "plt.plot(y)\n",
+ "b, a = PEQ(CENTER, BW, -GAIN_DB, DEFAULT_SF)\n",
+ "y = signal.lfilter(b, a, y)\n",
+ "plt.plot(y);"
]
},
{
"cell_type": "markdown",
- "metadata": {},
+ "metadata": {
+ "slideshow": {
+ "slide_type": "slide"
+ }
+ },
"source": [
"# Cascades of biquads\n",
"\n",
- "The performance of a single biquad filter may not be adequate for a given application: a second-order lowpass filter, for instance, may not provide a sufficent amount of rejection in the stopband because of its rather slow roll-off characteristic; or we may want to design an equalizer with multiple peaks and dips. In all cases, we usually want to implement the final design as a cascade of biquad sections, because of their inherent numerical robustness."
+ "Often the performance of a second-order filter is not enough.\n",
+ "\n",
+ " * higher-order filters can be decomposed into second-order stages\n",
+ " * \"canonical\" biquads can be repeated to improve or change overall performance\n",
+ "\n",
+ "Applications:\n",
+ " * increase rejection in stopband\n",
+ " * combine shelves to provide tone control\n",
+ " * design multiband equalizers"
]
},
{
"cell_type": "markdown",
- "metadata": {},
+ "metadata": {
+ "slideshow": {
+ "slide_type": "slide"
+ }
+ },
"source": [
"## Factorization of higher-order filters\n",
"\n",
- "The first solution if a filter does not meet the requires specifications is to design a higher-order filter, possibly using different filter \"recipes\"; in the case of bandpass filters, for instance, we could try a Chebyshev or elliptic design. The resulting high-order transfer function can be then factored into a cascade of second-order sections (or, in the case of an odd-order filter, a cascade of second order-sections followed by a first-order filter):\n",
- "The biquad elements returned by the factorization are not related to the \"cookbook\" prototypes of the previous section and therefore this method is simply an implementation strategy that focuses on second-order structures; the design algorithm, in other words, will be dependent on the particular type of filter. Nevertheless, both the design and the factorization are usually available in numerical packages such as Scipy, for instance and, in the following example, we illustrate the difference between a 6th-order elliptic lowpass and a single second-order butterworth. First we will use the full high-order realization and then we will show how a cascade of three second-order sections implements the same characteristic.\n",
- "\n",
- "Note that, when cascading transfer functions, the equivalent higher-order filter coefficients can be obtained simply by polynomial multiplication."
+ " analog_response(b, a, DEFAULT_SF, dB=-60, axis=ax, color=f'C{n}:')\n",
+ " cb = np.polymul(b, cb)\n",
+ " ca = np.polymul(a, ca)\n",
+ "analog_response(cb, ca, DEFAULT_SF, dB=-60, axis=ax, color='C3')"
]
},
{
"cell_type": "markdown",
- "metadata": {},
+ "metadata": {
+ "slideshow": {
+ "slide_type": "slide"
+ }
+ },
"source": [
"## Cascading lowpass and highpass biquads\n",
"\n",
- "A cascade of $N$ identical sections with transfer function $H(z)$ will yield the overall transfer function $H_c(z) = H^N(z)$ and thus the stopband attenuation in decibels will increase $N$-fold. For instance, the following example shows the cumulative magnitude responses obtained by cascading up to five identical second-order Butterworth lowpass sections:"
+ " * biquad's transfer function $H(z)$\n",
+ " * cascade $N$ times\n",
+ " * overall transfer function $H_c(z) = H^N(z)$\n",
+ " * magnitude response in dB multiplied by $N$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "slideshow": {
+ "slide_type": "slide"
+ }
+ },
+ "source": [
+ "Example: cascade of five second-order Butterworth lowpass"
- "As shown by the previous plot, a cascade of identical maximally flat lowpass sections yields a steeper roll-off and preserves the monotonicity of the response. However, since the passband of each filter is not perfectly flat, the $-3~\\mathrm{dB}$ cutoff frequency of the cascade becomes smaller with each added section and the effective bandwidth of the filter is reduced. In the previous example, the original $-3~\\mathrm{dB}$ cutoff frequency was $f_c = 1000~\\mathrm{Hz}$ but the magnitude response of the cascade at $f_c$ is $-15~\\mathrm{dB}$ whereas the actual $-3~\\mathrm{dB}$ point has shifted close to $600~\\mathrm{Hz}$.\n",
+ "Good things:\n",
+ " * still monotonic\n",
+ " * steeper transition to stopband\n",
+ "\n",
+ "Problems:\n",
+ " * cutoff frequency moves\n",
+ " * bandwidth decreases\n",
+ " \n",
+ "Can we (easily) obtain a cascade with controlled $f_c$ and Butterworth characteristic?\n",
+ " * we can always use Python and get a high-order filtered in factored form\n",
+ " * but there's a neat alternative way..."
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "slideshow": {
+ "slide_type": "slide"
+ }
+ },
+ "source": [
+ "### Poles of a Butterworth filter\n",
+ "\n",
+ "$H_N(s)$: transfer function of an order-$N$ Butterworth lowpass\n",
"\n",
- "If our goal is to obtain a cascade with a maximally flat (Butterworth) response with a given $f_c$, an obvious approach is simply to factorize the transfer function of a high-order Butterworth as explained in the previous section. There is however a clever and simpler design strategy that is based on the geometric arrangement of the poles of an analog Butterworth filter of order $N$:\n",
" * the $N$ complex-conjugate poles are equally spaced along a circular contour centered on the origin of the $s$-plane\n",
" * the angle between poles is equal to $\\pi/N$\n",
" \n",
- "With this, the pole angles in the upper $s$-plane are given by \n",
- "Now, a generic second-order analog filter will have a single pair of complex-conjugate poles at $p_{1,2} = \\rho e^{\\pm \\theta}$ on the $s$-plane and, by cascading identical sections, we will only manage to increase the poles' multiplicity but we will not be able to change their position. In order to achieve a Butterworth pole configuration we will thus need to adjust the pole angle for each section; this is a simple task because it turns out that a second-order filter's quality factor $Q$ is related to the pole angle as \n",
- "$$\n",
- " 1/Q = 2\\cos \\theta\n",
- "$$ \n",
- "which means that we can choose the suitable pole angle for each section simply by setting $Q_n = 1/(2\\cos \\theta_n)$. We can now design $N/2$ discrete-time biquads with the same $Q_n$ values to obtain the desired result.\n",
+ "### Poles of an analog second-order lowpass\n",
"\n",
- "Below is the example for a cascade of five lowpass sections (i.e. a 10th-order filter) compared to a single biquad, both with cutoff $f_c = 1000~\\mathrm{Hz}$; notice how the $-3~\\mathrm{dB}$ point has not moved in spite of the much steeper rolloff."
+ " * one pair of complex-conjugate poles at $\\rho e^{\\pm \\theta}$\n",
+ " * fundamental relationship:\n",
+ " $$ 1/Q = 2\\cos \\theta $$\n",
+ " \n",
+ "Solution: simply choose $Q_n = 1/(2\\cos \\theta_n)$ for stage number $n$ in the cascade!"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "slideshow": {
+ "slide_type": "slide"
+ }
+ },
+ "source": [
+ "Example: five lowpass sections (i.e. a 10th-order filter)"
"The resulting digital filter has its poles arranged on a circular contour centered in $z=1$ if the cutoff frequency is less than $\\pi/2$ and centered on $z=-1$ otherwise. "
- "Finally, the following plot shows the individual magnitude responses of the five sections. You can observe that the required $Q_n$ values lead to some biquad sections with a clear peak at the cutoff frequency, although the overall response is monotonic:"
+ "Individual magnitude responses of the five sections. Some $Q_n$ values lead to peaks!"
+ "# rescale response to have unit gain at band edges\n",
+ "cb *= 10 ** (-GAIN / 20)\n",
+ "analog_response(cb, ca, DEFAULT_SF, dB=-50, points=10001);"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "slideshow": {
+ "slide_type": "slide"
+ }
+ },
"source": [
"## Parametric equalization\n",
"\n",
"Peaking equalizers with distinct bandwidths can be cascaded to obtain an arbitrary equalization curve for the entire range of input frequencies; indeed, this is the technique behind so-called _parametric equalizers_ where a bank of logarithmically spaced peaking eq's with independent gain controls allow the user to easily define a global equalization response."
" * $x[n] = \\sin(\\omega_0 n + \\theta)$ with $0 \\le \\omega_0 < 2\\pi$. \n",
" * $\\eta[n] = \\eta(\\sin(\\omega_0 n + \\theta))$ and we are interested in computing its spectrum. \n",
"\n",
"\n",
" * using complex exponentials for the Fourier series:\n",
"$$\n",
" \\eta(x) = \\sum_{k \\neq 0} \\frac{(-1)^{kM}}{j2\\pi k}e^{j\\pi k M x}.\n",
"$$\n",
"\n",
"Now we need to replace $x$ by $\\sin(\\omega_0 n + \\theta)$ and we end up with terms of the form $e^{j \\alpha \\sin \\beta}$; these can be expanded in terms of Bessel functions using the so-called Jacobi-Anger formula:\n",
" * $\\Omega(\\omega_0) = \\{(2m+1)\\omega_0 \\mod 2\\pi\\}_{m \\in \\mathbb{Z}}$, i.e., all the odd multiples of the fundamental frequency aliased over the $[0, 2\\pi]$ interval;\n",
" \n",
" \n",
" * for each frequency $\\varphi \\in \\Omega(\\omega_0)$:\n",
"Assume $\\omega_0 = 2\\pi(A/B)$, with $A$ and $B$ coprime, as in the numerical experiments\n",
"\n",
" * the set $\\Omega(\\omega_0)$ is finite: <br /> $\\displaystyle\\Omega\\left(2\\pi\\frac{A}{B}\\right) = \\left\\{\\frac{2i\\pi}{B}\\right\\}_i, \\quad \\begin{cases}\n",
" i = 0, 1, 2, \\ldots, B-1 & \\mbox{if $A$ or $B$ even} \\\\\n",
" i = 1, 3, 5, \\ldots, B-1 & \\mbox{if $A$ and $B$ odd} \n",
" * contains a finite number of spectral lines at multiples of $2\\pi/B$ \n",
" * the power associated to each line $|b(2i\\pi/B)|^2$ should correspond to the square magnitude of the $i$-th coefficient of the $B$-point DFS of the error signal.\n",
"\n",
"\n",
"\n",
"The following function computes an approximation of the coefficients $|b(2i\\pi/B)|^2$ for $\\omega_0 = 2\\pi(A/B)$, scaled to represent the non-normalized quantization error:"
]
},
{
"cell_type": "code",
"execution_count": 28,
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"outputs": [],
"source": [
"def nqe_sin_psd(A, B, M, phase=1):\n",
" s = [1, -1 if M % 2 == 1 else 1]\n",
" b = np.zeros(B, dtype=complex)\n",
" m_lim, k_lim = max(1500, 2 * B), 600\n",
" for m in range(-m_lim, m_lim):\n",
" c = 0\n",
" for k in range(1, k_lim):\n",
" c += s[k % 2] * ss.jv(2 * m + 1, np.pi * k * M) / k\n",
" c /= 1j * np.pi\n",
" b[((2 * m + 1) * A) % B] += c * np.exp(1j * phase * (2 * m + 1))\n",
" # undo error normalization to obtain the real error PSD\n",
" * the max NHD is dominated by the first aliased component:<br />max NHD is $|c_M(m_0)|^2$ where $m_0$ is the minimum integer for which $(2m_0+1)\\nu > 1/2$. \n",
"\n",
"\n",
"\n",
" * for $\\nu > 1/6$, NHD $\\approx |c_M(1)|^2$\n",