0
$\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$
  • 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 '17 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 '17 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$ – Olli Niemitalo Mar 4 '17 at 8:49
5
$\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$
  • $\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 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$ – a concerned citizen Aug 20 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 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$ – robert bristow-johnson Aug 21 at 23:07
5
$\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$
  • $\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 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$ – robert bristow-johnson Aug 17 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 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 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$ – a concerned citizen Aug 20 at 6:49

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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