Calculating a Digital Linear Filter#

This is an example for the add-on fdesign. The example is taken from the article Werthmüller et al., 2019. Have a look at the article repository on empymod/article-fdesign for many more examples.

Reference

Werthmüller, D., K. Key, and E. Slob, 2019, A tool for designing digital filters for the Hankel and Fourier transforms in potential, diffusive, and wavefield modeling: Geophysics, 84(2), F47-F56; DOI: 10.1190/geo2018-0069.1.

import empymod
import numpy as np
import matplotlib.pyplot as plt
plt.style.use('ggplot')

inp = {'r': np.logspace(0, 10, 1000),
       'r_def': (1, 1, 2),
       'n': 201,
       'name': 'test',
       'full_output': True,
       'fI': (empymod.fdesign.j0_1(5), empymod.fdesign.j1_1(5))}

1. Rough overview over a wide range#

filt1, out1 = empymod.fdesign.design(
        spacing=(0.01, 0.2, 10), shift=(-4, 0, 10), save=False, **inp)

Minimal recovered fields, Filter values of best filter

   brute fct calls : 1/100
   brute fct calls : 2/100
   brute fct calls : 3/100
   brute fct calls : 4/100
   brute fct calls : 5/100
   brute fct calls : 6/100
   brute fct calls : 7/100
   brute fct calls : 8/100
   brute fct calls : 9/100
   brute fct calls : 10/100
   brute fct calls : 11/100
   brute fct calls : 12/100
   brute fct calls : 13/100
   brute fct calls : 14/100
   brute fct calls : 15/100
   brute fct calls : 16/100
   brute fct calls : 17/100
   brute fct calls : 18/100
   brute fct calls : 19/100
   brute fct calls : 20/100
   brute fct calls : 21/100
   brute fct calls : 22/100
   brute fct calls : 23/100
   brute fct calls : 24/100
   brute fct calls : 25/100
   brute fct calls : 26/100
   brute fct calls : 27/100
   brute fct calls : 28/100
   brute fct calls : 30/100
   brute fct calls : 31/100
   brute fct calls : 32/100
   brute fct calls : 33/100
   brute fct calls : 34/100
   brute fct calls : 35/100
   brute fct calls : 36/100
   brute fct calls : 37/100
   brute fct calls : 38/100
   brute fct calls : 39/100
   brute fct calls : 40/100
   brute fct calls : 41/100
   brute fct calls : 42/100
   brute fct calls : 43/100
   brute fct calls : 44/100
   brute fct calls : 45/100
   brute fct calls : 46/100
   brute fct calls : 47/100
   brute fct calls : 48/100
   brute fct calls : 49/100
   brute fct calls : 50/100
   brute fct calls : 51/100
   brute fct calls : 52/100
   brute fct calls : 53/100
   brute fct calls : 54/100
   brute fct calls : 55/100
   brute fct calls : 56/100
   brute fct calls : 58/100
   brute fct calls : 59/100
   brute fct calls : 60/100
   brute fct calls : 61/100
   brute fct calls : 62/100
   brute fct calls : 63/100
   brute fct calls : 64/100
   brute fct calls : 65/100
   brute fct calls : 66/100
   brute fct calls : 67/100
   brute fct calls : 68/100
   brute fct calls : 69/100
   brute fct calls : 70/100
   brute fct calls : 71/100
   brute fct calls : 72/100
   brute fct calls : 73/100
   brute fct calls : 74/100
   brute fct calls : 75/100
   brute fct calls : 76/100
   brute fct calls : 77/100
   brute fct calls : 78/100
   brute fct calls : 79/100
   brute fct calls : 80/100
   brute fct calls : 81/100
   brute fct calls : 82/100
   brute fct calls : 83/100
   brute fct calls : 84/100
   brute fct calls : 85/100
   brute fct calls : 86/100
   brute fct calls : 87/100
   brute fct calls : 88/100
   brute fct calls : 89/100
   brute fct calls : 90/100
   brute fct calls : 91/100
   brute fct calls : 92/100
   brute fct calls : 93/100
   brute fct calls : 94/100
   brute fct calls : 95/100
   brute fct calls : 96/100
   brute fct calls : 97/100
   brute fct calls : 98/100
   brute fct calls : 99/100
   Filter length   : 201
   Best filter
   > Min field     : 1.67419e-14
   > Spacing       : 0.07333333333
   > Shift         : -2.222222222
   > Base min/max  : 7.080680e-05 / 1.658545e+02

