{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib notebook"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# Digital Linear Filters\n\nGraphical explanation of the differences between standard DLF, lagged\nconvolution DLF, and splined DLF.\n\nThe comments here apply generally to the digital linear filter method. Having\n``empymod`` in mind, they are particularly meant for the Hankel\n(Bessel-Fourier) transform from the wavenumber-frequency domain ($k-f$)\nto the space-frequency domain ($x-f$), and for the Fourier transform from\nthe space-frequency domain ($x-f$) to the space-time domain\n($x-t$).\n\n## 1. Introduction\n\n*This introduction is taken from Werthm\u00fcller et al. (2018), which can be found\nin the repo* [empymod/article-fdesign](https://github.com/emsig/article-fdesign).\n\nIn electromagnetics we often have to evaluate integrals of the form\n\n\\begin{align}F(r) = \\int^\\infty_0 f(l)K(l r)\\,dl \\ ,\\end{align}\n\nwhere $l$ and $r$ denote input and output evaluation values,\nrespectively, and $K$ is the kernel function. In the specific case of the\nHankel transform $l$ corresponds to wavenumber, $r$ to offset, and\n$K$ to Bessel functions; in the case of the Fourier transform $l$\ncorresponds to frequency, $r$ to time, and $K$ to sine or cosine\nfunctions. In both cases it is an infinite integral which numerical integration\nis very time-consuming because of the slow decay of the kernel function and its\noscillatory behaviour.\n\nBy substituting $r = e^x$ and $l = e^{-y}$ we get\n\n\\begin{align}e^x F(e^x) = \\int^\\infty_{-\\infty} f(e^{-y})K(e^{x-y})e^{x-y}\\,dy\\ .\\end{align}\n\nThis can be re-written as a convolution integral and be approximated for an\n$N$-point filter by\n\n\\begin{align}F(r) \\approx \\sum^N_{n=1} \\frac{f(b_n/r) h_n}{r}\\ ,\\end{align}\n\nwhere $h$ is the digital linear filter, and the logarithmically spaced\nfilter abscissae is a function of the spacing $\\Delta$ and the shift\n$\\delta$,\n\n\\begin{align}b_n = \\exp\\left\\{\\Delta(-N/2+n) + \\delta\\right\\} \\ .\\end{align}\n\nFrom the penultimate equation it can be seen that the filter method requires\n$N$ evaluations at each $r$. For example, to compute the frequency\ndomain result for 100 offsets with a 201 pt filter requires 20,100 evaluations\nin the wavenumber domain. This is why the DLF often uses interpolation to\nminimize the required evaluations, either for $F(r)$ in what is referred\nto as *lagged convolution DLF*, or for $f(l)$, which we call here\n*splined DLF*.\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import empymod\nimport numpy as np\nimport matplotlib.pyplot as plt\nfrom copy import deepcopy as dc\nplt.style.use('ggplot')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## 2. How the DLF works\n\n### Design a very short, digital linear filter\n\nFor this we use ``empymod.fdesign``. This is outside the scope of this\nnotebook. If you are interested have a look at the [article-fdesign](https://github.com/emsig/article-fdesign)-repo for more information\nregarding the design of digital linear filters.\n\nWe design a 5pt filter using the theoretical transform pair\n\n\\begin{align}\\int^\\infty_0 l \\exp\\left(-l^2\\right) J_0(lr) dl =\n  \\exp\\left(\\frac{-r^2}{4}\\right)/2 \\ .\\end{align}\n\nA 5 pt filter is very short for this problem, so we expect a considerable\nerror level. In designing the filter we force the filter to be better than a\nrelative error of 5 %.\n\nIf you want to play around with this example I recommend to set ``verb=2``\nand ``plot=2`` to get some feedback from the minimization procedure.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "filt = empymod.fdesign.design(\n        n=5,                        # 5 point filter\n        spacing=(0.55, 0.65, 101),\n        shift=(0.6, 0.7, 101),\n        fI=empymod.fdesign.j0_1(),\n        r=np.logspace(0, 1, 100),\n        r_def=(1, 1, 10),\n        error=0.05,   # 5 % error level. If you set this too low you will\n        verb=1,       # #                not find a filter with only 5 points.\n        plot=0,\n)\n\nprint('Filter base    ::', filt.base)\nprint('Filter weights ::', filt.j0)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Now we carry out the DLF and check how good it is.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Desired x-f-domain points (output domain)\nx = np.array([0.5, 1, 2, 3])\n\n# Required k-f-domain points (input domain)\nk = filt.base/x[:, None]\n\n# Get the theoretical transform pair\ntp = empymod.fdesign.j0_1()\n\n# Compute the value at the five required wavenumbers\nk_val = tp.lhs(k)\n\n# Weight the values and sum them up\nx_val_filt = np.dot(k_val, filt.j0)/x\n\n# Compute the theoretical value for comparison\nx_val_theo = tp.rhs(x)\n\n# Compute relative error\nprint('A DLF for this problem with only a 5 pt filter is difficult. We used')\nprint('an error-limit of 0.05 in the filter design, so we expect the result')\nprint('to have a relative error of less than 5 %.\\n')\nprint('Theoretical value   ::', '; '.join([f'{i:G}' for i in x_val_theo]))\nprint('DLF value           ::', '; '.join([f'{i:G}' for i in x_val_filt]))\nrelerror = np.abs((x_val_theo-x_val_filt)/x_val_theo)\nprint('Rel. error 5 pt DLF ::', ' %   ; '.join(\n      [f'{i:G}' for i in np.round(relerror*100, 1)]), '%')\n\n# Figure\nplt.figure(figsize=(10, 4))\nplt.suptitle(r'DLF example for $J_0$ Hankel transform using 5 pt filter',\n             y=1.05)\n\n# x-axis values for the theoretical plots\nx_k = np.logspace(-1, np.log10(13))\nx_x = np.logspace(-0.5, np.log10(4))\n\n# k-f-domain\nplt.subplot(121)\nplt.title(r'$k-f$-domain')\nplt.loglog(x_k, tp.lhs(x_k), label='Theoretical')\nfor i, val in enumerate(x):\n    plt.loglog(k[i, :], k_val[i, ], 'o', label='5 pt DLF input x ='+str(val))\nplt.legend()\nplt.xlabel(r'$k$')\nplt.xlim([x_k.min(), x_k.max()])\n\n# x-f-domain\nplt.subplot(122)\nplt.title(r'$x-f$-domain')\nplt.loglog(x_x, tp.rhs(x_x), label='Theoretical')\nfor i, val in enumerate(x):\n    plt.loglog(val, x_val_filt[i], 'o', label='5 pt DLF output x = '+str(val))\nplt.legend()\nplt.xlabel(r'$x$')\nplt.xlim([x_x.min(), x_x.max()])\n\nplt.tight_layout()\nplt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## 3. Difference between standard, lagged convolution, and splined DLF\n\nFilter weights and the actual DLF are ignored in the explanation, we only\nlook at the required data points in the input domain given our desired points\nin the output domain.\n\n**General parameters**\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Wanted points in the output domain (x-f or x-t domain)\nd_out = np.array([0.2, 0.7, 5, 25, 100])\n\n# Made-up filter base for explanatory purpose\nbase = np.array([1e-2, 1e-1, 1e0, 1e1, 1e2])"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## 3.1 Standard DLF\n\nFor each point in the output domain you have to compute $n$ points in\nthe input domain, where $n$ is the filter length.\n\n**Implementation in ``empymod``**\n\nThis is the most precise one, as no interpolation is used, but generally the\nslowest one. It is the default method for the Hankel transform.\n\nFor the **Hankel transform**, use these parameters in ``empymod.dipole`` or\n``empymod.bipole`` (from version ``v1.6.0`` onwards):\n\n|   ht = 'dlf'                  # Default\n|   htarg = {'pts_per_dec': 0}  # Default\n\nFor the **Fourier transform**, use these parameters in ``empymod.dipole`` or\n``empymod.bipole`` (from version ``v1.6.0`` onwards):\n\n|   ft = 'sin' or 'cos'            # Default is 'sin'\n|   ftarg = {'pts_per_dec': 0}\n\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Required points in the input domain (k-f or x-f domain)\nd_in = base/d_out[:, None]\n\n# Print information\nprint('Points in output domain  ::', d_out.size)\nprint('Filter length            ::', base.size)\nprint('Req. pts in input domain ::', d_in.size)\n\n# Figure\nplt.figure()\nplt.title('Standard DLF')\nplt.hlines(1, 1e-5, 1e5)\nplt.hlines(0, 1e-5, 1e5)\n\n# Print scheme\nfor i, val in enumerate(d_out):\n    for ii, ival in enumerate(d_in[i, :]):\n        plt.plot(ival, 0, 'C'+str(i)+'x')\n        plt.plot([ival, val], [0, 1], 'C'+str(i))\n    plt.plot(val, 1, 'C'+str(i)+'o')\n\nplt.xscale('log')\nplt.xlim([5e-5, 5e3])\n\n# Annotations\nplt.text(2e3, 0.5, '--- DLF --->', rotation=90, va='center')\nplt.text(1e-4, 0.5, '--- DLF --->', rotation=90, va='center')\n\nplt.yticks([0, 1], (r'$k-f$', r'$x-f$'))\nplt.ylabel('Hankel transform')\nplt.ylim([-0.1, 1.1])\n\nplt.gca().twinx()\nplt.yticks([0, 1], (r'$x-f$', r'$x-t$'))\nplt.ylabel('Fourier transform')\nplt.ylim([-0.1, 1.1])\n\nplt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## 3.2 Lagged Convolution DLF\n\nThe spacing of the filter base is used to get from minimum to maximum\nrequired input-domain point ($k$ in the case of the Hankel transform,\n$f$ in the case of the Fourier transform); for each complete set the\nDLF is executed to compute the output-domain response ($f$ in the case\nof the Hankel transform, $t$ in the case of the Fourier transform), and\ninterpolation is carried out in the output-domain.\n\n**Implementation in ``empymod``**\n\nThis is usually the fastest option, and generally still more than\nsufficiently precise. It is the default method for the Fourier transform.\n\nFor the **Hankel transform**, use these parameters in ``empymod.dipole`` or\n``empymod.bipole`` (from version ``v1.6.0`` onwards):\n\n|   ht = 'dlf'                      # Default\n|   htarg = {'pts_per_dec': int<0}\n\nFor the **Fourier transform**, use these parameters in ``empymod.dipole`` or\n``empymod.bipole`` (from version ``v1.6.0`` onwards):\n\n|   ft = 'sin' or 'cos'             # Default is 'sin'\n|   ftarg = {'pts_per_dec': int<0}  # Default\n\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Required points in the k-domain\nd_in2 = np.array([1e-4, 1e-3, 1e-2, 1e-1, 1e0, 1e1, 1e2, 1e3])\n\n# Intermediat values in the f-domain\nd_out2 = np.array([1e-1, 1e0, 1e1, 1e2])\n\n# Print information\nprint('Points in output-domain  ::', d_out.size)\nprint('Filter length            ::', base.size)\nprint('Req. pts in input-domain ::', d_in2.size)\n\n# Figure\nplt.figure()\nplt.title('Lagged convolution DLF')\nplt.hlines(1, 1e-5, 1e5)\nplt.hlines(0, 1e-5, 1e5)\nplt.hlines(0.5, 1e-5, 1e5)\n\n# Plot scheme\nfor i, val in enumerate(d_out2):\n    for ii in range(base.size):\n        plt.plot([d_in2[-1-ii-i], val], [0, 0.5], str(0.7-0.2*i))\n\nfor iii, val2 in enumerate(d_out):\n    plt.plot(val2, 1, 'C'+str(iii)+'o')\n    plt.plot([val2, val2], [0.5, 1], 'C'+str(iii))\n\nplt.plot(d_in2, d_in2*0, 'ko', ms=8)\nplt.plot(d_out2, d_out2*0+0.5, 'ks', ms=8)\nplt.xscale('log')\nplt.xlim([5e-5, 5e3])\n\n# Annotations\nplt.text(2e3, 0.75, '- Interpolation ->', rotation=90, va='center')\nplt.text(2e3, 0.25, '--- DLF --->', rotation=90, va='center')\nplt.text(1e-4, 0.75, '- Interpolation ->', rotation=90, va='center')\nplt.text(1e-4, 0.25, '--- DLF --->', rotation=90, va='center')\n\nplt.yticks([0, 0.5, 1], (r'$k-f$', r'$x-f$', r'$x-f$'))\nplt.ylabel('Hankel transform')\nplt.ylim([-0.1, 1.1])\n\nplt.gca().twinx()\nplt.yticks([0, 0.5, 1], (r'$x-f$', r'$x-t$', r'$x-t$'))\nplt.ylabel('Fourier transform')\nplt.ylim([-0.1, 1.1])\n\nplt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## 3.3 Splined DLF\n\nIn the splined DLF $m$ points per decade are used from minimum to\nmaximum required input-domain point ($k$ in the case of the Hankel\ntransform, $f$ in the case of the Fourier transform); then the required\ninput-domain responses are interpolated in the input domain, and the DLF is\nexececuted subsequently.\n\n**Implementation in ``empymod``**\n\nThis option can, at times, yield more precise results than the lagged\nconvolution DLF, while being slower than the lagged convolution DLF but\nfaster than the standard DLF. However, you have to carefully choose (or\nbetter, test) the parameter ``pts_per_dec``.\n\nFor the **Hankel transform**, use these parameters in ``empymod.dipole`` or\n``empymod.bipole`` (from version ``v1.6.0`` onwards):\n\n|   ht = 'dlf'                     # Default\n|   htarg = {'pts_per_dec': int>0}\n\nFor the **Fourier transform**, use these parameters in ``empymod.dipole`` or\n``empymod.bipole`` (from version ``v1.6.0`` onwards):\n\n|   ft = 'sin' or 'cos'            # Default is 'sin'\n|   ftarg = {'pts_per_dec': int>0}\n\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Required points in the k-domain\nd_in_min = np.log10(d_in).min()\nd_in_max = np.ceil(np.log10(d_in).max())\npts_per_dec = 3\nd_in2 = np.logspace(d_in_min, d_in_max, int((d_in_max-d_in_min)*pts_per_dec+1))\n\n# Print information\nprint('Points in input-domain    ::', d_out.size)\nprint('Filter length             ::', base.size)\nprint('Points per decade         ::', pts_per_dec)\nprint('Req. pts in output-domain ::', d_in2.size)\n\n# Figure\nplt.figure()\nplt.title('Splined DLF')\nplt.hlines(1, 1e-5, 1e5)\nplt.hlines(0.5, 1e-5, 1e5)\nplt.hlines(0, 1e-5, 1e5)\n\n# Plot scheme\nfor i, val in enumerate(d_out):\n    for ii, ival in enumerate(d_in[i, :]):\n        plt.plot(ival, 0.5, 'C'+str(i)+'x')\n        plt.plot([ival, ival], [0, 0.5], str(0.6-0.1*i))\n        plt.plot([ival, val], [0.5, 1], 'C'+str(i))\n    plt.plot(val, 1, 'C'+str(i)+'o')\n\nplt.plot(d_in2, d_in2*0, 'ks')\nplt.xscale('log')\nplt.xlim([5e-5, 5e3])\n\n# Annotations\nplt.text(2e3, 0.25, '- Interpolation ->', rotation=90, va='center')\nplt.text(2e3, 0.75, '--- DLF --->', rotation=90, va='center')\nplt.text(1.5e-4, 0.25, '- Interpolation ->', rotation=90, va='center')\nplt.text(1.5e-4, 0.75, '--- DLF --->', rotation=90, va='center')\n\nplt.yticks([0, 0.5, 1], (r'$k-f$', r'$k-f$', r'$x-f$'))\nplt.ylabel('Hankel transform')\nplt.ylim([-0.1, 1.1])\n\nplt.gca().twinx()\nplt.yticks([0, 0.5, 1], (r'$x-f$', r'$x-f$', r'$x-t$'))\nplt.ylabel('Fourier transform')\nplt.ylim([-0.1, 1.1])\n\nplt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## 4. Example for the Hankel transform\n\nThe following is an example for the Hankel transform. **Be aware that the\nactual differences in time and accuracy depend highly on the model.** If time\nor accuracy is a critical issue in your computation I suggest to run some\npreliminary tests. It also depends heavily if you have many offsets, or many\nfrequencies, or many layers, as one method might be better for many\nfrequencies but few offsets, but the other method might be better for many\noffsets but few frequencies.\n\nAs general rules we can state that\n\n- the longer the used filter is, or\n- the more offsets you have\n\nthe higher is the time gain you get by using the lagged convolution or\nsplined version of the DLF.\n\nHere we compare the analytical halfspace solution to the numerical result,\nusing the standard DLF, the lagged convolution DLF, and the splined DLF. Note\nthe oscillating behaviour of the error of the lagged convolution and the\nsplined versions, which comes from the interpolation and is not present in\nthe standard version.\n\n**Define model, compute analytical solution**\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "x = (np.arange(1, 1001))*10\nparams = {\n    'src': [0, 0, 150],\n    'rec': [x, x*0, 200],\n    'depth': 0,\n    'res': [2e14, 1],\n    'freqtime': 1,\n    'ab': 11,\n    'aniso': [1, 2],\n    'xdirect': False,\n    'verb': 0,\n}\n\n# Used Hankel filter\nhfilt = empymod.filters.key_201_2009()\n\n# Compute analytical solution\nresp = empymod.analytical(\n        params['src'], params['rec'], params['res'][1], params['freqtime'],\n        solution='dhs', aniso=params['aniso'][1], ab=params['ab'],\n        verb=params['verb']\n)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "### Compute numerically the model using different Hankel options\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "standard = empymod.dipole(**params, htarg={'dlf': hfilt, 'pts_per_dec': 0})\nlaggedco = empymod.dipole(**params, htarg={'dlf': hfilt, 'pts_per_dec': -1})\nspline10 = empymod.dipole(**params, htarg={'dlf': hfilt, 'pts_per_dec': 10})\nspline30 = empymod.dipole(**params, htarg={'dlf': hfilt, 'pts_per_dec': 30})\nsplin100 = empymod.dipole(**params, htarg={'dlf': hfilt, 'pts_per_dec': 100})"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "### Results\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "plt.figure(figsize=(10, 4))\nplt.suptitle('Hankel transform example; frequency = ' +\n             str(params['freqtime'])+' Hz', y=1.05, fontsize=15)\n\nplt.subplot(121)\nplt.title('Amplitude (V/m)')\nplt.semilogy(x/1000, np.abs(resp), 'k', label='Analytical')\nplt.semilogy(x/1000, np.abs(standard), label='Standard')\nplt.semilogy(x/1000, np.abs(laggedco), label='Lagged')\nplt.semilogy(x/1000, np.abs(spline10), label='Splined 10/dec')\nplt.semilogy(x/1000, np.abs(spline30), label='Splined 30/dec')\nplt.semilogy(x/1000, np.abs(splin100), label='Splined 100/dec')\nplt.xlabel('Offset (km)')\nplt.legend()\n\nplt.subplot(122)\nplt.title('Relative Error (-); Frequency = '+str(params['freqtime'])+' Hz.')\nplt.semilogy(x/1000, np.abs((standard-resp)/resp), label='Standard')\nplt.semilogy(x/1000, np.abs((laggedco-resp)/resp), label='Lagged')\nplt.semilogy(x/1000, np.abs((spline10-resp)/resp), label='Splined 10/dec')\nplt.semilogy(x/1000, np.abs((spline30-resp)/resp), label='Splined 30/dec')\nplt.semilogy(x/1000, np.abs((splin100-resp)/resp), label='Splined 100/dec')\nplt.xlabel('Offset (km)')\n\nplt.tight_layout()\nplt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Runtimes and number of required wavenumbers for each method:\n\n+-------------------+-----------+------------------+\n|Hankel DLF Method  | Time (ms) | # of wavenumbers |\n+===================+===========+==================+\n|Standard           | 169       | 201000           |\n+-------------------+-----------+------------------+\n|Lagged Convolution |   4       | 295              |\n+-------------------+-----------+------------------+\n|Splined 10/dec     |  95       | 96               |\n+-------------------+-----------+------------------+\n|Splined 30/dec     |  98       | 284              |\n+-------------------+-----------+------------------+\n|Splined 100/dec    | 111       | 944              |\n+-------------------+-----------+------------------+\n\nSo the lagged convolution has a relative error between roughly 1e-6 to 1e-4,\nhence 0.0001 % to 0.01 %, which is more then enough for real-world\napplications.\n\nIf you want to measure the runtime on your machine set ``params['verb'] =\n2``.\n\nNote: If you use the splined version with about 500 samples per decade you\nget about the same accuracy as in the standard version. However, you get also\nabout the same runtime.\n\n## 5. Example for the Fourier transform\n\nThe same now for the Fourier transform. Obviously, when we use the Fourier\ntransform we also use the Hankel transform. However, in this example we use\nonly one offset. If there is only one offset then the lagged convolution or\nsplined DLF for the Hankel transform do not make sense, and the standard is\nthe fastest. So we use the standard Hankel DLF here. If you have many offsets\nthen that would be different.\n\nAs general rules we can state that\n\n- the longer the used filter is, or\n- the more times you have\n\nthe higher is the time gain you get by using the lagged convolution or\nsplined version of the DLF.\n\n### Define model, compute analytical solution\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "t = np.logspace(0, 2, 100)\nxt = 2000\n\ntparam = dc(params)\ntparam['rec'] = [xt, 0, 200]\ntparam['freqtime'] = t\ntparam['signal'] = 0  # Impulse response\n\n# Used Fourier filter\nffilt = empymod.filters.key_81_CosSin_2009()\n\n# Compute analytical solution\ntresp = empymod.analytical(\n        tparam['src'], tparam['rec'], tparam['res'][1], tparam['freqtime'],\n        signal=tparam['signal'], solution='dhs', aniso=tparam['aniso'][1],\n        ab=tparam['ab'], verb=tparam['verb']\n)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "### Compute numerically the model using different Fourier options\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "htarg = {'dlf': hfilt}\n\ntstandard = empymod.dipole(\n        **tparam, htarg=htarg, ftarg={'dlf': ffilt, 'pts_per_dec': 0})\ntlaggedco = empymod.dipole(\n        **tparam, htarg=htarg, ftarg={'dlf': ffilt, 'pts_per_dec': -1})\ntsplined4 = empymod.dipole(\n        **tparam, htarg=htarg, ftarg={'dlf': ffilt, 'pts_per_dec': 4})\ntspline10 = empymod.dipole(\n        **tparam, htarg=htarg, ftarg={'dlf': ffilt, 'pts_per_dec': 10})"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "### Results\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "plt.figure(figsize=(10, 4))\nplt.suptitle('Fourier transform example: impulse response; offset = ' +\n             str(xt/1000)+' km', y=1.05, fontsize=15)\n\nplt.subplot(121)\nplt.title('Amplitude (V/[m.s])')\nplt.semilogy(t, np.abs(tresp), 'k', label='Analytical')\nplt.semilogy(t, np.abs(tstandard), label='Standard')\nplt.semilogy(t, np.abs(tlaggedco), label='Lagged')\nplt.semilogy(t, np.abs(tsplined4), label='Splined 4/dec')\nplt.semilogy(t, np.abs(tspline10), label='Splined 10/dec')\nplt.xlabel('Time (s)')\nplt.legend()\n\nplt.subplot(122)\nplt.title('Relative Error (-)')\nplt.semilogy(t, np.abs((tstandard-tresp)/tresp), label='Standard')\nplt.semilogy(t, np.abs((tlaggedco-tresp)/tresp), label='Lagged')\nplt.semilogy(t, np.abs((tsplined4-tresp)/tresp), label='Splined 4/dec')\nplt.semilogy(t, np.abs((tspline10-tresp)/tresp), label='Splined 10/dec')\nplt.xlabel('Time (s)')\n\nplt.tight_layout()\nplt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Runtimes and number of required frequencies for each method:\n\n+-------------------+-----------+------------------+\n|Fourier DLF Method | Time (ms) | # of frequencies |\n+===================+===========+==================+\n|Standard           | 1442      | 8100             |\n+-------------------+-----------+------------------+\n|Lagged Convolution |   17      |  105             |\n+-------------------+-----------+------------------+\n|Splined 4/dec      |   14      |   37             |\n+-------------------+-----------+------------------+\n|Splined 10/dec     |   32      |   91             |\n+-------------------+-----------+------------------+\n\nAll methods require 201 wavenumbers (1 offset, filter length is 201).\n\nSo the lagged convolution has a relative error of roughly 1e-5, hence 0.001\n%, which is more then enough for real-world applications.\n\nIf you want to measure the runtime on your machine set ``tparam['verb'] =\n2``.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "empymod.Report()"
      ]
    }
  ],
  "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.8.6"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}