% Example 4: Examine how polynomial data behaves in the filter bank
% when the lowpass filter has p zeros at pi.

p = 3;              % number of zeros at pi.
%a = [1 1]/128;
a = [1 3 3 1]/128;    % Coefficients of the polynomial (low order to high order.)
n = 0:127;
q = length(a) - 1;  % Degree of the polynomial.
len = length(n);
x = zeros(1,len);
nq = ones(1,len);
for k = 1:q+1
  x = x + a(k) * nq;
  nq = nq .* n;
end

% Compute the DWT.
N = 2 * p;
h0 = daub(N);
h1 = (-1).^[0:N-1]' .* flipud(h0);
[y0,y1] = dwt(x,h0,h1);

% Plot the results.
clf
subplot(3,1,1);
plot(x)
axis([0 len-1 min(x) max(x)])
title(sprintf('Wavelet transform of degree %d polynomial data.  H0(w) has %d zeros at pi.', q, p))
ylabel('x')
subplot(3,1,2);
plot(y0)
axis([0 len-1 min(y0) max(y0)])
ylabel('y0')
subplot(3,1,3);
plot(y1)
minval = min(0,min(y1(p:length(y1)-p+1)));
maxval = max(0,max(y1(p:length(y1)-p+1)));
s = maxval - minval;
minval = minval - 0.2*s;
maxval = maxval + 0.2*s;
axis([0 len-1 minval maxval])
ylabel('y1')

fprintf('Maximum value of y1 (excluding boundary effects) = %d.\n',max(abs(y1(p:length(y1)-p+1))))