:: empymod END; runtime = 0:00:01.729699 ::

* QC: Overview of brute-force inversion:

2. First focus#

filt2, out2 = empymod.fdesign.design(
        spacing=(0.04, 0.1, 10), shift=(-3, -0.5, 10), save=False, **inp)

   brute fct calls : 1/100
   brute fct calls : 2/100
   brute fct calls : 3/100
   brute fct calls : 4/100
   brute fct calls : 5/100
   brute fct calls : 6/100
   brute fct calls : 7/100
   brute fct calls : 8/100
   brute fct calls : 9/100
   brute fct calls : 10/100
   brute fct calls : 11/100
   brute fct calls : 12/100
   brute fct calls : 13/100
   brute fct calls : 14/100
   brute fct calls : 15/100
   brute fct calls : 16/100
   brute fct calls : 17/100
   brute fct calls : 18/100
   brute fct calls : 19/100
   brute fct calls : 20/100
   brute fct calls : 21/100
   brute fct calls : 22/100
   brute fct calls : 23/100
   brute fct calls : 24/100
   brute fct calls : 25/100
   brute fct calls : 26/100
   brute fct calls : 27/100
   brute fct calls : 28/100
   brute fct calls : 30/100
   brute fct calls : 31/100
   brute fct calls : 32/100
   brute fct calls : 33/100
   brute fct calls : 34/100
   brute fct calls : 35/100
   brute fct calls : 36/100
   brute fct calls : 37/100
   brute fct calls : 38/100
   brute fct calls : 39/100
   brute fct calls : 40/100
   brute fct calls : 41/100
   brute fct calls : 42/100
   brute fct calls : 43/100
   brute fct calls : 44/100
   brute fct calls : 45/100
   brute fct calls : 46/100
   brute fct calls : 47/100
   brute fct calls : 48/100
   brute fct calls : 49/100
   brute fct calls : 50/100
   brute fct calls : 51/100
   brute fct calls : 52/100
   brute fct calls : 53/100
   brute fct calls : 54/100
   brute fct calls : 55/100
   brute fct calls : 56/100
   brute fct calls : 58/100
   brute fct calls : 59/100
   brute fct calls : 60/100
   brute fct calls : 61/100
   brute fct calls : 62/100
   brute fct calls : 63/100
   brute fct calls : 64/100
   brute fct calls : 65/100
   brute fct calls : 66/100
   brute fct calls : 67/100
   brute fct calls : 68/100
   brute fct calls : 69/100
   brute fct calls : 70/100
   brute fct calls : 71/100
   brute fct calls : 72/100
   brute fct calls : 73/100
   brute fct calls : 74/100
   brute fct calls : 75/100
   brute fct calls : 76/100
   brute fct calls : 77/100
   brute fct calls : 78/100
   brute fct calls : 79/100
   brute fct calls : 80/100
   brute fct calls : 81/100
   brute fct calls : 82/100
   brute fct calls : 83/100
   brute fct calls : 84/100
   brute fct calls : 85/100
   brute fct calls : 86/100
   brute fct calls : 87/100
   brute fct calls : 88/100
   brute fct calls : 89/100
   brute fct calls : 90/100
   brute fct calls : 91/100
   brute fct calls : 92/100
   brute fct calls : 93/100
   brute fct calls : 94/100
   brute fct calls : 95/100
   brute fct calls : 96/100
   brute fct calls : 97/100
   brute fct calls : 98/100
   brute fct calls : 99/100
   Filter length   : 201
   Best filter
   > Min field     : 3.65947e-16
   > Spacing       : 0.06
   > Shift         : -1.333333333
   > Base min/max  : 6.533920e-04 / 1.063427e+02

