% Write the discrete Haar wavelet transform 
% in matrix form. Here we include several scales.
%
% Samuli Siltanen Dec 2014

% Load a sample signal
load rowdata
signal = rowdata;
%signal = [1:8]; % For debugging and simple analysis
N = length(signal);

% Check that the signal length is a power of two
if (log2(N)-floor(log2(N)))>0
    error('Signal length should be a power of two')
end

% Build the Haar transform matrix
W = eye(N);
for iii=1:log2(N)
    tmp = zeros(N);
    % Calculate the size of the convolution part of the matrix
    cN = N/(2^(iii-1));
    % Calculate the size of the unit block matrix
    uN = N-cN;
    % Construct the block matrix
    tmp(1:uN,1:uN) = eye(uN);
    tmp((end-cN+1):end,(end-cN+1):end) = HaarTransformMatrix(cN);
    % Update transform matrix
    W = tmp*W;
end

% Compute Haar transform
HT = W*signal(:);

% Compute inverse Haar transform
recon = (W.')*HT;

% Take a look at the original signal 
% and reconstruction
figure(1)
clf
plot(signal,'r')
hold on
plot(recon,'b')

% Take a look at the original signal 
% and the transform
figure(2)
clf
subplot(2,1,1)
plot(signal,'r')
subplot(2,1,2)
plot(HT,'b')