{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"# Exercises from Linear algebra, signal processing, and wavelets. A unified approach.\\\\ Python version\n",
"\n",
" **Øyvind Ryan**\n",
"\n",
"Date: **Oct 24, 2018**\n",
"\n",
"\n",
"\n",
"\n",
"# Motivation for wavelets and some simple examples\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"## Example 4.9: Computing the DWT by hand\n",
"\n",
""
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"import os, sys\n",
"sys.path.append(os.path.join(os.getcwd(), 'python'))"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"%matplotlib inline\n",
"\n",
"from numpy import *\n",
"import matplotlib.pyplot as plt\n",
"from scipy import integrate\n",
"from dwt import *\n",
"from sound import *\n",
"from forward_compress_reverse import *"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In some cases, the DWT can be computed by hand, keeping in mind its definition as a change of coordinates. \n",
"Consider the simple vector $\\boldsymbol{x}$ of length $2^{10}=1024$ defined by"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"x_n=\\begin{cases} 1 & \\text{for } 0\\leq n < 512 \\\\ 0 & \\text{for } 512\\leq n < 1024,\\end{cases}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"and let us compute the $10$-level DWT of $\\boldsymbol{x}$ by first visualizing the function with these coordinates. \n",
"Since $m=10$, $\\boldsymbol{x}$ are the coordinates in the basis $\\boldsymbol{\\phi}_{10}$ of a function $f\\in V_{10}$. \n",
"More precisely, $f(t)=\\sum_{n=0}^{511} \\phi_{10,n}(t)$, and since $\\phi_{10,n}$ is supported on\n",
"$[2^{-10}n,2^{-10}(n+1))$, the support of $f$ has width $512\\times 2^{-10}=1/2$ (512 translates, each of width $2^{-10}$). \n",
"Moreover, since $\\phi_{10,n}$ is $2^{10/2}=2^5=32$ on $[2^{-10}n,2^{-10}(n+1))$ and $0$ elsewhere, it is clear that"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"f(t) = \\begin{cases} 32 & \\text{for } 0\\leq t < 1/2 \\\\ 0 & \\text{for } t\\geq 1/2. \\end{cases}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This is by definition a function in $V_1$: $f$ must in fact be a multiple of $\\phi_{1,0}$, since this also is supported on $[0,1/2)$. \n",
"We can thus write $f(t)=c\\phi_{1,0}(t)$ for some $c$. We can find $c$ by setting $t=0$. This gives that\n",
"$32=2^{1/2}c$ (since $f(0)=32$, $\\phi_{1,0}(0)=2^{1/2}$), so that $c=32/\\sqrt{2}$. \n",
"This means that $f(t)=\\frac{32}{\\sqrt{2}}\\phi_{1,0}(t)$, $f$ is in $V_1$, and with coordinates $(32/\\sqrt{2},0,\\ldots,0)$ in $\\boldsymbol{\\phi}_1$.\n",
"\n",
"When we run a $10$-level DWT we make a change of coordinates from $\\boldsymbol{\\phi}_{10}$ to $(\\boldsymbol{\\phi}_0,\\boldsymbol{\\psi}_0,\\cdots,\\boldsymbol{\\psi}_{9})$. \n",
"The first $9$ levels give us the coordinates in $(\\boldsymbol{\\phi}_1,\\boldsymbol{\\psi}_1,\\boldsymbol{\\psi}_2,\\ldots,\\boldsymbol{\\psi}_9)$, and these are $(32/\\sqrt{2},0,\\ldots,0)$ from what we showed.\n",
"It remains thus only to perform the last level in the DWT, i.e. perform the change of coordinates from $\\boldsymbol{\\phi}_1$ to $(\\boldsymbol{\\phi}_0,\\boldsymbol{\\psi}_0)$. \n",
"Since $\\phi_{1,0}(t)=\\frac{1}{\\sqrt{2}}(\\phi_{0,0}(t)+\\psi_{0,0}(t))$ we get"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"f(t)=\\frac{32}{\\sqrt{2}}\\phi_{1,0}(t)=\\frac{32}{\\sqrt{2}}\\frac{1}{\\sqrt{2}}(\\phi_{0,0}(t)+\\psi_{0,0}(t))=16\\phi_{0,0}(t)+16\\psi_{0,0}(t).\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"From this we see that the coordinate vector of $f$ in $(\\boldsymbol{\\phi}_0,\\boldsymbol{\\psi}_0,\\cdots,\\boldsymbol{\\psi}_{9})$, i.e. the 10-level DWT of $\\boldsymbol{x}$, is \n",
"$(16,16,0,0,\\ldots,0)$. Note that here $V_0$ and $W_0$ are both $1$-dimensional, since $V_{10}$ was assumed to\n",
"be of dimension $2^{10}$ (in particular, $N=1$).\n",
"\n",
"It is straightforward to verify what we just found as follows:"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"x = concatenate([ones(512), zeros(512)])\n",
"dwt_impl(x, 'Haar', m=10)\n",
"print(x)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The reason why the method from this example worked was that the vector we started with had a simple representation in the wavelet basis, actually it equaled\n",
"the coordinates of a basis function in $\\boldsymbol{\\phi}_1$. Usually this is not the case, and our only possibility is to run the DWT on a computer.\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"## Example 4.10: DWT on sound\n",
"\n",
"\n",
"\n",
"Let us plot the samples of our audio sample file, and compare them with the first order DWT."
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"x, fs = audioread('sounds/castanets.wav')\n",
"nvals = arange(0, 2**17)\n",
"plt.plot(nvals[::100],x[:2**17:100,0],'k-')"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"newx=(x[:2**17,0]).copy()\n",
"dwt_impl(newx, 'Haar')\n",
"plt.plot(nvals[::100],newx[::100],'k-')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The first part of the DWT plot represents the low resolution part, the second the detail. \n",
"Since $\\phi(2^mt-n)\\in V_m$ oscillates more quickly than $\\phi(t-n)\\in V_0$,\n",
"one is lead to believe that coefficients from lower order resolution spaces correspond to lower frequencies. \n",
"The functions $\\phi_{m,n}$ do not correspond to pure tones in the setting of wavelets, however, \n",
"but let us nevertheless listen to sound from the different resolution spaces. \n",
"The library includes a function `forw_comp_rev_dwt1` which runs an $m$-level DWT on the first samples \n",
"of the audio sample file, extracts the detail or the low-resolution approximation, and runs an IDWT to reconstruct the sound. \n",
"Since the returned values may lie outside the legal range $[-1,1]$, the values are normalized at the end.\n",
"To listen to the low-resolution approximation, write"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"m = 1"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"x, fs = forw_comp_rev_dwt1(m, 'Haar')\n",
"play(x, fs)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"For $m=2$ we clearly hear a degradation in the sound. For $m=4$ and above most of the sound is unrecognizable, as too much of the detail is omitted. \n",
"To be more precise, when listening to the sound by throwing away detail from $W_0$, $W_1$,...,$W_{m-1}$, we are left with a $2^{-m}$ share of the data. \n",
"\n",
"For $m=1$ and $m=2$ the detail can be played as follows"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"x, fs = forw_comp_rev_dwt1(1, 'Haar', False)\n",
"play(x, fs)"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"x, fs = forw_comp_rev_dwt1(2, 'Haar', False)\n",
"play(x, fs)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"and plotted as follows"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"vals=arange(0,2**17,100)\n",
"plt.plot(vals, x[1::100],'k-')"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"vals=arange(0,2**17,100)\n",
"plt.plot(vals, x[1::100],'k-')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"It is seen that the detail is larger in the part of the sound where there are bigger variations. \n",
"The detail is clearly a very degraded version of the sound, but if you increase $m$, the sound will improve in quality.\n",
"The reason is that more and more information is contained in the detail components for large $m$, which can be seen when we compare the plot for $m=1$ and $m=2$. \n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"The previous example illustrates that wavelets as well may be used to perform operations on sound. In this book the main application for wavelets\n",
"will be images, however, where they have found a more important role. Images typically display variations which are less abrupt than the ones found in sound.\n",
"The main motivation behind wavelets comes from the fact that the detail components often are very small, and in less abrupt data such as images, the detail components will be even smaller.\n",
"After a DWT one is therefore often left with a couple of significant coefficients, while most of the detail is small. This is a very good starting point for compression methods. \n",
"When we later look at how wavelets are applied to images, we will need to handle one final hurdle, namely that images are two-dimensional.\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"## Example 4.11: DWT on the samples of a mathematical function\n",
"\n",
"\n",
"\n",
"Above we plotted the DWT coefficients of a sound, as well as the detail/error. We can also experiment with samples generated from a mathematical function.\n",
"The following code computes and returns the detail for $m=6$ and $m=8$ for the vector `x`:"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"def findetailm6m8(x, wave_name):\n",
" N = len(x)\n",
" m6zeroout = N\n",
" for mres in range(6):\n",
" m6zeroout = ceil(m6zeroout/2)\n",
" m8zeroout = m6zeroout\n",
" m8zeroout = ceil(m8zeroout/2)\n",
" m8zeroout = ceil(m8zeroout/2)\n",
" m6zeroout = int(m6zeroout)\n",
" m8zeroout = int(m8zeroout)\n",
" \n",
" m=6\n",
" detail1=x.copy()\n",
" dwt_impl(detail1, wave_name, m=m)\n",
" detail1[:m6zeroout] = 0\n",
" idwt_impl(detail1, wave_name, m=m)\n",
" m=8\n",
" detail2=x.copy()\n",
" dwt_impl(detail2, wave_name, m=m)\n",
" detail2[:m8zeroout] = 0\n",
" idwt_impl(detail2, wave_name, m=m)\n",
" return detail1, detail2"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The following code computes input vectors corresponding to a piecewise constant, a trigonometric, and a piecewise linear function."
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"N = 2**10\n",
"nvals = linspace(0, N-1, N)\n",
"\n",
"fpwc = ones(N)\n",
"fpwc[:(int(N/2))] *= 0.75\n",
"fpwc[(int(N/2)):] *= 0.25\n",
"\n",
"ftrig = 0.5+0.5*cos(2*pi*nvals/float(N))\n",
"\n",
"flin = nvals*2/float(N-1) \n",
"flin = 1-abs(1-flin)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Finally, the following code compares the computed details with the entries in the vectors."
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"detail1, detail2 = findetailm6m8(fpwc, 'Haar')\n",
"plt.plot(nvals, abs(fpwc),'r.', nvals, abs(detail1), 'g.',nvals,abs(detail2),'b.')"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"detail1, detail2 = findetailm6m8(ftrig, 'Haar')\n",
"plt.plot(nvals,abs(ftrig),'r.',nvals,abs(detail1),'g.',nvals,abs(detail2),'b.')"
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"detail1, detail2 = findetailm6m8(flin, 'Haar') \n",
"plt.plot(nvals,abs(flin),'r.',nvals,abs(detail1),'g.',nvals,abs(detail2),'b.')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In these cases, we see that we require large $m$ before the detail/error becomes significant. We see also that there is no error for the square wave. The reason is\n",
"that the square wave is a piecewise constant function, so that it can be represented exactly by the basis functions. For the other functions this is not\n",
"the case, however, so we get an error.\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"## Exercise 4.17: Computing the DWT of a simple vector\n",
"\n",
"\n",
"\n",
"Suppose that we have the vector $\\boldsymbol{x}$ with length $2^{10}=1024$, defined by $x_n=1$ for $n$ even, $x_n=-1$ for $n$ odd. \n",
"What will be the result if you run a 10-level DWT on $\\boldsymbol{x}$?\n",
"Use the function `dwt_impl` to verify what you have found.\n",
"\n",
"\n",
"\n",
"**Hint.**\n",
"We defined $\\psi$ by $\\psi(t)=(\\phi_{1,0}(t)-\\phi_{1,1}(t))/\\sqrt{2}$. From this connection it follows that \n",
"$\\psi_{9,n}=(\\phi_{10,2n}-\\phi_{10,2n+1})/\\sqrt{2}$, and thus $\\phi_{10,2n}-\\phi_{10,2n+1}=\\sqrt{2}\\psi_{9,n}$. \n",
"Try to couple this identity with the alternating sign you see in $\\boldsymbol{x}$.\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"**Solution.**\n",
"The vector $\\boldsymbol{x}$ is the coordinate vector of the function $f(t)=\\sum_{n=0}^{1023} (-1)^n \\phi_{10,n}$ in the basis $\\boldsymbol{\\phi}_{10}$ for $V_{10}$. \n",
"Since $\\phi_{10,2n}-\\phi_{10,2n+1}=\\sqrt{2}\\psi_{9,n}$, we can write $f(t)=\\sum_{n=0}^{511} \\sqrt{2}\\psi_{9,n}$. \n",
"Since a 10-level-DWT gives as a result the coordinate vector of $f$ in"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"(\\boldsymbol{\\phi}_0,\\boldsymbol{\\psi}_0,\\boldsymbol{\\psi}_1,\\boldsymbol{\\psi}_2,\\boldsymbol{\\psi}_3,\\boldsymbol{\\psi}_4,\\boldsymbol{\\psi}_5,\\boldsymbol{\\psi}_6,\\boldsymbol{\\psi}_7,\\boldsymbol{\\psi}_8,\\boldsymbol{\\psi}_9),\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"(the DWT is nothing but the change of coordinates from $\\boldsymbol{\\phi}_{10}$ to this basis), and since $f(t)=\\sum_{n=0}^{511} \\sqrt{2}\\psi_{9,n}$, it is clear that the\n",
"coordinate vector of $f$ in this basis has $\\sqrt{2}$ in the second part (the $\\boldsymbol{\\psi}_9$-coordinates), and $0$ elsewhere.\n",
"The 10-level DWT of $\\boldsymbol{x}$ therefore gives the vector of length 1024 which is $0$ on the first half, and equal to $\\sqrt{2}$ on the second half. \n",
"$m=10$ is here arbitrarily chosen: The result would have been the same for $m=1,m=2,$ and so on.\n",
"The following code verifies the result:"
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"x = tile([1.,-1.], 512)\n",
"dwt_impl(x, 'Haar', m=10)\n",
"print(x)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"## Exercise 4.19: In-place DWT and partial bit-reversal\n",
"\n",
"\n",
"\n",
"The kernel transformations in this book can be computed in-place. A DWT is also required to reorder coordinates as in the $(\\boldsymbol{\\phi}_0,\\boldsymbol{\\psi}_0,...,\\boldsymbol{\\psi}_{m-1})$-basis, however, and this exercise addresses \n",
"how this reordering also can be computed in-place. As a result the entire DWT can be computed in-place.\n",
"\n",
"\n",
"**a)**\n",
"Show that the coordinates in $\\boldsymbol{\\phi}_0$ after $m$ iterations of a kernel transformation end up at indices $k2^m$, $k=0,1,2,\\ldots$, and that the coordinates in $\\boldsymbol{\\psi}_0$ end up at indices $2^{m-1} + k2^m$, $k=0,1,2,\\ldots$.\n",
"\n",
"\n",
"\n",
"**Solution.**\n",
"This is easily shown by induction. \n",
"For $m=1$ the statement says that the $\\boldsymbol{\\phi}_0$-coordinates are at the even indices, and the $\\psi_0$-coordinates are at the odd indices. Both these statements are clearly true. \n",
"For the induction step, if the coordinates in $\\boldsymbol{\\phi}_{1}$ after $m-1$ iterations of the kernel are at the indices $k2^{m-1}$, \n",
"after the next iteration every second element of these (i.e. the ones at indices $2^{m-1}\\cdot 2k =k2^{m}$) will be coordinates in $\\boldsymbol{\\phi}_0$. \n",
"while the ones at indices $2^{m-1}(2k+1)=2^{m-1}+k2^m$ will be coordinates in $\\psi_0$.\n",
"\n",
"\n",
"\n",
"**b)**\n",
"After $m$ iterations of a kernel transformation, show that if the indices are rearranged so that the last $m$ bits are placed in front in reverse order, \n",
"then coordinates from $\\boldsymbol{\\phi}_0$ are placed first, while coordinates from $\\boldsymbol{\\psi}_i$ are placed before those from $\\{\\boldsymbol{\\psi}_j\\}_{j>i}$. \n",
"This is also called a *partial bit-reversal*. When an $m$-level DWT on a vector of length $2^m$ is computed, the partial bit-reversal is actually a full bit-reversal. \n",
"\n",
"After a partial bit-reversal, what can you say about the internal ordering in $\\{\\boldsymbol{\\psi}_j\\}_{j>i}$?\n",
"\n",
"\n",
"\n",
"\n",
"**Solution.**\n",
"As proved in a), the $\\boldsymbol{\\psi}_{m-n}$-coordinates end up at indices $2^{n-1} + k2^{n}$. After a partial bit-reversal, their bit representations have $0\\cdots 01$ as the first $n$ bits. \n",
"Lower $n$ thus implies higher numbers, so that coordinates from $\\boldsymbol{\\psi}_i$ are placed before those from $\\{\\boldsymbol{\\psi}_j\\}_{j>i}$. \n",
"Since $\\boldsymbol{\\phi}_0$-coordinates end up at indices $k2^m$, their first $m$ bits are $0\\cdots 0$, so that these are placed first. \n",
"\n",
"The internal ordering of $\\boldsymbol{\\psi}_j$ will be violated by a partial bit-reversal, however. If one requires an in-place ordering that produces the same order as $(\\boldsymbol{\\phi}_0,\\boldsymbol{\\psi}_0,...,\\boldsymbol{\\psi}_{m-1})$, one should secure that, \n",
"for the $\\boldsymbol{\\psi}_{m-j}$-coordinates, where the last $j$ bits are $0\\cdots 01$, only those $j$ bits should be placed in front in reverse order. Thus, partial bit-reversal reverses too many bits, violating the internal order in the $\\boldsymbol{\\psi}_{m-j}$. \n",
"Reversing the last bits and placing them in front should thus be replaced with a more refined version which adapts to the different resolutions.\n",
"\n",
"\n",
"\n",
"**c)**\n",
"Write a function `partial_bit_reversal(x, m)` which computes a partial bit-reversal of the vector $\\boldsymbol{x}$ in-place. Assume that $\\boldsymbol{x}$ has length $2^n$ with $n\\geq m$.\n",
"\n",
"\n",
"\n",
"**Hint.**\n",
"Partial bit-reversal can be split in two: First a full bit-reversal, and then reversing the last $n-m$ bits again. It turns out that reversal of the last bits can be computed easily using bit-reversal of $n-m$ bits. \n",
"Due to this we can easily implement partial bit-reversal by adapting our previous function for bit-reversal.\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"**Solution.**\n",
"The general code can be found in the solution to d).\n",
"\n",
"\n",
"\n",
"**d)**\n",
"Contrary to full bit-reversal, partial bit-reversal is not its own inverse. It is straightforward, however, to compute the inverse following the same lines as in c). Extend the function `partial_bit_reversal` so that \n",
"it takes takes a third parameter `forward`, which indicates whether a forward- or reverse partial bit-reversal should be computed.\n",
"\n",
"\n",
"\n",
"**Solution.**\n",
"The code can look as follows. The code calls a function which is a slight modification of `bit_reversal`, adapated to the current setting."
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"def partial_bit_reversal(x, m, forward):\n",
" N = shape(x)[0]\n",
" if forward:\n",
" bit_reversal_last_bits(x, int(log2(N)))\n",
" bit_reversal_last_bits(x, m)\n",
" else:\n",
" bit_reversal_last_bits(x, m)\n",
" bit_reversal_last_bits(x, int(log2(N)))\n",
"\n",
"def bit_reversal_last_bits(x, num_bits):\n",
" N = 2**num_bits\n",
" j = 0\n",
" temp=zeros(int(len(x)/N))\n",
" for i in range(0, int(N/2), 2):\n",
" if (j > i): \n",
" temp[:] = x[j::2**num_bits]\n",
" x[j::2**num_bits] = x[i::2**num_bits]\n",
" x[i::2**num_bits] = temp\n",
" \n",
" temp[:] = x[(N - j - 1)::2**num_bits]\n",
" x[(N - j - 1)::2**num_bits] = x[(N - i - 1)::2**num_bits]\n",
" x[(N - i - 1)::2**num_bits] = temp\n",
" temp[:] = x[(i+1)::2**num_bits]\n",
" x[(i+1)::2**num_bits] = x[(j + int(N/2))::2**num_bits]\n",
" x[(j + int(N/2))::2**num_bits] = temp\n",
" \n",
" m = int(N/4)\n",
" while (m >= 1 and j >= m):\n",
" j -= m\n",
" m =int(m/2)\n",
" j += m"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"## Example 4.20: DWT on sound\n",
"\n",
"\n",
"\n",
"With the new wavelet, let us plot and listen to the new low resolution approximations and detail, as in Example 4.10. First we listen to the low-resolution approximation."
]
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"m = 1"
]
},
{
"cell_type": "code",
"execution_count": 20,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"x, fs = forw_comp_rev_dwt1(m, 'pwl0')\n",
"play(x, fs)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"There is a new and undesired effect when we increase $m$ here: The castanet sound seems to grow strange. \n",
"The sounds from the castanets are perhaps the sound with the highest frequencies. \n",
"\n",
"For $m=1$ and $m=2$ the detail can be played as follows"
]
},
{
"cell_type": "code",
"execution_count": 21,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"x, fs = forw_comp_rev_dwt1(1, 'pwl0', False)\n",
"play(x, fs)"
]
},
{
"cell_type": "code",
"execution_count": 22,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"x, fs = forw_comp_rev_dwt1(2, 'pwl0', False)\n",
"play(x, fs)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The detail can be plotted as follows"
]
},
{
"cell_type": "code",
"execution_count": 23,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"vals=arange(0,2**17,100)\n",
"plt.plot(vals,x[::100],'k-')"
]
},
{
"cell_type": "code",
"execution_count": 24,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"vals=arange(0,2**17,100)\n",
"plt.plot(vals,x[1::100],'k-')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"When comparing with Example 4.10 we see much of the same, but it seems here that the error is bigger than before. \n",
"In the next section we will try to explain why this is the case, and attempt to modify the definition of $\\psi$ to remedy this.\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"## Example 4.21: DWT on the samples of a mathematical function\n",
"\n",
"\n",
"\n",
"Let us also repeat Example 4.11, where we plotted the detail/error at different resolutions for the samples of a mathematical function. \n",
"\n",
"We first compute new input vectors corresponding to a piecewise constant, a trigonometric, and a piecewise linear function."
]
},
{
"cell_type": "code",
"execution_count": 25,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"N = 2**10 + 1\n",
"nvals = linspace(0, N-1, N)\n",
"\n",
"fpwl = ones(N)\n",
"fpwl[:int((N+1)/2)] *= 0.75\n",
"fpwl[int((N+1)/2):] *= 0.25\n",
"\n",
"ftrig = 0.5+0.5*cos(2*pi*nvals/float(N))\n",
"\n",
"flin = nvals*2/float(N-1) \n",
"flin = 1-abs(1-flin)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The following code compares the computed details with the entries in the new vectors."
]
},
{
"cell_type": "code",
"execution_count": 26,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"detail1, detail2 = findetailm6m8(fpwl, 'pwl0')\n",
"plt.plot(nvals,abs(fpwl),'r.',nvals,abs(detail1),'g.',nvals,abs(detail2),'b.')"
]
},
{
"cell_type": "code",
"execution_count": 27,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"detail1, detail2 = findetailm6m8(ftrig, 'pwl0')\n",
"plt.plot(nvals,abs(ftrig),'r.',nvals,abs(detail1),'g.',nvals,abs(detail2),'b.')"
]
},
{
"cell_type": "code",
"execution_count": 28,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"f = nvals*2/float(N-1) \n",
"f = 1-abs(1-f)\n",
"detail1, detail2 = findetailm6m8(flin, 'pwl0')\n",
"plt.plot(nvals,abs(flin),'r.',nvals,abs(detail1),'g.',nvals,abs(detail2),'b.')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"With the square wave we see now that there is an error. The reason is that a piecewise constant function can not be represented exactly by piecewise linear\n",
"functions, due to discontinuity. For the second function we see that there is no error, since this function is piecewise linear, so there is no error when we represent the function from the space $V_0$. \n",
"With the third function we see an error as before: a trigonometric function can not be represented exactly by piecewise constant nor piecewise linear functions.\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"## Exercise 4.25: Computing projections\n",
"\n",
"\n",
"\n",
"In this exercise we will show how the projection of $\\phi_{1,1}$ onto $V_0$ can be computed. We will see from this that it is nonzero, and that its support is the entire $[0,N]$. \n",
"Let $f=\\text{proj}_{V_0}\\phi_{1,1}$, and let $x_n=f(n)$ for $0\\leq n < N$. This means that, on $(n,n+1)$, $f(t)=x_n+(x_{n+1}-x_n)(t-n)$.\n",
"\n",
"\n",
"**a)**\n",
"Show that $\\int_n^{n+1} f(t)^2dt=(x_n^2+x_nx_{n+1}+x_{n+1}^2)/3$.\n",
"\n",
"\n",
"\n",
"**Solution.**\n",
"We have that"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\begin{align*}\n",
"\\int_n^{n+1} f(t)^2dt &= \\int_n^{n+1} (x_n+(x_{n+1}-x_n)(t-n))^2 dt = \\int_0^1 (x_n+(x_{n+1}-x_n)t)^2 dt \\\\ \n",
"&= \\int_0^1 (x_n^2 + 2x_n(x_{n+1}-x_n)t + (x_{n+1}-x_n)^2t^2) dt \\\\ \n",
"&= \\left[ x_n^2t + x_n(x_{n+1}-x_n)t^2 + (x_{n+1}-x_n)^2t^3/3 \\right]_0^1 \\\\ \n",
"&= x_n^2 + x_n(x_{n+1}-x_n) + (x_{n+1}-x_n)^2/3 = \\frac{1}{3}(x_n^2+x_nx_{n+1}+x_{n+1}^2).\n",
"\\end{align*}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"**b)**\n",
"Show that"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\begin{align*}\n",
"\\int_0^{1/2} (x_0+(x_1-x_0)t)\\phi_{1,1}(t)dt &= 2\\sqrt{2}\\left(\\frac{1}{12}x_0 + \\frac{1}{24}x_1\\right) \\\\ \n",
"\\int_{1/2}^1 (x_0+(x_1-x_0)t)\\phi_{1,1}(t)dt &= 2\\sqrt{2}\\left(\\frac{1}{24}x_0 + \\frac{1}{12}x_1\\right).\n",
"\\end{align*}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"**Solution.**\n",
"We have that"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\begin{align*}\n",
"\\lefteqn{\\int_0^{1/2} (x_0+(x_1-x_0)t)\\phi_{1,1}(t)dt} \\\\ \n",
"&= \\int_0^{1/2} (x_0+(x_1-x_0)t)2\\sqrt{2}tdt = 2\\sqrt{2}\\int_0^{1/2} (x_0t+(x_1-x_0)t^2)dt \\\\ \n",
"&= 2\\sqrt{2}\\left[ \\frac{1}{2}x_0t^2 + \\frac{1}{3}(x_1-x_0)t^3 \\right]_0^{1/2} = 2\\sqrt{2}\\left( \\frac{1}{8}x_0 + \\frac{1}{24}(x_1-x_0))\\right) \\\\ \n",
"&= 2\\sqrt{2}\\left(\\frac{1}{12}x_0 + \\frac{1}{24}x_1\\right).\n",
"\\end{align*}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In the same way"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\begin{align*}\n",
"& \\int_{1/2}^1 (x_0+(x_1-x_0)t)\\phi_{1,1}(t)dt \\\\ \n",
"&= \\int_{1/2}^1 (x_0+(x_1-x_0)t)2\\sqrt{2}(1-t)tdt = 2\\sqrt{2}\\int_{1/2}^1 (x_0 + (x_1-2x_0)t - (x_1-x_0)t^2)tdt \\\\ \n",
"&= 2\\sqrt{2}\\left[x_0t + \\frac{1}{2}(x_1-2x_0)t^2 - \\frac{1}{3}(x_1-x_0)t^3 \\right]_{1/2}^1 = 2\\sqrt{2}\\left(\\frac{1}{2}x_0 + \\frac{3}{8}(x_1-2x_0) - \\frac{7}{24}(x_1-x_0) \\right) \\\\ \n",
"&= 2\\sqrt{2}\\left(\\frac{1}{24}x_0 + \\frac{1}{12}x_1\\right).\n",
"\\end{align*}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"**c)**\n",
"Use the fact that"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\begin{align*}\n",
"& \\int_0^N (\\phi_{1,1}(t)-\\sum_{n=0}^{N-1} x_n\\phi_{0,n}(t))^2 dt \\\\ \n",
"&= \\int_0^1 \\phi_{1,1}(t)^2dt - 2\\int_0^{1/2}(x_0+(x_1-x_0)t)\\phi_{1,1}(t)dt -2\\int_{1/2}^1 (x_0+(x_1-x_0)t)\\phi_{1,1}(t)dt \\\\ \n",
"&+ \\sum_{n=0}^{N-1}\\int_n^{n+1} (x_n+(x_{n-1}-x_n)t)^2dt\n",
"\\end{align*}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"and a) and b) to find an expression for $\\|\\phi_{1,1}(t)-\\sum_{n=0}^{N-1} x_n\\phi_{0,n}(t)\\|^2$.\n",
"\n",
"\n",
"\n",
"**Solution.**\n",
"Using a) and b) we see that the above can be written as"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\begin{align*}\n",
"& \\frac{2}{3} + \\sum_{n=0}^{N-1} \\frac{1}{3}(x_n^2+x_nx_{n+1}+x_{n+1}^2)- 2\\left( 2\\sqrt{2}\\left(\\frac{1}{12}x_0 + \\frac{1}{24}x_1\\right) - 2\\sqrt{2}\\left(\\frac{1}{24}x_0 +\n",
"\\frac{1}{12}x_1\\right)\\right)\n",
"\\\\ \n",
"&= \\frac{2}{3} + \\frac{2}{3}\\sum_{n=0}^{N-1} x_n^2 + \\frac{1}{3}\\sum_{n=0}^{N-1} x_nx_{n+1} - \\frac{\\sqrt{2}}{2}(x_0 + x_1).\n",
"\\end{align*}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"**d)**\n",
"To find the minimum least squares error, we can set the gradient of the expression in c) to zero, and thus find the expression for the projection of $\\phi_{1,1}$ onto $V_0$. \n",
"Show that the values $\\{x_n\\}_{n=0}^{N-1}$ can be found by solving the equation $S\\boldsymbol{x}=\\boldsymbol{b}$, where $S=\\frac{1}{3}\\{1,\\underline{4},1\\}$ is an $N\\times N$ symmetric filter, and $\\boldsymbol{b}$ is\n",
"the vector with components $b_0=b_1=\\sqrt{2}/2$, and $b_k=0$ for $k\\geq 2$.\n",
"\n",
"\n",
"\n",
"**Solution.**\n",
"We see that the partial derivatives of the function in c) are"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\begin{align*}\n",
"\\frac{\\partial f}{\\partial x_0} &= \\frac{1}{3}x_{N-1} + \\frac{4}{3}x_0 + \\frac{1}{3}x_1 - \\frac{\\sqrt{2}}{2} \\\\ \n",
"\\frac{\\partial f}{\\partial x_1} &= \\frac{1}{3}x_0 + \\frac{4}{3}x_1 + \\frac{1}{3}x_1 - \\frac{\\sqrt{2}}{2} \\\\ \n",
"\\frac{\\partial f}{\\partial x_i} &= \\frac{1}{3}x_{i-1} + \\frac{4}{3}x_i + \\frac{1}{3}x_{i+1} \\text{ $2\\leq i < N-1$}\\\\ \n",
"\\frac{\\partial f}{\\partial x_{N-1}} &= \\frac{1}{3}x_{N-2} + \\frac{4}{3}x_{N-1} + \\frac{1}{3}x_0.\n",
"\\end{align*}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Moving the two terms $\\frac{\\sqrt{2}}{2}$ over to the right hand side, setting the gradient equal to zero is the same as solving the system $S\\boldsymbol{x}=\\boldsymbol{b}$ which we stated.\n",
"\n",
"\n",
"\n",
"**e)**\n",
"Solve the system in d. for some values of $N$ to verify that the projection of $\\phi_{1,1}$ onto $V_0$ is nonzero, and that its support covers the entire $[0,N]$.\n",
"\n",
"\n",
"\n",
"**Solution.**\n",
"The following code can be used"
]
},
{
"cell_type": "code",
"execution_count": 29,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"N = 16\n",
"S = zeros((N, N))\n",
"S[0,N-1] = 1/3.; S[0,0] = 4/3.; S[0,1] = 1/3.; # First row\n",
"for k in range(1,N-1):\n",
" S[k,(k-1):(k+2)] = [1/3., 4/3., 1/3.]\n",
"S[N-1,N-2] = 1/3.; S[N-1,N-1]=4/3.; S[N-1,0]=1/3.; # Last row\n",
"b=zeros(N); b[0]=sqrt(2)/2; b[1]=sqrt(2)/2;\n",
"plt.plot(range(0,N),linalg.solve(S,b)) # Plots the projection"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"## Example 4.27: DWT on sound\n",
"\n",
"\n",
"\n",
"Using the new kernels, let us again listen to the low resolution approximations and the detail. First the low-resolution approximation:"
]
},
{
"cell_type": "code",
"execution_count": 30,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"m = 1"
]
},
{
"cell_type": "code",
"execution_count": 31,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"x, fs = forw_comp_rev_dwt1(m, 'pwl2')\n",
"play(x, fs)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The undesired effect in the castanets from Example 4.20 seems to be gone. \n",
"The detail for $m=1$ and $m=2$ can be played as follows"
]
},
{
"cell_type": "code",
"execution_count": 32,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"x, fs = forw_comp_rev_dwt1(1, 'pwl2', False)\n",
"play(x, fs)"
]
},
{
"cell_type": "code",
"execution_count": 33,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"x, fs = forw_comp_rev_dwt1(2, 'pwl2', False)\n",
"play(x, fs)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"and plotted as follows"
]
},
{
"cell_type": "code",
"execution_count": 34,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"vals=arange(0,2**17,100)\n",
"plt.plot(vals,x[1::100],'k-')"
]
},
{
"cell_type": "code",
"execution_count": 35,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"vals=arange(0,2**17,100)\n",
"plt.plot(vals,x[1::100],'k-')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Again, when comparing with Example 4.10 we see much of\n",
"the same. It is difficult to see an improvement from this figure. However, this figure also clearly shows a smaller error than the piecewise linear wavelet. \n",
"A partial explanation is that the wavelet we now have constructed has two vanishing moments, while the other had not.\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"## Example 4.28: DWT on the samples of a mathematical function\n",
"\n",
"\n",
"\n",
"Let us also repeat Example 4.11 for our alternative wavelet, where we plotted the detail/error at different resolutions for the samples of a\n",
"mathematical function. \n",
"The following code compares the computed details with the entries in the vectors."
]
},
{
"cell_type": "code",
"execution_count": 36,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"detail1, detail2 = findetailm6m8(fpwl, 'pwl2')\n",
"plt.plot(nvals,abs(fpwl),'r.', nvals,abs(detail1),'g.', nvals,abs(detail2),'b.')"
]
},
{
"cell_type": "code",
"execution_count": 37,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"detail1, detail2 = findetailm6m8(ftrig, 'pwl2')\n",
"plt.plot(nvals,abs(ftrig),'r.',nvals,abs(detail1),'g.',nvals,abs(detail2),'b.')"
]
},
{
"cell_type": "code",
"execution_count": 38,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"detail1, detail2 = findetailm6m8(flin, 'pwl2')\n",
"plt.plot(nvals,abs(flin),'r.',nvals, abs(detail1),'g.',nvals, abs(detail2),'b.')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Again for the square wave there is an error, which seems to be slightly lower than for the previous wavelet. For the second function we see that there is no error,\n",
"as before. The reason is the same as before, since the function is piecewise linear. With the third function there is an error. The error seems to be slightly\n",
"lower than for the previous wavelet, which fits well with the increased number of vanishing moments.\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"## Exercise 4.31: More than two vanishing moments\n",
"\n",
"\n",
"\n",
"In the previous exercise we ended up with a lot of calculations to find $\\alpha,\\beta$ in \n",
"equation (4.30). Let us try to make a program which does\n",
"this for us, and which also makes us able to generalize the result.\n",
"\n",
"\n",
"**a)**\n",
"Define"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\begin{align*}\n",
"a_k &= \\int_{-1}^1 t^k(1-|t|)dt, & b_k &= \\int_0^2 t^k(1-|t-1|)dt, & e_k &= \\int_0^1 t^k(1-2|t-1/2|)dt,\n",
"\\end{align*}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"for $k\\geq 0$. Explain why finding $\\alpha,\\beta$ so that we have two vanishing moments in equation \n",
"(4.30) \n",
"is equivalent to solving the following equation:"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\begin{pmatrix}\n",
"a_0 & b_0 \\\\ \n",
"a_1 & b_1\n",
"\\end{pmatrix}\n",
"\\begin{pmatrix} \\alpha \\\\ \\beta \\end{pmatrix}\n",
"=\n",
"\\begin{pmatrix}\n",
"e_0 \\\\ e_1\n",
"\\end{pmatrix}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Write a program which sets up and solves this system of equations, and use this program to verify the values for $\\alpha,\\beta$ we previously have found.\n",
"\n",
"\n",
"\n",
"**Hint.**\n",
"you can integrate functions in Python \n",
"with the function `quad` in the package\n",
"`scipy.integrate`.\n",
"As an example, the function $\\phi(t)$, which is nonzero only on $[-1,1]$, can be integrated as follows:"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
" res, err = quad(lambda t: t**k*(1-abs(t)), -1, 1)\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"\n",
"\n",
"**Solution.**\n",
"In order for $\\psi$ to have vanishing moments we must have that $\\int\\hat\\psi(t)dt=\\int t\\hat\\psi(t)dt=0$\n",
"Substituting $\\hat\\psi=\\psi-\\alpha \\phi_{0,0}-\\beta \\phi_{0,1}$ we see that, for $k=0,1$,"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\int t^k\\left(\\alpha \\phi_{0,0}+\\beta \\phi_{0,1}\\right)dt = \\int t^k\\psi(t)dt.\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The left hand side can here be written"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\begin{align*}\n",
"\\int t^k\\left(\\alpha \\phi_{0,0}+\\beta \\phi_{0,1}\\right)dt &= \\alpha\\int t^k\\phi_{0,0}dt + \\beta\\int t^k\\phi_{0,1}(t)dt \\\\ \n",
"&=\\alpha\\int_{-1}^1 t^k(1-|t|)dt + \\beta\\int_0^2 t^k(1-|t-1|)dt =\\alpha a_k + \\beta b_k.\n",
"\\end{align*}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The right hand side is"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\int t^k\\psi(t)dt = \\int t^k\\phi_{1,1}(t)dt = \\int_0^1 (1-2|t-1/2|)dt=e_k.\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The following program sets up the corresponding equation systems, and solves it by finding $\\alpha,\\beta$."
]
},
{
"cell_type": "code",
"execution_count": 39,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"A = zeros((2, 2))\n",
"b = zeros((2, 1))\n",
"for k in range(2):\n",
" res1, err1 = integrate.quad(lambda t: t**k*(1-abs(t)), -1, 1)\n",
" res2, err2 = integrate.quad(lambda t: t**k*(1-abs(t-1)), 0, 2)\n",
" res3, err3 = integrate.quad(lambda t: t**k*(1-2*abs(t-1/2.)), 0, 1)\n",
" A[k,:] = [res1, res2]\n",
" b[k] = res3\n",
"linalg.solve(A,b)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"**b)**\n",
"The procedure where we set up a matrix equation in a) allows for generalization to more vanishing moments. Define"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\hat\\psi=\\psi_{0,0}-\\alpha \\phi_{0,0}-\\beta \\phi_{0,1} -\\gamma\\phi_{0,-1}-\\delta\\phi_{0,2}.\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We would like to choose $\\alpha,\\beta,\\gamma,\\delta$ so that we have $4$ vanishing moments. Define also"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\begin{align*}\n",
"g_k &= \\int_{-2}^0 t^k(1-|t+1|)dt, & d_k &= \\int_1^3 t^k(1-|t-2|)dt\n",
"\\end{align*}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"for $k\\geq 0$. Show that $\\alpha,\\beta,\\gamma,\\delta$ must solve the equation"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\begin{pmatrix}\n",
"a_0 & b_0 & g_0 & d_0 \\\\ \n",
"a_1 & b_1 & g_1 & d_1 \\\\ \n",
"a_2 & b_2 & g_2 & d_2 \\\\ \n",
"a_3 & b_3 & g_3 & d_3\n",
"\\end{pmatrix}\n",
"\\begin{pmatrix}\n",
"\\alpha \\\\ \\beta \\\\ \\gamma\\\\ \\delta\n",
"\\end{pmatrix}\n",
"=\n",
"\\begin{pmatrix}\n",
"e_0 \\\\ e_1 \\\\ e_2 \\\\ e_3\n",
"\\end{pmatrix},\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"and solve this with your computer.\n",
"\n",
"\n",
"\n",
"**Solution.**\n",
"Similarly to a), the equation in b) gives that"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\int t^k\\left(\\alpha \\phi_{0,0}+\\beta \\phi_{0,1} +\\gamma\\phi_{0,-1}+\\delta\\phi_{0,2}\\right)dt = \\int t^k\\psi(t)dt.\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The corresponding equation system is deduced exactly as in a). \n",
"The following program sets up the corresponding equation systems, and solves it by finding $\\alpha,\\beta,\\gamma,\\delta$."
]
},
{
"cell_type": "code",
"execution_count": 40,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"A=zeros((4, 4))\n",
"b=zeros((4, 1))\n",
"for k in range(4):\n",
" res1, err1 = integrate.quad(lambda t: t**k*(1-abs(t)), -1, 1)\n",
" res2, err2 = integrate.quad(lambda t: t**k*(1-abs(t-1)), 0, 2)\n",
" res3, err3 = integrate.quad(lambda t: t**k*(1-abs(t+1)), -2, 0)\n",
" res4, err4 = integrate.quad(lambda t: t**k*(1-abs(t-2)), 1, 3)\n",
" res5, err5 = integrate.quad(lambda t: t**k*(1-2*abs(t-1/2.)), 0, 1)\n",
" A[k,:] = [res1, res2, res3, res4]\n",
" b[k] = res5\n",
"coeffs = linalg.solve(A,b)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"**c)**\n",
"Plot the function defined by the equation from b).\n",
"\n",
"\n",
"\n",
"**Hint.**\n",
"If `t` is the vector of $t$-values, and you write"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
" (t >= 0)*(t <= 1)*(1-2*abs(t-0.5))\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"you get the points $\\phi_{1,1}(t)$.\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"**Solution.**\n",
"The function $\\hat\\psi$ now is supported on $[-2,3]$, and can be plotted as follows:"
]
},
{
"cell_type": "code",
"execution_count": 41,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"t=linspace(-2,3,100)\n",
"plt.plot(t, (t >= 0)*(t <= 1)*(1-2*abs(t - 0.5)) \\\n",
" -coeffs[0]*(t >= -1)*(t <= 1)*(1 - abs(t)) \\\n",
" -coeffs[1]*(t >= 0)*(t <= 2)*(1 - abs(t - 1)) \\\n",
" -coeffs[2]*(t >= -2)*(t <= 0)*(1 - abs(t + 1)) \\\n",
" -coeffs[3]*(t >= 1)*(t <= 3)*(1 - abs(t - 2)))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"**d)**\n",
"Explain why the coordinate vector of $\\hat{\\psi}$ in the basis $(\\boldsymbol{\\phi}_0,\\boldsymbol{\\psi}_0)$ is"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"[\\hat{\\psi}]_{(\\boldsymbol{\\phi}_0,\\boldsymbol{\\psi}_0)}=(-\\alpha,-\\beta,-\\delta,0,\\ldots,0-\\gamma)\\oplus(1,0,\\ldots,0).\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"**Hint.**\n",
"The placement of $-\\gamma$ may seem a bit strange here, and has to do with that $\\phi_{0,-1}$\n",
"is not one of the basis functions $\\{\\phi_{0,n}\\}_{n=0}^{N-1}$. However, we have that $\\phi_{0,-1}=\\phi_{0,N-1}$, i.e. $\\phi(t+1)=\\phi(t-N+1)$, since we always\n",
"assume that the functions we work with have period $N$.\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"**Solution.**\n",
"The coordinates of $\\phi_{0,0}$, $\\phi_{0,1}$, and $\\phi_{0,2}$ in $\\hat\\psi$ are $-\\alpha$, $-\\beta$, and $-\\delta$, respectively. Since these three basis functions occur first in $(\\boldsymbol{\\phi}_0,\\boldsymbol{\\psi}_0)$, they appear first in $[\\hat{\\psi}]_{(\\boldsymbol{\\phi}_0,\\boldsymbol{\\psi}_0)}$. \n",
"Since $\\phi_{0,-1}=\\phi_{0,N-1}$, and since this is the last basis function in the first half of $(\\boldsymbol{\\phi}_0,\\boldsymbol{\\psi}_0)$, $-\\gamma$ (the corresponding coordinate) is placed accordingly. \n",
"Finally, $\\psi_{0,0}$ is th first basis function in the second half of $(\\boldsymbol{\\phi}_0,\\boldsymbol{\\psi}_0)$, so that a $1$ (the coordinate) should be placed at that location.\n",
"\n",
"\n",
"\n",
"**e)**\n",
"Compute the coordinates of $\\hat{\\psi}$ in the basis $\\boldsymbol{\\phi}_1$ (i.e. $[\\hat\\psi]_{\\boldsymbol{\\phi}_1}$) with $N=8$, i.e. compute the IDWT of"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"[\\hat{\\psi}]_{(\\boldsymbol{\\phi}_0,\\boldsymbol{\\psi}_0)}=(-\\alpha,-\\beta,-\\delta,0,0,0,0,-\\gamma)\\oplus(1,0,0,0,0,0,0,0),\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"which is the coordinate vector you computed in d).\n",
"For this, you should use the function `idwt_impl`, with the kernel of the piecewise linear wavelet without symmetric extension as input.\n",
"\n",
"\n",
"\n",
"**Solution.**\n",
"The code which can be used looks like this:"
]
},
{
"cell_type": "code",
"execution_count": 42,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"g1 = array([-coeffs[0],-coeffs[1],-coeffs[3],0,0,0,0,-coeffs[2],\\\n",
" 1, 0, 0, 0, 0, 0, 0, 0])\n",
"idwt_impl(g1, 'pwl0', bd_mode='per')\n",
"g1 = hstack([g1[13:16],g1[0:6]]) # Compact filter notation"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Note that we have used a kernel which does not make symmetric extensions.\n",
"\n",
"\n",
"\n",
"**f)**\n",
"Sketch a more general procedure than the one you found in b), which can be used to find wavelet bases where we have even more vanishing moments.\n",
"\n",
"\n",
"\n",
"**Solution.**\n",
"If we define"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\hat\\psi=\\psi_{0,0} - \\sum_{k=0}^K \\left( \\alpha_k \\phi_{0,-k} -\\beta_k\\phi_{0,k+1}\\right),\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"we have $2k$ unknowns. These can be determined if we require $2k$ vanishing moments.\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"## Exercise 4.32: Two vanishing moments for the Haar wavelet\n",
"\n",
"\n",
"\n",
"Let $\\phi$ be the scaling function of the Haar wavelet.\n",
"\n",
"\n",
"**a)**\n",
"Compute $\\text{proj}_{V_0}f$, where $f(t)=t^2$, and where $f$ is defined on $[0,N)$.\n",
"\n",
"\n",
"\n",
"**Solution.**\n",
"We have that"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\begin{align*}\n",
"\\text{proj}_{V_0}f &= \\sum_{n=0}^{N-1} \\langle t^2, \\phi_{0,n}\\rangle \\phi_{0,n} = \\sum_{n=0}^{N-1} \\int_n^{n+1} t^2 dt \\phi_{0,n} \\\\ \n",
"&= \\sum_{n=0}^{N-1} \\frac{1}{3}((n+1)^3-n^3)\\phi_{0,n} = \\sum_{n=0}^{N-1} (n^2+n+1/3)\\phi_{0,n}. \n",
"\\end{align*}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"**b)**\n",
"Find constants $\\alpha,\\beta$ so that $\\hat\\psi=\\psi-\\alpha\\phi_{0,0}-\\beta\\phi_{0,1}$ has two vanishing moments, i.e. so that \n",
"$\\langle \\hat\\psi,1\\rangle=\\langle \\hat\\psi,t\\rangle=0$. Plot also the function $\\hat\\psi$.\n",
"\n",
"\n",
"\n",
"**Hint.**\n",
"Start with computing the integrals $\\int\\psi(t)dt$, $\\int t\\psi(t)dt$, $\\int\\phi_{0,0}(t)dt$, $\\int \\phi_{0,1}(t)dt$, and \n",
"$\\int t\\phi_{0,0}(t)dt$, $\\int t\\phi_{0,1}(t)dt$.\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"**Solution.**\n",
"We get that"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\begin{align*}\n",
"\\int\\psi(t)dt &= 0 \\\\ \n",
"\\int t\\psi(t)dt &= \\int_{0}^{1/2} tdt - \\int_{1/2}^{1} tdt = [t^2/2]_0^{1/2} - [t^2/2]_{1/2}^{1} = 1/8 - 1/2 + 1/8 = -1/4 \\\\ \n",
"\\int \\phi_{0,0}(t)dt &= \\int \\phi_{0,1}(t)dt = 1 \\\\ \n",
"\\int t\\phi_{0,0}(t)dt &= \\int_0^1 tdt = 1/2 \\\\ \n",
"\\int t\\phi_{0,1}(t)dt &= \\int_1^2 tdt = 3/2.\n",
"\\end{align*}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$\\langle \\hat\\psi,1\\rangle=\\langle \\hat\\psi,t\\rangle=0$ can be written as the system"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\begin{align*}\n",
"\\alpha\\langle\\phi_{0,0},1\\rangle + \\beta\\langle\\phi_{0,1},1\\rangle &= \\langle \\psi,1\\rangle \\\\ \n",
"\\alpha\\langle\\phi_{0,0},t\\rangle + \\beta\\langle\\phi_{0,1},t\\rangle &= \\langle \\psi,t\\rangle\n",
"\\end{align*}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"so that"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\begin{align*}\n",
"\\alpha + \\beta &= 0 & \\alpha/2 + 3\\beta/2 &= -1/4 \n",
"\\end{align*}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"which has the solution $\\alpha=1/4$, $\\beta=-1/4$. $\\hat\\psi=\\psi-\\phi_{0,0}/4+\\beta\\phi_{0,1}/4$ can be plotted as follows."
]
},
{
"cell_type": "code",
"execution_count": 43,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"t=linspace(0,4,100)\n",
"phi=lambda t: ((0<=t)*(t<1)).astype(float)\n",
"psi=lambda t: ((0<=t)*(t<1/2)).astype(float) - ((1/2<=t)*(t<1)).astype(float)\n",
"plt.plot(t,psi(t)-phi(t)/4+phi(t-1)/4)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"**c)**\n",
"Express $\\phi$ and $\\hat{\\psi}$ with the help of functions from $\\boldsymbol{\\phi}_1$, and use this to write down the change of coordinate matrix from $(\\boldsymbol{\\phi}_0,\\hat{\\boldsymbol{\\psi}}_0)$\n",
"to $\\boldsymbol{\\phi}_1$.\n",
"\n",
"\n",
"\n",
"**Solution.**\n",
"We have that"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\begin{align*}\n",
"\\hat\\psi &= \\psi-\\phi_{0,0}/4+\\beta\\phi_{0,1}/4 \\\\ \n",
"&= \\frac{1}{\\sqrt{2}}(\\phi_{1,0}-\\phi_{1,1}) - \\frac{1}{4\\sqrt{2}}(\\phi_{1,0}+\\phi_{1,1}) + \\frac{1}{4\\sqrt{2}}(\\phi_{1,2}+\\phi_{1,3}) \\\\ \n",
"&= \\frac{3}{4\\sqrt{2}}\\phi_{1,0} -\\frac{5}{4\\sqrt{2}}\\phi_{1,1} + \\frac{1}{4\\sqrt{2}}\\phi_{1,2}+\\frac{1}{4\\sqrt{2}}\\phi_{1,3}\n",
"\\end{align*}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This means that the last half of the columns in the change of coordinates matrix are the translates of $\\frac{1}{4\\sqrt{2}}(3,-5,1,1)$. The remaining columns coincide with those of the matrix for the Haar wavelet.\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"## Exercise 4.33: More than two vanishing moments for the Haar wavelet\n",
"\n",
"\n",
"\n",
"It is also possible to add more vanishing moments to the Haar wavelet. Define"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\hat\\psi=\\psi_{0,0} - a_0 \\phi_{0,0} - \\cdots - a_{k-1} \\phi_{0,k-1}.\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Define also $c_{r,l}=\\int_l^{l+1} t^rdt$, and $e_{r}=\\int_0^1 t^r\\psi(t)dt$.\n",
"\n",
"\n",
"**a)**\n",
"Show that $\\hat\\psi$ has $k$ vanishing moments if and only if $a_0,\\ldots,a_{k-1}$ solves the system"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\begin{pmatrix}\n",
"c_{0,0} & c_{0,1} & \\cdots & c_{0,k-1} \\\\ \n",
"c_{1,0} & c_{1,1} & \\cdots & c_{1,k-1} \\\\ \n",
"\\vdots & \\vdots & \\vdots & \\vdots \\\\ \n",
"c_{k-1,0} & c_{k-1,1} & \\cdots & c_{k-1,k-1}\n",
"\\end{pmatrix}\n",
"\\begin{pmatrix}\n",
"a_0 \\\\ a_1 \\\\ \\vdots \\\\ a_{k-1}\n",
"\\end{pmatrix}\n",
"=\n",
"\\begin{pmatrix}\n",
"e_0 \\\\ e_1 \\\\ \\vdots \\\\ e_{k-1}\n",
"\\end{pmatrix}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"**Solution.**\n",
"Multiplying with $t^r$ and computing the integral we obtain"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\langle\\hat\\psi,t^r\\rangle = \\langle\\psi_{0,0},t^r\\rangle - a_0 \\langle\\phi_{0,0},t^r\\rangle - \\cdots - a_{k-1} \\langle\\phi_{0,k-1},t^r\\rangle.\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Setting those to zero we obtain"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"0 = e_r - a_0c_{r,0} - \\cdots - a_{k-1}c_{r,k-1},\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"which is row $r$ in the stated equation system.\n",
"\n",
"\n",
"\n",
"**b)**\n",
"Write a function `vanishing_moms_haar` which takes $k$ as input, solves the system in a), and returns the vector $\\boldsymbol{a}=(a_0,a_1,\\ldots,a_{k-1})$.\n",
"\n",
"\n",
"\n",
"**Solution.**\n",
"We first find that $c_{r,l}=\\frac{1}{r+1}((l+1)^{r+1}-l^{r+1})$, and"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"e_r = [t^{r+1}/(r+1)]_0^{1/2} - [t^{r+1}/(r+1)]_{1/2}^{1} = \\frac{1}{r+1}(2^{-r}-1).\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The following code can be used. If you set $k=2$ the values from the previous exercise are obtained."
]
},
{
"cell_type": "code",
"execution_count": 44,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"def vanishing_moms_haar(k):\n",
" l = arange(0,k).reshape((1,k))\n",
" r=transpose(l)\n",
" C = ((l+1)**(r+1)-l**(r+1))/(r+1)\n",
" e = (2**(-r)-1)/(r+1)\n",
" return linalg.solve(C, e)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"## Exercise 4.34: Listening experiments\n",
"\n",
"\n",
"\n",
"Run the function `forw_comp_rev_dwt1` for different $m$ for the Haar wavelet, \n",
"the piecewise linear wavelet, and the alternative piecewise linear wavelet, but listen to the detail components instead. \n",
"Describe the sounds you hear for different $m$, and try to explain why the sound seems to get louder when you increase $m$.\n",
"\n",
"\n",
"\n",
"**Solution.**\n",
"The following code can be used:"
]
},
{
"cell_type": "code",
"execution_count": 45,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"x, fs = forw_comp_rev_dwt1(m, 'Haar')\n",
"play(x, fs)"
]
},
{
"cell_type": "code",
"execution_count": 46,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"x, fs = forw_comp_rev_dwt1(m, 'pwl0')\n",
"play(x, fs)"
]
},
{
"cell_type": "code",
"execution_count": 47,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"x, fs = forw_comp_rev_dwt1(m, 'pwl2')\n",
"play(x, fs)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"## Exercise 5.8: Computing filters and frequency responses\n",
"\n",
"\n",
"\n",
"\n",
"**a)**\n",
"Write down the filter coefficients for the corresponding filter $G_1$ obtained in Exercise 4.31, and plot its frequency response.\n",
"\n",
"\n",
"\n",
"**Solution.**\n",
"The code can look as follows:"
]
},
{
"cell_type": "code",
"execution_count": 48,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"omega = linspace(0,2*pi,100)\n",
"plt.plot(omega, g1[4] + g1[5]*2*cos(omega) + g1[6]*2*cos(2*omega) \\\n",
" + g1[7]*2*cos(3*omega) + g1[8]*2*cos(4*omega))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"**b)**\n",
"Write down the corresponding filters $G_0$ and $G_1$ for Exercise 4.32. \n",
"Plot their frequency responses, and characterize the filters as low-pass- or high-pass.\n",
"\n",
"\n",
"\n",
"**Solution.**\n",
"The filter $G_0$ is the same as for the Haar wavelet, and this frequency response has already been plotted. In Exercise 4.32 we obtained $G_1=\\frac{1}{4\\sqrt{2}}(3,-5,1,1)$, and the magnitude of this frequency response can be plotted as follows."
]
},
{
"cell_type": "code",
"execution_count": 49,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"omega=linspace(0,2*pi,100)\n",
"x= concatenate([ array([3, -5, 1, 1])/(4*sqrt(2)), zeros(96)])\n",
"plt.plot(omega, abs(fft.fft(x)))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"**c)**\n",
"Repeat b) for the Haar wavelet as in Exercise 4.33, and plot the corresponding frequency responses for $k=2,4,6$.\n",
"\n",
"\n",
"\n",
"**Solution.**\n",
"Using the function we implemented in Exercise 4.33, the following code can be used."
]
},
{
"cell_type": "code",
"execution_count": 50,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"for k in [2, 4, 6]:\n",
" plt.figure()\n",
" x = vanishing_moms_haar(k)\n",
" x = x.reshape(len(x))\n",
" x = concatenate([-x, zeros(50-len(x)), [1], zeros(49)])\n",
" idwt_impl(x, 'haar', m=1, bd_mode='per')\n",
" plt.plot(omega, abs(fft.fft(x)))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"## Exercise 5.19: Using the cascade algorithm\n",
"\n",
"\n",
"\n",
"In Exercise 4.31 we constructed a new mother wavelet $\\hat{\\psi}$ for piecewise linear functions by finding constants\n",
"$\\alpha,\\beta,\\gamma,\\delta$ so that"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\hat\\psi=\\psi-\\alpha\\phi_{0,0}-\\beta\\phi_{0,1}-\\delta\\phi_{0,2}-\\gamma\\phi_{0,N-1}.\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Use the cascade algorithm to plot $\\hat\\psi$. Do this by using the wavelet kernel for the piecewise linear wavelet \n",
"(do not use the code above, since we have not implemented kernels for this wavelet yet).\n",
"\n",
"\n",
"\n",
"**Solution.**\n",
"Assuming that the vector `coeffs` has been set as in Exercise 4.31, the code can look as follows"
]
},
{
"cell_type": "code",
"execution_count": 51,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"m = 10\n",
"t = linspace(-2, 6, 8*2**m)\n",
"coordsvm = hstack([[-coeffs[0], -coeffs[1], -coeffs[3], 0, 0, 0, 0,\\\n",
" -coeffs[2], 1, 0, 0, 0, 0, 0, 0, 0], \\\n",
" zeros(8*2**m-16)])\n",
"idwt_impl(coordsvm, 'pwl0', m=m, bd_mode='per')\n",
"coordsvm *= 2**(m/2.)\n",
"plt.plot(t,hstack([coordsvm[(6*2**m):(8*2**m+1)], coordsvm[0:(6*2**m)]]))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
""
]
}
],
"metadata": {},
"nbformat": 4,
"nbformat_minor": 2
}