{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# Convolution\nThis example shows how to use the :py:class:`pylops.signalprocessing.Convolve1D`,\n:py:class:`pylops.signalprocessing.Convolve2D` and\n:py:class:`pylops.signalprocessing.ConvolveND` operators to perform convolution\nbetween two signals.\n\nSuch operators can be used in the forward model of several common application\nin signal processing that require filtering of an input signal for the\ninstrument response. Similarly, removing the effect of the instrument\nresponse from signal is equivalent to solving linear system of equations\nbased on Convolve1D, Convolve2D or ConvolveND operators.\nThis problem is generally referred to as *Deconvolution*.\n\nA very practical example of deconvolution can be found in the geophysical\nprocessing of seismic data where the effect of the source response\n(i.e., airgun or vibroseis) should be removed from the recorded signal\nto be able to better interpret the response of the subsurface. Similar examples\ncan be found in telecommunication and speech analysis.\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import matplotlib.pyplot as plt\nimport numpy as np\nfrom scipy.sparse.linalg import lsqr\n\nimport pylops\nfrom pylops.utils.wavelets import ricker\n\nplt.close(\"all\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We will start by creating a zero signal of lenght $nt$ and we will\nplace a unitary spike at its center. We also create our filter to be\napplied by means of :py:class:`pylops.signalprocessing.Convolve1D` operator.\nFollowing the seismic example mentioned above, the filter is a\n[Ricker wavelet](http://subsurfwiki.org/wiki/Ricker_wavelet)\nwith dominant frequency $f_0 = 30 Hz$.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "nt = 1001\ndt = 0.004\nt = np.arange(nt) * dt\nx = np.zeros(nt)\nx[int(nt / 2)] = 1\nh, th, hcenter = ricker(t[:101], f0=30)\n\nCop = pylops.signalprocessing.Convolve1D(nt, h=h, offset=hcenter, dtype=\"float32\")\ny = Cop * x\n\nxinv = Cop / y\n\nfig, ax = plt.subplots(1, 1, figsize=(10, 3))\nax.plot(t, x, \"k\", lw=2, label=r\"$x$\")\nax.plot(t, y, \"r\", lw=2, label=r\"$y=Ax$\")\nax.plot(t, xinv, \"--g\", lw=2, label=r\"$x_{ext}$\")\nax.set_title(\"Convolve 1d data\", fontsize=14, fontweight=\"bold\")\nax.legend()\nax.set_xlim(1.9, 2.1)\nplt.tight_layout()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We show now that also a filter with mixed phase (i.e., not centered\naround zero) can be applied and inverted for using the\n:py:class:`pylops.signalprocessing.Convolve1D`\noperator.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "Cop = pylops.signalprocessing.Convolve1D(nt, h=h, offset=hcenter - 3, dtype=\"float32\")\ny = Cop * x\ny1 = Cop.H * x\nxinv = Cop / y\n\nfig, ax = plt.subplots(1, 1, figsize=(10, 3))\nax.plot(t, x, \"k\", lw=2, label=r\"$x$\")\nax.plot(t, y, \"r\", lw=2, label=r\"$y=Ax$\")\nax.plot(t, y1, \"b\", lw=2, label=r\"$y=A^Hx$\")\nax.plot(t, xinv, \"--g\", lw=2, label=r\"$x_{ext}$\")\nax.set_title(\n    \"Convolve 1d data with non-zero phase filter\", fontsize=14, fontweight=\"bold\"\n)\nax.set_xlim(1.9, 2.1)\nax.legend()\nplt.tight_layout()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We repeat a similar exercise but using two dimensional signals and\nfilters taking advantage of the\n:py:class:`pylops.signalprocessing.Convolve2D` operator.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "nt = 51\nnx = 81\ndt = 0.004\nt = np.arange(nt) * dt\nx = np.zeros((nt, nx))\nx[int(nt / 2), int(nx / 2)] = 1\n\nnh = [11, 5]\nh = np.ones((nh[0], nh[1]))\n\nCop = pylops.signalprocessing.Convolve2D(\n    (nt, nx),\n    h=h,\n    offset=(int(nh[0]) / 2, int(nh[1]) / 2),\n    dtype=\"float32\",\n)\ny = Cop * x\nxinv = (Cop / y.ravel()).reshape(Cop.dims)\n\nfig, axs = plt.subplots(1, 3, figsize=(10, 3))\nfig.suptitle(\"Convolve 2d data\", fontsize=14, fontweight=\"bold\", y=0.95)\naxs[0].imshow(x, cmap=\"gray\", vmin=-1, vmax=1)\naxs[1].imshow(y, cmap=\"gray\", vmin=-1, vmax=1)\naxs[2].imshow(xinv, cmap=\"gray\", vmin=-1, vmax=1)\naxs[0].set_title(\"x\")\naxs[0].axis(\"tight\")\naxs[1].set_title(\"y\")\naxs[1].axis(\"tight\")\naxs[2].set_title(\"xlsqr\")\naxs[2].axis(\"tight\")\nplt.tight_layout()\nplt.subplots_adjust(top=0.8)\n\nfig, ax = plt.subplots(1, 2, figsize=(10, 3))\nfig.suptitle(\"Convolve in 2d data - traces\", fontsize=14, fontweight=\"bold\", y=0.95)\nax[0].plot(x[int(nt / 2), :], \"k\", lw=2, label=r\"$x$\")\nax[0].plot(y[int(nt / 2), :], \"r\", lw=2, label=r\"$y=Ax$\")\nax[0].plot(xinv[int(nt / 2), :], \"--g\", lw=2, label=r\"$x_{ext}$\")\nax[1].plot(x[:, int(nx / 2)], \"k\", lw=2, label=r\"$x$\")\nax[1].plot(y[:, int(nx / 2)], \"r\", lw=2, label=r\"$y=Ax$\")\nax[1].plot(xinv[:, int(nx / 2)], \"--g\", lw=2, label=r\"$x_{ext}$\")\nax[0].legend()\nax[0].set_xlim(30, 50)\nax[1].legend()\nax[1].set_xlim(10, 40)\nplt.tight_layout()\nplt.subplots_adjust(top=0.8)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Finally we do the same using three dimensional signals and\nfilters taking advantage of the\n:py:class:`pylops.signalprocessing.ConvolveND` operator.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "ny, nx, nz = 13, 10, 7\nx = np.zeros((ny, nx, nz))\nx[ny // 3, nx // 2, nz // 4] = 1\nh = np.ones((3, 5, 3))\noffset = [1, 2, 1]\n\nCop = pylops.signalprocessing.ConvolveND(\n    dims=(ny, nx, nz), h=h, offset=offset, axes=(0, 1, 2), dtype=\"float32\"\n)\ny = Cop * x\nxlsqr = lsqr(Cop, y.ravel(), damp=0, iter_lim=300, show=0)[0]\nxlsqr = xlsqr.reshape(Cop.dims)\n\nfig, axs = plt.subplots(3, 3, figsize=(10, 12))\nfig.suptitle(\"Convolve 3d data\", y=0.95, fontsize=14, fontweight=\"bold\")\naxs[0][0].imshow(x[ny // 3], cmap=\"gray\", vmin=-1, vmax=1)\naxs[0][1].imshow(y[ny // 3], cmap=\"gray\", vmin=-1, vmax=1)\naxs[0][2].imshow(xlsqr[ny // 3], cmap=\"gray\", vmin=-1, vmax=1)\naxs[0][0].set_title(\"x\")\naxs[0][0].axis(\"tight\")\naxs[0][1].set_title(\"y\")\naxs[0][1].axis(\"tight\")\naxs[0][2].set_title(\"xlsqr\")\naxs[0][2].axis(\"tight\")\naxs[1][0].imshow(x[:, nx // 2], cmap=\"gray\", vmin=-1, vmax=1)\naxs[1][1].imshow(y[:, nx // 2], cmap=\"gray\", vmin=-1, vmax=1)\naxs[1][2].imshow(xlsqr[:, nx // 2], cmap=\"gray\", vmin=-1, vmax=1)\naxs[1][0].axis(\"tight\")\naxs[1][1].axis(\"tight\")\naxs[1][2].axis(\"tight\")\naxs[2][0].imshow(x[..., nz // 4], cmap=\"gray\", vmin=-1, vmax=1)\naxs[2][1].imshow(y[..., nz // 4], cmap=\"gray\", vmin=-1, vmax=1)\naxs[2][2].imshow(xlsqr[..., nz // 4], cmap=\"gray\", vmin=-1, vmax=1)\naxs[2][0].axis(\"tight\")\naxs[2][1].axis(\"tight\")\naxs[2][2].axis(\"tight\")\nplt.tight_layout()"
      ]
    }
  ],
  "metadata": {
    "kernelspec": {
      "display_name": "Python 3",
      "language": "python",
      "name": "python3"
    },
    "language_info": {
      "codemirror_mode": {
        "name": "ipython",
        "version": 3
      },
      "file_extension": ".py",
      "mimetype": "text/x-python",
      "name": "python",
      "nbconvert_exporter": "python",
      "pygments_lexer": "ipython3",
      "version": "3.9.15"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}