Search code examples
numpymatlabsignal-processingwindowing

Range of Beta in Kaiser Window

Why is the biggest value for Beta in the Kaiser Window function 709? I get errors when I use values of 710 or more (in Matlab and NumPy). Seems like a stupid question but I'm still not able to find an answer... Thanks in advance!

Kaiser Window definition in NumPy

This is the Kaiser window with a length of M=64 and Beta = 710. Kaiser Window M=64 and Beta=710

Solution

  • For large x, i0(x) is approximately exp(x), and exp(x) overflows at x near 709.79.

    If you are programming in Python and you don't mind the dependency on SciPy, you can use the function scipy.special.i0e to implement the function in a way that avoids the overflow:

    In [47]: from scipy.special import i0e

In [48]: def i0_ratio(x, y):
    ...:     return i0e(x)*np.exp(x - y)/i0e(y)
    ...:

    Verify that the function returns the same value as np.i0(x)/np.i0(y):

    In [49]: np.i0(3)/np.i0(4)
Out[49]: 0.4318550956673735

In [50]: i0_ratio(3, 4)
Out[50]: 0.43185509566737346

    An example where the naive implementation overflows, but i0_ratio does not:

    In [51]: np.i0(650)/np.i0(720)
/Library/Frameworks/Python.framework/Versions/3.10/lib/python3.10/site-packages/numpy/lib/function_base.py:3359: RuntimeWarning: overflow encountered in exp
  return exp(x) * _chbevl(32.0/x - 2.0, _i0B) / sqrt(x)
Out[51]: 0.0

In [52]: i0_ratio(650, 720)
Out[52]: 4.184118428217628e-31

    In Matlab, to get the same scaled Bessel function as i0e(x), you can use besseli(0, x, 1).