function [w, idx] = wkaiser(n, b)

% function [w, idx] = wkaiser(n, b)

k1  = 1.0/besselizero(b);
k2  = 1 - mod(n, 2);
ende = fix((n + 1)/2);

idx = zeros(n, 1);
% Calculate window coefficients
for k=0:ende-1,
    tmp = (2*k + k2) / (n - 1.0);
    tmp2 = k1 * besselizero(b*sqrt(1.0 - tmp*tmp));
    mm = ende-(mod(not(k2), 2))+k+1;
    nn = ende-k;
    w(mm) = tmp2;
    w(nn) = tmp2;
    idx(nn) = idx(nn) + 1;
    idx(mm) = idx(mm) + 1;
end;

function sum = besselizero(x)
    BIZ_EPSILON = 1E-21;    % Max error acceptable 
    sum   = 1.0;
    u     = 1.0;
    halfx = x/2.0;
    n     = 1;
    while(1)
        temp = halfx/n;
        u = u * temp * temp;
        sum = sum + u;
        n = n + 1;
        if (u < (BIZ_EPSILON * sum))
            break;
        end;
    end;