% Eigenvalues of the transition matrix, T = down2(2 H H')

h = [1 1]' / 2; p = 1;  % smax = 0.5;
h = daub(8) / 2; p = 2; % smax = 1.0 
%h = [1 2 1]'/4; p = 2;  % smax = 1.5
h = [1 4 6 4 1]' / 16; p = 4; 
%[x,phi,psi] = phivals(h,5);

%[f,h] = biorwavf('bior2.2'); p = 2; h = h';  % 5/3
%[h,f] = biorwavf('bior2.2'); p = 2; h = h';  % 3/5
%[f,h] = biorwavf('bior4.4'); p = 4; h = h'; % 9/7
%[h,f] = biorwavf('bior4.4'); p = 4; h = h'; % 7/9
%h = [-1 0 9 16 9 0 -1]' / 32; f = 1; p = 4; % Halfband filter
%[h0,h1,f0,f1] = biorfilt(h',f);
%[x,phi,psi,phitilde,psitilde] = biphivals(h0,h1,f0,f1,5);

%plot(x,phi);
%pause

a = conv(h, flipud(h));

% Method 1: Use the function given in Ch 7, p221.
T = down(a');

% Method 2: Use basic Matlab commands.
N = length(a);
v = [a; zeros(N,1)];
V = toeplitz(v, [v(1) zeros(1, N-1)]);
TT = 2 * dyaddown(V, 'r', 1);

if norm(T - TT) ~= 0
  error('Something is wrong!')
end

lambda = flipud(sort(eig(T)))


pause

hs = 1;
for i = 1:p
  hs = conv(hs, [1,1]/2);
end
hq = deconv(h, hs);
aq = conv(hq, flipud(hq));
TQ = down(aq');
lambdaQ = flipud(sort(eig(TQ)));
lambdaMax = max(lambdaQ) / 4^p;

disp(sprintf('Largest nonspecial eigenvalue = %0.14f', lambdaMax));
disp(sprintf('Smallest special eigenvalue = %0.14f', 1/2^(2*p-1)));
disp(sprintf('Smoothness, s_max, = %0.14f', -log2(lambdaMax)/2));

