{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# Fourier Transform\nThis example shows how to use the :py:class:`pylops.signalprocessing.FFT`,\n:py:class:`pylops.signalprocessing.FFT2D`\nand :py:class:`pylops.signalprocessing.FFTND` operators to apply the Fourier\nTransform to the model and the inverse Fourier Transform to the data.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import matplotlib.pyplot as plt\nimport numpy as np\n\nimport pylops\n\nplt.close(\"all\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Let's start by applying the one dimensional FFT to a one dimensional\nsinusoidal signal $d(t)=sin(2 \\pi f_0t)$ using a time axis of\nlenght $nt$ and sampling $dt$\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "dt = 0.005\nnt = 100\nt = np.arange(nt) * dt\nf0 = 10\nnfft = 2 ** 10\nd = np.sin(2 * np.pi * f0 * t)\n\nFFTop = pylops.signalprocessing.FFT(dims=nt, nfft=nfft, sampling=dt, engine=\"numpy\")\nD = FFTop * d\n\n# Adjoint = inverse for FFT\ndinv = FFTop.H * D\ndinv = FFTop / D\n\nfig, axs = plt.subplots(1, 2, figsize=(10, 4))\naxs[0].plot(t, d, \"k\", lw=2, label=\"True\")\naxs[0].plot(t, dinv, \"--r\", lw=2, label=\"Inverted\")\naxs[0].legend()\naxs[0].set_title(\"Signal\")\naxs[1].plot(FFTop.f[: int(FFTop.nfft / 2)], np.abs(D[: int(FFTop.nfft / 2)]), \"k\", lw=2)\naxs[1].set_title(\"Fourier Transform\")\naxs[1].set_xlim([0, 3 * f0])"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "In this example we used numpy as our engine for the ``fft`` and ``ifft``.\nPyLops implements a second engine (``engine='fftw'``) which uses the\nwell-known `FFTW <http://www.fftw.org>`_ via the python wrapper\n:py:class:`pyfftw.FFTW`. This optimized fft tends to outperform the one from\nnumpy in many cases but it is not inserted in the mandatory requirements of\nPyLops. If interested to use ``FFTW`` backend, read the `fft routines`\nsection at `performance`.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "FFTop = pylops.signalprocessing.FFT(dims=nt, nfft=nfft, sampling=dt, engine=\"fftw\")\nD = FFTop * d\n\n# Adjoint = inverse for FFT\ndinv = FFTop.H * D\ndinv = FFTop / D\n\nfig, axs = plt.subplots(1, 2, figsize=(10, 4))\naxs[0].plot(t, d, \"k\", lw=2, label=\"True\")\naxs[0].plot(t, dinv, \"--r\", lw=2, label=\"Inverted\")\naxs[0].legend()\naxs[0].set_title(\"Signal\")\naxs[1].plot(FFTop.f[: int(FFTop.nfft / 2)], np.abs(D[: int(FFTop.nfft / 2)]), \"k\", lw=2)\naxs[1].set_title(\"Fourier Transform with fftw\")\naxs[1].set_xlim([0, 3 * f0])"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We can also apply the one dimensional FFT to to a two-dimensional\nsignal (along one of the first axis)\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "dt = 0.005\nnt, nx = 100, 20\nt = np.arange(nt) * dt\nf0 = 10\nnfft = 2 ** 10\nd = np.outer(np.sin(2 * np.pi * f0 * t), np.arange(nx) + 1)\n\nFFTop = pylops.signalprocessing.FFT(dims=(nt, nx), dir=0, nfft=nfft, sampling=dt)\nD = FFTop * d.ravel()\n\n# Adjoint = inverse for FFT\ndinv = FFTop.H * D\ndinv = FFTop / D\ndinv = np.real(dinv).reshape(nt, nx)\n\nfig, axs = plt.subplots(2, 2, figsize=(10, 6))\naxs[0][0].imshow(d, vmin=-20, vmax=20, cmap=\"seismic\")\naxs[0][0].set_title(\"Signal\")\naxs[0][0].axis(\"tight\")\naxs[0][1].imshow(np.abs(D.reshape(nfft, nx)[:200, :]), cmap=\"seismic\")\naxs[0][1].set_title(\"Fourier Transform\")\naxs[0][1].axis(\"tight\")\naxs[1][0].imshow(dinv, vmin=-20, vmax=20, cmap=\"seismic\")\naxs[1][0].set_title(\"Inverted\")\naxs[1][0].axis(\"tight\")\naxs[1][1].imshow(d - dinv, vmin=-20, vmax=20, cmap=\"seismic\")\naxs[1][1].set_title(\"Error\")\naxs[1][1].axis(\"tight\")\nfig.tight_layout()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We can also apply the two dimensional FFT to to a two-dimensional signal\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "dt, dx = 0.005, 5\nnt, nx = 100, 201\nt = np.arange(nt) * dt\nx = np.arange(nx) * dx\nf0 = 10\nnfft = 2 ** 10\nd = np.outer(np.sin(2 * np.pi * f0 * t), np.arange(nx) + 1)\n\nFFTop = pylops.signalprocessing.FFT2D(\n    dims=(nt, nx), nffts=(nfft, nfft), sampling=(dt, dx)\n)\nD = FFTop * d.ravel()\n\ndinv = FFTop.H * D\ndinv = FFTop / D\ndinv = np.real(dinv).reshape(nt, nx)\n\nfig, axs = plt.subplots(2, 2, figsize=(10, 6))\naxs[0][0].imshow(d, vmin=-100, vmax=100, cmap=\"seismic\")\naxs[0][0].set_title(\"Signal\")\naxs[0][0].axis(\"tight\")\naxs[0][1].imshow(\n    np.abs(np.fft.fftshift(D.reshape(nfft, nfft), axes=1)[:200, :]), cmap=\"seismic\"\n)\naxs[0][1].set_title(\"Fourier Transform\")\naxs[0][1].axis(\"tight\")\naxs[1][0].imshow(dinv, vmin=-100, vmax=100, cmap=\"seismic\")\naxs[1][0].set_title(\"Inverted\")\naxs[1][0].axis(\"tight\")\naxs[1][1].imshow(d - dinv, vmin=-100, vmax=100, cmap=\"seismic\")\naxs[1][1].set_title(\"Error\")\naxs[1][1].axis(\"tight\")\nfig.tight_layout()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Finally can apply the three dimensional FFT to to a three-dimensional signal\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "dt, dx, dy = 0.005, 5, 3\nnt, nx, ny = 30, 21, 11\nt = np.arange(nt) * dt\nx = np.arange(nx) * dx\ny = np.arange(nx) * dy\nf0 = 10\nnfft = 2 ** 6\nnfftk = 2 ** 5\n\nd = np.outer(np.sin(2 * np.pi * f0 * t), np.arange(nx) + 1)\nd = np.tile(d[:, :, np.newaxis], [1, 1, ny])\n\nFFTop = pylops.signalprocessing.FFTND(\n    dims=(nt, nx, ny), nffts=(nfft, nfftk, nfftk), sampling=(dt, dx, dy)\n)\nD = FFTop * d.ravel()\n\ndinv = FFTop.H * D\ndinv = FFTop / D\ndinv = np.real(dinv).reshape(nt, nx, ny)\n\nfig, axs = plt.subplots(2, 2, figsize=(10, 6))\naxs[0][0].imshow(d[:, :, ny // 2], vmin=-20, vmax=20, cmap=\"seismic\")\naxs[0][0].set_title(\"Signal\")\naxs[0][0].axis(\"tight\")\naxs[0][1].imshow(\n    np.abs(np.fft.fftshift(D.reshape(nfft, nfftk, nfftk), axes=1)[:20, :, nfftk // 2]),\n    cmap=\"seismic\",\n)\naxs[0][1].set_title(\"Fourier Transform\")\naxs[0][1].axis(\"tight\")\naxs[1][0].imshow(dinv[:, :, ny // 2], vmin=-20, vmax=20, cmap=\"seismic\")\naxs[1][0].set_title(\"Inverted\")\naxs[1][0].axis(\"tight\")\naxs[1][1].imshow(\n    d[:, :, ny // 2] - dinv[:, :, ny // 2], vmin=-20, vmax=20, cmap=\"seismic\"\n)\naxs[1][1].set_title(\"Error\")\naxs[1][1].axis(\"tight\")\nfig.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.6.12"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}