{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- dom:TITLE: Exercises from Linear algebra, signal processing, and wavelets. A unified approach.\\\\ Python version -->\n",
    "# Exercises from Linear algebra, signal processing, and wavelets. A unified approach.\\\\ Python version\n",
    "<!-- dom:AUTHOR: Øyvind Ryan -->\n",
    "<!-- Author: --> **Øyvind Ryan**\n",
    "\n",
    "Date: **Oct 24, 2018**\n",
    "\n",
    "<!-- Externaldocuments: applinalg -->\n",
    "<!-- Mapping from exercise labels to numbers: label2numbers = {'example:haarfreqresp': '5.1', 'example:freqresp3first': '5.18', 'example:cascade_pwc_pwl': '5.17', 'example:filtersaltplw': '5.3', 'exercise:plot_freqresp_mp3_wavfilters': '5.29'} -->\n",
    "\n",
    "# The filter representation of wavelets\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "<!-- --- begin exercise --- -->\n",
    "\n",
    "## Example 5.1: The Haar wavelet\n",
    "<div id=\"example:haarfreqresp\"></div>\n",
    "<!-- keywords = wavfilters -->"
   ]
  },
  {
   "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 dwt import *\n",
    "from sound import *\n",
    "from forward_compress_reverse import *\n",
    "from scipy import integrate\n",
    "from mp3funcs import *"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "N = 128\n",
    "n = arange(0, N)\n",
    "omega = 2*pi*n/float(N)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For the Haar wavelet we saw that, in $G$, the matrix"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\begin{pmatrix} \\frac{1}{\\sqrt{2}} & \\frac{1}{\\sqrt{2}} \\\\ \\frac{1}{\\sqrt{2}} & -\\frac{1}{\\sqrt{2}} \\end{pmatrix}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "repeated along the diagonal. The filters $G_0$ and $G_1$ can be found directly from this as"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\begin{align*}\n",
    "G_0 &=\\{\\underline{1/\\sqrt{2}},1/\\sqrt{2}\\} \\\\ \n",
    "G_1 &=\\{1/\\sqrt{2},\\underline{-1/\\sqrt{2}}\\}.\n",
    "\\end{align*}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We have seen these filters previously: $G_0$ is a moving average filter of two elements (up to multiplication with a constant), i.e. a low-pass filter, \n",
    "and $G_1$ is the corresponding high-pass filter (up to a delay), which also can be seen as an approximation to the derivative. \n",
    "The frequency responses are"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\begin{align*}\n",
    "\\lambda_{G_0}(\\omega) &= \\frac{1}{\\sqrt{2}}+\\frac{1}{\\sqrt{2}}e^{-i\\omega}=\\sqrt{2}e^{-i\\omega/2}\\cos(\\omega/2) \\\\ \n",
    "\\lambda_{G_1}(\\omega) &= \\frac{1}{\\sqrt{2}}e^{i\\omega}-\\frac{1}{\\sqrt{2}}=\\sqrt{2}ie^{i\\omega/2}\\sin(\\omega/2).\n",
    "\\end{align*}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "By considering filters which repeat with the same rows, it is clear that"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\begin{align*}\n",
    "H_0 &= \\{1/\\sqrt{2},\\underline{1/\\sqrt{2}}\\} \\\\ \n",
    "H_1 &=\\{\\underline{-1/\\sqrt{2}},1/\\sqrt{2}\\},\n",
    "\\end{align*}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "and these frequency responses have the same low-pass/high-pass characteristics. The frequency responses can be plotted as follows."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "plt.plot(omega, abs(1+exp(-1j*omega))/sqrt(2))\n",
    "plt.figure()\n",
    "plt.plot(omega, abs(exp(1j*omega)-1)/sqrt(2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- --- end exercise --- -->\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "<!-- --- begin exercise --- -->\n",
    "\n",
    "## Example 5.3: The alternative piecewise linear wavelet\n",
    "<div id=\"example:filtersaltplw\"></div>\n",
    "<!-- keywords = wavfilters -->\n",
    "\n",
    "Previously we found the first two columns in $P_{\\boldsymbol{\\phi}_m\\leftarrow{\\cal C}_m}$ for the alternative piecewise linear wavelet. This\n",
    "gives us that the filters $G_0$ and $G_1$ are"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\begin{align*}\n",
    "G_0 &= \\frac{1}{\\sqrt{2}}\\{1/2,\\underline{1},1/2\\} \\\\ \n",
    "G_1 &= \\frac{1}{\\sqrt{2}}\\{-1/8,-1/4,\\underline{3/4},-1/4,-1/8\\}.\n",
    "\\end{align*}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here $G_0$ is as for the wavelet of piecewise linear functions since we use the same scaling function. $G_1$ is different, however. \n",
    "It tums out that $G_1$ now has high-pass characteristics.\n",
    "\n",
    "We can perform the matrix multiplication in equation (4.37) to compute the DWT filter components:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\begin{align*}\n",
    "H_0 &= \\sqrt{2}\\{ -1/8,1/4,\\underline{3/4},1/4,-1/8\\} \\\\ \n",
    "H_1 &= \\sqrt{2}\\{ -1/2,\\underline{1},-1/2\\}.\n",
    "\\end{align*}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The frequency responses can be plotted as follows"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "plt.plot(omega, abs(1+cos(omega))/sqrt(2))\n",
    "plt.figure()\n",
    "plt.plot(omega, abs(3/4. - cos(omega)/2. - cos(2*omega)/4.)/sqrt(2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The filters $G_0,G_1,H_0,H_1$ are particularly important in applications: Apart from the scaling factors $1/\\sqrt{2}$, $\\sqrt{2}$ in front, we see that the filter\n",
    "coefficients are all dyadic fractions, i.e. they are on the form $\\beta/2^j$. Arithmetic operations with dyadic fractions can be carried out exactly on a\n",
    "computer, due to representations as binary numbers. These filters are thus important in applications, since they can be used for\n",
    "lossless coding. The same argument can be made for the Haar wavelet, but this wavelet had one less vanishing moment. \n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "Note that the role of $H_1$ as the high-pass filter corresponding to $G_0$ is the case in both previous examples. We will prove in the next chapter that this is a much more general result.\n",
    "\n",
    "<!-- --- end exercise --- -->\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "<!-- --- begin exercise --- -->\n",
    "\n",
    "## Example 5.17: Plotting scaling functions and mother wavelets\n",
    "<div id=\"example:cascade_pwc_pwl\"></div>\n",
    "<!-- keywords = wavfilters -->\n",
    "\n",
    "It turns out that the interval $[-2,6]$ contains the supports of all scaling functions and mother wavelets we have considered up to now. \n",
    "We can therefore use $a=-2$ and $b=6$ as parameters to `cascade_alg`, in order to plot all these functions.\n",
    "We will also use $m=10$ levels. \n",
    "The following code then runs the cascade algorithm for the three wavelets we have considered, to reproduce all previous scaling functions and mother wavelets."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "cascade_alg(10, -2, 6, 'Haar', True, False)\n",
    "cascade_alg(10, -2, 6, 'Haar', False, False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "cascade_alg(10, -2, 6, 'pwl0', True, False)\n",
    "cascade_alg(10, -2, 6, 'pwl0', False, False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "cascade_alg(10, -2, 6, 'pwl2', True, False)\n",
    "cascade_alg(10, -2, 6, 'pwl2', False, False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the above code you are encouraged to also try the other wavelets we encounter later in these notes. You may then also have to change the support interval $[a,b]$ accordingly.\n",
    "\n",
    "<!-- --- end exercise --- -->\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "<!-- --- begin exercise --- -->\n",
    "\n",
    "## Example 5.18: Plotting frequency responses\n",
    "<div id=\"example:freqresp3first\"></div>\n",
    "<!-- keywords = wavfilters -->\n",
    "\n",
    "In order to verify the low-pass/high-pass characteristics of $G_0$ and $G_1$, \n",
    "let us plot the frequency responses of the wavelets we have considered. \n",
    "To plot $\\lambda_{G_0}(\\omega)$ and $\\lambda_{G_1}(\\omega)$ for the Haar wavelet, we can write"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "freqresp_alg('Haar', True, False)\n",
    "freqresp_alg('Haar', False, False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To plot the same frequency response for the alternative piecewise linear wavelet, we can write"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "freqresp_alg('pwl2', True, False)\n",
    "freqresp_alg('pwl2', False, False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- --- end exercise --- -->\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "<!-- --- begin exercise --- -->\n",
    "\n",
    "## Exercise 5.29: Verifying a connection between the tables in the MP3 standard\n",
    "<div id=\"exercise:plot_freqresp_mp3_wavfilters\"></div>\n",
    "<!-- keywords = wavfilters -->\n",
    "\n",
    "In addition to the function `mp3_ctable()`, the library has a function `mp3_dtable()` which returns the table $D$ used by the reverse transform. \n",
    "Verify the connection we have stated between the tables $C$ and $D$, i.e. that $D_i=32C_i$ for all $i$.\n",
    "\n",
    "\n",
    "<!-- --- begin solution of exercise --- -->\n",
    "**Solution.**\n",
    "This can be verified as follows."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "C = mp3_ctable()\n",
    "D = mp3_dtable()\n",
    "print(abs(D-32*C).max())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- --- end solution of exercise --- -->\n",
    "\n",
    "<!-- --- end exercise --- -->"
   ]
  }
 ],
 "metadata": {},
 "nbformat": 4,
 "nbformat_minor": 2
}