{
  "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_gpu.signalprocessing.Convolve1D`\noperator to perform convolution between two signals.\n\nThis example closely follow the equivalent\n\n`PyLops example <https://pylops.readthedocs.io/en/latest/gallery/plot_convolve.html#sphx-glr-gallery-plot-convolve-py>`_.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import torch\nimport numpy as np\nimport matplotlib.pyplot as plt\nimport pylops_gpu\n\nfrom pylops.utils.wavelets import ricker\nfrom pylops_gpu.utils.backend import device\nfrom pylops_gpu.optimization.cg import cg\n\ndev = device()\nprint('PyLops-gpu working on %s...' % dev)\nplt.close('all')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We will start by creating a zero signal of length $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_gpu.signalprocessing.Convolve1D`\noperator.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "nt = 1001\ndt = 0.004\nt = np.arange(nt)*dt\n\nx = torch.zeros(nt, dtype=torch.float32)\nx[int(nt/2)] = 1\n\nh, th, hcenter = ricker(t[:101], f0=30)\nh = torch.from_numpy(h.astype(np.float32))\nCop = pylops_gpu.signalprocessing.Convolve1D(nt, h=h, offset=hcenter,\n                                             dtype=torch.float32)\ny = Cop*x\n\nxinv = Cop / y\n\nfig, ax = plt.subplots(1, 1, figsize=(10, 3))\nax.plot(t, x.cpu().numpy(), 'k', lw=2, label=r'$x$')\nax.plot(t, y.cpu().numpy(), 'r', lw=2, label=r'$y=Ax$')\nax.plot(t, xinv.cpu().numpy(), '--g', lw=2, label=r'$x_{ext}$')\nax.set_title('Convolve in 1st direction', fontsize=14, fontweight='bold')\nax.legend()\nax.set_xlim(1.9, 2.1)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We show now that also a filter with mixed phase (i.e., not centered around zero)\ncan be applied and inverted for using the :py:class:`pylops.signalprocessing.Convolve1D`\noperator.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "Cop = pylops_gpu.signalprocessing.Convolve1D(nt, h=h, offset=hcenter - 3,\n                                             dtype=torch.float32)\ny = Cop * x\ny1 = Cop.H * x\nxinv = cg(Cop.H*Cop, Cop.H*y, niter=100)[0]\n\nfig, ax = plt.subplots(1, 1, figsize=(10, 3))\nax.plot(t, x.cpu().numpy(), 'k', lw=2, label=r'$x$')\nax.plot(t, y.cpu().numpy(), 'r', lw=2, label=r'$y=Ax$')\nax.plot(t, y1.cpu().numpy(), 'b', lw=2, label=r'$y=A^Hx$')\nax.plot(t, xinv.cpu().numpy(), '--g', lw=2, label=r'$x_{ext}$')\nax.set_title('Convolve in 1st direction', fontsize=14, fontweight='bold')\nax.set_xlim(1.9, 2.1)\nax.legend()"
      ]
    }
  ],
  "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.6.12"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}