{
  "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 numpy as np\nimport matplotlib.pyplot as plt\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,\n                                    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)],\n            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,\n                                    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)],\n            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,\n                                    nfft=nfft, sampling=dt)\nD = FFTop*d.flatten()\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(dims=(nt, nx),\n                                      nffts=(nfft, nfft), sampling=(dt, dx))\nD = FFTop*d.flatten()\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(np.abs(np.fft.fftshift(D.reshape(nfft, nfft),\n                                        axes=1)[:200, :]), cmap='seismic')\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(dims=(nt, nx, ny),\n                                      nffts=(nfft, nfftk, nfftk),\n                                      sampling=(dt, dx, dy))\nD = FFTop*d.flatten()\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(np.abs(np.fft.fftshift(D.reshape(nfft, nfftk, nfftk),\n                                        axes=1)[:20, :, nfftk//2]),\n                 cmap='seismic')\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(d[:, :, ny//2]-dinv[:, :, ny//2],\n                 vmin=-20, vmax=20, cmap='seismic')\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
}