Page MenuHomec4science
Files F105172466

2_iir.tex
No OneTemporary
Actions

File Metadata

Created
Sat, Mar 15, 04:14

2_iir.tex
View Options

\documentclass[aspectratio=169]{beamer}
\def\stylepath{../styles}
\usepackage{\stylepath/com202}
\begin{document}
\begin{frame} \frametitle{IIR: conversion of analog design}
Filter design was an established art long before digital processing appeared
\begin{itemize}
\item lots of nice analog filters exist
\item methods exist to ``translate'' the analog design into a rational transfer function
\item most numerical packages (Matlab, Numpy, etc.) provide ready-made routines
\item design involves specifying some parameters and testing that the specs are fulfilled
\end{itemize}
\end{frame}
\begin{frame} \frametitle{Three classic filter families to be aware of}
\begin{itemize}
\item Butterworth (smooth monotonic frequency response)
\item Chebyshev (monotonic/equiripple)
\item Elliptic (equiripple)
\end{itemize}
\end{frame}
\begin{frame} \frametitle{Butterworth lowpass}
\begin{columns}
\begin{column}{0.4\paperwidth}
Magnitude response:
\begin{itemize}
\item maximally flat
\item monotonic over $[0,\pi]$
\end{itemize}
\vspace{2em}
Design parameters:
\begin{itemize}
\item order $N$ ($N$ poles and $N$ zeros)
\item cutoff frequency
\end{itemize}
\end{column}
\pause
\begin{column}{0.4\paperwidth}
Design test criterion:
\begin{itemize}
\item width of transition band
\item passband error
\end{itemize}
\end{column}
\end{columns}
\end{frame}
\begin{frame}[fragile] \frametitle{Butterworth lowpass design with SciPy}
\begin{verbatim}
import scipy.signal as sp
b, a = sp.butter(4, 0.25)
wb, Hb = sp.freqz(b, a, 1024);
plt.plot(wb/np.pi, np.abs(Hb));
\end{verbatim}
\end{frame}
\begin{frame} \frametitle{Butterworth lowpass example}
\centering
$N=4, \omega_c = \pi/4$
\begin{figure}
\begin{dspPlot}[xtype=freq,xticks=4,ylabel={$|H(e^{j\omega})|$}]{-1,1}{0,1.1}
\moocStyle
% coeffs: [b a]=butter(4,0.25)
\dspFunc{x \dspTFM{ 0.0102 0.0408 0.0613 0.0408 0.0102}{ 1.0000 -1.9684 1.7359 -0.7245 0.1204}}
\end{dspPlot}
\end{figure}
\end{frame}
\begin{frame} \frametitle{Butterworth lowpass example}
\centering
$N=8, \omega_c = \pi/4$
\begin{figure}
\begin{dspPlot}[xtype=freq,xticks=4,ylabel={$|H(e^{j\omega})|$}]{-1,1}{0,1.1}
\moocStyle
% coeffs: [b a]=butter(8,0.25)
\dspFunc{x \dspTFM{0.000107911284731 0.000863290277849 0.003021515972471 0.006043031944942 0.007553789931178 0.006043031944942 0.003021515972471 0.000863290277849 0.000107911284731
}{ 1.000000000000000 -3.983784273174195 7.536234110120903 -8.599815064801408 6.400154060347649 -3.156025260730575 1.001696579551288 -0.186342477677486 0.015507615254987}}
\end{dspPlot}
\end{figure}
\end{frame}
\begin{frame} \frametitle{Chebyshev lowpass}
\begin{columns}
\begin{column}{0.4\paperwidth}
Magnitude response:
\begin{itemize}
\item equiripple in passband
\item monotonic in stopband
\item (or vice-versa)
\end{itemize}
\vspace{2em}
Design parameters:
\begin{itemize}
\item order $N$ ($N$ poles and $N$ zeros)
\item passband max error
\item cutoff frequency
\end{itemize}
\end{column}
\pause
\begin{column}{0.4\paperwidth}
Design test criterion:
\begin{itemize}
\item width of transition band
\item stopband error
\end{itemize}
\end{column}
\end{columns}
\end{frame}
\begin{frame}[fragile] \frametitle{Chebyshev lowpass design with SciPy}
\begin{verbatim}
b, a = sp.cheby1(4, .12, 0.25)
\end{verbatim}
\end{frame}
\begin{frame} \frametitle{Chebyshev lowpass example}
\centering
$N=4, \omega_c = \pi/4, e_{\max} = 12\%$
\begin{figure}
\begin{dspPlot}[xtype=freq,xticks=4,ylabel={$|H(e^{j\omega})|$}]{-1,1}{0,1.1}
\moocStyle
% coeffs: [b a]=cheby1(4,0.5,0.25)
\dspFunc{x \dspTFM{0.0056 0.0225 0.0337 0.0225 0.0056}{1.0000 -2.5614 2.9222 -1.6586 0.3931}}
\end{dspPlot}
\end{figure}
\end{frame}
\begin{frame} \frametitle{Chebyshev lowpass example}
\centering
$N=8, \omega_c = \pi/4, e_{\max} = 12\%$
\begin{figure}
\begin{dspPlot}[xtype=freq,xticks=4,ylabel={$|H(e^{j\omega})|$}]{-1,1}{0,1.1}
\moocStyle
% coeffs: [b a]=cheby1(8,0.5,0.25)
\dspFunc{x \dspTFM{0.000008952611389 0.000071620891113 0.000250673118897 0.000501346237795 0.000626682797244 0.000501346237795 0.000250673118897 0.000071620891113 0.000008952611389}{ 1.000000000000000 -5.975292291885454 16.581223292021008 -27.714232735429224 30.395097583553124 -22.347296704268793 10.745098004349103 -3.089246336974975 0.407076858898017}}
\end{dspPlot}
\end{figure}
\end{frame}
\begin{frame} \frametitle{Elliptic lowpass}
\begin{columns}
\begin{column}{0.4\paperwidth}
Magnitude response:
\begin{itemize}
\item equiripple in passband and stopband
\end{itemize}
\vspace{2em}
Design parameters:
\begin{itemize}
\item order $N$
\item cutoff frequency
\item passband max error
\item stopband min attenuation
\end{itemize}
\end{column}
\pause
\begin{column}{0.4\paperwidth}
Design test criterion:
\begin{itemize}
\item width of transition band
\end{itemize}
\end{column}
\end{columns}
\end{frame}
\begin{frame}[fragile] \frametitle{Elliptic lowpass design with SciPy}
\begin{verbatim}
b, a = sp.ellip(4, .1, 50, 0.25)
\end{verbatim}
\end{frame}
\begin{frame} \frametitle{Elliptic lowpass example}
\centering
$N=4, \omega_c = \pi/4, e_{\max} = 12\%, \mbox{att}_{\min} = 0.03$
\begin{figure}
\begin{dspPlot}[xtype=freq,xticks=4,ylabel={$|H(e^{j\omega})|$}]{-1,1}{0,1.1}
\moocStyle
% coeffs: [b a]=ellip(4,0.5,15,0.25)
\dspFunc{x \dspTFM{0.1866 -0.3248 0.4782 -0.3248 0.1866}{1.0000 -2.4512 2.8891 -1.6623 0.4380}}
\end{dspPlot}
\end{figure}
\end{frame}
\begin{frame} \frametitle{Elliptic lowpass example}
\centering
$N=6, \omega_c = \pi/4, e_{\max} = 12\%, \mbox{att}_{\min} = 0.03$
\begin{figure}
\begin{dspPlot}[xtype=freq,xticks=4,ylabel={$|H(e^{j\omega})|$}]{-1,1}{0,1.1}
\moocStyle
% coeffs: [b a]=ellip(6,0.5,15,0.25)
\dspFunc{x \dspTFM{ 0.183242079035887 -0.594902462167290 1.145062255460029 -1.360125010039982 1.145062255460030 -0.594902462167291 0.183242079035887}{1.000000000000000 -3.910500497053838 7.473515172378240 -8.353747075193558 5.781357181637871 -2.318511593591088 0.440886658862916}}
\end{dspPlot}
\end{figure}
\end{frame}
\begin{frame} \frametitle{Elliptic lowpass example: numerical errors for high-order}
\centering
$N=8, \omega_c = \pi/4, e_{\max} = 12\%, \mbox{att}_{\min} = 0.03$
\begin{figure}
\begin{dspPlot}[xtype=freq,xticks=4,ylabel={$|H(e^{j\omega})|$}]{-1,1}{0,1.1}
\moocStyle
% coeffs: [b a]=ellip(8,0.5,15,0.25)
\dspFunc{x \dspTFM{ 0.182945412317538 -0.854706836398519 2.173711958822330 -3.584001416594707 4.225696280882447 -3.584001416594708 2.173711958822331 -0.854706836398519 0.182945412317538
}{ 1.000000000000000 -5.331625887296879 14.031349749371518 -22.887835522028467 25.132537021945247 -18.895472499481155 9.522145636275305 -2.947011314805802 0.441157037789125}}
\end{dspPlot}
\end{figure}
\end{frame}
\begin{frame} \frametitle{Let's compare}
\begin{itemize}
\item compare magnitude response of 4th-order lowpass filters
\item same cutoff frequency and transition band width
\item plot the magnitude response in dB
\end{itemize}
\end{frame}
\begin{frame} \frametitle{The decibel for amplitude ratios}
Relative measure of amplitude in log scale:
\[
|H(e^{j\omega})|_{\text{dB}} = 20\log_{10}\frac{|H(e^{j\omega})|}{H_{\text{ref}}}
\]
Here we choose $H_{\text{ref}} = 1$, target value in passband.
\vspace{2em}
\begin{itemize}
\item -6 dB = half the amplitude
\item -20 dB = one tenth of the amplitude
\end{itemize}
\end{frame}
\begin{frame} \frametitle{4-th order IIR lowpass comparison}
\begin{figure}
\begin{dspPlot}[xtype=freq,xticks=4,yticks=20,ylabel={$|H(e^{j\omega})|$ (dB)},xout=true]{-1,1}{-130,0}
\moocStyle
\dspFunc{x \dspTFM{ 0.0102 0.0408 0.0613 0.0408 0.0102}{ 1.0000 -1.9684 1.7359 -0.7245 0.1204} log 20 mul}
\dspFunc[linecolor=blue!60]{x \dspTFM{0.0056 0.0225 0.0337 0.0225 0.0056}{1.0000 -2.5614 2.9222 -1.6586 0.3931} log 20 mul}
\dspFunc[linecolor=green!60]{x \dspTFM{0.1866 -0.3248 0.4782 -0.3248 0.1866}{1.0000 -2.4512 2.8891 -1.6623 0.4380} log 20 mul}
\dspLegend(-.3,-80){darkred {Butterworth} blue!60 {Chebyshev} green!60 {elliptic}}
\end{dspPlot}
\end{figure}
\centering
\small
all filters require 9 multiplications per output sample
\end{frame}
\begin{frame} \frametitle{Qualitative comparison}
For a given order $N$
\begin{itemize}
\item sharpness of transition band: Elliptic $>$ Chebyshev $>$ Butterworth
\item phase distortion: Butterworth $<$ Chebyshev $<$ Elliptic
\item passband ripples Butterworth $<$ Chebyshev $<$ Elliptic
\item stopband attenuation: Elliptic $>$ Chebyshev $>$ Butterworth
\end{itemize}
\end{frame}
\begin{frame} \frametitle{Elliptic lowpass example: numerical errors for high-order}
\centering
$N=8, \omega_c = \pi/4, e_{\max} = 12\%, \mbox{att}_{\min} = 0.03$
\begin{figure}
\begin{dspPlot}[xtype=freq,xticks=4,ylabel={$|H(e^{j\omega})|$}]{-1,1}{0,1.1}
\moocStyle
% coeffs: [b a]=ellip(8,0.5,15,0.25)
\dspFunc{x \dspTFM{ 0.182945412317538 -0.854706836398519 2.173711958822330 -3.584001416594707 4.225696280882447 -3.584001416594708 2.173711958822331 -0.854706836398519 0.182945412317538
}{ 1.000000000000000 -5.331625887296879 14.031349749371518 -22.887835522028467 25.132537021945247 -18.895472499481155 9.522145636275305 -2.947011314805802 0.441157037789125}}
\end{dspPlot}
\end{figure}
\end{frame}
\begin{frame} \frametitle{Numerical precision issues}
\begin{itemize}
\item all digital devices represent numbers using finite precision
\item poles are the roots of the denominator of the transfer function
\item filter algorithms store the value of the coefficients, not of the poles
\item the value of a pole is a nonlinear function of the filter coefficients
\item insufficient numerical precision may cause poles to drift out of unit circle
\end{itemize}
\end{frame}
\begin{frame} \frametitle{Pole drifting: example}
\begin{itemize}
\item nominal pole: $p = \rho e^{j\theta}$
\item second-order transfer function: $P(z) = (1 - p z^{-1})(1 - p^* z^{-1})$
\item $P(z) = 1 + a_1 z^{-1} + a_2 z^{-2}$, with $a_1 = - 2\rho\cos\theta$ and $a_2 = \rho^2$
\item $p' = (\sqrt{a_1^2 - 4a_2} - a_1)/2$
\end{itemize}
\vspace{1em}
\centering
\small
\begin{tabular}{c||l}
\# digits & $|p| - |p|'$ \\ \hline
8 & 2.220446049250313e-16 \\
7 & 5.000500014062936e-09 \\
4 & 5.000499903040634e-09 \\
3 & 0.00040012506253905844 \\
2 & 0.004912562893379935
\end{tabular}
\end{frame}
\begin{frame} \frametitle{Poles of the 8th order elliptic lowpass}
\centering
\begin{dspPZPlot}[width=6cm,clabel={$1$}]{1.2}
\moocStyle
\dspPZ[label=none,linecolor=blue]{ 0.7068075307693427 , 0.706979289475587 }
\dspPZ[label=none,linecolor=blue]{ 0.7068075307693427 , -0.706979289475587 }
\dspPZ[label=none,linecolor=blue]{ 0.7050364061542383 , 0.704114656343972 }
\dspPZ[label=none,linecolor=blue]{ 0.7050364061542383 , -0.704114656343972 }
\dspPZ[label=none,linecolor=blue]{ 0.6870469975049976 , 0.6738034200692913 }
\dspPZ[label=none,linecolor=blue]{ 0.6870469975049976 , -0.6738034200692913 }
\dspPZ[label=none,linecolor=blue]{ 0.5669885494214224 , 0.39833927419126536 }
\dspPZ[label=none,linecolor=blue]{ 0.5669885494214224 , -0.39833927419126536 }
\end{dspPZPlot}
\end{frame}
\begin{frame} \frametitle{Pole magnitude}
\centering
Magnitude of poles as a function of the number of digits used to store coefficients
\vspace{2em}
\begin{tabular}{c||c|c|c|c|}
\# digits \\ \hline
9 & 0.99969893 & 0.99641971 & 0.96231223 & 0.6929287 \\ \hline
8 & 0.99970234 & 0.99641583 & 0.96231266 & 0.69292873 \\ \hline
7 & 0.99987231 & 0.99622669 & 0.96233196 & 0.69292855 \\ \hline
6 & \color{red} 1.0027213 & 0.99267273 & 0.96304264 & 0.69292212 \\ \hline
5 & \color{red} 1.00418091 & 0.99647046 & 0.95797945 & 0.69292331 \\ \hline
\end{tabular}
\end{frame}
\begin{frame} \frametitle{Numerical precision: how to mitigate}
\begin{itemize}
\item design filter in factored form
\item use a cascade of second-order sections
\item in Python: \begin{tt}
b, a = sp.ellip(4, .1, 50, 0.25, output='sos')
\end{tt}
\end{itemize}
\end{frame}
\end{document}

Event Timeline

c4science · Help