2
$\begingroup$

Are there any "good" low-computational-cost approximations for parameterized Kaiser window generation suitable for small systems/languages that do not include Bessel functions in their standard library?

(maybe only to 8 to 12-bit accuracy, and something that can be run in Basic on an Apple II, et.al.)

$\endgroup$
3
  • 1
    $\begingroup$ I don't know how good this technique is, but CORDIC has been adapted to compute the Kaiser window: link.springer.com/article/10.1007%2Fs11265-013-0781-z $\endgroup$
    – MBaz
    Feb 18, 2017 at 21:59
  • $\begingroup$ Apple II? Computational wise LUT will probably be the lowest costing approach. Only the storage is the concern. Or may be you can truncate the polynomial Bessel series with enough terms for 12 bit accuracy. $\endgroup$
    – Fat32
    Feb 18, 2017 at 23:46
  • 2
    $\begingroup$ Kaiser is an approximation of the DPSS (discrete prolate spheroidal sequence) a.k.a. Slepian window. So do you need Kaiser in itself or would you accept something with similar performance? $\endgroup$ Mar 4, 2017 at 8:49

2 Answers 2

7
$\begingroup$

just implement the Modified Bessel function. it's easy.

i always like my window definitions centered about zero, since pretty much all of them are even symmetry.

i'll do this in discrete-time, but it's essentially the same thing in continuous-time.

Kaiser window:

$$ w[n] \triangleq \begin{cases} \frac{1}{I_0(\beta)} I_0\left(\beta \sqrt{1 - \left(\frac{n}{M/2}\right)^2} \right) \quad \quad & |n| \le M/2 \\ 0 & |n|>M/2 \\ \end{cases} $$

$I_0(x)$ is the 0th-order modified Bessel function of the 1st kind. $M+1$ is the number of non-zero samples or FIR taps (the FIR filter order is $M$ and, in my centered and symmetrical case, must be even). $\beta$ is a "shape parameter" and O&S recommend this heuristic:

$$ \beta = \begin{cases} 0.1102 \cdot (A-8.7) & A>50 \\ 0.5842 \cdot (A-21)^{2/5} + 0.07886 \cdot (A-21) \quad & 21 \le A \le 50 \\ 0.0 & A<21 \\ \end{cases}$$

$$ M = 2 \left\lceil \frac{A-8}{4.57 \cdot \Delta\omega} \right\rceil $$

$A$ is the desired stopband attention in dB and $\Delta\omega$ is the desired width of the transition band in normalized angular frequency.

finally, the Bessel is evaluated as:

$$ I_0(x) = 1 \ + \ \lim_{K \to \infty} \ \sum\limits_{k=1}^{K} \left(\frac{x^2}{4}\right)^{k} (k!)^{-2} $$

when you evaluate this with a computer, pick a $K$ decently large (my guess is that $K=32$ is good enough) and evaluate the summation starting with $k=K$ and work it backwards to $k=1$ to keep numerical accuracy. you might want to use Horner's method.

$$ I_0(x) \approx 1 + x^2\left( \tfrac{1}{(1!)^2 \, 4^1} + x^2\left(\tfrac{1}{(2!)^2 \, 4^2} + x^2\left(... + \, x^2\left(\tfrac{1}{((K-1)!)^2 \, 4^{K-1}} + x^2 \tfrac{1}{(K!)^2 \, 4^K} \right) \right) \right) \right) $$

you can evaluate all of the $(k!)^{-2}$ in advance with a short table.


Alright, someone made me do some work. This took about 45 minutes to code up and debug. So here is my MATLAB code for implementing the 0th-order Bessel function of the first find (which is $I_0(x)$ above):

function y = mybessel(x)
%
%   Computes the 0th-order Modified Bessel function of the first kind
%

    K = 32;

    bessel_coef =  zeros(1,K);

    kfac = 1;
    two_to_the_k = 1;

    for k = 1:K
        kfac = kfac * k;
        two_to_the_k = two_to_the_k * 2;
        bessel_coef(k) = 1/(kfac*two_to_the_k)^2;   % compute power series coefficients in advance
    end

    x = x.^2;

        y = x .* bessel_coef(K);
    for k = K-1:-1:1
        y = x .* (bessel_coef(k) + y);              % Horner's method
    end

    y = 1 + y;

end

and here is the test code:

x = linspace(-16, 16, 32*4096+1);

I_0 = besseli(0, x);            % MATLAB's modified bessel

y = mybessel(x);                % my bessel

figure(1)
plot(x, I_0)
hold on
plot(x, y)
hold off

figure(2)
plot(x, y - I_0)                % plot error

with results:

enter image description here

and error:

enter image description here

for $|x| \le 16$, the relative error is less than $10^{-15}$. the error increases with increasing $|x|$.

