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

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

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

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

* 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     : 1.98558e-16
> Spacing       : 0.06166666667
> Shift         : -1.3
> Base min/max  : 5.718312e-04 / 1.298872e+02

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

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

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
```
```* 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()
```
 Sat Oct 15 19:19:46 2022 UTC OS Linux CPU(s) 2 Machine x86_64 Architecture 64bit RAM 7.6 GiB Environment Python File system ext4 Python 3.8.6 (default, Oct 19 2020, 15:10:29) [GCC 7.5.0] numpy 1.23.4 scipy 1.9.2 numba 0.56.3 empymod 2.2.1 IPython 8.5.0 matplotlib 3.6.1

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

Estimated memory usage: 10 MB

Gallery generated by Sphinx-Gallery