function coef = jackson_dct(n,r)

% Compute the coefficients for the generalized Jackson kernel
%    K_{m,2r}, see (4.2).
% n:    the size of the Toeplitz matrix T;
% r:    2, 3 or 4;
% coef: the coefficients {b_k^(m,r):|k|<=r(m-1)} of
%           the generalized Jackson kernel K_{m,2r}.
% The code obtain the coefficients by interpoltating on K_{m,2r}
%   which is not the same as the method presented in the book.

m = floor(n/r);
k1 = r*(m-1)+1;
k = 2^nextpow2(r *(m-1)+1);

a = zeros(k,1);

a(1:m) = 2:-2/m:2/m;  a(1) = 1; % K_{m,2}(xi) = a(1) + a(2)*cos(xi) + ... + a(m-1)*cos((m-1)*xi);
z = exp([0:pi/(2*k):(k-1)*pi/(2*k)]'*sqrt(-1)); 

coef = (k*real(ifft(a.*z))).^r; % evaluate K_{m,2r}((2j-1)pi/(2k)) = [K_{m,2}((2j-1)pi/(2k))]^r, j = 1,2,...,k

% coef also can be obtained as follows by using inverse fast cosine transform.
% a(1) = a(1)*sqrt(k); a(2:m) = a(2:m)*sqrt(k/2);
% coef2 = idct(a).^r;
% coef(1:k/2) = coef2(1:2:k-1); coef(k/2+1:k) = coef2(k:-2:2);

coef = real(conj(z).*fft(coef)); % compute k b_j^{m,2r}, j = 0,1,...,r(m-1)
coef = [coef(1:k1)/k; zeros(n-k1,1)];

% coef also can be obtained as follows by using fast cosine transform.
% coef = dct(coef2);
% coef = [coef(1)/sqrt(k); coef(2:k1)/sqrt(2*k); zeros(n-k1,1)];
coef = coef';
