{
 "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:secondstrategy': '6.11', 'exercise:view_freq_resp_spline': '6.9', 'exercise:oblig2wavelet': '6.8'} -->\n",
    "\n",
    "# Constructing interesting wavelets\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "<!-- --- begin exercise --- -->\n",
    "\n",
    "## Exercise 6.8: Computing filters\n",
    "<div id=\"exercise:oblig2wavelet\"></div>\n",
    "<!-- keywords = wavinteresting; student -->"
   ]
  },
  {
   "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": [
    "from sympy import * \n",
    "from forward_compress_reverse import *\n",
    "from sound import *\n",
    "from dwt import *"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compute the filters $H_0$, $G_0$ in Theorem 8 \n",
    "when $N=N_1=N_2=4$, and $Q_1=Q^{(4)}$, $Q_2=1$. \n",
    "Compute also filters $H_1$, $G_1$ so that we have perfect\n",
    "reconstruction (note that these are not unique).\n",
    "\n",
    "\n",
    "<!-- --- begin solution of exercise --- -->\n",
    "**Solution.**\n",
    "We have that"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\begin{align*}\n",
    "  \\lambda_{H_0}(\\omega) &= \\left(\\frac{1}{2}(1+\\cos\\omega)\\right)^{N_1/2} Q_1\\left(\\frac{1}{2}(1-\\cos\\omega)\\right) \\\\ \n",
    "  &= \\left(\\frac{1}{2}(1+\\cos\\omega)\\right)^2 Q^{(4)}\\left(\\frac{1}{2}(1-\\cos\\omega)\\right) \\\\ \n",
    "  \\lambda_{G_0}(\\omega) &= \\left(\\frac{1}{2}(1+\\cos\\omega)\\right)^{N_2/2} Q_2\\left(\\frac{1}{2}(1-\\cos\\omega)\\right) = \\left(\\frac{1}{2}(1+\\cos\\omega)\\right)^2 \\\\ \n",
    "                        &= \\frac{1}{4}\\left(1+\\frac{1}{2}e^{i\\omega}+\\frac{1}{2}e^{-i\\omega}\\right)^2 = \\frac{1}{16}(e^{2i\\omega}+4e^{i\\omega}+6+4e^{-i\\omega}+e^{-2i\\omega}).\n",
    "\\end{align*}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Therefore $G_0=\\frac{1}{16}\\{1,4,\\underline{6},4,1\\}$. We do not recommend to compute $H_0$ by hand. \n",
    "With \n",
    "the package `sympy` in Python\n",
    "you can do as follows to compute $H_0$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "x = Symbol('x')\n",
    "z = expand( ((1+x/2+1/(2*x))/2)**2* \\\n",
    "        (2+8*((1-x/2-1/(2*x))/2)+20*((1-x/2-1/(2*x))/2)**2\\\n",
    "        +40*((1-x/2-1/(2*x))/2)**3) )\n",
    "print(z)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we have substituted `x` for $e^{i\\omega}$, `1/x` for $e^{-i\\omega}$. The first part represents $\\left(\\frac{1}{2}(1+\\cos\\omega)\\right)^2$, the second part represents\n",
    "$Q^{(4)}(u)=2+8u+20u^2+40u^3$ with $u=\\frac{1}{2}(1-\\cos\\omega)=\\frac{1}{2}\\left(1-\\frac{1}{2}e^{i\\omega}-\\frac{1}{2}e^{-i\\omega}\\right)$. This gives"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "H_0=\\frac{1}{128}\\{-5,20,-1,-96,70,\\underline{280},70,-96,-1,20,-5\\}.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using Exercise 5.10 with $\\alpha=1$, $d=0$, we get"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\begin{align*}\n",
    "H_1 &= \\frac{1}{16}\\{1,-4,\\underline{6},-4,1\\} \\\\ \n",
    "G_1 &= \\frac{1}{128}\\{5,20,1,-96,-70,\\underline{280},-70,-96,1,20,5\\}\n",
    "\\end{align*}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- --- end solution of exercise --- -->\n",
    "\n",
    "<!-- --- end exercise --- -->\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "<!-- --- begin exercise --- -->\n",
    "\n",
    "## Exercise 6.9: Plotting frequency responses of Spline wavelets\n",
    "<div id=\"exercise:view_freq_resp_spline\"></div>\n",
    "<!-- keywords = wavinteresting; student -->\n",
    "\n",
    "The library represents a Spline wavelet with `x` and `y` vanishing moments (i.e. $N_1$ and $N_2$) with the name `\"splinex.y\"`.\n",
    "\n",
    "\n",
    "**a)**\n",
    "Listen to the low-resolution approximations and detail components in sound for the `\"spline4.4\"` wavelet.\n",
    "\n",
    "\n",
    "<!-- --- begin solution of exercise --- -->\n",
    "**Solution.**\n",
    "The following code can be used for listening to the low-resolution approximations for a given value of $m$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "m = 1\n",
    "x, fs = forw_comp_rev_dwt1(m, 'spline4.4')\n",
    "play(x, fs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "x, fs = forw_comp_rev_dwt1(m, 'spline4.4', False)\n",
    "play(x, fs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- --- end solution of exercise --- -->\n",
    "\n",
    "**b)**\n",
    "Plot all scaling functions and mother wavelets (using the cascade algorithm), and frequency responses for the `\"spline4.4\"` wavelet.\n",
    "\n",
    "\n",
    "<!-- --- begin solution of exercise --- -->\n",
    "**Solution.**\n",
    "The code to plot the functions can look as follows."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "m=10;\n",
    "cascade_alg(m, -4, 4, 'spline4.4', True, False)\n",
    "cascade_alg(m, -4, 4, 'spline4.4', False, False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The frequency response can be plotted as follows"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "        freqresp_alg('spline4.4', lowpass, dual)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- --- end solution of exercise --- -->\n",
    "\n",
    "<!-- --- end exercise --- -->\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "<!-- --- begin exercise --- -->\n",
    "\n",
    "## Example 6.11: The CDF 9/7 wavelet\n",
    "<div id=\"example:secondstrategy\"></div>\n",
    "<!-- keywords = wavinteresting -->\n",
    "\n",
    "We now choose $N_1=N_2=4$. In equation (6.20) we pair inverse terms to obtain"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\begin{align*}\n",
    "\\lefteqn{Q^{(3)}\\left(\\frac{1}{2}(1-\\cos\\omega)\\right)} \\\\ \n",
    "    &= \\frac{5}{8}\\frac{1}{3.0407}\\frac{1}{7.1495}(e^{i\\omega}-3.0407)(e^{-i\\omega}-3.0407) \\\\ \n",
    "    & \\times (e^{2i\\omega}-4.0623e^{i\\omega}+7.1495)(e^{-2i\\omega}-4.0623e^{-i\\omega}+7.1495) \\\\ \n",
    "    &= \\frac{5}{8}\\frac{1}{3.0407}\\frac{1}{7.1495}(-3.0407e^{i\\omega} + 10.2456 - 3.0407e^{-i\\omega}) \\\\ \n",
    "    & \\times (7.1495e^{2i\\omega} - 33.1053e^{i\\omega} + 68.6168 - 33.1053e^{-i\\omega} + 7.1495e^{-2i\\omega}).\n",
    "\\end{align*}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can write this as $Q_1Q_2$ with $Q_1(0)=Q_2(0)$ when"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\begin{align*}\n",
    "Q_1(\\omega) &= -1.0326e^{i\\omega} + 3.4795 - 1.0326e^{-i\\omega} \\\\ \n",
    "Q_2(\\omega) &= 0.6053e^{2i\\omega}   -2.8026e^{i\\omega}    + 5.8089   -2.8026e^{-i\\omega}    + 0.6053e^{-2i\\omega},\n",
    "\\end{align*}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "from which we obtain"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\begin{align*}\n",
    "  \\lambda_{G_0}(\\omega) &= \\left(\\frac{1}{2}(1+\\cos\\omega)\\right)^2    Q_1(\\omega) \\\\ \n",
    "         &= -0.0645e^{3i\\omega}   -0.0407e^{2i\\omega}    +0.4181e^{i\\omega}    + 0.7885    \\\\ \n",
    "         &  +0.4181e^{-i\\omega}   -0.0407e^{-2i\\omega}   -0.0645e^{-3i\\omega} \\\\ \n",
    "  \\lambda_{H_0}(\\omega) &= \\left(\\frac{1}{2}(1+\\cos\\omega)\\right)^2 40 Q_2(\\omega) \\\\ \n",
    "         &= 0.0378e^{4i\\omega}   -0.0238e^{3i\\omega}   -0.1106e^{2i\\omega}    +0.3774e^{i\\omega}    +0.8527    \\\\ \n",
    "         & +0.3774e^{-i\\omega}   -0.1106e^{-2i\\omega}  -0.0238e^{-3i\\omega} + 0.0378e^{-4i\\omega}.\n",
    "\\end{align*}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The filters $G_0$, $H_0$ are thus"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\begin{align*}\n",
    "G_0 &= \\{0.0645,0.0407,-0.4181,\\underline{-0.7885},-0.4181,0.0407,0.0645\\} \\\\ \n",
    "H_0 &= \\{-0.0378,0.0238,0.1106,-0.3774,\\underline{-0.8527},-0.3774,0.1106,0.0238,-0.0378\\}.\n",
    "\\end{align*}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The corresponding frequency responses are plotted \n",
    "with the following code. \n",
    "It is seen that both filters are low-pass also here, and that they are closer to an ideal band-pass filter. The frequency responses seem very flat near $\\pi$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "freqresp_alg('cdf97', True, False)\n",
    "freqresp_alg('cdf97', True, True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We also get"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\begin{align*}\n",
    "G_1 &= \\{-0.0378,-0.0238,0.1106,0.3774,\\underline{-0.8527},0.3774,0.1106,-0.0238,-0.0378\\} \\\\ \n",
    "H_1 &= \\{-0.0645,0.0407,0.4181,\\underline{-0.7885},0.4181,0.0407,-0.0645\\}.\n",
    "\\end{align*}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The length of the filters are $9$ and $7$ in this case, which is why this wavelet is called the *CDF 9/7 wavelet*. \n",
    "\n",
    "\n",
    "In Example 5.3 we saw that we had analytical expressions for the scaling functions and the mother wavelet, but that we could not obtain this for\n",
    "the dual functions. For the CDF 9/7 wavelet it turns out that none of the four functions have analytical expressions. Let us therefore use the cascade algorithm \n",
    "to plot them. Note first that since $G_0$ has $7$ filter coefficients, and $G_1$ has $9$ filter coefficients, it follows from\n",
    "Exercise 5.16 that $\\text{Supp}(\\phi)=[-3,3]$, $\\text{Supp}(\\psi)=[-3,4]$, $\\text{Supp}(\\tilde\\phi)=[-4,4]$, and $\\text{Supp}(\\tilde\\psi)=[-3,4]$. \n",
    "The scaling functions and mother wavelets over these supports are \n",
    "plotted with the following code. \n",
    "Again they have irregular shapes, but now at least the functions and dual functions resemble each other more. The functions would be equal to their dual counterparts if the wavelet was orthogonal, so this suggests that this wavelet \n",
    "is close to being orthogonal. This may also explain partially why the CDF 9/7 wavelet works well for compression purposes: It is in fact used for lossy\n",
    "compression with the JPEG2000 standard."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "cascade_alg(10, -4, 4, 'cdf97', True, False)\n",
    "cascade_alg(10, -4, 4, 'cdf97', False, False)\n",
    "cascade_alg(10, -4, 4, 'cdf97', True, True)\n",
    "cascade_alg(10, -4, 4, 'cdf97', False, True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the above example there was a unique way of factoring $Q$ into a product of real polynomials. For higher degree polynomials there is no unique way to\n",
    "distribute the factors, and we will not go into what strategy can be used for this. In general, we must go through the following steps:\n",
    "\n",
    "* Compute the polynomial $Q$, and find its roots.\n",
    "\n",
    "* Pair complex conjugate roots into real second degree polynomials, and form polynomials $Q_1$, $Q_2$.\n",
    "\n",
    "* Compute the coefficients in equations (6.16) and (6.17).\n",
    "\n",
    "<!-- --- end exercise --- -->"
   ]
  }
 ],
 "metadata": {},
 "nbformat": 4,
 "nbformat_minor": 2
}