{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- dom:TITLE: Exercises from Linear algebra, signal processing, and wavelets. A unified approach.\\\\ MATLAB version -->\n",
    "# Exercises from Linear algebra, signal processing, and wavelets. A unified approach.\\\\ MATLAB 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 = {'exercise:verify_statement_filters': '7.14', 'exercise:5397many': '7.8', 'exercise:ex5kernel': '7.9', 'exercise:verify_near_perfect': '7.15', 'exercise:080301': '7.13'} -->\n",
    "\n",
    "# The polyphase representation of filter bank transforms\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "<!-- --- begin exercise --- -->\n",
    "\n",
    "## Exercise 7.8: Lifting-based implementations of the Spline 5/3 and CDF 9/7 wavelets\n",
    "<div id=\"exercise:5397many\"></div>\n",
    "<!-- keywords = wavlifting; student -->\n",
    "\n",
    "Let us use the different lifting factorizations obtained in this chapter to implement the corresponding wavelet kernels.\n",
    "Your functions should call the functions from Exercise 4.24 and 4.29. \n",
    "You will need equations (7.8)-(7.11) here, in order to complete the kernels."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "% Initialization code"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**a)**\n",
    "Write functions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "        dwt_kernel_53(x, bd_mode)\n",
    "        idwt_kernel_53(x, bd_mode)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "which implement the DWT and IDWT kernels for the Spline 5/3 wavelet.\n",
    "Use the lifting factorization obtained in Example 7.2.3.\n",
    "\n",
    "\n",
    "<!-- --- begin solution of exercise --- -->\n",
    "**Solution.**\n",
    "The code can look like this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%%file dwt_kernel_53.m\n",
    "function x=dwt_kernel_53(x, bd_mode)\n",
    "    lambda1 = -1;\n",
    "    lambda2 = 0.125;\n",
    "    alpha = 2;\n",
    "    beta = 0.5;\n",
    "  \n",
    "    x(1:2:end, :) = x(1:2:end, :)*alpha;\n",
    "    x(2:2:end, :) = x(2:2:end, :)*beta;\n",
    "    x = lifting_odd_symm(-lambda2, x, bd_mode);\n",
    "    x = lifting_even_symm(-lambda1, x, bd_mode);\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%%file idwt_kernel_53.m\n",
    "function x=idwt_kernel_53(x, bd_mode)\n",
    "    lambda1 = -1;\n",
    "    lambda2 = 0.125;\n",
    "    alpha = 2;\n",
    "    beta = 0.5;\n",
    "  \n",
    "    x = lifting_even_symm(lambda1, x, bd_mode);\n",
    "    x = lifting_odd_symm(lambda2, x, bd_mode);\n",
    "    x(1:2:end, :) = x(1:2:end, :)/alpha;\n",
    "    x(2:2:end, :) = x(2:2:end, :)/beta;\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- --- end solution of exercise --- -->\n",
    "\n",
    "\n",
    "**b)**\n",
    "Write functions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "        dwt_kernel_97(x, bd_mode)\n",
    "        idwt_kernel_97(x, bd_mode)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "which implement the DWT and IDWT kernels for the CDF 9/7 wavelet. \n",
    "Use the lifting factorization obtained in Example 7.2.4.\n",
    "\n",
    "\n",
    "<!-- --- begin solution of exercise --- -->\n",
    "**Solution.**\n",
    "The code can look like this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%%file dwt_kernel_97.m\n",
    "function x=dwt_kernel_97(x, bd_mode)\n",
    "    lambda1 = -0.586134342059950;\n",
    "    lambda2 = -0.668067171029734;\n",
    "    lambda3 = 0.070018009414994;\n",
    "    lambda4 = 1.200171016244178;\n",
    "    alpha = -1.149604398860250;\n",
    "    beta = -0.869864451624777;\n",
    "  \n",
    "    x(1:2:end, :) = x(1:2:end, :)*alpha;\n",
    "    x(2:2:end, :) = x(2:2:end, :)*beta;\n",
    "    x = lifting_odd_symm(-lambda4, x, bd_mode);\n",
    "    x = lifting_even_symm(-lambda3, x, bd_mode);\n",
    "    x = lifting_odd_symm(-lambda2, x, bd_mode);\n",
    "    x = lifting_even_symm(-lambda1, x, bd_mode);\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%%file idwt_kernel_97.m\n",
    "function x=idwt_kernel_97(x, bd_mode)\n",
    "    lambda1 = -0.586134342059950;\n",
    "    lambda2 = -0.668067171029734;\n",
    "    lambda3 = 0.070018009414994;\n",
    "    lambda4 = 1.200171016244178;\n",
    "    alpha = -1.149604398860250;\n",
    "    beta = -0.869864451624777;\n",
    "  \n",
    "    x = lifting_even_symm(lambda1, x, bd_mode);\n",
    "    x = lifting_odd_symm(lambda2, x, bd_mode);\n",
    "    x = lifting_even_symm(lambda3, x, bd_mode);\n",
    "    x = lifting_odd_symm(lambda4, x, bd_mode);\n",
    "    x(1:2:end, :) = x(1:2:end, :)/alpha;\n",
    "    x(2:2:end, :) = x(2:2:end, :)/beta;\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- --- end solution of exercise --- -->\n",
    "\n",
    "**c)**\n",
    "In Chapter 4, \n",
    "we listened to the low-resolution approximations and detail components in sound for three different wavelets. \n",
    "Repeat these experiments with the Spline 5/3 and the CDF 9/7 wavelet, using the new kernels you implemented\n",
    "in a) and b).\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": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "m = 1;\n",
    "[x, fs] = forw_comp_rev_dwt1(m, 'spline53');\n",
    "playerobj = audioplayer(x, fs);\n",
    "playblocking(playerobj)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "m = 1;\n",
    "[x, fs] = forw_comp_rev_dwt1(m, 'cdf97');\n",
    "playerobj = audioplayer(x, fs);\n",
    "playblocking(playerobj)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- --- end solution of exercise --- -->\n",
    "\n",
    "**d)**\n",
    "Plot all scaling functions and mother wavelets for the Spline 5/3\n",
    "and the CDF 9/7 wavelets, using the cascade algorithm and the kernels you have implemented.\n",
    "\n",
    "\n",
    "<!-- --- begin solution of exercise --- -->\n",
    "**Solution.**\n",
    "The code can look as follows."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "m=10;\n",
    "cascade_alg(m, -4, 4, 'spline53', 1, 0)\n",
    "cascade_alg(m, -4, 4, 'spline53', 0, 0)\n",
    "\n",
    "cascade_alg(m, -4, 4, 'cdf97', 1, 0)\n",
    "cascade_alg(m, -4, 4, 'cdf97', 0, 0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the plot for the CDF 9/7 wavelet, it is seen that the functions and their dual counterparts are close to being equal. \n",
    "This reflects the fact that this wavelet is close to being orthogonal.\n",
    "\n",
    "<!-- --- end solution of exercise --- -->\n",
    "\n",
    "<!-- --- end exercise --- -->\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "<!-- --- begin exercise --- -->\n",
    "\n",
    "## Exercise 7.9: Lifting-based implementation of orthonormal wavelets\n",
    "<div id=\"exercise:ex5kernel\"></div>\n",
    "<!-- keywords = wavlifting -->\n",
    "\n",
    "\n",
    "\n",
    "**a)**\n",
    "Write functions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "        lifting_even(lambda1, lambda2, x, bd_mode)\n",
    "        lifting_odd(lambda1, lambda2, x, bd_mode)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "which apply the elementary\n",
    "lifting matrices (7.14) to $\\boldsymbol{x}$. Assume that $N$ is even.\n",
    "\n",
    "\n",
    "<!-- --- begin solution of exercise --- -->\n",
    "**Solution.**\n",
    "The library versions of these functions are as follows."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "        function x=lifting_even(lambda1, lambda2, x, bd_mode)\n",
    "            N = size(x, 1);\n",
    "            if strcmpi(bd_mode, 'per')\n",
    "                x(1, :) = lambda1*x(2, :) + x(1, :) + lambda2*x(N, :);\n",
    "            elseif strcmpi(bd_mode, 'bd') || strcmpi(bd_mode, 'none')\n",
    "                x(1, :) = lambda1*x(2, :) + x(1, :);\n",
    "            end\n",
    "            x(3:2:(N-1), :) = ...\n",
    "                lambda1*x(4:2:N, :) + x(3:2:(N-1), :) + lambda2*x(2:2:(N-2), :);\n",
    "            if strcmpi(bd_mode, 'bd') || strcmpi(bd_mode, 'none')\n",
    "                if mod(N,2) == 1 % last must also be included\n",
    "                    x(N, :) = x(N, :) + lambda2*x(N-1, :);\n",
    "                end\n",
    "            end\n",
    "        end\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "        function x=lifting_odd(lambda1, lambda2, x, bd_mode)\n",
    "            N = size(x, 1);\n",
    "            x(2:2:(N-1), :) = ...\n",
    "                lambda1*x(3:2:N, :) + x(2:2:(N-1), :) + lambda2*x(1:2:(N-2), :);\n",
    "            if mod(N,2) == 0 % last must also be included\n",
    "                if strcmpi(bd_mode, 'per')\n",
    "                    x(N, :) = lambda1*x(1, :) + x(N, :) + lambda2*x(N-1, :);\n",
    "                elseif strcmpi(bd_mode, 'bd') || strcmpi(bd_mode, 'none')\n",
    "                    x(N, :) = x(N, :) + lambda2*x(N-1, :);\n",
    "                end\n",
    "            end\n",
    "        end\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- --- end solution of exercise --- -->\n",
    "\n",
    "\n",
    "**b)**\n",
    "Write functions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "        dwt_kernel_ortho(x, bd_mode, wav_props)\n",
    "        idwt_kernel_ortho(x, bd_mode, wav_props)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "which apply the DWT and IDWT kernels for orthonormal wavelets, using the functions `lifting_even` and `lifting_odd`. \n",
    "Assume that you can access lifting steps so that the lifting factorization (7.7) holds through the object `wav_props` by means of writing \n",
    "`wav_props.lambdas`, `wav_props.alpha`, and `wav_props.beta`, and that `wav_props.last_even` indicates whether the last lifting step is even.\n",
    "`wav_props.lambdas` is an $n\\times 2$-matrix so that the filter coefficients \n",
    "$\\{\\underline{\\lambda_1},\\lambda_2\\}$ or $\\{\\lambda_1,\\underline{\\lambda_2}\\}$ in the $i$'th\n",
    "lifting step are found in row $i$. \n",
    "\n",
    "You can now define standard DWT and IDWT kernels in the following way, once the `wav_props` object has been computed for a given $N$:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "        dwt_kernel = @(x, bd_mode) dwt_kernel_ortho(x, bd_mode, wav_props);\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- --- begin solution of exercise --- -->\n",
    "**Solution.**\n",
    "The library versions of these functions are as follows, and again support more general boundary modes:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "        function x=dwt_kernel_ortho(x, bd_mode, dual_wav_props)\n",
    "            if strcmpi(bd_mode, 'bd')\n",
    "                y1 = dual_wav_props.A_L'*x(1:size(dual_wav_props.A_L,1), :);\n",
    "                y2 = dual_wav_props.A_R'*x((end-size(dual_wav_props.A_R,1)+1):end, :);\n",
    "            end\n",
    "            x(1:2:end, :) = x(1:2:end, :)/dual_wav_props.alpha;\n",
    "            x(2:2:end, :) = x(2:2:end, :)/dual_wav_props.beta;\n",
    "            iseven = ~dual_wav_props.last_even;\n",
    "            for stepnr = (size(dual_wav_props.lambdas, 1)):(-1):1\n",
    "                if iseven\n",
    "                    x = lifting_even(dual_wav_props.lambdas(stepnr,2), ...\n",
    "                                     dual_wav_props.lambdas(stepnr,1), x, bd_mode);\n",
    "                else\n",
    "                    x = lifting_odd(dual_wav_props.lambdas(stepnr,2), ...\n",
    "                                    dual_wav_props.lambdas(stepnr,1), x, bd_mode);\n",
    "                end\n",
    "                iseven = ~iseven;\n",
    "            end\n",
    "            if strcmpi(bd_mode, 'bd')\n",
    "                x(1:size(dual_wav_props.A_L,2), :) = ...\n",
    "                    x(1:size(dual_wav_props.A_L,2), :) + y1;\n",
    "                x((end-size(dual_wav_props.A_R,2)+1):end, :) = ...\n",
    "                    x((end-size(dual_wav_props.A_R,2)+1):end, :) + y2;\n",
    "            end\n",
    "        end\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "        function x=idwt_kernel_ortho(x, bd_mode, wav_props)    \n",
    "            if strcmpi(bd_mode, 'bd')\n",
    "               y1 = wav_props.A_L*x(1:size(wav_props.A_L,2), :);\n",
    "               y2 = wav_props.A_R*x((end-size(wav_props.A_R,2)+1):end, :);\n",
    "            end\n",
    "            iseven = ( mod(size(wav_props.lambdas, 1), 2) == wav_props.last_even );\n",
    "            for stepnr = 1:(size(wav_props.lambdas, 1))\n",
    "                if iseven\n",
    "                    x = lifting_even(wav_props.lambdas(stepnr, 1), ...\n",
    "                                     wav_props.lambdas(stepnr, 2), x, bd_mode);\n",
    "                else\n",
    "                    x = lifting_odd(wav_props.lambdas(stepnr, 1), ...\n",
    "                                    wav_props.lambdas(stepnr, 2), x, bd_mode);\n",
    "                end\n",
    "                iseven = ~iseven;\n",
    "            end\n",
    "            x(1:2:end, :)=x(1:2:end, :)/wav_props.alpha;\n",
    "            x(2:2:end, :)=x(2:2:end, :)/wav_props.beta;\n",
    "        \n",
    "            if strcmpi(bd_mode, 'bd')\n",
    "                x(1:size(wav_props.A_L,1), :) = ...\n",
    "                    x(1:size(wav_props.A_L,1), :) + y1;\n",
    "                x((end-size(wav_props.A_R,1)+1):end, :) = ...\n",
    "                    x((end-size(wav_props.A_R,1)+1):end, :) + y2;\n",
    "            end\n",
    "        end\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- --- end solution of exercise --- -->\n",
    "\n",
    "**c)**\n",
    "Listen to the low-resolution approximations and detail components in sound for orthonormal wavelets for $N=1,2,3,4$.\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": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "[x, fs] = forw_comp_rev_dwt1(m, 'db2');\n",
    "playerobj = audioplayer(x, fs);\n",
    "playblocking(playerobj)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "[x, fs] = forw_comp_rev_dwt1(m, 'db3');\n",
    "playerobj = audioplayer(x, fs);\n",
    "playblocking(playerobj)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "[x, fs] = forw_comp_rev_dwt1(m, 'db4');\n",
    "playerobj = audioplayer(x, fs);\n",
    "playblocking(playerobj)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- --- end solution of exercise --- -->\n",
    "\n",
    "**d)**\n",
    "Plot all scaling functions and mother wavelets for the orthonormal\n",
    "wavelets for $N=1,2,3,4$, using the cascade algorithm. \n",
    "Since the wavelets are orthonormal, we should have that $\\phi=\\tilde{\\phi}$, and $\\psi=\\tilde{\\psi}$. \n",
    "In other words, you should see that the bottom plots\n",
    "equal the upper plots.\n",
    "\n",
    "\n",
    "<!-- --- begin solution of exercise --- -->\n",
    "**Solution.**\n",
    "The code can look as follows."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "m=10;\n",
    "cascade_alg(m, -4, 4, 'db2', 1, 0)\n",
    "cascade_alg(m, -4, 4, 'db2', 0, 0)\n",
    "\n",
    "cascade_alg(m, -4, 4, 'db3', 1, 0)\n",
    "cascade_alg(m, -4, 4, 'db3', 0, 0)\n",
    "\n",
    "cascade_alg(m, -4, 4, 'db4', 1, 0)\n",
    "cascade_alg(m, -4, 4, 'db4', 0, 0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- --- end solution of exercise --- -->\n",
    "\n",
    "<!-- --- end exercise --- -->\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "<!-- --- begin exercise --- -->\n",
    "\n",
    "## Exercise 7.13: Run forward and reverse transform\n",
    "<div id=\"exercise:080301\"></div>\n",
    "<!-- keywords = wavlifting; student -->\n",
    "\n",
    "Run the forward and then the reverse transform from Exercise 5.30 on the vector $(1,2,3,\\ldots,8192)$.\n",
    "Verify that there seems to be a delay on $481$ elements, as promised in Section 7.3.4. \n",
    "Do you get the exact same result back?\n",
    "\n",
    "\n",
    "<!-- --- begin solution of exercise --- -->\n",
    "**Solution.**\n",
    "The following code can be used:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "x = (1:8192)';\n",
    "x = mp3_reverse_fbt(mp3_forward_fbt(x));\n",
    "plot(x)\n",
    "\n",
    "figure()\n",
    "\n",
    "x = (1:8192)';\n",
    "x = mp3_reverse_fbt(x);\n",
    "x(1:size(x)-481) = x(482:size(x));\n",
    "x = mp3_forward_fbt(x);\n",
    "plot(x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are some small errors from the original vector in the resulting vector, when one compensates for the delay of $481$ elements.\n",
    "\n",
    "<!-- --- end solution of exercise --- -->\n",
    "\n",
    "<!-- --- end exercise --- -->\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "<!-- --- begin exercise --- -->\n",
    "\n",
    "## Exercise 7.14: Verify statement of filters\n",
    "<div id=\"exercise:verify_statement_filters\"></div>\n",
    "<!-- keywords = wavlifting; student -->\n",
    "\n",
    "Use your computer to verify the following symmetries for the prototype filters:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "C_i = \\begin{cases} -C_{512-i} & i\\neq 64,128,\\ldots,448 \\\\ C_{512-i} & i=64,128,\\ldots,448. \\end{cases}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Explain also that this implies that $h_i=h_{512-i}$ for $i=1,\\ldots,511$. Recall that the modified version had the symmetry $h_i=h_{511-i}$.  \n",
    "When the filters $V^{(m)}$ are defined as in this section, explain why (7.25) should be replaced by"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "V^{(64-k)} =-E_{14}(V^{(k)})^T\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "in the MP3 standard\n",
    "\n",
    "\n",
    "<!-- --- begin solution of exercise --- -->\n",
    "**Solution.**\n",
    "The following code can be used"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "ctable = mp3_ctable();\n",
    "ctable = ctable(2:end);\n",
    "ctable(64:64:255) = -ctable(64:64:255); \n",
    "res = ctable + ctable(end:-1:1);\n",
    "res(256)=0;\n",
    "max(abs(res))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that the middle element of the table is set to $0$ in this verification. The reason is that, for $i=256$, there is nothing to check. \n",
    "So, since the verification works by producing zero when there is symmetry, the element corresponding to $i=256$ is set to zero. \n",
    "Also, every $64$'th element in one half of the table changed sign, in order to establish symmetry for the corresponding entries, \n",
    "rather than anti-symmetry as for the other entries.\n",
    "\n",
    "<!-- --- end solution of exercise --- -->\n",
    "\n",
    "<!-- --- end exercise --- -->\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "<!-- --- begin exercise --- -->\n",
    "\n",
    "## Exercise 7.15: Verify near perfect reconstruction in the MP3 standard\n",
    "<div id=\"exercise:verify_near_perfect\"></div>\n",
    "<!-- keywords = wavlifting -->\n",
    "\n",
    "In Section 7.3.3 we saw that the polyphase representation of $GH$ is"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "16\\begin{pmatrix} S_0 & \\cdots & \\boldsymbol{0} \\\\ \\vdots & \\ddots & \\vdots \\\\ \\boldsymbol{0} & \\cdots & S_{31} \\end{pmatrix}E_{480},\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "where $S_i = (V^{(16+i)})^TV^{(16+i)} + (V^{(48+i)})^TV^{(48+i)}$. As explained in Section 7.3.4, the MP3 standard chooses the prototype filter so that the $S_i$ are equal\n",
    "(so that the block matrix above is a filter), but only near perfect reconstruction is possible.\n",
    "\n",
    "\n",
    "**a)**\n",
    "With $S_i=\\{...,t_{-1},\\underline{t_0},t_1,...\\}$, explain why the (magnitude of the) frequency response of (the filter) $GH$ equals"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "16\\left|\\sum_{k} t_k e^{-iMk\\omega}\\right|.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- --- begin solution of exercise --- -->\n",
    "**Solution.**\n",
    "Disregarding the delay $E_{480}$, which does not affect the absolute value of the frequency response, \n",
    "the filter coefficients of $GH$ come at indices $0$, $M$, $2M$, and so on, and are $16t_0$, $16t_1$, $16t_2$, and so on. This gives"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\begin{align*}\n",
    "|\\lambda_{GH}(\\omega)| &= 16\\left|\\sum_k t_ke^{-iMk\\omega}\\right|.\n",
    "\\end{align*}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- --- end solution of exercise --- -->\n",
    "\n",
    "**b)**\n",
    "Use the values in the table $C$ to find the filter coefficients $t_k$ for"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "S_1= (V^{(17)})^TV^{(17)} + (V^{(49)})^TV^{(49)}.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(use (7.19), even though this filter bank does not give perfect reconstruction), and plot this frequency response. \n",
    "If you get an almost flat frequency response, you have verified near perfect reconstruction in the MP3 standard. It is important that you scale the vertical axis in order to see that the frequency response is close to constant.\n",
    "\n",
    "\n",
    "<!-- --- begin solution of exercise --- -->\n",
    "**Solution.**\n",
    "The following code can be used."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "C = mp3_ctable();\n",
    "\n",
    "V1=zeros(15,1);\n",
    "V1(1:2:end)=C(18:64:end);\n",
    "\n",
    "V2=zeros(15,1); \n",
    "V2(1:2:end)=C(50:64:end);\n",
    "diagentry = conv(flip(V1),V1)+conv(flip(V2),V2);\n",
    "\n",
    "S = zeros(32*length(diagentry),1);\n",
    "S(1:32:end) = diagentry; % upsampled\n",
    "S = [S; zeros(16384,1)];\n",
    "N = length(S);\n",
    "n = (0:(N-1))';\n",
    "\n",
    "plot(2*pi*n/N, abs(fft(S)), 'k-')\n",
    "axis([0 2*pi 0 1.5/512])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- --- end solution of exercise --- -->\n",
    "\n",
    "**c)**\n",
    "Verify that the frequency response of $S_1$ in b) should be close to $1/(16\\cdot 32) = 1/512$. \n",
    "The factor $32$ has to do with the relative scaling of the $C$ and $D$ tables.\n",
    "\n",
    "\n",
    "<!-- --- begin solution of exercise --- -->\n",
    "**Solution.**\n",
    "Since the values in the table D are $32$ times those in the table C, we will instead have"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\lambda_S(\\omega)=(16\\cdot 32)\\left|\\sum_{k} t_k e^{-iMk\\omega}\\right| = 512\\left|\\sum_{k} t_k e^{-iMk\\omega}\\right|\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(recall that $S_i$ represents values from the table C). For this to be close to one we must clearly have $\\left|\\sum_{k} t_k e^{-iMk\\omega}\\right|\\approx 1/512$.\n",
    "\n",
    "<!-- --- end solution of exercise --- -->\n",
    "\n",
    "\n",
    "\n",
    "<!-- --- end exercise --- -->"
   ]
  }
 ],
 "metadata": {},
 "nbformat": 4,
 "nbformat_minor": 2
}