% 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;