:: empymod END; runtime = 0:00:01.732921 ::

* QC: Overview of brute-force inversion:

3. Final focus#

filt, out = empymod.fdesign.design(
        spacing=(0.047, 0.08, 10), shift=(-2.4, -0.75, 10), finish=False,
        save=False, **inp)

   brute fct calls : 1/100
   brute fct calls : 2/100
   brute fct calls : 3/100
   brute fct calls : 4/100
   brute fct calls : 5/100
   brute fct calls : 6/100
   brute fct calls : 7/100
   brute fct calls : 8/100
   brute fct calls : 9/100
   brute fct calls : 10/100
   brute fct calls : 11/100
   brute fct calls : 12/100
   brute fct calls : 13/100
   brute fct calls : 14/100
   brute fct calls : 15/100
   brute fct calls : 16/100
   brute fct calls : 17/100
   brute fct calls : 18/100
   brute fct calls : 19/100
   brute fct calls : 20/100
   brute fct calls : 21/100
   brute fct calls : 22/100
   brute fct calls : 23/100
   brute fct calls : 24/100
   brute fct calls : 25/100
   brute fct calls : 26/100
   brute fct calls : 27/100
   brute fct calls : 28/100
   brute fct calls : 30/100
   brute fct calls : 31/100
   brute fct calls : 32/100
   brute fct calls : 33/100
   brute fct calls : 34/100
   brute fct calls : 35/100
   brute fct calls : 36/100
   brute fct calls : 37/100
   brute fct calls : 38/100
   brute fct calls : 39/100
   brute fct calls : 40/100
   brute fct calls : 41/100
   brute fct calls : 42/100
   brute fct calls : 43/100
   brute fct calls : 44/100
   brute fct calls : 45/100
   brute fct calls : 46/100
   brute fct calls : 47/100
   brute fct calls : 48/100
   brute fct calls : 49/100
   brute fct calls : 50/100
   brute fct calls : 51/100
   brute fct calls : 52/100
   brute fct calls : 53/100
   brute fct calls : 54/100
   brute fct calls : 55/100
   brute fct calls : 56/100
   brute fct calls : 58/100
   brute fct calls : 59/100
   brute fct calls : 60/100
   brute fct calls : 61/100
   brute fct calls : 62/100
   brute fct calls : 63/100
   brute fct calls : 64/100
   brute fct calls : 65/100
   brute fct calls : 66/100
   brute fct calls : 67/100
   brute fct calls : 68/100
   brute fct calls : 69/100
   brute fct calls : 70/100
   brute fct calls : 71/100
   brute fct calls : 72/100
   brute fct calls : 73/100
   brute fct calls : 74/100
   brute fct calls : 75/100
   brute fct calls : 76/100
   brute fct calls : 77/100
   brute fct calls : 78/100
   brute fct calls : 79/100
   brute fct calls : 80/100
   brute fct calls : 81/100
   brute fct calls : 82/100
   brute fct calls : 83/100
   brute fct calls : 84/100
   brute fct calls : 85/100
   brute fct calls : 86/100
   brute fct calls : 87/100
   brute fct calls : 88/100
   brute fct calls : 89/100
   brute fct calls : 90/100
   brute fct calls : 91/100
   brute fct calls : 92/100
   brute fct calls : 93/100
   brute fct calls : 94/100
   brute fct calls : 95/100
   brute fct calls : 96/100
   brute fct calls : 97/100
   brute fct calls : 98/100
   brute fct calls : 99/100
   Filter length   : 201
   Best filter
   > Min field     : 3.61036e-16
   > Spacing       : 0.06166666667
   > Shift         : -1.3
   > Base min/max  : 5.718312e-04 / 1.298872e+02

:: empymod END; runtime = 0:00:01.694686 ::

