% function [p0,b,q] = prodfilt(p)
%
% Generate the halfband product filter of degree 4p-2.
%
%              Kevin Amaratunga
%              19 February, 2003
%  p0 = coefficients of product filter of degree 4p-2.
%  b = coefficients of binomial (spline) filter of degree 2p
%  q = coefficients of filter of degree 2p-2 that produces the halfband
%      filter p0 when convolved with b.

function [p0,b,q] = prodfilt(p)

% Binomial filter (1 + z^-1)^2p
tmp1 = [1 1];
b = 1;
for k = 0:2*p-1
  b = conv(b, tmp1);
end

%  Q(z)
tmp2 = [-1 2 -1] / 4;
q = zeros(1,2*p-1);
vec = zeros(1,2*p-1);
vec(p) = 1;
for k=0:p-1
  q = q + vec;
  vec = conv(vec, tmp2) * (p + k) / (k + 1);
  vec = wkeep(vec, 2*p-1);
end
q = q / 2^(2*p-1);

% Halfband filter, P0(z).  
p0 = conv(b, q);