Kaiser window approximation
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 &amp; |n| \le M/2 \\
0 &amp; |n|&gt;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&amp;S recommend this heuristic:

$$ \beta = \begin{cases} 0.1102 \cdot (A-8.7) &amp; A&gt;50 \\
0.5842 \cdot (A-21)^{2/5} + 0.07886 \cdot (A-21) \quad &amp; 21 \le A \le 50 \\
0.0 &amp; A&lt;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|$.