This file is larger than 256 KB, so syntax highlighting was skipped.
diff --git a/teaching/COM303/2017/lectures/11 - z-trans and more filters/0_ztransform.tex b/teaching/COM303/2017/lectures/11 - z-trans and more filters/0_ztransform.tex
index d425c54..5c80481 100644
--- a/teaching/COM303/2017/lectures/11 - z-trans and more filters/0_ztransform.tex
+++ b/teaching/COM303/2017/lectures/11 - z-trans and more filters/0_ztransform.tex
@@ -1,556 +1,649 @@
\documentclass[aspectratio=169]{beamer}
\def\stylepath{../styles}
\usepackage{pst-3dplot}
\usepackage{\stylepath/com303}
\begin{document}
\begin{frame}
\frametitle{Overview:}
\begin{itemize}[<+->]
\item Constant-Coefficient Difference Equations
\item The $z$-transform
\item System transfer function
\item Region of convergence and system stability
\end{itemize}
\end{frame}
\begin{frame} \frametitle{The most general LTI systems}
\begin{itemize}[<+->]
\item ideal filters cannot be implemented
\item what is the most general, realizable LTI transformation?
\begin{itemize}
\item linearity: only sums and multiplications
\item time-invariance: only multiplications by constants
\item realizability: only finite number of past and future samples
diff --git a/teaching/COM303/2017/lectures/11 - z-trans and more filters/1_structures.tex b/teaching/COM303/2017/lectures/11 - z-trans and more filters/1_structures.tex
index 91d8ccc..8708d80 100644
--- a/teaching/COM303/2017/lectures/11 - z-trans and more filters/1_structures.tex
+++ b/teaching/COM303/2017/lectures/11 - z-trans and more filters/1_structures.tex
\item use Chebyshev polynomials to transform $H_d(e^{j\omega})$ into $P(x)$
\item rewrite specifications to obtain global error function $E(x)$ over passband + stopband
\item iteratively adjust extrema of the error function
\item use alternation theorem as stopping condition
\item verify stopband and passband error; increase $M$ and repeat if necessary
\end{itemize}
\end{frame}
\begin{frame}
\frametitle{Alternation Theorem}
Ingredients:
\begin{itemize}
\item $\{I_k \}$ set of disjoint intervals on the real line, $I = \bigcup_k I_k$
\item $P(x)$ polynomial of degree $L$,
\item $D(x)$ target function
\item $W(x)$ positive weighting function
\item $E(x) = W(x)[D(x) - P(x)]$ approximation error function
\item $E_{\max} = \max_{x \in I} \{|E(x)|\}$
\end{itemize}
\pause
\vspace{1em}
$P(x)$ is the \emph{unique} order-$L$ polynomial which minimizes $E_{\max}$ if and only if there exist \emph{at least} $L+2$ successive values $x_i$ in $I$ such that $|E(x_i)| = E_{\max}$ and $E(x_i) = -E(x_{i+1})$.
\item use Chebyshev polynomials to transform $H_d(e^{j\omega})$ into $P(x)$
\item rewrite specifications to obtain global error function $E(x)$ over passband + stopband
\item iteratively adjust extrema of the error function
\item use alternation theorem as stopping condition
\item verify stopband and passband error; increase $M$ and repeat if necessary
\end{itemize}
\end{frame}
\begin{frame}
\frametitle{Alternation Theorem}
Ingredients:
\begin{itemize}
\item $\{I_k \}$ set of disjoint intervals on the real line, $I = \bigcup_k I_k$
\item $P(x)$ polynomial of degree $L$,
\item $D(x)$ target function
\item $W(x)$ positive weighting function
\item $E(x) = W(x)[D(x) - P(x)]$ approximation error function
\item $E_{\max} = \max_{x \in I} \{|E(x)|\}$
\end{itemize}
\pause
\vspace{1em}
$P(x)$ is the \emph{unique} order-$L$ polynomial which minimizes $E_{\max}$ if and only if there exist \emph{at least} $L+2$ successive values $x_i$ in $I$ such that $|E(x_i)| = E_{\max}$ and $E(x_i) = -E(x_{i+1})$.
"The Parks-McClellan (PM) algorithm is a clever application of advanced polynomial fitting techniques to the problem of FIR filter design. In this notebook we will explore the key ideas behind the method by considering the design of a simple Type-I lowpass filter."
]
},
{
"cell_type": "code",
- "execution_count": 1,
+ "execution_count": 13,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# standard bookkeeping\n",
"%matplotlib inline\n",
"import matplotlib.pyplot as plt\n",
"import numpy as np\n",
"from IPython.display import Audio, display\n",
"from ipywidgets import interactive, fixed\n",
"plt.rcParams[\"figure.figsize\"] = (14,4)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The first intuition is to recognize that finding the filter's coefficients is equivalent to solving a polynomial fitting problem. Since we want a linear-phase filter, the filter will be symmetric and we can easily show that in this case its frequency response is the real-valued function:\n",
"In the above expression, the $N = 2M+1$ nonzero taps of the impulse response $h[n]$ are $h[-M]$ to $h[M]$. Given a (positive) passband $[0, \\omega_p]$ and a stopband $[\\omega_s, \\pi]$, we need to fit $H(e^{j\\omega})$ to one in the passband and zero in the stopband as in the following figure. \n",
"Enter [Chebishev polynomials](https://en.wikipedia.org/wiki/Chebyshev_polynomials) (CP): they have a ton of interesting properties but the relevant one here is that $T_n(\\cos\\omega) = \\cos n\\omega$, where $T_n(x)$ is the CP of order $n$. Since the coefficients for the CP are easy to determine, with a simple substitution we can write:\n",
"While the relation between the $M+1$ $p_n$ coefficients and the $2M+1$ filter coefficients is nonlinear, it is easily invertible and therefore we can just concentrate on fitting the degree-$M$ polynomial over the desired response."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"For the sake of simplicity, let's skip some of the details associated to the Chebishev substitution (for instance, the mapping $x=\\cos\\omega$ reverses the abscissae) and let's examine the following equivalent problem:\n",
"\n",
"> consider the desired function: $$ D(x) = \\begin{cases} 1 & \\mbox{for $x\\in [0, A]$} \\\\ 0 & \\mbox{for $x\\in [B,1]$} \\end{cases}$$ <br> find the coefficients of a degree-$M$ polynomial $P(x)$ so that $P(x)\\approx D(x)$ over $[0,A]\\cup [B,1]$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Standard polynomial fitting\n",
"\n",
"The simplest approach is to use standard polynomial fitting: choose $a$ points in $[0, A]$ and $b$ points in $[B, 1]$ so that $a+b = M+1$ and find an interpolator over the coordinate pairs:\n",
"The result will minimize the mean square error between the interpolator and the piecewise characteristic we are trying to approximate. \n",
"\n",
"We can write a simple Python function to test this approach; you can play with the sliders below and see the result of the plain interpolation as the order and the size of the intervals change:"
]
},
{
"cell_type": "code",
- "execution_count": 2,
+ "execution_count": 14,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"def MSE_fit(A, B, order):\n",
" if order < 3:\n",
" raise ValueError(\"order is too small, we need at least 3\")\n",
" \n",
" # interpolation points always one more than the order of the interpolator\n",
" pts = order+1\n",
" # split number of interpolation points across intervals proportionately \n",
" # with the length of each interval \n",
" ptsA = int(pts * A / (A+(1-B)))\n",
" if ptsA < 2:\n",
" ptsA = 2\n",
" ptsB = pts - ptsA\n",
" \n",
" # for the MSE fit, place a point at each interval edge and distribute the rest\n",
" # (if any) evenly over the interval\n",
" x = np.concatenate((\n",
" np.arange(0, ptsA) * (A / (ptsA-1)),\n",
" B + np.arange(0, ptsB) * ((1-B) / (ptsB-1)) \n",
" ))\n",
" y = np.concatenate((\n",
" np.ones(ptsA),\n",
" np.zeros(ptsB) \n",
" ))\n",
" \n",
" # now just solve the linear interpolation problem\n",
"As you can see, simple polynomial interpolation, while minimizing the MSE has two problems:\n",
"\n",
" * it becomes numerically unstable as soon as the order of the interpolation exceeds 16 or 17\n",
" * although the MSE is minimized, the **maximum** error can become very large\n",
" \n",
"Because of these problems, direct interpolation is rarely used in numerical analysis and filter design is no exception. "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Minimax fitting\n",
"\n",
"As we said, the first clever intuition behind the the Parks-McClellan algorithm is rephrasing the filter design problem as a polynomial interpolation. The real stroke of genius, however, is the use of a different kind of polynomial fitting called *minimax* approximation. In this kind of fitting the goal is to minimixe the *maximum* error between the polynomial and the desired function over the intervals of interest:\n",
"\n",
"$$\n",
" E = \\min\\max_{[0,A]\\cup [B,1]} | P(x) - D(x) |\n",
"$$\n",
"\n",
"This is brilliant for two reasons: \n",
"\n",
"* the minimax criterion gives us a guarantee on the worst-case error for the filter response\n",
"* an obscure yet powerful theorem, called the **alternation theorem**, gives us a remarkably straightforward recipe to build a robust numerical algorithm that finds the solution."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### The Alternation Theorem\n",
"\n",
"Because of the fundamental theorem of algebra, a non-constant polynomial cannot be constant over an interval. Since our polynomial cannot be a constant (it needs to move from 1 in the passband to zero in the stopband), it will necessarily oscillate over the approximation intervals. As you could see from the previous demo, MSE minimization tries to keep the oscillations small over the approximation intervals, but the price to pay is potentially wide oscillations at the band edges; on the other hand, the minimax approximation will allow for oscillations that are larger overall but that do not swing wildly. Since the polynomial oscillates around the target value, the error will oscillate between positive and negative peaks; the alternation theorem states that \n",
"\n",
"> $P(x)$ is the minimax approximation to $D(x)$ if and only if $P(x) - D(x)$ alternates $M+2$ times between $+E$ and $-E$ over $[0,A]\\cup [B,1]$\n",
"\n",
"The alternation theorem gives one thing right away: the ability to recognize if a polynomial is the minimax solution. All we need to do is look at the extrema of the error and check that \n",
"\n",
"* they are $M+2$\n",
"* they alternate in sign\n",
"* their magnitude is exactly the same\n",
"\n",
"We can write a simple Python function to find the extrema of the error:"
]
},
{
"cell_type": "code",
- "execution_count": 5,
+ "execution_count": 17,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"def find_error_extrema(p, A, B):\n",
" intervals = {\n",
" (0, A): 1, \n",
" (B, 1): 0\n",
" }\n",
" loc = [] # locations of the extrema\n",
" err = [] # values of the extrema\n",
" for rng, val in intervals.iteritems():\n",
" # we don't need enormous precision, 100 points per interval will do\n",
" t = np.linspace(rng[0], rng[1], 100)\n",
" y = val - p(t) # error values\n",
" # this finds all the points where the error changes sign:\n",
" ix = np.diff(np.sign(np.diff(y))).nonzero()[0] + 1 # local min+max\n",
" loc = np.hstack((loc, t[0], t[ix], t[-1]))\n",
" err = np.hstack((err, y[0], y[ix], y[-1]))\n",
" return loc, err"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"With this, it's easy to verify that the MSE fit does not satisfy the alternation theorem: the magnitude of the peaks is not constant."
"The alternation theorem provides us also with a very efficient way to find the coefficients of $P(x)$. Suppose we knew the exact ordinates $x_i$ of the $M+2$ alternations; in this case we could solve the following system of linear equations\n",
"and find in one go both the $M+1$ polynomial coefficients *and* the value of the minimax error $E=|\\epsilon|$ (we use $\\epsilon$ instead of $E$ in the linear system because we don't know the sign of the first alternation, it could be positive or negative). Although the above is a non-standard linear system of equations, it can be shown rather easily that, as long as the $x_i$ are distinct, it does have a solution. \n",
"\n",
"Of course we don't know the $x_i$ in advance but we can start with a guess, and solve the system anyway. Once we find the polynomial coefficients from the guess, we use Remez's exchange algorithm:\n",
"\n",
"* find the locations for the maxima and minima of the error for the $P(x)$ we just found\n",
"* if the extrema satisfy the alternation theorem, we are done\n",
"* otherwise, use the new locations as a new guess, solve the system again and repeat.\n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {
"collapsed": true
},
"source": [
"The Remez algorithm is remarkably fast and robust. Here is an implementation you can play with. First, we need an auxiliary function to solve the system of equations above; we will use standard linear algebra functions."
]
},
{
"cell_type": "code",
- "execution_count": 22,
+ "execution_count": 20,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"def solve(x, y):\n",
" # simple solver for the extended interpolation problem\n",
" # first build a Vandermonde matrix using the interpolation locations\n",
" # There are N+2 locations, so the matrix will be (N+2)x(N+2) but we \n",
" # will replace the last column with the error sign\n",
" V = np.vander(x, increasing=True)\n",
" # replace last column\n",
" V[:,-1] = pow(-1, np.arange(0, len(x)))\n",
" # just solve Ax = y\n",
" v = np.linalg.solve(V, y)\n",
" # need to reverse the vector because poly1d starts with the highest degree\n",
" p = np.poly1d(v[-2::-1])\n",
" e = np.abs(v[-1])\n",
" return p, e"
]
},
{
"cell_type": "markdown",
"metadata": {
"collapsed": true
},
"source": [
"And here comes the main course: the Remez routine. The code is quite straightforward; it doesn't have a termination conditions since the number of iterations is passed as a parameter (we want to be able to show intermediate results)."
]
},
{
"cell_type": "code",
- "execution_count": 23,
+ "execution_count": 21,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"def remez_fit(A, B, order, iterations):\n",
" if order < 3:\n",
" raise ValueError(\"order is too small, we need at least 3\")\n",
" pts = order+2\n",
" \n",
" # initial choice of interpolation points: distribute them evenly\n",
" # across the two regions as a proportion of each region's width\n",
" ptsA = int(pts * A / (A-B+1))\n",
" if ptsA < 2:\n",
" ptsA = 2\n",
" ptsB = pts - ptsA\n",
" \n",
" x = np.concatenate((\n",
" np.arange(1, ptsA+1) * (A / (ptsA+1)),\n",
" B + np.arange(1, ptsB+1) * ((1-B) / (ptsB+1)) \n",
" ))\n",
" y = np.concatenate((\n",
" np.ones(ptsA),\n",
" np.zeros(ptsB) \n",
" ))\n",
" \n",
" # the \"data\" dictionary only holds values that we will use in plotting\n",
" data = {}\n",
" \n",
" for n in range(0, iterations):\n",
" # previous interpolation points\n",
" data['prev_x'] = x\n",
" data['prev_y'] = y\n",
" \n",
" # solve the interpolation problem \n",
" p, e = solve(x, y)\n",
" data['err'] = e\n",
" # find the extrema of the error\n",
" loc, err = find_error_extrema(p, A, B) \n",
" \n",
" # find the alternations \n",
" alt = []\n",
" for n in range(0, len(loc)):\n",
" # each extremum is a new candidate for an alternation\n",
" c = {\n",
" 'loc': loc[n],\n",
" 'sign': np.sign(err[n]),\n",
" 'err_mag': np.abs(err[n])\n",
" }\n",
" # only keep extrema that are larger or equal than the minimum\n",
" # error returned by the interpolation solution\n",
" if c['err_mag'] >= e - 1e-3:\n",
" # ensure that the error alternates; if the candidate has the \n",
" # same sign, replace the last alternation with the candidate\n",
" # if its error value is larger\n",
" if alt == [] or alt[-1]['sign'] != c['sign']:\n",
" alt.append(c)\n",
" elif alt[-1]['err_mag'] < c['err_mag']:\n",
" alt.pop()\n",
" alt.append(c)\n",
" \n",
" # if there are more than the necessary number of alternations, trim \n",
" # from the left or the right keeping the largest errors\n",
" while len(alt) > order + 2:\n",
" if alt[0]['err_mag'] > alt[-1]['err_mag']:\n",
" alt.pop(-1)\n",
" else:\n",
" alt.pop(0)\n",
" \n",
" # the new set of interpolation points\n",
" x = [c['loc'] for c in alt]\n",
" y = [1 if c <= A else 0 for c in x]\n",
" data['new_x'] = x\n",
"\n",
" return p, data"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Finally, a simple auxiliary function to plot the results; the yellow dots indicate the guess used for the current interpolation while the blue stars show the new maxima that will be used as the new guess. As you can see the algorithm converges very rapidly. The error in passband and stopband is shown magnified in the bottom panels."
"In the example above we used a very simple piecewise constant target function $D(x)$ but the Alternation Theorem is actually much more general and the Remez algorithm can be used to approximate more complicated characteristics. \n",
"\n",
"In particular, a *weighting function* can be used in the minimax problem formulation in order to tune the error in the approximation regions. Think of minimax filter design as a budgeting problem: for a given number of coefficients, we need to make sure that the bandpass and bandstop reqirements are fulfilled first, while the error is a consequence of the budget. With a weighting function we can decide to allocate a larger part of the budget to either passband or stopband. \n",
"\n",
"Mathematically, the minimax problem becomes:\n",
"\n",
"$$\n",
" E = \\min\\max_{[0,A]\\cup [B,1]} | W(x)[P(x) - D(x)] |\n",
"$$\n",
"\n",
"but the alternation theorem still holds:\n",
"\n",
"> $P(x)$ is the minimax approximation to $D(x)$ *weighted by $W(x)$* if and only if $W(x)[P(x) - D(x)]$ alternates $M+2$ times between $+E$ and $-E$ over $[0,A]\\cup [B,1]$\n",
"\n",
"For instance, suppose we want the error in the first interval to be 10 times smaller than in the second interval. In this case the weighting function will be equal to 0.1 over $[0,A]$ and just one over $[B,1]$. We can rewrite the extended interpolation problem as\n",
" Your browser does not support the audio element.\n",
" </audio>\n",
" "
],
"text/plain": [
"<IPython.lib.display.Audio object>"
]
},
- "execution_count": 2,
+ "execution_count": 4,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"SF, s = wavfile.read('sm.wav')\n",
"IPython.display.Audio(s, rate=SF)"
]
},
{
"cell_type": "markdown",
"metadata": {
"collapsed": false
},
"source": [
"Suppose I'm interested in the bass part; I know the bass will have most of its energy in the low frequency range. I could try to filter the audio with a lowpass with cutoff, say, 200Hz and see what comes out of it (you will need headphones to hear properly).\n",
"\n",
"Remember that the clock of the system is given by the sampling rate of the audio clip, so we will need to normalize the frequencies by $SF/2$ so that the highest frequency corresponds to 1. \n",
" Your browser does not support the audio element.\n",
" </audio>\n",
" "
],
"text/plain": [
"<IPython.lib.display.Audio object>"
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# let's filter the signal and hear the result\n",
"y = sp.lfilter(b, a, s)\n",
"IPython.display.Audio(data = y, rate = SF)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Well, it's certainly not a perfect source separation but, if you're trying to learn the bass part, this would be a good start. \n",
"\n",
"We can try to achieve the same effect using FIR filters, and you'll see it's possible but we will need a lot of taps since the bandwidth is small. We can use the built-in normalization facility of Scipy's `remez()` function:"
"As you can see, the FIR filter will be very expensive! Is it worth it? Well, you can judge for yourself: arguably linear phase preserves the instrumental attack more, but the high frequency rejection is not as good."
" Your browser does not support the audio element.\n",
" </audio>\n",
" "
],
"text/plain": [
"<IPython.lib.display.Audio object>"
]
},
"execution_count": 7,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"y = sp.lfilter(h, 1, s)\n",
"IPython.display.Audio(y, rate=SF)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Just like we tried to extract the bass, we can try to extract parts of the drum patterd. Usually, we get a good feel for the hi-hat and cymbals by keeping frequencies above 7KHz.\n",
"\n",
"Let's use an FIR highpass; note that, to design a highpass, we choose a type II filter, so the length will be odd"