$\endgroup$
18
  • $\begingroup$ I tried the Horner's method from above to calculate a Bessel approximation. It works great for small numbers < 1.0. But my beta at 90 dB is 8.959. When I calculate the Bessel approximation for numbers that large then it gets pretty wild. Here are Bessel values for K = 7,8,9. 7, -4100.13715 8, 1284.70023 9, -318.22047 $\endgroup$
    – philburk
    Aug 17, 2019 at 17:42
  • $\begingroup$ well @philburk , of course the values of $\left(-\frac{x^2}{4}\right)^{k}$ will be increasing with increasing $k$, but they are divided by $(k!)^2$ which i think is increasing even faster. so the terms should start to converge. $\endgroup$ Aug 17, 2019 at 20:40
  • $\begingroup$ If x >= k then the numerator is significantly larger than the denominator. My x is 8.959, so a K=8 is not close to converging. I found that a K of 16 was needed for the direct summation method, and a K of 18 was needed for the Horner's method. $\endgroup$
    – philburk
    Aug 19, 2019 at 16:41
  • $\begingroup$ Another thing I find puzzling. This Bessel function has three zero crossings before it reaches x=8.959. That means the Kaiser window will also have three zero crossings. Also what happens if i0(beta) is zero, or very close to zero? $\endgroup$
    – philburk
    Aug 19, 2019 at 17:15
  • $\begingroup$ @philburk Something must be wrong, because Io(0)=1 and the function keeps on going up, there should not be any zero crossings. $\endgroup$ Aug 20, 2019 at 6:49
6
$\begingroup$

You could try the exponential window:

$$w_n=\frac{ \exp \left[\alpha \sqrt{1-\left(\frac{n-M}{M}\right)^2}\right]}{\exp(\alpha)}$$ $$\alpha=-427.5*10^{-6} A_s^2+0.1808*A_s-3.516$$

or the hyperbolic cosine window:

$$w_n=\frac{\cosh \left[ \alpha \sqrt{1-\left(\frac{n-M}{M}\right)^2} \right]}{\cosh(\alpha)}$$ $$\alpha=-325.1*10^{-6}*A_s^2+0.1677*A_s-3.149$$

Both have similar results with themselves and the Kaiser window, while being less computational intensive.

If you want to go with Robert's answer, you could simplify the Bessel function by using Horner's method, or by applying a Gauss-Kronrod (or whatever other) integration. For example, I got this formula (after some reductions&co):

$$f(x)=\frac{\left[ ~~ \cosh(x) ~+~ 2*(\cosh(0.970941817426052*x) ~+~ \cosh(0.8854560256532099*x) ~+~ \cosh(0.7485107481711011*x) ~+~ \cosh(0.5680647467311558*x) ~+~ \cosh(0.3546048870425356*x) ~+~ \cosh(0.120536680255323*x)) \right]}{13}$$

For $\beta=12.26526$ ($A_s=120dB$), results in a these values:

$I_0(\beta)=24430.40185694905$ $f(\beta)=24430.4018627702$

which, IMHO, is pretty good and it's faster than using giga-km long numbers you'd deal with in the case of the classic formula, either as a sum or as Horner -- either way, I think even long double would not be enough for 20 terms (which is how I've truncated it).

$\endgroup$
7
  • $\begingroup$ I implemented your Bessel function f(x) and get values that match the numpy.i0() function in Python. But they seem unrelated to the Bessel values in Robert's code. From Robert's code with K=16 I get values that start at 1 and wiggle down, like: mathworld.wolfram.com/BesselFunctionoftheFirstKind.html But your function looks more like a parabola than a Bessel function. oreilly.com/library/view/numpy-15/9781849515306/ch07s32.html But maybe that results in a better Kaiser window that does not have zero crossings... $\endgroup$
    – philburk
    Aug 19, 2019 at 17:30
  • $\begingroup$ What I did was to first express the Bessel Io(x) as its integral representation, $\frac{1}{\pi}\int_0^{\pi}\cosh(x\cos\theta)d\theta$, then used that with Gauss-Kronrod (or Gauss-Legendre? don't remember). If you expand Robert's series, you get monstruous numbers, even with Horner, which makes it difficult to preserve accuracy. Also, IIRC, you need some 20 terms for the series expansion to have decent approximation for $\beta\ge 120dB$. $\endgroup$ Aug 20, 2019 at 6:46
  • $\begingroup$ I'm curious why your Bessel function looks so different than Robert's. Yours starts at one and goes up like a parabola. Robert's goes down and has zero crossings. Are they different kinds of Bessel functions? Is one better for the Kaiser window. (BTW, I used your hyperbolic cosine window from above and it worked great. Thanks.) $\endgroup$
    – philburk
    Aug 20, 2019 at 18:54
  • $\begingroup$ @philburk, the Kaiser window is based on the 0th-order Bessel function of the first kind. i think that's well-defined. $\endgroup$ Aug 21, 2019 at 23:07
  • $\begingroup$ The exponential and cosh windows have nearly identical Fourier spectra. The cosh might be slightly better but I'm unsure. The same author published both, and the papers are nearly identical but somehow avoid comparisons to each other, even indirectly through known windows. They are marginally worse than the Kaiser but the ease of implementation makes the tradeoff worth it. $\endgroup$
    – Kevin Yin
    Jun 29, 2022 at 16:39

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.