RC series, from analog to digital filters
In this post I invite you to follow me along the analysis of an RC series analog filter and its derivation towards a digital approximation. We will cover the subject from electrical equations, s-domain, transfer functions to bilinear transforms, get onboard !
Electronic circuit setup
First, let's introduce the subject of experimentation, the RC series circuit.
Let's derive the equation of voltage for this circuit. $$ \begin{aligned} V_T &= V_R + V_C \cr & = Ri(t)+V_C \cr & = V_C\left(RC\frac{dV_C(t)}{dt}+1\right)\cr \end{aligned} $$ This is our starting point, we can go further by switching our point of view.
Analog analysis : Laplace transform
For analysis purposes, we apply the Laplace transform on this equation. Here, we only need to know that derivation (of order $n$) is the same as multiplying by $s^n$.
$$ \begin{aligned} \mathcal{L}[V_T] &= \mathcal{L}\left[V_C\left(RC\frac{dV_C(t)}{dt}+1\right)\right]\cr Y(s) &= X(s)(sRC + 1) \cr \end{aligned} $$
From which we obtain the transfer function in Laplace domain. $$H(s) = \frac{Y(s)}{X(s)} = \frac{1}{sRC+1}$$
This is single pole transfer function, the pole makes the denominator equal to $0$.
$$p_0 = -\frac{1}{RC}$$
Rewriting our transfer function with $p_0$:
$$H(s) =\frac{k}{s-p_0}, k=-p_0$$
We can plot the module surface of $H(s)$ in the complex $s$-domain along $\sigma$ and $\omega$, here is a contour view of it with our pole $p_0$.
$p_0$ is real, which is a special type of pole because it means that the inverse Laplace transform (which one can use to find to associated time expression) is simply a real exponential.
$$ \mathcal{L}^{-1}\left[\frac{1}{sRC+1}\right]=e^{-\frac{t}{RC}}$$
From it, we can derive the time response of our circuit to a voltage step $\gamma(t)$.
$$ \begin{aligned} Y(s)&=\mathcal{L}[\gamma(t)]\cdot H(s) \cr &= \Gamma(s)H(s) \cr &= \frac{1}{s}H(s) \cr &= \frac{1}{s}\frac{k}{s-p_0} \end{aligned} $$
From partial fraction expansion, we know we can express this as:
$$Y(s)=\frac{a}{s}+\frac{b}{s-p_0}$$
And with a bit of equation scrambling, we find $a=1$ and $b=-1$. Finding the temporal response of our circuit to $\gamma(t)$ is now as simple as finding the inverse Laplace transform of $Y(s)$.
$$\mathcal{L}^{-1}\left[\frac{1}{s}-\frac{1}{s-p_0}\right] = 1 - e^{p_0t} = 1 - e^{-\frac{t}{RC}}$$
Familiar ? Of course !
If we remember the definition of the Laplace transform, we can plot the frequency response of this filter by estimating the transfer function on the y-axis of the Laplace complex domain.
So, we replace $s=j\omega$ in $H(s)$:
$$H(j\omega) = \frac{1}{j\omega RC+1}$$
And we get the response in amplitude and phase from the module and angle of $H(j\omega)$.
Our analog filter is well defined, but what if we wanted to bring it in the digital world ? Introducing, z-transform !
Digital analysis : z-transform
There are several ways to transpose an analog system into a digital one, but today we'll use the bilinear transform.
The bilinear transform consists in replacing the Laplace variable $s$ by a fraction doubly linear in $z$.
$$s = \frac{2}{T}\frac{1-z^{-1}}{1+z^{-1}} $$
With the sampling period, $T=1/f_s=2\pi/\omega_s$. The transform has the property to map the infinite $j\omega_{analog}$ axis to the finite unit circle $e^{jw_{digital}T}$. The mapping is however non-linear and presents frequency warping around Nyquist. Below is plotted the mapping of $s$ to $z$ for several frequencies, notice the warp around $-.5$ and $.5$ normalized frequencies (i.e. same spacing on the $j\omega$ axis equals smaller and smaller spacings on the unit circle).
The objective of this transform is to find the difference equation associated with this filter. A difference equation links current output as a weighted sum of past inputs and outputs, like so :
$$y[n] + a_1y[n-1]+...+a_ny[n-m]= b_0x[n]+b_1x[n-1]+...+b_nx[n-m]$$
The digital transfer function can be found by applying the z-transform :
$$H(z)=\frac{b_0+b_1z^{-1} + ...+ b_nz^{-n}}{1+a_1z^{-1}+... +a_nz^{-n}}$$
One way of understanding why delays of $n$ become $z^{-n}$ is because, on the unit circle, we are dealing with complex exponential sampled signals in the form of $e^{jwkT}$, so powers are actually a way to go back and forth in discrete time.
Now, let's derive the expression for our RC filter.
$$ H(z)=\frac{1}{\frac{2RC}{T}\frac{1-z^{-1}}{1+z^{-1}}+1} = \frac{T+Tz^{-1}}{T+2RC+ (T-2RC)z^{-1}} $$
$a_0$ should be equal to $1$, so we normalize by $T+2RC$. $$H(z)=\frac{\frac{T}{T+2RC}+\frac{T}{T+2RC}z^{-1}}{1+\frac{T-2RC}{T+2RC}z^{-1}}$$
Which gives us our coefficients $a_i$ and $b_i$. $$a_0 = 1, a_1=\frac{T-2RC}{T+2RC}$$ $$b_0 = b_1 =\frac{T}{T+2RC}$$
We can now start to filter sampled signals and test if our digital filter works! I will generate a random gaussian noise at $f_s=$, which is supposed to have a flat frequency response, and filter it using our newly found coefficients. I used the function lfilter from scipy.signalmodule. Next, I plot the simulated frequency response by computing an average Discrete Fourier Transform next to theoretical analog gain.
The response looks perfect, but it isn't the whole picture, if I expand the $x$ axis and zoom toward the Nyquist frequency, we see the impact of the non-linear frequency warping...
The deep dive happens because we are actually mapping the whole range of frequency up to $+\inf$ inside the range $[0, f_s/2]$. Fortunately, one can tweak the non-linear mapping to preserve frequency responses profiles at certain key frequencies by pre-warping the $j\omega$ axis before the bilinear transform, but that's for another post.
Also as a quick note, no one really design digital filters like that for real applications under strict specifications, the method we used in this post is mainly useful for emulating existing continuous system responses (think about audio processing/modeling).
I hope you liked this post, and see you next time.
