[Jupyter notebook]

1  Table of Contents

•  
  • Synthesis of IIR filters by the bilinear transformation method
    • The bilinear transform
    • Synthesis of low-pass filters -- procedure
    • Synthesis of other type of filters
    • Numerical results

$$\begin{equation}\require{color} \require{cancel} \def\tf#1{{\mathrm{FT}\left\{ #1 \right\}}} \def\flecheTF{\rightleftharpoons } \def\TFI#1#2#3{{\displaystyle{\int_{-\infty}^{+\infty} #1 ~e^{j2\pi #2 #3} ~\dr{#2}}}} \def\TF#1#2#3{{\displaystyle{\int_{-\infty}^{+\infty} #1 ~e^{-j2\pi #3 #2} ~\dr{#2}}}} \def\sha{ш} \def\dr#1{\mathrm{d}#1} \def\egalpardef{\mathop{=}\limits^\triangle} \def\sinc#1{{\mathrm{sinc}\left( #1 \right)}} \def\rect{\mathrm{rect}} \definecolor{lightred}{rgb}{1,0.1,0} \def\myblueeqbox#1{{\fcolorbox{blue}{lightblue}{$ extcolor{blue}{ #1}$}}} \def\myeqbox#1#2{{\fcolorbox{#1}{light#1}{$ extcolor{#1}{ #2}$}}} \def\eqbox#1#2#3#4{{\fcolorbox{#1}{#2}{$\textcolor{#3}{ #4}$}}} \def\eqboxa#1{{\boxed{#1}}} \def\eqboxb#1{{\eqbox{green}{white}{green}{#1}}} \def\eqboxc#1{{\eqbox{blue}{white}{blue}{#1}}} \def\eqboxd#1{{\eqbox{blue}{lightblue}{blue}{#1}}} \def\E#1{\mathbb{E}\left[#1\right]} \def\ta#1{\left<#1\right>} \def\egalparerg{{\mathop{=}\limits_\mathrm{erg}}} \def\expo#1{\exp\left(#1\right)} \def\d#1{\mathrm{d}#1} \def\wb{\mathbf{w}} \def\sb{\mathbf{s}} \def\xb{\mathbf{x}} \def\Rb{\mathbf{R}} \def\rb{\mathbf{r}} \def\mystar{{*}} \def\ub{\mathbf{u}} \def\wbopt{\mathop{\mathbf{w}}\limits^\triangle} \def\deriv#1#2{\frac{\mathrm{d}#1}{\mathrm{d}#2}} \def\Ub{\mathbf{U}} \def\db{\mathbf{d}} \def\eb{\mathbf{e}} \def\vb{\mathbf{v}} \def\Ib{\mathbf{I}} \def\Vb{\mathbf{V}} \def\Lambdab{\mathbf{\Lambda}} \def\Ab{\mathbf{A}} \def\Bb{\mathbf{B}} \def\Cb{\mathbf{C}} \def\Db{\mathbf{D}} \def\Kb{\mathbf{K}} \def\sinc#1{\mathrm{sinc\left(#1\right)}} \end{equation}$$

import warnings
warnings.filterwarnings("ignore", category=DeprecationWarning) 

import seaborn as sns                 #this can be removed safely  if seaborn is not available
sns.mpl.rc("figure", figsize=(7,4))   #idem 

1.1  Synthesis of IIR filters by the bilinear transformation method

A simple and effective method for designing IIR digital filters with prescribed magnitude response specifications is the bilinear transformation method. The point is that already exists well-known optimal methods for the design of analog filters, such as Butterworth, Chebyshev, or elliptic filter designs. Then, the idea is to map the digital filter into an equivalent analog filter, which can be designed optimally, and map back the design to the digital domain. The key for this procedure is to dispose of a reversible mapping from the analog domain to the digital domain.

1.1.1  The bilinear transform

Recall that in the analog domain, the equivalent of the Z-transform is the Laplace transform which associates a signal $s(t)$ with a function $S(p)$ of the complex variable $p$. When $p$ lives on the imaginary axis of the complex plane, then the Laplace transform reduces to the Fourier transform (for causal signals). For transfer functions, stability is linked to the positions of poles in the complex plane. They must have a negative real part (that is belong to the left half plane) to ensure the stability of the underlying system.

The formula for the bilinear transform comes from a formal analogy between the derivation operator in the Laplace and Z domains.

