% Matt Stone under the guidance of Mladen Victor Wickerhauser
% This program replicates the graphics obtained in Strang & Shen's "The
% Zeros of the Daubechies Polynomials"
% I attempted to not just generate pictures for fixed length polynomials but
% dynamically plot increasing filters to "show" convergence
% We first generate the coefficients for each polynomial from
% p = 3 to p = maxP
% According to the paper, after 80, bifurcation occurs due to
% computer roundoff.
% I thought I might be able to increase it, but this still holds today,
% ~20 years after the paper.
maxP = 80;
% Each column represents the coefficients of B_p(y) where the diagonal has
% the coefficients of the constant term and the 1st row of column p has
% the coefficients of the term y^(p-1)
% This indexing is to correspond with MATLAB's roots method
% columns are ordered from p = 1 to p = 80 from left to right
% Coefficients are obtained for B_p(y) passed through the transform
% y -> 4y in order to obtain a higher number of accurate roots
coeff = zeros(maxP);
for p = 1:maxP
b = zeros(maxP, 1);
b(p) = 1;
for i = p - 1:-1:1
b(i) = b(i+1)*(2*p - i - 1)/(4*(p-i));
end
coeff(:, p) = b;
end
% With this, we plot the roots
% Note that we plot not roots(b) but roots(b)/4 to account for our
% earlier transformation
% uncommenting the hold on and hold off command
% plots all roots of consecutive Daubechies polynomials on one graph
% hold on;
for p = 3:maxP
rDaub = roots(coeff(1:p, p))/4;
plot(rDaub, '.');
pause(0.5);
end
% hold off;