%% Haar Transform 1D
% David Ruch and Patrick J. Van Fleet
%
% Minicourse #4 January 2008 Joint Mathematics Meetings
%
% San Diego, CA

%% Objectives
% In this m-file, we will explore the one-dimensional discrete Haar
% wavelet transform.

%% Conventions
% This m-file uses the DiscreteWavelets Toolbox.  Help is available on all
% functions in the toolbox by accessing the Help menu and clicking on the
% DiscreteWavelets Toolbox.  
%
% Demos are also available.  Click the Demos tab under Help to access all
% demos.

%% The One-Dimensional Haar Transform
% Here is the m-file for computing the one-dimensional Haar wavelet
% transformation.  Note the routine checks to see if the input vector is of
% even length.

type HaarWT

%% Using the Haar Wavelet Transform
% It is easy to use the HaarWT function

x=1:8
y=HaarWT(x);
y'

%%
% What should HaarWT do to a constant vector?

x=ones(1,10);
y=HaarWT(x);
y'

%%
% Here is an example that introduces a command we can use to plot the
% wavelet transform.

x=1:100;
v=x.^2;
y=HaarWT(v);
plot(v);
figure;
WaveletVectorPlot(y,1,'UseColors','True');

%%
% Plot the differences only

WaveletVectorPlot(y,1,'Region','HighPass','Iteration',1,'UseColors','True');

%% One-Dimensional Inverse Haar Wavelet Transformation
% The function IHaarWT computes the one-dimensional inverse discrete Haar
% wavelet transform.  It takes as input a vector of length n.  The number
% of elements must be even or the function returns an error message.
%
% The routine starts with an error check on the length of the input.  Then
% the input is partitioned into a 2 x n/2 matrix.  The matrix is multiplied
% by [1 -1] to obtain the odd elements of the inverse transform vector and
% (1,1) to form the even elements of the inverse transform.  The two parts
% are then intertwined and returned.

type IHaarWT

%% Using IHaarWT
% The cell below contain illustrates the use of IHWT1D.   Feel free to
% change the input vector x.  For something different, I've formed x as a 
% list of 20 random integers whose values range from 0 to 100.

x = round(100*rand(1,20))
y = HaarWT(x)'
newx = IHaarWT(y)'

%% An Illuminating Example for Students
% I use this example every time I introduce the Haar wavelet transform.  
% The student can create a basic function from Calculus I (some examples 
% are provided).  The function is then sampled 100 times from a to b.  
% This creates a list x that we plot.  Note the x-axis on the plot is in 
% terms of the list elements while the y-axis is in terms of the list 
% values.
% 
% Feel free to uncomment other functions or define your own.

n=100;
t=0:1/n:1-1/n;
v=sin(2*pi*t);
%v=cos(2*pi*t);
%v=t.^2;
%v=t;
%v=1;
%v=exp(t);
%v=log(1+t);
plot(v);

%% 
% Next compute the HWT of x.    The first plot is the entire transform.  
% The second plot shows only the approximation portion of the transform 
% while the third plot shows only the differences portion of the transform.  
%
% How are the second and third plots related?  Try a couple of other 
% functions.

y = HaarWT(v);

figure;
WaveletVectorPlot(y,1);
title('HWT');
figure;
WaveletVectorPlot(y,1,'Region','LowPass');
title('Averages');
figure;
WaveletVectorPlot(y,1,'Region','HighPass')
title('Differences');

%% Iterating the Process
% In many applications, scientists will iterate the wavelet transform.
% Suppose you are given n-vector x and you compute the HWT y.    The next
% iteration of transform is applied to the approximation portion of the y.
% Here is some code (bulky) to do the job.  Note that to iterate a second 
% time requires n to be divisible by 4.  

x = [2 4 6 8 10 12 14 16];
n = length(x);
y = HaarWT(x)
approx = y(1:n/2);
diff = y(n/2+1:end);

z = HaarWT(approx);
iterated = [z; diff]'

%%
% In general, if the length of x contains the factor 2^i, we can iterate i
% times.  There is a routine in DiscreteWavelets for computing the 
% iterated Haar Wavelet Transform.  It is called HWT1D.  It needs two 
% arguments.  The first is the input vector.  The second input is the 
% number of iterations.  
%
% Here is an example.  Note that the length of the input vector is 
% divisible by 32 so we can perform as many as 5 iterations.

n = 96;
x = (0:n-1)/n;
y = HWT1D(x,2);

close all;
WaveletVectorPlot(y,2)

%%
% Here is three iterations of the transform.  Note the approximation
% portion contains 96/2^3 = 12 elements.

y = HWT1D(x,3);
WaveletVectorPlot(y,3);

%%
% You can look at particular portions of the transform.

WaveletVectorPlot(y,3,'Iteration',2,'Region','HighPass');
figure;
WaveletVectorPlot(y,3,'Region','LowPass');

%%
% The process is invertible using the IHWT1D routine from DiscreteWavelets.
% Like WT1D, the routine requires two arguments.  The first is the input 
% vector and the second is the number of iterations.

x = round(100*rand(1,32))
y = HWT1D(x,5);
orig = IHWT1D(y,5)'

%% Orthogonalizing the Transform
% In the early development of wavelets, the transforms were designed to be
% orthogonal matrices.  Consider the Haar Wavelet transform matrix of size
% 8 x 8 and its inverse:

W=[.5 .5 0 0 0 0 0 0; 0 0 .5 .5 0 0 0 0; 0 0 0 0 .5 .5 0 0; 0 0 0 0 0 0 ...
    .5 .5; -.5 .5 0 0 0 0 0 0; 0 0 -.5 .5 0 0 0 0; 0 0 0 0 -.5 .5 0 0;...
    0 0 0 0 0 0 -.5 .5]

inv(W)

disp('How can we modify the matrix so that it is orthogonal?');

%%
% When we use orthogonal transformations, we can use the same list as the
% second arguments for WT1D and IWT1D.  There is a command in 
% DiscreteWavelets for generating this list:

Haar()

%%
close all;
displayEndOfDemoMessage(mfilename)