The bilinear transform makes an association between analog and digital frequencies, as follows:

$$\begin{equation} \eqboxd{ \displaystyle{ p=k \frac{1-z^{-1}}{1+z^{-1}}}}, \label{eq:BLT} \end{equation}$$

where $k$ is an arbitrary constant. The usual derivation leads to $k=2/T_s$, where $T_s$ is the sampling period. However, using a general parameter $k$ does not change the methodology and offers a free parameter that enables to simplify the procedure.

The point is that this transform presents some interesting and useful features:

1. It preserves stability and minimum phase property (the zeros of the transfer function are with negative real part (analog case) or are inside the unit circle (discrete case).
2. It maps the infinite analog axis into a periodic frequency axis in the frequency domain for discrete signals. That mapping is highly non linear and warp the frequency components, but it recovers the well-known property of periodicity of the Fourier transform of discrete signals.

The corresponding mapping of frequencies is obtained as follows. Letting $p=j\omega_a=j2\pi f_a$ and $z=\exp(\omega_d)=\exp(j2\pi f_d)$. Plugging this in ($\ref{eq:BLT}$), we readily obtain

$$\begin{equation} \eqboxc{\displaystyle{ \omega_a=k \tan\left( \frac{\omega_d }{2}\right)}}, \label{eq:prewarping} \end{equation}$$

or

$$\begin{equation} \eqboxc{\displaystyle{ \omega_d=2 \arctan\left( \frac{\omega_a }{k}\right)}}. \label{eq:prewarping_b} \end{equation}$$

The transformation ($\ref{eq:prewarping}$) corresponds to the initial transform of the specifications in the digital domain into analog domain specifications. It is often called a pre-warping . Figure shows the mapping of pulsations from one domain into the other one.

k = 2
xmin = -5 * pi
xmax = -xmin
omegaa = np.arange(xmin, xmax, 0.1)
omegad = 2 * np.arctan(omegaa / k)
plt.plot(omegaa, omegad)
plt.plot([xmin, xmax], [-pi, -pi], '--', color='lightblue')
plt.plot([xmin, xmax], [pi, pi], '--', color='lightblue')
#plt.text(-3.7,0.4,'Fs/2', color='blue',fontsize=14)
plt.xlabel("Analog pulsations $\omega_a$")
plt.ylabel("Digital pulsations $\omega_d$")
_ = plt.xlim([xmin, xmax])
plt.title("Frequency mapping of the bilinear transform")
figcaption("Frequency mapping of the bilinear transform", label="fig:BLT")

Caption: Frequency mapping of the bilinear transform

When designing a digital filter using an analog approximation and the bilinear transform, we follow these steps: Pre-warp the cutoff frequencies Design the necessary analog filter apply the bilinear transform to the transfer function Normalize the resultant transfer function to be monotonic and have a unity passband gain (0dB).

1.1.2  Synthesis of low-pass filters -- procedure

Let $\omega_p$ and $\omega_s$ denote the edges of the pass-band and of the stop-band.

1. For the synthesis of the analog filter, it is convenient to work with a normalized filter such that $\Omega_p=1$. Therefore, as a first step, we set $$\begin{equation} k = \arctan(\omega_p/2) \end{equation}$$ which ensures that $\Omega_p=1$. Then, we compute $\Omega_s=2 \arctan\left( \omega_s /k\right)$.

2. Synthetize the optimum filter in the analog domain, given the type of filter, the frequency and gain constraints. This usually consists in determining the order of the filter such that the gain constraints (ripples, minimum attenuation, etc) are satisfied, and then select the corresponding polynomial. This yields a transfer function $H_a(p)$.

3. Map back the results to the digital domain, using the bilinear transform ($\ref{eq:BLT}$), that is compute $$\begin{equation} \eqboxc{\displaystyle{ H(z) = \left. H_a(p) \right|_{p=k \frac{1-z^{-1}}{1+z^{-1}}}}} \end{equation}$$

Exercise 1

We want to synthetize a digital filter with $f_p=6kHz$, with a maximum attenuation of -3dB, and a stop-band frequency of $f_s=9kHz$, with a minimum attenuation of -9dB. The Nyquist rate (sampling frequency) is $F_s=36kHz$. - Represent the template for this sythesis, - Compute the pulsations in the analog domain (fix $\Omega_p=1$). - if we choose a Butterworth filter, the best order is $n=2$ and the corresponding polynomial is $D(p)=p^2+\sqrt{2}p+1$, and the transfer function is $H_a(p)=1/D(p)$. Compute $H(z)$. - Plot $H(f)$. Use `sig.freqz` for computing the transfer function. We also provide a function `plt_LPtemplate(omega, A, Abounds=None)` which displays the template of the filter. Import it using `from plt_LPtemplate import *`.

Elements of solution

• $k=1/\tan(\pi/6)=\sqrt{3}$
• $\Omega_s = k \tan(\pi/4) = \sqrt{3}$
• $$\begin{equation} H(z) = \frac{1+2z^{-1}+z^{-2}}{(4+\sqrt{6})-4 z^{-1} + (4-\sqrt{6})z^{-2}} \end{equation}$$

# compute the transfer function using freqz 
w, H = sig.freqz([1, 2, 1], [4 + sqrt(6), -4, 4 - sqrt(6)], whole=True)
# plot the result --w-pi corresponds to a shift of the pulsation 
#axis associated with the fftshift of the transfer function. 
plt.plot(w - pi, 20 * np.log10(fftshift(np.abs(H))))
# plot the filter template
from plt_LPtemplate import *
plt_LPtemplate([pi / 3, pi / 2], [-3, -9], Abounds=[5, -35])
plt.title("Realized filter and its template")
figcaption("Realized filter and its template")

/home/bercherj/.local/lib/python3.5/site-packages/ipykernel_launcher.py:5: RuntimeWarning: divide by zero encountered in log10
  """

Caption: Realized filter and its template

In practice, the transfer function can be generated by transforming the poles and zeros of the analog filter into the poles and zeros of the digital filter. This is simply done using the transformation $$\begin{equation} z=\frac{1+p/k}{1-p/k}. \end{equation}$$ It remains a global gain that can be fixed using a value at $H(1)$ ($\omega=0)$.This dramatically simplifies the synthesis in practice.

1.1.3  Synthesis of other type of filters

For other kind of filters (namely band-pass, high-pass, band-stop), we can either:

• use the bilinear transform and the pre-warping to obtain a filter of the same type in the analog domain; then transferring the problem to the analog designer.
• use an auxiliary transform that converts a digital low-pass filter into another type of filter. Using the second option, we see that the low-pass synthesis procedure has to be completed with

1. Transform the specifications of the digital filter into specifications for a low-pass digital filter.

1.1.3.1  Transposition from a low-pass filter ($\omega_p$) to another type

Let $\omega_1$ and $\omega_2$ denote the low and high cut-off frequencies (only $\omega_1$ for the transposition of a low-pass into a high-pass). These transformations preserve the unit circle. That is, $z=e^{j\theta}$ is transformed into $z=e^{j\theta'}$. There is an additional frequency warping, but the notion of frequency is preserved.

1. low-pass $\omega_p$-- low-pass $\omega_1$ $$\begin{equation} z^{-1} \rightarrow \frac{z^{-1}-\alpha}{1-\alpha z^{-1}} \end{equation}$$ with $$\begin{equation} \alpha=\frac{\sin\left(\frac{\omega_p-\omega_1}{2}\right)}{\sin\left(\frac{\omega_p+\omega_1}{2}\right)} \end{equation}$$

2. low-pass $\omega_p$-- band-pass $\omega_1, \omega_2$ $$\begin{equation} z^{-1} \rightarrow - \frac {z^{-2}-\frac{2\alpha \beta}{\beta+1} z^{-1} + \frac{\beta-1}{\beta+1}} {\frac{\beta-1}{\beta+1}z^{-2} -\frac{2\alpha \beta}{\beta+1} z^{-1} + 1} \end{equation}$$ with $$\begin{equation} \alpha=\frac{\cos\left(\frac{\omega_p+\omega_1}{2}\right)}{\cos\left(\frac{\omega_p-\omega_1}{2}\right)} \end{equation}$$ and $$\begin{equation} \beta=\tan\left(\frac{\omega_p}{2}\right)\tan\left(\frac{\omega_2-\omega_1}{2}\right) \end{equation}$$

Exercise 2

We want to synthetize an high-pass digital filter with edge frequencies $2, 4, 12$ and $14$ kHz, with a maximum attenuation of 3dB in the band-pass, and a minimum attenuation of -12dB in the stop-band. The Nyquist rate (sampling frequency) is $F_s=32kHz$. - Represent the template for this sythesis, - Compute the pulsations in the analog domain (fix $\Omega_p=1$). - if we choose a Butterworth filter, the best order is $n=2$ and the corresponding polynomial is $D(p)=p^2+\sqrt{2}p+1$, and the transfer function is $H_a(p)=1/D(p)$. Compute $H(z)$. - Plot $H(f)$. Use `sig.freqz` for computing the transfer function.

Sketch of solution

• normalized pulsations: $\omega_0=\frac{\pi}{8}, \omega_1=\frac{\pi}{4}, \omega_2=\frac{3\pi}{4}, \omega_3=\frac{7\pi}{8}$
• transposition into a digital low-pass $\alpha=0, \beta=\tan(\omega_p/2)$. Choosing $\beta=1$, we get $\omega_p=\pi/2$
• the transform is thus $z^{-1} \rightarrow -z^{-2}$
• In the BLT, we take $k=1$; thus $\omega_s=1$ and $\omega_s=\tan(3\pi/8)$
• Finally, we obtain $$\begin{equation} H(z) = \frac{1+2z^{-1}+z^{-2}}{(2+\sqrt{2}) + (2-\sqrt{2})z^{-2}} \end{equation}$$ for the digital low-pass, which after the transform $z^{-1} \rightarrow -z^{-2}$ gives $$\begin{equation} H(z) = \frac{1-2z^{-2}+z^{-4}}{(2+\sqrt{2}) + (2-\sqrt{2})z^{-4}} \end{equation}$$

# compute the transfer function using freqz 
w, H = sig.freqz(
    [1, 0, -2, 0, 1], [2 + sqrt(2), 0, 0, 0, 2 - sqrt(2)], whole=True)
# plot the result --w-pi corresponds to a shift of the pulsation 
#axis associated with the fftshift of the transfer function. 
plt.plot((w - pi) / (2 * pi) * 32, 20 * np.log10(fftshift(np.abs(H))))
plt.xlabel("Frequencies")
_ = plt.ylim([-15, 0])

/home/bercherj/.local/lib/python3.5/site-packages/ipykernel_launcher.py:6: RuntimeWarning: divide by zero encountered in log10
  

1.1.4  Numerical results

Finally, we point out that these procedures have been systematized into computer programs. Two functions are available in scipy to design an IIR filter: sig.iirfilter that computes the coefficients of a filter for a given order, and sig.iirdesign that even determine a correct order given constraint on maximum and/or minimum attenuation. It is instructive to consult the help of these functions (e.g. help(sig.iirdesign)) and try to reproduce the results we obtained analytically above. Possible solutions are provided below.

b, a = sig.iirfilter(
    2, [1 / (pi)],
    rp=None,
    rs=None,
    btype='lowpass',
    analog=False,
    ftype='butter',
    output='ba')
# compute the ttransfer function using freqz 
w, H = sig.freqz(b, a, whole=True)
# plot the result --w-pi corresponds to a shift of the pulsation 
#axis associated with the fftshift of the transfer function. 
plt.plot(w - pi, 20 * np.log10(fftshift(np.abs(H))))
# plot the filter template
from plt_LPtemplate import *
plt_LPtemplate([pi / 3, pi / 2], [-3, -9], Abounds=[12, -35])
plt.title("Realized filter and its template")
figcaption("Realized filter and its template")

/home/bercherj/.local/lib/python3.5/site-packages/ipykernel_launcher.py:13: RuntimeWarning: divide by zero encountered in log10
  del sys.path[0]

Caption: Realized filter and its template

b, a = sig.iirdesign(
    [4 / 16, 12 / 16], [2 / 16, 14 / 16],
    3,
    12,
    analog=False,
    ftype='butter',
    output='ba')
# compute the ttransfer function using freqz 
w, H = sig.freqz(b, a, whole=True)
# plot the result --w-pi corresponds to a shift of the pulsation 
#axis associated with the fftshift of the transfer function. 
plt.plot((w - pi) / (2 * pi) * 32, 20 * np.log10(fftshift(np.abs(H))))
_ = plt.ylim([-20, 0])

/home/bercherj/.local/lib/python3.5/site-packages/ipykernel_launcher.py:12: RuntimeWarning: divide by zero encountered in log10
  if sys.path[0] == '':

---

Index - Back - Next