  function P4Grade()
% Checks out HaarThreshApprox
Y = imread('Moon.jpg');
tau = 0.05;
% My implementation...
tic
[Ytilde_CVL,CF_CVL] = HaarThreshApproxCVL(double(Y),tau);
tCVL = toc;
% Your implementation...
tic
[Ytilde_YOU,CF_YOU] = HaarThreshApprox(double(Y),tau);
tYOU = toc;
err = max(max(abs(Ytilde_YOU-Ytilde_CVL)));
fprintf('tYOU/tCVL = %5.2f \n',tYOU/tCVL)
fprintf('Ytilde_Error = %10.2e\n',err)
fprintf('CF_CVL = %6.3f\n',CF_CVL)
fprintf('CF_YOU = %6.3f\n',CF_YOU)

% Show the compressed image produced by your code
imshow(uint8(Ytilde_YOU))
title(sprintf('tau = %4.2f   CF = %4.2f',tau,CF_YOU),'fontsize',16)
shg

  function [Ytilde,CF] = HaarThreshApproxCVL(Y,tau)
% Compressed Haar representation of the n1xn2 matrix Y
% where n1 and n2 are powers of two.
% tau satisfies 0<=tau<=1
% Ytilde is n1xn2 and is given by W_n1*Xtilde*W_n2' where X_tilde is
% obtained from the inverse transform of Y with thresholding.
% CF is the ratio of n1*n2 to the number of nonzeros in Xtilde
[n1,n2] = size(Y);

% Inverse transform of Y...
X = IFHTvec(Y);
X = IFHTvec(X')';
% Thresholding...
mX = max(max(abs(X)));
Xtilde = X.*(abs(X)>=tau*mX);
CF = n1*n2/sum(sum(abs(Xtilde)>0));
% Forward transform of Xtilde...
Ytilde = FHTvec(Xtilde);
Ytilde = FHTvec(Ytilde')';

  function Y = FHTvec(X)
% Fast Haar transform on multiple vectors.
% X is nxp matrix and n is a power of two.
% Y is nxp matrix and Y(:,k) is the Haar transform of X(:,k), k=1:p
[n,p] = size(X);
if n==1
    Y = X;
else
    m = n/2;
    Y = zeros(n,p);
    Z = FHTvec(X(1:m,:));
    Y(1:2:n,:) = Z + X(m+1:n,:);
    Y(2:2:n,:) = Z - X(m+1:n,:);
end

  function X = IFHTvec(Y)
% Inverse Fast Haar transform on multiple vectors.
% Y is nxp matrix and n is a power of two.
% X is nxp matrix and X(:,k) is the inverse Haar transform of Y(:,k), k=1:p

[n,p] = size(Y);
if n==1
    X = Y;
else
    YO = Y(1:2:n,:);
    YE = Y(2:2:n,:);
    X = [IFHTvec((YO+YE)/2) ; (YO-YE)/2];
end

% -2 for no vectorization
% -1 for thresholding errors
