% Product filter examples

%p = 2

switch p

case 1,
  % Degree 2
  b = [1 2 1];  % (1 + z^-1)^2
  q = 1 / 2;
  p0 = [1 2 1] / 2;   % conv(b, q)

case 2,
  % Degree 6
  b = [1 4 6 4 1];  % (1 + z^-1)^4
  q = [-1 4 -1] / 16;
  p0 = [-1 0 9 16 9 0 -1] / 16;  % conv(b, q)

case 3,
  % Degree 10
  b = [1 6 15 20 15 6 1];   % (1 + z^-1)^6
  q = [3 -18 38 -18 3] / 256;
  p0 = [3 0 -25 0 150 256 150 0 -25 0 3] / 256;  % conv(b,q)

case 4,
  % Degree 14
  b = [1 8 28 56 70 56 28 8 1];  % (1 + z^-1)^8
  q = [-5 40 -131 208 -131 40 -5] / 2048;
  p0 = [-5 0 49 0 -245 0 1225 2048 1225 0 -245 0 49 0 -5] / 2048;  % conv(b,q)

otherwise,
  % Degree 4p-2
  [p0,b,q] = prodfilt(p);
end

zplane(p0);
title(sprintf('Zeros of the product filter with degree %d', 4*p-2))
pause
[P,W] = dtft(p0,512);
plot(W/pi, abs(P))
xlabel('Angular frequency (normalized by pi)')
ylabel('Frequency response magnitude')
title(sprintf('Frequency response of the product filter with degree %d', 4*p-2))