{
  "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\n:py:class:`pylops_distributed.signalprocessing.Convolve1D` operator to perform\nconvolution between two signals.\n\nNote that when the input dataset is distributed across multiple nodes,\nadditional care should be taken when applying convolution.\nIn PyLops-distributed, we leverage the :func:`dask.array.map_block` and\n:func:`dask.array.map_overlap` functionalities of dask when chunking is\nperformed along the dimension of application of the convolution and along any\nother dimension, respectively. This is however handled by the inner working of\n:py:class:`pylops_distributed.signalprocessing.Convolve1D`.\n\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import numpy as np\nimport dask.array as da\nimport matplotlib.pyplot as plt\nfrom scipy.sparse.linalg import lsqr\n\nimport pylops_distributed\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\n:py:class:`pylops-distributed.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\nx = da.from_array(x, chunks=nt // 2 + 1)\nh, th, hcenter = ricker(t[:101], f0=30)\n\nCop = pylops_distributed.signalprocessing.Convolve1D(nt, h=h, offset=hcenter,\n                                                     chunks=(nt // 2 + 1,\n                                                             nt // 2 + 1),\n                                                     dtype='float32')\ny = Cop * 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, 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)"
      ]
    },
    {
      "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_distributed.signalprocessing.Convolve1D(nt, h=h, offset=hcenter - 3,\n                                                     chunks=(nt // 2 + 1,\n                                                             nt // 2 + 1),\n                                                     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('Convolve 1d data with non-zero phase filter', fontsize=14,\n             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
}