* QC: Overview of brute-force inversion:

To reproduce exactly the same filter as wer_201_2018:#

filt_orig, out_orig = fdesign.load_filter('wer201', True)
fdesign.plot_result(filt_orig, out_orig)
filt_orig, out_orig = fdesign.design(
        spacing=out_orig[0][0], shift=out_orig[0][1], **inp)

Plot the result#

Plot function#

def plotresult(depth, res, zsrc, zrec):
    x = np.arange(1, 101)*200
    inp = {
        'src': [0, 0, depth[1]-zsrc],
        'rec': [x, x*0, depth[1]-zrec],
        'depth': depth,
        'res': res,
        'ab': 11,
        'freqtime': 1,
        'verb': 1,
    }

    kong241 = empymod.dipole(htarg={'dlf': 'kong_241_2007'}, **inp)
    key201 = empymod.dipole(htarg={'dlf': 'key_201_2012'}, **inp)
    and801 = empymod.dipole(htarg={'dlf': 'anderson_801_1982'}, **inp)
    test = empymod.dipole(htarg={'dlf': filt}, **inp)
    wer201 = empymod.dipole(htarg={'dlf': 'wer_201_2018'}, **inp)
    qwe = empymod.dipole(ht='qwe', **inp)

    plt.figure(figsize=(8, 3.5))
    plt.subplot(121)
    plt.semilogy(x, qwe.amp(), c='0.5', label='QWE')
    plt.semilogy(x, kong241.amp(), 'k--', label='Kong241')
    plt.semilogy(x, key201.amp(), 'k:', label='Key201')
    plt.semilogy(x, and801.amp(), 'k-.', label='And801')
    plt.semilogy(x, test.amp(), 'r', label='This filter')
    plt.semilogy(x, wer201.amp(), 'b', label='Wer201')
    plt.legend()
    plt.xticks([0, 5e3, 10e3, 15e3, 20e3])
    plt.xlim([0, 20e3])

    plt.subplot(122)
    plt.semilogy(x, np.abs((kong241-qwe)/qwe), 'k--', label='Kong241')
    plt.semilogy(x, np.abs((key201-qwe)/qwe), 'k:', label='Key201')
    plt.semilogy(x, np.abs((and801-qwe)/qwe), 'k-.', label='And801')
    plt.semilogy(x, np.abs((test-qwe)/qwe), 'r', label='This filter')
    plt.semilogy(x, np.abs((wer201-qwe)/qwe), 'b', label='Wer201')
    plt.legend()
    plt.xticks([0, 5e3, 10e3, 15e3, 20e3])
    plt.xlim([0, 20e3])
    plt.ylim([1e-14, 1])

Plot the individual models#

plotresult([-1e50, 2000], [2e14, 1/3.2, 1], 50, 0)  # KONG model
plotresult([0, 1000, 2000, 2100], [1/1e-12, 1/3.3, 1, 100, 1], 10, 0)  # KEY m.
plotresult([0, 1, 1000, 1100], [2e14, 10, 10, 500, 10], 0.5, 0.2)  # LAND model

* WARNING :: Hankel-quadrature did not converge at least once;
             => desired `atol` and `rtol` might not be achieved.
* WARNING :: Hankel-quadrature did not converge at least once;
             => desired `atol` and `rtol` might not be achieved.

empymod.Report()

Fri Oct 25 09:12:19 2024 UTC
OS	Linux (Ubuntu 22.04)	CPU(s)	2	Machine	x86_64
Architecture	64bit	RAM	7.5 GiB	Environment	Python
File system	ext4
Python 3.11.9 (main, Jun 18 2024, 09:40:25) [GCC 11.4.0]
numpy	2.0.2	scipy	1.14.1	numba	0.60.0
empymod	2.4.0	libdlf	0.3.0	IPython	8.28.0
matplotlib	3.9.2

Total running time of the script: (0 minutes 20.383 seconds)

Estimated memory usage: 188 MB

Gallery generated by Sphinx-Gallery