{"worksheets": [{"cells": [{"cell_type": "markdown", "metadata": {}, "source": ["1-D Daubechies Wavelets", "=======================", "", "*Important:* Please read the [installation page](http://gpeyre.github.io/numerical-tours/installation_python/) for details about how to install the toolboxes.", "$\\newcommand{\\dotp}[2]{\\langle #1, #2 \\rangle}$", "$\\newcommand{\\enscond}[2]{\\lbrace #1, #2 \\rbrace}$", "$\\newcommand{\\pd}[2]{ \\frac{ \\partial #1}{\\partial #2} }$", "$\\newcommand{\\umin}[1]{\\underset{#1}{\\min}\\;}$", "$\\newcommand{\\umax}[1]{\\underset{#1}{\\max}\\;}$", "$\\newcommand{\\umin}[1]{\\underset{#1}{\\min}\\;}$", "$\\newcommand{\\uargmin}[1]{\\underset{#1}{argmin}\\;}$", "$\\newcommand{\\norm}[1]{\\|#1\\|}$", "$\\newcommand{\\abs}[1]{\\left|#1\\right|}$", "$\\newcommand{\\choice}[1]{ \\left\\{  \\begin{array}{l} #1 \\end{array} \\right. }$", "$\\newcommand{\\pa}[1]{\\left(#1\\right)}$", "$\\newcommand{\\diag}[1]{{diag}\\left( #1 \\right)}$", "$\\newcommand{\\qandq}{\\quad\\text{and}\\quad}$", "$\\newcommand{\\qwhereq}{\\quad\\text{where}\\quad}$", "$\\newcommand{\\qifq}{ \\quad \\text{if} \\quad }$", "$\\newcommand{\\qarrq}{ \\quad \\Longrightarrow \\quad }$", "$\\newcommand{\\ZZ}{\\mathbb{Z}}$", "$\\newcommand{\\CC}{\\mathbb{C}}$", "$\\newcommand{\\RR}{\\mathbb{R}}$", "$\\newcommand{\\EE}{\\mathbb{E}}$", "$\\newcommand{\\Zz}{\\mathcal{Z}}$", "$\\newcommand{\\Ww}{\\mathcal{W}}$", "$\\newcommand{\\Vv}{\\mathcal{V}}$", "$\\newcommand{\\Nn}{\\mathcal{N}}$", "$\\newcommand{\\NN}{\\mathcal{N}}$", "$\\newcommand{\\Hh}{\\mathcal{H}}$", "$\\newcommand{\\Bb}{\\mathcal{B}}$", "$\\newcommand{\\Ee}{\\mathcal{E}}$", "$\\newcommand{\\Cc}{\\mathcal{C}}$", "$\\newcommand{\\Gg}{\\mathcal{G}}$", "$\\newcommand{\\Ss}{\\mathcal{S}}$", "$\\newcommand{\\Pp}{\\mathcal{P}}$", "$\\newcommand{\\Ff}{\\mathcal{F}}$", "$\\newcommand{\\Xx}{\\mathcal{X}}$", "$\\newcommand{\\Mm}{\\mathcal{M}}$", "$\\newcommand{\\Ii}{\\mathcal{I}}$", "$\\newcommand{\\Dd}{\\mathcal{D}}$", "$\\newcommand{\\Ll}{\\mathcal{L}}$", "$\\newcommand{\\Tt}{\\mathcal{T}}$", "$\\newcommand{\\si}{\\sigma}$", "$\\newcommand{\\al}{\\alpha}$", "$\\newcommand{\\la}{\\lambda}$", "$\\newcommand{\\ga}{\\gamma}$", "$\\newcommand{\\Ga}{\\Gamma}$", "$\\newcommand{\\La}{\\Lambda}$", "$\\newcommand{\\si}{\\sigma}$", "$\\newcommand{\\Si}{\\Sigma}$", "$\\newcommand{\\be}{\\beta}$", "$\\newcommand{\\de}{\\delta}$", "$\\newcommand{\\De}{\\Delta}$", "$\\newcommand{\\phi}{\\varphi}$", "$\\newcommand{\\th}{\\theta}$", "$\\newcommand{\\om}{\\omega}$", "$\\newcommand{\\Om}{\\Omega}$"]}, {"cell_type": "markdown", "metadata": {}, "source": ["This numerical tour explores 1-D multiresolution analysis", "with Daubechies wavelets with a varying number of vanishing moments", "(varying order)."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["from __future__ import division", "import nt_toolbox as nt", "from nt_solutions import wavelet_3_daubechies1d as solutions", "%matplotlib inline", "%load_ext autoreload", "%autoreload 2"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Wavelets Filters", "----------------", "The 2-D wavelet transform of a continuous signal $f(x)$ computes the set", "of inner products", "$$ d_j[n] = \\dotp{f}{\\psi_{j,n}} $$", "for scales $ j \\in \\ZZ", "$ and position $ n \\in \\ZZ $.", "", "", "The wavelet atoms are defined by scaling and translating a mother", "atom $ \\psi(x) $:", "$$ \\psi_{j,n}(x) = \\frac{1}{2^j}\\psi\\pa{\\frac{x-2^j n}{2^j}}.  $$", "", "", "Associated to this oscillating (high pass) wavelet function $\\psi$", "is a non-oscillating (low pass) scaling function $\\phi$.", "", "", "The fast wavelet transform algorithm does not make use of the wavelet and scaling functions,", "but of the filters $h$ and $g$ that caracterize their interaction:", "$$ g[n] = \\frac{1}{\\sqrt{2}}\\dotp{\\psi(x/2)}{\\phi(x-n)}", "\\qandq h[n] = \\frac{1}{\\sqrt{2}}\\dotp{\\phi(x/2)}{\\phi(x-n)}. $$", "", "", "The simplest filters are the Haar filters", "$$ h = [1, 1]/\\sqrt{2} \\qandq g = [-1, 1]/\\sqrt{2}. $$", "", "", "Daubechies wavelets extends the haar wavelets by using longer", "filters, that produce smoother scaling functions and wavelets.", "Furthermore, the larger the size $p=2k$ of the filter, the higher is the number", "$k$ of vanishing moment.", "", "", "A high number of vanishing moments allows to better compress regular", "parts of the signal. However, increasing the number of vanishing moments", "also inceases the size of the support of the wavelets, wich can be", "problematic in part where the signal is singular (for instance", "discontinuous).", "", "", "Choosing the _best_ wavelet, and thus choosing $k$, that is adapted to a", "given class of signals, thus corresponds to", "a tradeoff between efficiency in regular and singular parts.", "", "", "* The filter with $k=1$ vanishing moments corresponds to the Haar filter.", "* The filter with $k=2$ vanishing moments corresponds to the famous |D4| wavelet, which compresses perfectly linear signals.", "* The filter with $k=3$ vanishing moments compresses perfectly quadratic signals.", "", "", "Set the support size.", "To begin, we select the D4 filter."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["p = 4"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Create the low pass filter $h$ and the high pass $g$. We add a zero to ensure that it has a odd", "length. Note that the central value of $h$ corresponds to the 0 position."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["[h, g] = compute_wavelet_filter('Daubechies', p)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Note that the high pass filter $g$ is computed directly from the low", "pass filter as:", "$$g[n] = (-1)^{1-n}h[1-n]$$", "", "", "Display."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["disp(['h filter = [' num2str(h) ']'])", "disp(['g filter = [' num2str(g) ']'])"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Up and Down Filtering", "---------------------", "The basic wavelet operation is low/high filtering, followed by down", "sampling.", "", "", "Starting from some 1-D signal $f \\in \\RR^N$, one thus computes the", "low pass signal $a \\in \\RR^{N/2}$ and the high pass", "signal $d \\in \\RR^{N/2}$ as", "$$ a = (f \\star h) \\downarrow 2 \\qandq", "d = (f \\star g) \\downarrow 2$$", "where the sub-sampling is defined as", "$$ (u \\downarrow 2)[k] = u[2k]. $$", "", "", "Create a random signal $f \\in \\RR^N$."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["N = 256", "f = rand(N, 1)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Low/High pass filtering followed by sub-sampling."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["a = subsampling(cconvol(f, h))", "d = subsampling(cconvol(f, g))"]}, {"cell_type": "markdown", "metadata": {}, "source": ["For orthogonal filters, the reverse of this process is its dual", "(aka its transpose), which is upsampling followed by low/high pass", "filtering with the reversed filters and summing:", "$$ (a \\uparrow h) \\star \\tilde h + (d \\uparrow g) \\star \\tilde g = f $$", "where $\\tilde h[n]=h[-n]$ (computed modulo $N$) and", "$ (u \\uparrow 2)[2n]=u[n] $ and  $ (u \\uparrow 2)[2n+1]=0 $.", "", "", "Up-sampling followed by filtering."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["f1 =  cconvol(upsampling(a), reverse(h)) + cconvol(upsampling(d), reverse(g))"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Check that we really recover the same signal."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["disp(strcat((['Error |f-f1|/ |f| = ' num2str(norm(f-f1)/ norm(f))])))"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Forward Wavelet Transform", "-------------------------", "The set of wavelet coefficients are computed with a fast algorithm that", "exploit the embedding of the approximation spaces $V_j$ spanned by the", "scaling function $ \\{ \\phi_{j,n} \\}_n $.", "", "", "First we load a 1-D signal."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["name = 'piece-regular'", "N = 512", "f = rescale(load_signal(name, N))"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Display it."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["plot(f)", "axis('tight')"]}, {"cell_type": "markdown", "metadata": {}, "source": ["We will store all the transformed coefficients $d_j$ in a single vector", "|fw|. This vector is initialized as |f| and the left sub-part", "|fw(1:1^j)| of |fw| will be retransformed at each iteration for a decreasing scale", "index |j|.", "", "", "Initialize the result vector."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["fw = f"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Initialize the scale index $j$ as $ j = J = \\log_2(N)-1 $."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["j = log2(N)-1"]}, {"cell_type": "markdown", "metadata": {}, "source": ["The wavelet transform of $f$ is computed by using intermediate discretized low", "resolution images obtained by projection on the spaces $V_j$:", "$$ a_j[n] = \\dotp{f}{\\phi_{j,n}}. $$", "", "", "The algorithm processes by moving from scale $j$ to the coarser scale", "$j-1$ using the filtering+sub-sampling:", "$$ a_{j-1} = (a_j \\star h) \\downarrow 2 \\qandq", "d_{j-1} = (a_j \\star g) \\downarrow 2$$", "", "", "Retrieve the coefficients $a_j$ from the variable |fw| and", "store them in the variable |a1|"]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["a1 = fw(1: 2^(j + 1))"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Apply high and low filtering+subsampling", "to obtain $a_{j-1}$ and $d_{j-1}$ (stored in |a| and |d|)."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["a = subsampling(cconvol(a1, h))", "d = subsampling(cconvol(a1, g))"]}, {"cell_type": "markdown", "metadata": {}, "source": ["_Note:_ |subsampling(A)| is equivalent to |A(1:2:end)|.", "", "", "Concatenate them to get the result and store it in |fw|."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["fw(1: 2^(j + 1)) = cat(1, a, d)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Display the result of the first step of the transform."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["subplot(2, 1, 1)", "plot(f); axis('tight'); title('Signal')", "subplot(2, 1, 2)", "plot(fw); axis('tight'); title('Transformed')"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Display the signal and its coarse coefficients."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["s = 400; t = 40", "", "subplot(2, 1, 1)", "plot(f, '.-'); axis([s-t s + t 0 1]); title('Signal (zoom)')", "subplot(2, 1, 2)", "plot(a, '.-'); axis([(s-t)/ 2 (s + t)/ 2 min(a) max(a)]); title('Averages (zoom)')"]}, {"cell_type": "markdown", "metadata": {}, "source": ["__Exercise 1__", "", "Implement a full wavelet transform that extract iteratively wavelet", "coefficients, by repeating these steps. Take care of choosing the", "correct number of steps."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["solutions.exo1()"]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["## Insert your code here."]}, {"cell_type": "markdown", "metadata": {}, "source": ["Check that the transform is", "orthogonal, which means that the energy of the coefficient is the same", "as the energy of the signal."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["disp(strcat(['Energy of the signal       = ' num2str(norm(f).^2, 3)]))", "disp(strcat(['Energy of the coefficients = ' num2str(norm(fw).^2, 3)]))"]}, {"cell_type": "markdown", "metadata": {}, "source": ["We display the whole set of coefficients |fw|, with red vertical", "separator between the scales.", "Can you recognize where are the low frequencies and the high frequencies", "? You can use the function |plot_wavelet| to help you."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["plot_wavelet(fw)", "axis([1 N -1 1])"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Backward Wavelet Transform", "--------------------------", "The backward wavelet transform reconstructs a signal $ f $ from a set of", "wavelet coeffcients $ \\{ d_j[n] \\}_{j,n} $. For continuous functions,", "it corresponds to the following reconstruction formula:", "$$ f(x) = \\sum_{j,n} d_j[n] \\psi_{j,n}(x).  $$", "", "", "For discrete signal, it reconstructs a signal $ f \\in \\RR^N $ by", "inverting the wavelet filtering/sub-sampling steps.", "", "", "It starts from the coarsest scale $ j=0 $, where $ a_0 \\in \\RR $ is", "the single remaining coefficient.", "", "", "The algorithm processes by moving from scale $j$ to the finer scale", "$j+1$ using the up-sampling/filtering:", "$$ a_{j+1} = (a_j \\uparrow 2) \\star \\tilde h", "   + (d_j \\uparrow 2) \\star \\tilde g $$", "", "", "Initialize the signal to recover |f1| as the transformed coefficient, and", "select the smallezt possible scale."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["f1 = fw", "j = 0"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Retrieve coarse and detail coefficients in the vertical direction."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["a = f1(1: 2^j)", "d = f1(2^j + 1: 2^(j + 1))"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Perform the up-sampling/filtering and summation:"]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["a = cconvol(upsampling(a, 1), reverse(h), 1) + cconvol(upsampling(d, 1), reverse(g), 1)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Replace the coefficients at the correct locations."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["f1(1: 2^(j + 1)) = a"]}, {"cell_type": "markdown", "metadata": {}, "source": ["__Exercise 2__", "", "Write the inverse wavelet transform that computes |f1| from the", "coefficients |fw|."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["solutions.exo2()"]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["## Insert your code here."]}, {"cell_type": "markdown", "metadata": {}, "source": ["Check that we have correctly recovered the signal."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["disp(strcat((['Error |f-f1|/ |f| = ' num2str(norm(f-f1)/ norm(f))])))"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Daubechies Wavelets Approximation", "---------------------------------", "Non-linear approximation is obtained by thresholding low amplitude", "wavelet coefficients.", "", "", "This defines the best $M$-terms approximation $f_M$ of $f$:", "", "", "$$ f_M = \\sum_{ \\abs{\\dotp{f}{\\psi_{j,n}}}>T } \\dotp{f}{\\psi_{j,n}}\\psi_{j,n}. $$", "", "", "Set the threshold value."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["T = .5"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Coefficients |fw(i)| smaller in magnitude than |T| are set to zero."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["fwT = fw .* (abs(fw) >T)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Display the coefficients before and after thresholding."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["subplot(2, 1, 1)", "plot_wavelet(fw); axis([1 N -1 1]); title('Original coefficients')", "subplot(2, 1, 2)", "plot_wavelet(fwT); axis([1 N -1 1]); title('Thresholded coefficients')"]}, {"cell_type": "markdown", "metadata": {}, "source": ["__Exercise 3__", "", "Find the threshold $T$ to obtained a given number $M$ of", "non thresholded coefficients.", "Try for an increasing number $M$ of coeffiients.", "ompute the threshold T"]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["solutions.exo3()"]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["## Insert your code here."]}, {"cell_type": "markdown", "metadata": {}, "source": ["__Exercise 4__", "", "Try with", "Different kind of wavelets, with an increasing number of vanishing", "moments.", "ompute the threshold T"]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["solutions.exo4()"]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["## Insert your code here."]}, {"cell_type": "markdown", "metadata": {}, "source": ["The Shape of a Wavelet", "----------------------", "A wavelet coefficient is an inner product", "$ d_j[n] = \\dotp{f}{\\psi_{j,n}} $ with a wavelet atom $\\psi_{j,n}$.", "", "", "A wavelet atom $\\psi_{j_0,n_0}$ can be computed by", "applying the inverse wavele transform to coefficients $ \\{d_j[n]\\}_{j,n} $", "such that", "$$ d_{j}[n]=\\choice{ 1 \\qifq j=j_0 \\qandq n=n_0, \\\\ 0 \\quad\\text{otherwise.} } $$"]}, {"cell_type": "markdown", "metadata": {}, "source": ["__Exercise 5__", "", "Compute wavelets at several positions and scales."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["solutions.exo5()"]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["## Insert your code here."]}, {"cell_type": "markdown", "metadata": {}, "source": ["__Exercise 6__", "", "Display Daubechies wavelets with an increasing number of vanishing", "moments."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["solutions.exo6()"]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["## Insert your code here."]}]}], "nbformat": 3, "metadata": {"name": ""}, "nbformat_minor": 0}