Kaiser window approximation

Question

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.)

Answer 1

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:

and error:

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

Answer 2

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).