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)

../../_images/sphx_glr_dlf_design_001.png

Out:

   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.67412e-14
   > Spacing       : 0.07333333333
   > Shift         : -2.222222222
   > Base min/max  : 7.080680e-05 / 1.658545e+02

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

* 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)

../../_images/sphx_glr_dlf_design_002.png

Out:

   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.62718e-16
   > Spacing       : 0.06
   > Shift         : -1.333333333
   > Base min/max  : 6.533920e-04 / 1.063427e+02

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

* 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)

../../_images/sphx_glr_dlf_design_003.png

Out:

   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.9984e-16
   > Spacing       : 0.06533333333
   > Shift         : -1.483333333
   > Base min/max  : 3.299179e-04 / 1.560225e+02

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

* 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='kong_241_2007', **inp)
    key201 = empymod.dipole(htarg='key_201_2012', **inp)
    and801 = empymod.dipole(htarg='anderson_801_1982', **inp)
    test = empymod.dipole(htarg=filt, **inp)
    wer201 = empymod.dipole(htarg='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])

    plt.show()

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

Out:

* 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()

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

Estimated memory usage: 28 MB

Gallery generated by Sphinx-Gallery