{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"# Exercises from Linear algebra, signal processing, and wavelets. A unified approach.\\\\ MATLAB version\n",
"\n",
" **Øyvind Ryan**\n",
"\n",
"Date: **Oct 24, 2018**\n",
"\n",
"\n",
"\n",
"\n",
"# Discrete time filters\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"## Exercise 3.6: Replacing a sparse circular convolution with shorter ones\n",
"\n",
""
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"[x, f_s] = audioread('sounds/castanets.wav');\n",
"k = 5;"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Let $\\boldsymbol{s},\\boldsymbol{x}\\in\\mathbb{R}^N$, and assume that $\\boldsymbol{s}$ has at most $k$ nonzero entries, gathered in one segment.\n",
"\n",
"\n",
"**a)**\n",
"Show that there exists a $\\tilde{\\boldsymbol{s}}\\in\\mathbb{R}^{M+k-1}$ and a number $a$ so that, for any $r$,"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"(\\boldsymbol{s}\\circledast\\boldsymbol{x})_{r,r+1...,r+M-1} = (\\tilde{\\boldsymbol{s}}\\circledast \\tilde{\\boldsymbol{x}})_{0,...,M-1},\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"where"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\tilde{\\boldsymbol{x}} = (x_{a+r},x_{a+r+1},...,x_{a+r+M+k-2}).\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In other words, any set of $M$ consecutive entries in $\\boldsymbol{s}\\circledast\\boldsymbol{x}$ can be obtained from a circular convolution of size $M+k-1$. \n",
"Thus, a circular convolution of size $N$ can be computed as $\\frac{N}{M}$ circular convolutions of size $M+k-1$.\n",
"\n",
"\n",
"\n",
"**Solution.**\n",
"Denote by $S$ the $N\\times N$ circulant Toeplitz matrix with first column equal to $\\boldsymbol{s}$.\n",
"Assume that the nonzero entries in the first row of $S$ occur in the interval $[a,a+k-1]$.\n",
"Then the nonzero entries in rows $r,...,r+M-1$ of $S$ are in the columns from $a+r$ to $a+r+M-1+k-1=a+r+M+k-2$, i.e. the interval $a+[r,r+M+k-2]$ (there are $M+k-1$ such columns). With"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\tilde{\\boldsymbol{x}} = (x_{a+r},x_{r+1},...,x_{a+r+M+k-2}),\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"and viewing the submatrix of $S$ from rows $r,...,r+M-1$ and columns from $a+[r,r+M+k-2]$ as the first $M$ rows of a $(M+k-1)\\times(M+k-1)$-circulant Toeplitz matrix $\\tilde{S}$ with first column $\\tilde{\\boldsymbol{s}}$, we obtain immediately that"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"(\\boldsymbol{s}\\circledast\\boldsymbol{x})_{r,r+1...,r+M-1} = (\\tilde{\\boldsymbol{s}}\\circledast \\tilde{\\boldsymbol{x}})_{0,...,M-1}.\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"**b)**\n",
"It is natural to implement the shorter circular convolutions in terms of the FFT. \n",
"Use this strategy with $M+k-1=2^r$ to find the number of operations needed to compute the $\\frac{N}{M}$ circular convolutions from a), \n",
"and write a program which finds the $r$ which gives the lowest number of operations. \n",
"Create a plot where $k$ is plotted against the number of operations corresponding to the optimal $r$, \n",
"and also against the number of operations required by a direct implementation of circular convolution.\n",
"\n",
"\n",
"\n",
"**Solution.**\n",
"Recall that an $N$-point split-radix FFT requires $4N\\log_2N-6N+8$ arithmetic operations. \n",
"With $M+k-1=2^r$ for $N$ this gives $(4r-6)2^r+8$ operations for each of the shorter convolutions. \n",
"A $2^r$-point convolution implemented in terms of the FFT (or rather FFT, multiplication in frequency, and IFFT) thus requires \n",
"\n",
"\\[\n",
"2((4r-6)2^r+8) + 6\\cdot 2^r = (8r-6)2^r+16\n",
"\\] \n",
"operations. The factor $2$ comes from computing both an FFT and an IFFT, and the summand $6\\cdot 2^r$ comes\n",
"from the multiplications in frequency (each complex multiplication requires 4 real multiplications and 2 real additions, a total of 6 operations \n",
"for each of the $2^r$ entries). \n",
"In total for all the $N/M$ (or rather $\\lceil N/M\\rceil$) shorter convolutions we thus need $((8r-6)2^r+16)\\left\\lceil\\frac{N}{2^r-k+1}\\right\\rceil$ \n",
"operations. In a direct implementation, for each of the $N$ entries we need to make $k$ complex multiplications and $k-1$ complex additions. \n",
"Since each complex addition requires $2$ real additions, this totals $6k+2(k-1)=8k-2$ operations for each entry, so that the overall total \n",
"is $(8k-2)N$ operations. \n",
"\n",
"The following code finds the $r$ which gives the minimum of $((8r-6)2^r+16)/(2^r-k+1)$ for a fixed $k$, \n",
"and then stores that minimum value. This is then repeated for all the $k$, and compared with $8k-2$, as we deduced for a direct implementation."
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"kvals = 1:50;\n",
"opt_vals = kvals;\n",
"for k=kvals\n",
" rvals = ceil(log2(k)):50;\n",
" num_ops = ((8*rvals-6).*2.^rvals + 16)./(2.^rvals-k+1);\n",
" opt_vals(k) = min(num_ops); \n",
"end\n",
"plot(kvals,8*kvals-2,kvals,opt_vals)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Although this code does not plot the optimal value for $r$, this value is very much needed as a guideline for implementations on how to split \n",
"convolutions into shorter convolutions.\n",
"The reason that we iterate $r$ from $\\log_2 k$ in the code is that $k \\leq 2^r$ secures that there is room for all the filter coefficients inside the vector $\\tilde{\\boldsymbol{s}}$ used in the smaller convolution. \n",
"From the resulting plot we see that when the number of filter coefficients gets above $5$, the FFT-based approach works better. \n",
"When the number of filter coefficients is below this, a direct implementation of the filter should be chosen. \n",
"\n",
"When the number of points in an FFT gets high, one usually not cares about the terms $-6N+8$ from the exact operation count of the Split radix FFT, \n",
"since they are not dominant terms. \n",
"In this exercise these terms really mattered, however, since short FFT's are in focus. \n",
"If you had dropped the terms $-6N+8$, you can verify that $((8r-6)2^r+16)/(2^r-k+1)$ should be replaced by $(8r+6)2^r/(2^r-k+1)$,\n",
"and that a new plot shows that at least 8 filter coefficients are needed in order for the FFT-based approach to outperform the direct implementation.\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"## Exercise 3.8: Finding sparse solutions to under-determined linear systems\n",
"\n",
"\n",
"\n",
"Let $\\boldsymbol{x}\\in\\mathbb{C}^N$ be a vector where we know that at most $s$ components are nonzero. When $s$ is small such vectors are called *sparse*.\n",
"Denote the indices of these nonzero components by $S$, and let $\\boldsymbol{y}=DFT_N\\boldsymbol{x}$. \n",
"In this exercise we will find a procedure for finding $\\boldsymbol{x}$ from the values $y_0,...,y_{2s-1}$, i.e. we can recover all components of a sparse vector from an under-determined set of measurements. \n",
"This procedure is also called a *Reed-Solomon code*. \n",
"In some sense this parallels the sampling theorem (which in a similar way fills in the values between the samples), with an assumption on the highest frequency replaced by an assumption on sparsity.\n",
"Many other results also exist on how to recover sparse vectors, and these results give rise to the field of *compressed sensing*. See [[foucart]](#foucart) for a review on this field.\n",
"\n",
"\n",
"**a)**\n",
"Let $\\boldsymbol{z}$ be the vector with components $z_k=\\frac{1}{N}\\prod_{n\\in S} \\left(1-e^{-2\\pi in/N}e^{2\\pi ik/N}\\right)$. \n",
"Show that $z_k=0$ for $k\\in S$, and also that $\\boldsymbol{x}\\circ\\boldsymbol{z}=\\boldsymbol{0}$.\n",
"\n",
"\n",
"\n",
"**Solution.**\n",
"This follows from that $1-e^{-2\\pi ik/N}e^{2\\pi ik/N}=1-1=0$ is a factor in $z_k$ when $k\\in S$. \n",
"Since $x_k=0$ when $k\\not\\in S$, and $z_k=0$ when $k\\in S$, we have for all $k$ that $(\\boldsymbol{x}\\circ\\boldsymbol{z})_k=x_kz_k=0$.\n",
"\n",
"\n",
"\n",
"**b)**\n",
"Let $\\boldsymbol{w}=DFT_N\\boldsymbol{z}$. Show that $w_0=1$ and that $w_n=0$ for $n>s$. \n",
"Conclude from this, that $\\boldsymbol{x}\\circ\\boldsymbol{z}=\\boldsymbol{0}$, and due the previous exercise that $w_1,w_2,...,w_s$ \n",
"must fulfill the equations"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\begin{array}{ccccccccc}\n",
"y_s &+& w_1y_{s-1} &+& \\cdots &+& w_sy_0 &=& 0 \\\\ \n",
"y_{s+1} &+& w_1y_s &+& \\cdots &+& w_sy_1 &=& 0 \\\\ \n",
"\\vdots & & \\vdots & & \\ddots & & \\vdots &=& 0 \\\\ \n",
"y_{2s-1} &+& w_1y_{2s-2} &+& \\cdots &+& w_sy_{s-1} &=& 0\n",
"\\end{array}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"**Solution.**\n",
"If we multiply out the factors in $z_k$ we obtain $1+\\sum_{r=1}^{s} c_r e^{2\\pi irk}$ for some values $c_r$. \n",
"Since also $\\boldsymbol{z}=\\sum_{n=1}^N w_ne^{2\\pi ikn/N}$ by the definition of the DFT, we see that $w_0=1$ and that \n",
"$w_n=0$ for $n>s$. From the previous exercise, multiplication in time corresponds to convolution in frequency, \n",
"so that $\\boldsymbol{y}\\circledast\\boldsymbol{w}=\\boldsymbol{0}$. The stated equation system corresponds to"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"(\\boldsymbol{y}\\circledast\\boldsymbol{w})_s=\\cdots=(\\boldsymbol{y}\\circledast\\boldsymbol{w})_{2s-1}=0.\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"where we inserted $w_0=1$, $w_n=0$ for $n>s$.\n",
"\n",
"\n",
"\n",
"We now have $s$ equations with $s$ unknowns. By construction this has a solution, \n",
"but there may not be a unique solution. In the last part of the exercise we will show how to find $\\boldsymbol{x}$, \n",
"regardless of which solution we choose. \n",
"\n",
"**c)**\n",
"Assume that $v_1,v_2,...,v_s$ is any solution to the system in b), and extend this to a vector $\\boldsymbol{v}\\in\\mathbb{R}^N$ by setting $v_0=1$, $v_{s+1}=\\cdots=v_{N-1}=0$. \n",
"Show that $(\\boldsymbol{v}\\circledast\\boldsymbol{y})_n=0$ for $s\\leq n\\leq 2s-1$, and conclude from the previous exercise that \n",
"$\\left(DFT_N\\left((IDFT_N\\boldsymbol{v})\\circ\\boldsymbol{x}\\right)\\right)_n=0$ for $s\\leq n\\leq 2s-1$.\n",
"\n",
"\n",
"\n",
"**Solution.**\n",
"For $s\\leq n\\leq 2s-1$ we have that"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"(\\boldsymbol{v}\\circledast\\boldsymbol{y})_n = \\sum_{k=0}^{N-1} v_k y_{(n-k)\\bmod N} = y_n + \\sum_{k=1}^{s} v_k y_{(n-k)\\bmod N} = 0.\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Since $DFT_N\\left((IDFT_N\\boldsymbol{v})\\circ\\boldsymbol{x}\\right) = \\frac{1}{N}\\boldsymbol{v}\\circledast\\boldsymbol{y}$ from the previous exercise, \n",
"the result follows.\n",
"\n",
"\n",
"\n",
"The vector $(IDFT_N\\boldsymbol{v})\\circ\\boldsymbol{x}$ has at most $s$ nonzero components (since $\\boldsymbol{x}$ has). \n",
"If we take the columns from $DFT_N$ with indices in $S$, and the rows with indices between $s$ and $2s-1$,\n",
"the resulting $s\\times s$-matrix was shown in Exercise 2.17 to be non-singular. \n",
"\n",
"\n",
"**d)**\n",
"Explain that $((IDFT_N\\boldsymbol{v})\\circ\\boldsymbol{x})_n=0$ for $n\\in S$, and conclude that $(IDFT_N\\boldsymbol{v})\\circ\\boldsymbol{x}=\\boldsymbol{0}$.\n",
"\n",
"\n",
"\n",
"**Solution.**\n",
"The stated non-singular matrix applied to the components of $(IDFT_N\\boldsymbol{v})\\circ\\boldsymbol{x}$ from $S$ produces zero. \n",
"These components must thus also be zero, so that $(IDFT_N\\boldsymbol{v})\\circ\\boldsymbol{x}$ have all components in $S$ equal to $0$. Since $\\boldsymbol{x}$ is zero outside $S$, \n",
"it follows that $(IDFT_N\\boldsymbol{v})\\circ\\boldsymbol{x}=\\boldsymbol{0}$.\n",
"\n",
"\n",
"\n",
"**e)**\n",
"In particular we must have that $(IDFT_N\\boldsymbol{v})_k=0$ for $k\\in S$, so that"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\frac{1}{N}\\sum_{n=0}^{N-1} v_ne^{2\\pi ikn/N} = \\frac{1}{N}\\sum_{n=0}^s v_ne^{2\\pi ikn/N} = 0\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"for $k\\in S$. This is a polynomial in $e^{2\\pi ik/N}$ of degree at most $s$. Explain how one can find the set $S$ from this polynomial.\n",
"\n",
"\n",
"\n",
"**Solution.**\n",
"Simply apply an IDFT, and choose $S$ as the set from the result with smallest coefficients. Due to roundoff errors \n",
"the coefficients will not be exactly zero, however.\n",
"\n",
"\n",
"\n",
"**f)**\n",
"Explain how you can recover $\\boldsymbol{x}$ from $y_0,...,y_{2s-1}$ when the set $S$ is known.\n",
"\n",
"\n",
"\n",
"**Solution.**\n",
"The submatrix $A$ from $\\text{DFT}_N$ with columns from $S$, and rows from $0,...,s-1$ is again non-singular. \n",
"The values of $\\boldsymbol{x}$ on $S$, denoted $\\boldsymbol{x}_s$ can be found by solving $A\\boldsymbol{x}_S = (y_0,...,y_{s-1})$.\n",
"\n",
"\n",
"\n",
"**g)**\n",
"Write a function"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
" recover_x(y, N)\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"which uses the procedure obtained by combining b), e), and f) to compute and return $\\boldsymbol{x}\\in\\mathbb{R}^N$ under the assumption that $\\boldsymbol{x}$ has at most $s$ nonzero components, \n",
"and that `y` contains the first $2s$ components in $\\boldsymbol{y}=DFT_N\\boldsymbol{x}$. \n",
"Test that the code works correctly on a vector $\\boldsymbol{x}$ of your own choosing.\n",
"\n",
"\n",
"\n",
"**Solution.**\n",
"The function can be implemented as follows."
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"%%file recover_x.m\n",
"function x=recover_x(y, N)\n",
" s= size(y,1)/2;\n",
" A=zeros(s);\n",
" for k=1:s\n",
" A(:,s-k+1) = y(k:(k+s-1));\n",
" end\n",
" \n",
" b=-y((s+1):(2*s));\n",
" w=A\\b;\n",
" v = [1; w; zeros(N-(s+1),1)];\n",
" inds = find(abs(ifft(v)) <= 0.00001);\n",
" x_S = exp(-2*pi*i*((0:(s-1))')*(inds-1)'/N)\\y(1:s);\n",
" x = zeros(N,1);\n",
" x(inds)=x_S;\n",
"end"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The following code test that the function manages to recover a sparse vector with only four nonzero integer entries, from the $8$ first Fourier coefficients."
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"N=32;\n",
"x = zeros(N,1);\n",
"x([3,5,9,12])=[1,4,3,2];\n",
"y=fft(x);\n",
"x=recover_x(y(1:8),N)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"## Exercise 3.9: Windows\n",
"\n",
"\n",
"\n",
"We mentioned in Section 2.2 that we obtain some undesirable effects in the frequency representation when we restrict to a block of the signal. Assume that $\\boldsymbol{x}\\in\\mathbb{R}^M$, and that"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\boldsymbol{w}=\\{\\underline{w_0},\\ldots,w_{N-1}\\} \\text{ with } N < M.\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We call $(w_0x_0,\\ldots,w_{N-1}x_{N-1},0,...,0)\\in\\mathbb{R}^M$ a *windowed signal* and $\\boldsymbol{w}\\in\\mathbb{R}^N$ a *window* of length $N$.\n",
"\n",
"\n",
"**a)**\n",
"Use Exercise 3.7 to show that the DFT of the windowed signal is"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\begin{align*}\n",
"\\frac{1}{M}(DFT_M(w_0,...,w_{N-1},0,...,0))\\circledast (DFT_M\\boldsymbol{x})\n",
"\\end{align*}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"**Solution.**\n",
"With $\\boldsymbol{s}_1=(w_0,...,w_{N-1},0,...,0)$, $\\boldsymbol{s}_2=\\boldsymbol{x}$ we have that $\\boldsymbol{s}_1 \\circ \\boldsymbol{s}_2$ is the windowed signal. The result now follows from the second property in Exercise 3.7.\n",
"\n",
"\n",
"\n",
"If $DFT_M(w_0,w_1,...,w_{N-1},0,...,0)$ is close to $(M,0,...,0)$, this will be close to $\\boldsymbol{y}=DFT_M\\boldsymbol{x}$. \n",
"In other words, a good window should satisfy"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"DFT_M(w_0,w_1,...,w_{N-1},0,...,0)\\approx(M,0,...,0).\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"There is a loss when we pass from the signal to the windowed signal, since we can only construct a DFT close to, not equal to, $(M,0,...,0)$. \n",
"We will not go into techniques for how to find values $w_i$ which are close to satisfying this, only evaluate three such in the rest of the exercise. \n",
"\n",
"\n",
"\n",
"**b)**\n",
"The *rectangular window* is defined by $\\boldsymbol{w}=\\{\\underline{1},1,\\ldots,1\\}$. Show that"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"(DFT_M(1,...,1,0,...,0))_n = \\sum_{k=0}^{N-1} e^{-2\\pi ikn/M}=\\frac{1-e^{-2\\pi inN/M}}{1-e^{-2\\pi in/M}},\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"and use this to check whether $DFT_M(w_0,w_1,...,w_{N-1},0,...,0)\\approx(M,0,...,0)$\n",
"\n",
"\n",
"\n",
"**Solution.**\n",
"This has absolute value $\\left|\\frac{\\sin(\\pi nN/M)}{\\sin(\\pi n/M)}\\right|$. The $y_n$ thus lie on the curve \n",
"$\\frac{\\sin(N\\omega/2)}{\\sin(\\omega/2)}$, from which one can see that $DFT_M(w_0,w_1,...,w_{N-1},0,...,0)\\approx(M,0,...,0)$ \n",
"is partially fulfilled. The frequency response can be plotted as follows"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"N=32;\n",
"omega = linspace(0, 2*pi, 1000);\n",
"wd = [ones(1, N) zeros(1, 1000-N)];\n",
"plot(omega, abs(fft(wd)), 'k-');"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"\n",
"**c)**\n",
"The *Hamming window* is defined by"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"w_n=2(0.54-0.46\\cos(2\\pi n/(N-1))).\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Make a plot of $DFT_M(w_0,w_1,...,w_{N-1},0,...,0)$, and compare with $(M,0,...,0)$ to see if the Hamming window is a good window.\n",
"\n",
"\n",
"\n",
"**Solution.**\n",
"The frequency response of the Hamming window can be plotted as follows"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"n = 0:(N-1);\n",
"wd = [ 2*(0.54-0.46*cos(2*pi*n/(N-1))) zeros(1, 1000-N)];\n",
"plot(omega, abs(fft(wd)), 'k-');"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The Hamming window is seen to differ from the rectangular window in two ways: \n",
"* It has much smaller values away from $0$,\n",
"\n",
"* the width of the \"main lobe\" (i.e. the main structure at the center) is bigger. \n",
"\n",
"As a consequence, windowing with the Hamming window reduces the contribution from higher frequencies, but increases the contribution from the smallest frequencies. \n",
"The window coefficients for both windows \n",
"can be plotted with the following code"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"N = 32;\n",
"hold on\n",
"for k = 0:(N-1)\n",
" plot([k k],[0 1],'k-')\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"hold on\n",
"for k = 0:(N-1)\n",
" plot([k k],[0 2*(0.54-0.46*cos(2*pi*k/(N-1)))], 'k-')\n",
"end"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"\n",
"**d)**\n",
"The *Hanning window* is defined by"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"w_n=1-\\cos(2\\pi n/(N-1)).\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Repeat the analysis you did above for the rectangular and Hamming windows, for a Hanning window for $N=32$. \n",
"\n",
"The Hanning window is used in the MP3 standard, where it is applied to blocks which overlap, contrary to the non-overlapping block pattern we have used.\n",
"After the Hanning window has been applied, the MP3 standard applies an FFT to the windowed signal in order to make a frequency analysis of that part of the sound.\n",
"\n",
"\n",
"\n",
"**Solution.**\n",
"The frequency response can be plotted as follows"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"wd = [1-cos(2*pi*n/(N-1)) zeros(1,1000-N)];\n",
"plot(omega, abs(fft(wd)), 'k-');"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"A good survey on the windows mentioned here, as well as other ones, can be found in [[harris01]](#harris01). \n",
"One lesson to be learned from this exercise is that the windowed signal may not give a good frequency representation. \n",
"One must therefore be careful when splitting a sound into blocks, as this alters the frequency representation. \n",
"You can use the function `forw_comp_rev_DFT` from Example 2.5 to experiment with this.\n",
"This function accepts a named parameter $N$, which can be used to split the DFT of a sound into blocks of length $N$, and eliminate low frequencies before taking an IDFT to reconstruct and listen to the new sound \n",
"(see also Example 3.31).\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"## Exercise 3.10: The multiplicative group of integers modulo $N$\n",
"\n",
"\n",
"\n",
"In number theory it is known that, when $N$ is a prime number, the numbers $1,...,N-1$ is a *group* under multiplication modulo $N$. \n",
"In particular this means that for each $1\\leq a\\leq N-1$ there exists a unique $1\\leq b\\leq N-1$ so that $ab=1\\bmod N$ \n",
"(the number $1$ acts as a *unit*). \n",
"$b$ us then called the inverse of $a$, and we write $b=a^{-1}$ (this is not the same as $1/a$!). \n",
"Also, it is known that there exists an integer $g$ so that the set $\\{g^q\\bmod N\\}_{q=0}^{N-2}$ constitutes all numbers $1,...,N-1$ \n",
"(the order these numbers appear in may have changed, though).\n",
"$g$ is called a *generator* for the group.\n",
"\n",
"\n",
"**a)**\n",
"Find generators for the multiplicative group of integers modulo $11$, $23$ and $41$. Are the generators for each group unique?\n",
"\n",
"\n",
"\n",
"**Solution.**\n",
"Even for small primes, one saves tedious computation here with an implementation. The following code returns the smallest generator."
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"%%file find_smallest_generator.m\n",
"function g=find_smallest_generator(N)\n",
" g = 0;\n",
" for k=2:(N-1)\n",
" res = zeros(N-1,1);\n",
" val = 1;\n",
" for s=1:(N-1)\n",
" if res(val) ~= 0\n",
" break;\n",
" end\n",
" res(val) = s;\n",
" val = mod(val*k,N);\n",
" if s == N-1\n",
" g = k; \n",
" return;\n",
" end\n",
" end\n",
" end\n",
"end"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The outer loop starts at $2$, because the unit $1$ is never a generator. If your run this function for the stated primes,"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"find_smallest_generator(11)\n",
"find_smallest_generator(23)\n",
"find_smallest_generator(41)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"the generators $2$, $5$, and $6$, respectively, are found. \n",
"\n",
"A generator may not be unique: In the next exercise we state that $g^{-1}$ is a generator whenever $g$ is.\n",
"\n",
"\n",
"\n",
"**b)**\n",
"Write a function"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
" reorder(x, g)\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"which, for a vector $\\boldsymbol{x}$ of length $N$, and a generator $g$ for the numbers $1,...,N-1$, \n",
"returns the vector $(x_{g^0},x_{g^1},...,x_{g^{N-2}})$.\n",
"\n",
"\n",
"In the next exercise you will see how the multiplicative group of integers modulo $N$ relates to circular convolution.\n",
"\n",
"\n",
"\n",
"**Solution.**\n",
"The following code can be used."
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"%%file reorder.m\n",
"function xnew=reorder(x, g)\n",
" N = length(x);\n",
" xnew = zeros(N-1,1);\n",
" val = 1;\n",
" for s=1:(N-1)\n",
" xnew(s) = x(val+1);\n",
" val = mod(val*g, N);\n",
" end\n",
"end"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"## Exercise 3.11: The FFT algorithm when $N$ is prime. Rader's algorithm\n",
"\n",
"\n",
"\n",
"It is possible to show that, if $g$ is a generator, then $g^{-1}$ is also a generator.\n",
"This means that $\\{g^{-p}\\bmod N\\}_{p=0}^{N-2}$ also constitute \n",
"the numbers $1,...,N-1$ (where $g^{-p}$ is defined as $(g^{-1})^p$). \n",
"\n",
"Replacing $n=g^{-p}$ and $k=g^q$ in the FFT algorithm with $N$ prime, \n",
"we can rephrase it as (we pull $k=0$ out of the sum, since the number $0$ is not an element in our group)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\begin{align*}\n",
"y_0 &= \\sum_{k=0}^{N-1} x_k \\\\ \n",
"y_{g^{-p}} &= x_0 + \\sum_{q=0}^{N-2} x_{g^q}e^{-2\\pi i g^{-(p-q)}/N} \\text{, } 0\\leq p\\leq N-2.\n",
"\\end{align*}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Define"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\begin{align*}\n",
"\\boldsymbol{a} &= (x_{g^0} ,x_{g^1} ,...,x_{g^{N-2}}) \\in\\mathbb{R}^{N-1}\\\\ \n",
"\\boldsymbol{b} &= (e^{-2\\pi i g^0/N}, e^{-2\\pi i g^{-1}/N},...,e^{-2\\pi i g^{-(N-2)}/N}) \\in\\mathbb{R}^{N-1}\\\\ \n",
"\\boldsymbol{c} &= (y_{g^0} ,y_{g^{-1}} ,...,y_{g^{-(N-2)}}) \\in\\mathbb{R}^{N-1}\n",
"\\end{align*}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Explain that $\\boldsymbol{c}=x_0+\\boldsymbol{a}\\circledast\\boldsymbol{b}$, where $x_0$ is added to every component in the vector. \n",
"\n",
"This explains how to compute an $N$-point DFT (with $N$ prime) using ($N-1$)-point circular convolution. \n",
"Since a circular convolution can be computed using a DFT/IDFT of the same size, so this method effectively reduces an $N$-point DFT to an $(N-1)$-point DFT. \n",
"Since $N-1$ is not prime, we can use the algorithm for the FFT for composite $N$, to reduce the problem to DFT's of lengths being smaller prime numbers.\n",
"Another possibility is to use Exercise 3.5 to replace the circular convolution with a longer one, of length a power of two, \n",
"effectively reducing the problem to a DFT of length a power of two.\n",
"\n",
"\n",
"\n",
"**Solution.**\n",
"Inserting the definitions of $\\boldsymbol{a}$, $\\boldsymbol{b}$, and $\\boldsymbol{c}$ in the DFT formula above gives"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"c_p = x_0 + \\sum_{q=0}^{N-2} a_q b_{p-q} = x_0 + \\sum_{q=0}^{N-2} a_q b_{(p-q)\\bmod(N-1)} = x_0 + (\\boldsymbol{a}\\circledast\\boldsymbol{b})_p.\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The formula follows. We used that $b_{p-q}=b_{(p-q)\\bmod(N-1)}$, which follows from that $g^{N-1}$ equals the unit $e$: If this was not the case, we would have that $g^{N-1}=g^i$ for some $i>0$. This gives that $g^{N-1-i}=e$, which would imply that $g$ does not generate all numbers from $1$ to $N-1$.\n",
"The following code tests that Rader's algorithm produces the same result as an $N$-point FFT when $N$ is prime."
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"N = 41;\n",
"x=rand(N,1);\n",
"g=find_smallest_generator(N);\n",
"a = reorder(x, g);\n",
"b = reorder(exp(-2*pi*i*(0:(N-1))/N), g); b(2:(N-1)) = flip(b(2:(N-1)));\n",
"c = cconv(a, b, N-1) + x(1);\n",
"c(2:(N-1)) = flip(c(2:(N-1)));\n",
"y = x;\n",
"y(1) = sum(x);\n",
"ind = 1;\n",
"for k = 1:(N-1)\n",
" y(ind+1) = c(k);\n",
" ind = mod(ind*g, N);\n",
"end\n",
"max(abs(fft(x)-y))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Here we \n",
"used the built-in function `cconv` for circular convolution in MATLAB. \n",
"Later we will consider our own implementation of circular convolution, but this will have a different assumption on the input. \n",
"In the code above, also note that we reversed the vectors $\\boldsymbol{b}$ and $\\boldsymbol{c}$. This is because these used $g^{-1}$ rather than $g$ to generate their values.\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"## Example 3.14: Plotting a frequency response\n",
"\n",
"\n",
"\n",
"The frequency response of the filter (3.13) is"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\begin{align*}\n",
"\\lambda_S(\\omega) &= \\frac{1}{4} e^{i\\omega} + \\frac{1}{2}e^{0} + \\frac{1}{4} e^{-i\\omega} = \\frac{1}{2} + \\frac{1}{4} (e^{i\\omega} + e^{-i\\omega}) = \\frac{1}{2} + \\frac{1}{2}\\cos(\\omega).\n",
"\\end{align*}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We can plot the frequency response over $[0,2\\pi]$ as follows,"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"omega=linspace(0, 2*pi, 100);\n",
"figure();\n",
"plot(omega,0.5+0.5*cos(omega), 'k-')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"and over $[-\\pi,\\pi]$ as follows."
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"omega=linspace(-pi, pi, 100);\n",
"plot(omega,0.5+0.5*cos(omega), 'k-')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The plots clearly show that the high frequencies are made smaller by the filter, and it is therefore called a *low-pass filter*. These are investigated further in the next section. \n",
"\n",
"A more general way to plot the frequency response of a FIR filter is to\n",
"combine Proposition 7 with Theorem 2: If $\\boldsymbol{t}$ is the set of filter coefficients, equation (3.3) gives that"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
" omega = 2*pi*(0:(N-1))/N;\n",
" s = [t; zeros(N - length(t), 1)];\n",
" plot(omega, abs(fft(s))); \n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Some comments are in order. First of all we have restricted to $[0,2\\pi)$. Secondly, $\\boldsymbol{t}$ was expanded to a longer vector $\\boldsymbol{s}$, to increase the number of plot points.\n",
"Finally, the expanded vector $\\boldsymbol{s}$ is correct only up to a delay, due to Proposition 7, \n",
"and as we have seen a delay of the filter does not affect the magnitude of the frequency response.\n",
"Here only the magnitude is plotted, as we often do.\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"## Example 3.16: Adding echo\n",
"\n",
""
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"[x, f_s] = audioread('sounds/castanets.wav');\n",
"d = 10000;\n",
"c = 0.5;\n",
"k = 5;"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"An echo is a delayed and softer copy of the original sound. If $\\boldsymbol{x}$ is this original sound, the sound $\\boldsymbol{z}$ with samples given by"
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"[N, nchannels] = size(x);\n",
"z = zeros(N, nchannels);\n",
"z(1:d,:) = x(1:d,:); % No echo at the start\n",
"z((d+1):N,:) = x((d+1):N,:)+c*x(1:(N-d),:); % Add echo\n",
"z = z/max(max(abs(z))); % Scale to within [-1,1]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"includes an echo of $\\boldsymbol{x}$. \n",
"$d$ is an integer which represents the delay in samples. If you need a delay in $t$ seconds, set $d=tf_s$ (and round to the nearest integer), with $f_s$ the sample rate. \n",
"$c$ is called the *damping factor*. An echo is usually weaker than the original sound, so that $c < 1$. \n",
"Let us listen to the sound after adding echo."
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"playerobj=audioplayer(z, f_s);\n",
"playblocking(playerobj);"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"You are encouraged to experiment and find the range of $d$ which makes the echo indistinguishable from the sound itself, \n",
"and how low can you choose $c$ in order to still hear the echo.\n",
"\n",
"The compact notation for a filter which adds echo is"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"S=\\{ \\underline{1},0,\\ldots,0,c \\},\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"where the damping factor $c$ appears at index $d$. \n",
"The frequency response of this is $\\lambda_S(\\omega)=1+ce^{-id\\omega}$, which is not real. \n",
"We can plot the frequency response as follows"
]
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"omega = linspace(0.01, 2*pi-0.01, 100);\n",
"plot(omega, abs(1+0.1*exp(-i*10*omega)), 'k-')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We see that the response varies between $1-c$ and $1+c$. The oscillation is controlled by the delay $d$.\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"## Exercise 3.19: Composing two filters\n",
"\n",
"\n",
"\n",
"Assume that the filters $S_1$ and $S_2$ have the frequency responses $\\lambda_{S_1}(\\omega)=2+4\\cos(\\omega)$, $\\lambda_{S_2}(\\omega)=3\\sin(2\\omega)$.\n",
"\n",
"\n",
"**a)**\n",
"Compute and plot the frequency response of the filter $S_1S_2$.\n",
"\n",
"\n",
"\n",
"**Solution.**\n",
"We have that"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\lambda_{S_1S_2}(\\omega) = \\lambda_{S_1}(\\omega) \\lambda_{S_2}(\\omega) = (2+4\\cos(\\omega))(3\\sin(2\\omega)).\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The frequency response can be plotted as follows."
]
},
{
"cell_type": "code",
"execution_count": 20,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"omega=linspace(-pi,pi,100);\n",
"plot(omega, (2+4*cos(omega)).*(3*sin(2*omega)));"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"**b)**\n",
"Write down the filter coefficients $t_k$ for the filter $S_1S_2$.\n",
"\n",
"\n",
"\n",
"**Solution.**\n",
"We have that"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\begin{align*}\n",
"\\lambda_{S_1S_2}(\\omega) &= (2+4\\cos(\\omega))(3\\sin(2\\omega)) = (2+2e^{i\\omega}+2e^{-i\\omega})\\frac{3}{2i}(e^{2i\\omega}-e^{-2i\\omega}) \\\\ \n",
" &= -3i (1+e^{i\\omega}+e^{-i\\omega})(e^{2i\\omega}-e^{-2i\\omega}) \\\\ \n",
" &= -3i(e^{2i\\omega}-e^{-2i\\omega} + e^{3i\\omega} + e^{i\\omega} - e^{-i\\omega} - e^{-3i\\omega}).\n",
"\\end{align*}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"From this we see that $S_1S_2=3i\\{-1,-1,-1,\\underline{0},1,1,1\\}$.\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"## Exercise 3.25: Execution times for convolution\n",
"\n",
"\n",
"\n",
"Implement code which computes $\\boldsymbol{t}\\ast\\boldsymbol{x}$ in the vectorized- and non-vectorized ways described in Section 3.2.1 (i.e. as a single loop over $k$ or $n$ with the other variable vectorized, or a double loop). \n",
"As your $\\boldsymbol{t}$, take $k$ randomly generated numbers. Compare execution times with the `conv` function, for different values of $k$. Present the result as a\n",
"plot where $k$ runs along the $x$-axis, and execution times run along the $y$-axis. Your result will depend on how vectorization is performed by the computing environment.\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"**Solution.**\n",
"The code can look as follows."
]
},
{
"cell_type": "code",
"execution_count": 21,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"[x, f_s] = audioread('sounds/castanets.wav');\n",
"x = x(:,1);\n",
"N= length(x);\n",
"\n",
"kmax=100;\n",
"vals1 = zeros(1, kmax/10);\n",
"vals2 = zeros(1, kmax/10);\n",
"vals3 = zeros(1, kmax/10);\n",
"ind = 1;\n",
"for k=10:10:kmax\n",
" t = rand(k, 1);\n",
" tic;\n",
" conv(t, x);\n",
" vals1(ind) = toc;\n",
" \n",
" z = zeros(N, 1);\n",
" tic;\n",
" for s = 1:k\n",
" z(k:N) = z(k:N) + t(s)*x((k-s+1):(N-s+1));\n",
" end\n",
" vals2(ind)=toc;\n",
"\n",
" z = zeros(N, 1);\n",
" tic;\n",
" for n=k:N\n",
" for s = 1:k\n",
" z(n) = z(n) + t(s)*x(n-s+1);\n",
" end\n",
" end\n",
" vals3(ind)=toc;\n",
" ind = ind+1;\n",
"end\n",
"plot(10:10:kmax,log(vals1),10:10:kmax,log(vals2), 10:10:kmax, log(vals3));\n",
"legend('conv','simple for','double for')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"## Example 3.30: Moving average filters\n",
"\n",
"\n",
"\n",
"A general way of reducing variations in sound is to replace one number by the average of itself and its neighbors. \n",
"If $\\boldsymbol{z}=(z_i)_{i=0}^{N-1}$ is the sound signal produced by taking the average of three successive samples, we have that"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"z_n = \\frac{1}{3}(x_{n+1} + x_n + x_{n-1}),\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"i.e. $S=\\{1/3,\\underline{1/3},1/3\\}$. This is also called a *moving average filter* (with three elements). \n",
"If we set $N=4$, the corresponding circulant Toeplitz matrix for the filter is"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"S\n",
"=\n",
"\\frac{1}{3}\n",
"\\begin{pmatrix}\n",
"1 & 1 & 0 & 1 \\\\ \n",
"1 & 1 & 1 & 0 \\\\ \n",
"0 & 1 & 1 & 1 \\\\ \n",
"1 & 0 & 1 & 1 \n",
"\\end{pmatrix}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The frequency response is"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\lambda_S(\\omega)=(e^{i\\omega}+1+e^{-i\\omega})/3=(1+2\\cos(\\omega))/3.\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"More generally we can construct the moving average filter of $2L+1$ elements, which is $S=\\{1,\\cdots,\\underline{1},\\cdots,1\\}/(2L+1)$ (there is symmetry around $t_0$). \n",
"Let us verify that these filters are low-pass filters.\n",
"\n",
"Clearly the first column of $S$ is $\\boldsymbol{s}=(\\underbrace{1,\\ldots,1}_{L+1\\text{ times}},0,\\ldots,0,\\underbrace{1,\\ldots,1}_{L\\text{ times}})/(2L+1)$.\n",
"In Example 2.2 we computed that"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"DFT_N\\left(\\underbrace{1,\\ldots,1}_{L+1\\text{ times}},0,\\ldots,0,\\underbrace{1,\\ldots,1}_{L\\text{ times}}\\right)=\\boldsymbol{y}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"where $\\boldsymbol{y}$ had components"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"y_n = \\frac{\\sin(\\pi n(2L+1)/N)}{\\sin(\\pi n/N)}.\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Since $\\boldsymbol{\\lambda}_S=\\text{DFT}_N\\boldsymbol{s}$, dividing by $2L+1$ and inserting $\\omega=2\\pi n/N$ gives that"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\lambda_S(\\omega)=\\frac{1}{2L+1}\\frac{\\sin((2L+1)\\omega/2)}{\\sin(\\omega/2)}.\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We clearly have"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"0\\leq\\frac{1}{2L+1}\\frac{\\sin((2L+1)\\omega/2)}{\\sin(\\omega/2)}\\leq 1,\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"and the frequency response approaches $1$ as $\\omega\\to 0$, so that it peaks at $0$. This peak gets narrower and narrower as $L$ increases. \n",
"This filter thus \"keeps\" only the lowest frequencies. It is also seen that the frequency response is small for $\\omega\\approx\\pi$. \n",
"In fact it is straightforward to see that $|\\lambda_S(\\pi)|=1/(2L+1)$. \n",
"The frequency responses for the moving average filters corresponding to $L=1$, $L=5$, and $L=20$\n",
"can be plotted as follows."
]
},
{
"cell_type": "code",
"execution_count": 22,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"omega = linspace(0.01, 2*pi-0.01, 100);\n",
"for L = [1 5 20]\n",
" figure()\n",
" plot(omega, sin((2*L+1)*omega/2)./(sin(omega/2)*(2*L+1)), 'k-')\n",
"end"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In conclusion, moving average filters are low-pass filters, but they are far from ideal such, since not all higher frequencies are annihilated, and since small frequencies also are changed.\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"## Example 3.31: Ideal low-pass filters\n",
"\n",
"\n",
"\n",
"By definition, the ideal low-pass filter keeps frequencies near $0$ unchanged, and completely removes frequencies near $\\pi$. \n",
"$S=(F_N)^HDF_N$ is an ideal low-pass filter when $D$ is diagonal with"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"(\\underbrace{1,\\ldots,1}_{L+1\\text{ times}},0,\\ldots,0,\\underbrace{1,\\ldots,1}_{L\\text{ times}})\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"on the diagonal. \n",
"If the filter should keep the angular frequencies $|\\omega|\\leq\\omega_c$ only, where $\\omega_c$ is the \"cutoff\" frequency, we should choose $L$ so that $\\omega_c=2\\pi L/N$. \n",
"An ideal high-pass filters similarly corresponds to a diagonal matrix with"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"(\\underbrace{0,\\ldots,0}_{N/2-L\\text{ times}},\\underbrace{1,\\ldots,1}_{2L+1\\text{ times}},\\underbrace{0,\\ldots,0}_{N/2-L-1\\text{ times}})\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"on the diagonal.\n",
"\n",
"Let us compute the filter coefficients for the ideal low-pass filter.\n",
"Again, in Example 2.2 we computed the DFT of the vector above, and it followed from Theorem 7 that the IDFT of this vector equals its DFT, up to a factor $1/N$. \n",
"This means that $\\boldsymbol{s}=IDFT_N\\boldsymbol{\\lambda}_S$ is"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\frac{1}{N}\\frac{\\sin(\\pi k(2L+1)/N)}{\\sin(\\pi k/N)}.\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The filter coefficients are thus $N$ points uniformly spaced between $0$ and $1$ on the curve $\\frac{1}{N}\\frac{\\sin(\\pi t(2L+1)/2)}{\\sin(\\pi t/2)}$. \n",
"This curve has been encountered many other places in the book. \n",
"Moving average filters and ideal low-pass filters thus have a duality between vectors which contain only zeros and ones on one side (i.e. windows), \n",
"and the vector $\\frac{1}{N}\\frac{\\sin(\\pi k(2L+1)/N)}{\\sin(\\pi k/N)}$ on the other side: \n",
"filters of the one type correspond to frequency responses of the other type, and vice versa. \n",
"\n",
"The extreme cases for $L$ are\n",
"* $L=1$: Only the lowest frequency is kept. All filter coefficients are equal to $\\frac{1}{N}$. \n",
"\n",
"* $L=N$: All frequencies are kept. The filter equals the identity matrix.\n",
"\n",
"Between these two extremes, $s_0$ is the biggest coefficient, while the others decrease towards $0$ along the curve we stated. \n",
"The bigger $L$ and $N$ are, the quicker they decrease to zero. All filter coefficients are usually nonzero for this filter, since this curve is zero only at certain points. \n",
"The filter is thus not a FIR filter. Many filters which are not FIR still have efficient implementations, but for this filter it turns out to be difficult to find one.\n",
"The best thing we can do is probably to use a DFT for computing the filter, followed by an inverse DFT.\n",
"\n",
"The function `forw_comp_rev_DFT` from Example 2.5 accepts named parameters `L` and `lower`, where `L` is as described above, and where `lower` states whether the lowest or the highest frequencies should be kept. \n",
"We can use the function to listen to the lower frequencies in the audio sample file.\n",
"We start with $L=13000$."
]
},
{
"cell_type": "code",
"execution_count": 23,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"[x, fs] = forw_comp_rev_DFT('L', 13000, 'lower', 1);\n",
"playerobj=audioplayer(x, fs);\n",
"playblocking(playerobj);"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Then with $L=5000$."
]
},
{
"cell_type": "code",
"execution_count": 24,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"[x, fs] = forw_comp_rev_DFT('L', 5000, 'lower', 1);\n",
"playerobj=audioplayer(x, fs);\n",
"playblocking(playerobj);"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"With $L=13000$ you can hear the disturbance in the sound, but we have not lost that much even if about 90\\% of the DFT coefficients are dropped. \n",
"The quality is much poorer when $L=5000$ (here we keep less than 5\\% of the DFT coefficients). \n",
"However we can still recognize the song, and this suggests that most of the frequency information is contained in the lower frequencies.\n",
"\n",
"Let us then listen to higher frequencies instead. \n",
"We start with $L=140000$."
]
},
{
"cell_type": "code",
"execution_count": 25,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"[x, fs] = forw_comp_rev_DFT('L', 140000, 'lower', 0);\n",
"playerobj=audioplayer(x, fs);\n",
"playblocking(playerobj);"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Then with $L=100000$."
]
},
{
"cell_type": "code",
"execution_count": 26,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"[x, fs] = forw_comp_rev_DFT('L', 100000, 'lower', 0);\n",
"playerobj=audioplayer(x, fs);\n",
"playblocking(playerobj);"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We find that we need very high values of $L$ to hear anything, suggesting again that most information is contained in the lowest frequencies.\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"## Example 3.32: Windowing an ideal low-pass filter\n",
"\n",
"\n",
"\n",
"In order to decrease the operation count of the ideal low-pass filter, one could apply a rectangular window to the filter coefficients (see Exercise 3.9), i.e. consider the filter"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\left\\{\\frac{1}{N}\\frac{\\sin(\\pi k(2L+1)/N)}{\\sin(\\pi k/N)}\\right\\}_{k=-N_0}^{N_0}.\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In light of that exercise, this may not be the best strategy - applying a window different from the rectangular one may better preserve the frequency response of the ideal low-pass filter. \n",
"\n",
"Consider the ideal low-pass filter with $N=128$, $L=32$ (i.e. the filter removes frequencies $\\omega>\\pi/2$). \n",
"The following code tests this for different values of $N_0$."
]
},
{
"cell_type": "code",
"execution_count": 27,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"N = 128;\n",
"L = 32;\n",
"k=1:(N-1);\n",
"filtercoeffs = sin(pi*k*(2*L+1)/N)./(sin(pi*k/N)*N);\n",
"omega=linspace(0, 2*pi, 100);\n",
" \n",
"for N0 = [2 4 16 64]\n",
" freqresp=((2*L+1)/N)*ones(1,length(omega));\n",
" for k = 1:N0\n",
" freqresp = freqresp + filtercoeffs(k)*2*cos(omega*k);\n",
" end\n",
" figure();\n",
" plot(omega, freqresp, 'k-')\n",
"end"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This shows that we should be careful when we omit filter coefficients: if we drop too many, the frequency response is far away from that of an ideal\n",
"low-pass filter. In particular, we see that the new frequency response oscillates wildly near the discontinuity of the ideal low-pass filter. \n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"## Example 3.34: Low-pass filters deduced from Pascal's triangle\n",
"\n",
"\n",
"\n",
"When computing an average, it is reasonable to let the middle sample count more than the neighbors. So, an alternative to that of moving averages is to compute"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"z_n = (x_{n-1} + 2x_n + x_{n+1})/4\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The coefficients $1,2,1$ are taken from row 2 in Pascal's triangle. It turns out that this is a good choice of coefficients. \n",
"In fact, other rows from Pascal's triangle are also good choices. \n",
"To explain why, let $S=\\{\\underline{1/2},1/2\\}$ be the moving average filter of two elements.\n",
"The frequency response of $S$ is"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\lambda_S(\\omega)=\\frac{1}{2}(1+e^{-i\\omega})=e^{-i\\omega/2}\\cos(\\omega/2).\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"If we apply this filter $k$ times, Theorem 6 states that"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\lambda_{S^k}(\\omega)=\\frac{1}{2^k}(1+e^{-i\\omega})^k=e^{-ik\\omega/2}\\cos^k(\\omega/2),\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"which is a polynomial in $e^{-i\\omega}$ with the coefficients taken from Pascal's triangle (the values in Pascal's triangle are the\n",
"coefficients of $x$ in the expression $(1+x)^k$, i.e. the binomial coefficients $\\binom{k}{r}$ for $0\\leq r\\leq k$). \n",
"Thus, the filter coefficients of $S^k$ are rows in Pascal's triangle. \n",
"The reason why the filters $S^k$ are more desirable than moving average filters is that,\n",
"since $(1+e^{-i\\omega})^k$ is a factor in its frequency response, it has a zero of multiplicity $k$ at $\\omega=\\pi$. \n",
"This implies that the frequency response is very flat for $\\omega\\approx\\pi$ when $k$ increases, i.e. the filter is good at removing the highest frequencies. This can be seen\n",
"with help from the following code, where we plot\n",
"the magnitude of the frequency response when $k=5$, and when $k=30$."
]
},
{
"cell_type": "code",
"execution_count": 28,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"omega = linspace(0, 2*pi, 100);\n",
"for k=[5 30]\n",
" figure()\n",
" plot(omega, abs((cos(omega/2)).^k), 'k-')\n",
" axis([0 2*pi 0 1.2])\n",
"end"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Clearly the latter frequency response is much flatter for $\\omega\\approx\\pi$. On the other side, the filters of Example 3.30 satisfied $|\\lambda_S(\\pi)|=1/(2L+1)$, with a frequency response not very flat near $\\pi$ \n",
"\n",
"While using $S^k$ gives a desirable behavior for $\\omega\\approx\\pi$, we see that the behavior is not so desirable for small frequencies $\\omega\\approx 0$: Only frequencies very close to $0$ are kept\n",
"unaltered. It should be possible to produce better low-pass filters than this also, as the frequency responses of the filters in the MP3 standard suggest.\n",
"\n",
"Filtering with $S^k$ can be computed by iterating the filter $\\{\\underline{1/2},1/2\\}$ $k$ times:"
]
},
{
"cell_type": "code",
"execution_count": 29,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"k = 5;\n",
"z = x(:, 1);\n",
"for kval=1:k\n",
" z = conv(z,[1/2 1/2]);\n",
"end"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This code disregards the circularity of $S$, and we introduce a time delay. These issues\n",
"will, however, not be audible when we listen to the output. \n",
"In Exercise 3.36 you will be asked to perform these steps our sample audio file.\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"## Example 3.35: High-pass filters deduced from Pascal's triangle\n",
"\n",
"\n",
"\n",
"If we apply Proposition 11 to the low-pass filter deduced from the fourth row in Pascal's triangle we obtain"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"z_n = \\frac{1}{16}(x_{n-2}-4x_{n-1}+6x_n-4x_{n+1}+x_{n+2})\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Clearly the high-pass filter arising from row $k$ in Pascal's triangle can be written as $S^k$, with $S=\\frac{1}{2}\\{\\underline{1},-1\\}$. \n",
"In other words, we can use convolution as in the previous example to compute the output from such filters. \n",
"In Exercise 3.36 you will be asked to apply these to the audio sample file. \n",
"The new sound will be difficult to hear for large $k$, and we will explain why later. \n",
"\n",
"The frequency response we obtain from using row $5$ of Pascal's triangle\n",
"can be plotted as follows"
]
},
{
"cell_type": "code",
"execution_count": 30,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"t=0:(1/4400):0.01;\n",
"vals=sin(2*pi*440*t);\n",
"arr2 = [ -1/4400 -2/4400 t 0.01+1/4400 0.01+2/4400];\n",
"plot(arr2,conv(vals,[1 -4 6 -4 1]/16),'-ko')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"## Exercise 3.36: Applying low-pass- and high-pass filters deduced from Pascal's triangle to the audio sample file\n",
"\n",
"\n",
"\n",
"\n",
"**a)**\n",
"Write code where you apply the low-pass and high-pass filters in examples 3.34 and 3.35 to the audio sample file, and verify that the sounds you get are the same as in these examples.\n",
"How high must $k$ be in order for you to hear difference from the actual sound? How high can you choose $k$ and still recognize the sound at all? \n",
"If you solved Exercise 3.28 you can also use the function `filter_impl` to perform the filtering, \n",
"rather than using convolution (which, as mentioned, discards circularity).\n",
"\n",
"\n",
"\n",
"**Solution.**\n",
"The code can look like this, when we apply low-pass filters."
]
},
{
"cell_type": "code",
"execution_count": 31,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"k = 5;\n",
"z = x(:, 1);\n",
"for kval=1:k\n",
" z = conv(z,[1/2 1/2]);\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": 32,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"playerobj=audioplayer(z, f_s);\n",
"playblocking(playerobj);"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"If we apply high-pass filters instead, the code can look like this."
]
},
{
"cell_type": "code",
"execution_count": 33,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"z = x(:, 1);\n",
"for kval=1:k\n",
" z = conv(z,[1/2 -1/2]);\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": 34,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"playerobj=audioplayer(z, f_s);\n",
"playblocking(playerobj);"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"**b)**\n",
"In your code, it is not necessary to scale the values after applying the low-pass or high-pass filters so that values fit inside $[-1,1]$. Explain why this is the case.\n",
"\n",
"\n",
"\n",
"**Solution.**\n",
"The low-pass filters compute weighted averages of the input samples. Therefore, as long as the input values are inside the legal range$ [-1,1]$, the output values will as well. \n",
"Since the filter coefficients sum to one it is easy to see that also the high-pass filters produce values in the legal range.\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"## Exercise 3.45: The IIR filter for the Karplus-Strong algorithm\n",
"\n",
"\n",
"\n",
"Let us rewrite the difference equation in Exercise 1.12 as"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"z_{n}-\\frac{1}{2}(z_{n-p}+z_{n-(p+1)}) = 0\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**a)**\n",
"Explain that the frequency response of this realization is 0. Therefore, this realization is only interesting if the initial conditions are chosen so that we don't obtain a filter.\n",
"\n",
"\n",
"\n",
"**Solution.**\n",
"The numerator in $\\lambda_S(\\omega)$ is zero, since all $b_k=0$. We therefore have that $\\lambda_S(\\omega)=0$ also.\n",
"\n",
"\n",
"\n",
"**b)**\n",
"Plot the zeros of the characteristic equation for $p=20$. Is the filter BIBO stable?\n",
"\n",
"\n",
"\n",
"**Solution.**\n",
"We have that $a_0=1$, $a_p=a_{p+1}=-1/2$, so that the characteristing equation is $\\sum_{k=0}^M a_k r^{-k} = 1 - \\frac{1}{2}r^{-p} - \\frac{1}{2}r^{-(p+1)}$. Its zeros can be plotted as follows."
]
},
{
"cell_type": "code",
"execution_count": 35,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"p=20;\n",
"rts = roots([1 zeros(1, p-1) -1/2 -1/2]);\n",
"plot(real(rts), imag(rts), 'kx')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We see from the resulting plot that the zeros are inside the unit circle, so that the filter is BIBO stable.\n",
"\n",
"\n",
"\n",
"**c)**\n",
"What happens with the zeros of the characteristic equation when you increase $p$? Attempt to explain why the resulting sound changes more slowly to a static sound when you increase $p$.\n",
"\n",
"\n",
"\n",
"**Solution.**\n",
"It seems that the absolute values of the zeros of the characteristic equation gets closer and closer to $1$ as $p$ increases. This also means that their powers decreases slower to zero, so that the resulting sound goes to zero somewhat more slowly.\n",
"\n",
"\n",
"\n",
"**d)**\n",
"A related IIR filter is"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"z_{n}-\\frac{1}{2}(z_{n-p}+z_{n-(p+1)}) = x_n.\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Plot the frequency response of this filter for $p=20$. Include only values in the response between 0 and 10.\n",
"\n",
"\n",
"\n",
"**Solution.**\n",
"The following code can be used"
]
},
{
"cell_type": "code",
"execution_count": 36,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"p=20;\n",
"omega = linspace(0,2*pi,500);\n",
"plot(omega, abs(1./(1-0.5*exp(-omega*p*i)-0.5*exp(-omega*(p+1)*i))) );\n",
"axis([0,2*pi,0,10]);"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"## Exercise 3.48: Computing eigenvalues\n",
"\n",
"\n",
"\n",
"Consider the matrix"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"S=\\frac{1}{3}\\begin{pmatrix}\n",
"2 & 1 & 0 & 0 & 0 & 0 \\\\ \n",
"1 & 1 & 1 & 0 & 0 & 0 \\\\ \n",
"0 & 1 & 1 & 1 & 0 & 0 \\\\ \n",
"0 & 0 & 1 & 1 & 1 & 0 \\\\ \n",
"0 & 0 & 0 & 1 & 1 & 1 \\\\ \n",
"0 & 0 & 0 & 0 & 1 & 2\n",
"\\end{pmatrix}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**a)**\n",
"Compute the eigenvalues and eigenvectors of $S$ using the results of this section. You should only need to perform one DFT or one DCT in order to achieve this.\n",
"\n",
"\n",
"\n",
"**Solution.**\n",
"In c) we show that this indeed corresponds to a symmetric filter. We then know that the DCT basis vectors are eigenvectors. The eigenvalues can be found on the diagonal in $\\text{DCT}_NS(\\text{DCT}_N)^T$. \n",
"Note that it is not true that the eigenvalues can be computed by taking a DCT on the first column, as is the case for filters and the DFT.\n",
"\n",
"\n",
"\n",
"**b)**\n",
"Use a computer to find the eigenvectors and eigenvalues of $S$ also. What are the differences from what you found in a)?\n",
"\n",
"\n",
"\n",
"**Solution.**\n",
"The following code compares what you found in a) with a direct computation."
]
},
{
"cell_type": "code",
"execution_count": 37,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"S = [2 1 0 0 0 0; ...\n",
" 1 1 1 0 0 0; ...\n",
" 0 1 1 1 0 0; ...\n",
" 0 0 1 1 1 0; ...\n",
" 0 0 0 1 1 1; ...\n",
" 0 0 0 0 1 2];\n",
"\n",
"idctmatr = idct(eye(6));\n",
"idctmatr % Eigenvectors\n",
"dct(S*idctmatr) % eigenvalues\n",
"\n",
"[V, D]=eig(S);\n",
"V % Eigenvectors\n",
"D % eigenvalues"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The two typically will list eigenvectors and eigenvalues in different order.\n",
"\n",
"\n",
"\n",
"**c)**\n",
"Find a filter $T$ so that $S=T_r$. What kind of filter is $T$?\n",
"\n",
"\n",
"\n",
"**Solution.**\n",
"Writing down the circulant $12\\times12$-matrix for the symmetric filter $T=\\frac{1}{3}\\{1,\\underline{1},1\\}$ we obtain the two upper blocks"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\begin{align*}\n",
"T_1 &= \\frac{1}{3}\n",
"\\begin{pmatrix}\n",
"1 & 1 & 0 & 0 & 0 & 0 \\\\ \n",
"1 & 1 & 1 & 0 & 0 & 0 \\\\ \n",
"0 & 1 & 1 & 1 & 0 & 0 \\\\ \n",
"0 & 0 & 1 & 1 & 1 & 0 \\\\ \n",
"0 & 0 & 0 & 1 & 1 & 1 \\\\ \n",
"0 & 0 & 0 & 0 & 1 & 1\n",
"\\end{pmatrix}\n",
"&\n",
"T_2 &= \\frac{1}{3}\n",
"\\begin{pmatrix}\n",
"0 & 0 & 0 & 0 & 0 & 1 \\\\ \n",
"0 & 0 & 0 & 0 & 0 & 0 \\\\ \n",
"0 & 0 & 0 & 0 & 0 & 0 \\\\ \n",
"0 & 0 & 0 & 0 & 0 & 0 \\\\ \n",
"0 & 0 & 0 & 0 & 0 & 0 \\\\ \n",
"1 & 0 & 0 & 0 & 0 & 0\n",
"\\end{pmatrix}.\n",
"\\end{align*}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"From these expressions it is clear that $T_r=T_1+T_2^{\\leftrightarrow}=S$.\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
""
]
}
],
"metadata": {},
"nbformat": 4,
"nbformat_minor": 2
}