Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
Implementation of Kaiser windowed sinc functions after Graham Hicks's 2002 Geophysics paper
HC_KAISER = {
1: 1.24,
2: 2.94,
3: 4.53,
4: 6.31,
5: 7.91,
6: 9.42,
7: 10.95,
8: 12.53,
9: 14.09,
10: 14.18,
}
def KaiserWindowedSinc(ireg, offset):
'''
Finds 2D source terms to approximate a band-limited point source, based on
Hicks, Graham J. (2002) Arbitrary source and receiver positioning in finite-difference
schemes using Kaiser windowed sinc functions. Geophysics (67) 1, 156-166.
KaiserWindowedSince(ireg, offset) --> 2D ndarray of size (2*ireg+1, 2*ireg+1)
Input offset is the 2D offsets in fractional gridpoints between the source location and
the nearest node on the modelling grid.
'''
from scipy.special import i0 as bessi0
import warnings
ireg = int(ireg)
try:
b = HC_KAISER.get(ireg)
except KeyError:
print('Kaiser windowed sinc function not implemented for half-width of %d!'%(ireg,))
raise
freg = 2*ireg+1
xOffset, zOffset = offset
# Grid from 0 to freg-1
Zi, Xi = numpy.mgrid[:freg,:freg]
# Distances from source point
dZi = (zOffset + ireg - Zi)
dXi = (xOffset + ireg - Xi)
# Taper terms for decay function
with warnings.catch_warnings():
warnings.simplefilter('ignore')
tZi = numpy.nan_to_num(numpy.sqrt(1 - (dZi / ireg)**2))
tXi = numpy.nan_to_num(numpy.sqrt(1 - (dXi / ireg)**2))
tZi[tZi == inf] = 0
tXi[tXi == inf] = 0
# Actual tapers for Kaiser window
taperZ = bessi0(b*tZi) / bessi0(b)
taperX = bessi0(b*tXi) / bessi0(b)
# Windowed sinc responses in Z and X
responseZ = numpy.sinc(dZi) * taperZ
responseX = numpy.sinc(dXi) * taperX
# Combined 2D source response
result = responseX * responseZ
return result
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment