Discrete Haar Wavelet Transform

See article at http://www.dsprelated.com/showabstract/3639.php for the theory.

function X = pm_haar(x,s)
%
% MATLAB function for pm Haar 1-D DWT algorithm. 
% This function receives N (a power of two) real values in x
% and the number of stages (scales) of decomposition required in s, 
% computes the Haar DWT, and returns the N DWT coefficients in X.
%
% For the theory behind this algorithm and example input/output,
% please refer the paper:
% Fundamentals of the discrete Haar wavelet transform
% Duraisamy Sundararajan
% dsprelated.com, 2011
% articles/paper section
%
N = length(x);  % length of the input vector 
b = N;sq2 = sqrt(2);fac = 1;
  for ns =1:s % outer loop stepping over stages
    p = 1;q = 1;
    while(p < b) % inner loop stepping over each butterfly in a stage
      r = p + 1;
      xmr(q) = x(p) - x(r); xpr(q) = x(p) + x(r);
      p = p + 2; q = q + 1;
    end
    fac = fac / sq2;
    x(1:b) = [xpr(1:b/2) xmr(1:b/2) * fac];
    b = b / 2;
  end
x(1:b) = fac * x(1:b);
X=x;

function x = pm_haar_inv(X,s)
%
% MATLAB function for pm Haar 1-D IDWT algorithm. 
% This function receives N (a power of two) Haar DWT coefficients in X
% and the number of stages (scales) decomposed in s, 
% computes the Haar IDWT, and returns the N real values in x.
%
% For the theory behind this algorithm and example input/output,
% please refer the paper:
% Fundamentals of the discrete Haar wavelet transform
% Duraisamy Sundararajan
% dsprelated.com, 2011
% articles/paper section
%
N = length(X); % length of the input vector 
b = N/2^s;sq2 = sqrt(2);fac=1/(sq2^s);
X(1:b) = X(1:b) * fac;
  for ns =1:s  % outer loop stepping over stages
    p = 1;q = 1;
    xpr(1:b) = X(1:b);
    xmr(1:b) = X(b+1:2*b) * fac;
    while(p <= b)  % inner loop stepping over each butterfly in a stage
      X(q + 1) = xpr(p) - xmr(p); X(q) = xpr(p) + xmr(p);
      p = p + 1; q = q + 2;
    end
    fac = sq2 * fac;
    b = b * 2;
  end
  x = X;
  

% Main program for calling functions to compute forward and inverse
% Haar 1-D DWT
% functions pm_haar.m and pm_haar_inv.m are required in the same directory
clear;
% Example input x = [0 1 2 3 4 5 6 7], s = 3
% Example output X = [9.8995, -5.6569, -2, -2,
% -0.7071, -0.7071, -0.7071, -0.7071]
%
x = [0 1 2 3 4 5 6 7];
s = 3; % number of stages of computation required, s = 1,2,...,log2(N),
% where N is the number of elements in the input vector x. N is an
% integral power of two
X = pm_haar(x,s); % Invoke pm_haar.m for
% computing the forward DWT
X % the Haar DWT of input vector x after s stages of computation
x = pm_haar_inv(X,s);  % Invoke pm_haar_inv.m 
% for computing the IDWT
x % the given input after s stages of inverse computation

% Main program for calling functions to compute forward and inverse
% 2-D Haar DWT
% functions pm_haar.m and pm_haar_inv.m are required in the same directory
clear;
% Example input x2 = [
% 1 2 -1 3
% 2 1 4 3
% 1 1 2 2
% 4 2 1 3
% ]
% scale = 2
% Example output X = [
% 7.7500   -0.7500         0   -1.5000
%-0.2500   -0.7500    1.0000   -1.0000
%      0   -2.5000   -1.0000   -2.5000
%-2.0000         0   -1.0000    1.0000
% ]
x2 = [1 2 -1 3;2 1 4 3;1 1 2 2;4 2 1 3];
N = size(x2,1); % size of each row, N must be an integral power of two 
scale = 2; % number of stages (scales) of computation required,
% s = 1,2,...,log2(N)
lim = N; % size of the row or column of the subimage in different stages 
% computation of the forward 2-D Haar DWT of x2 by row-column method
for k=1:scale % compute the 2-D DWT scale by scale
%    
for j=1:lim
    x = x2(j,1:lim); % copy each row of input
    X = pm_haar(x,1); % Invoke pm_haar.m to compute each row DWT
    x2(j,1:lim) = X;
end
%
for j = 1:lim
    x = x2(1:lim,j); % copy each column of the result of row DFTs
    X = pm_haar(x',1);  % Invoke pm_haar.m to compute each column DWT
    x2(1:lim,j) = X';
end
lim = lim / 2; % transform size reduced for each scale
end
x2 % 2-D DWT of x2 

% computation of the 2-D Haar IDWT of x2 by row-column method
lim = N/(2^(scale-1)); % size of the row or column of the subimage in
% different stages computation of the 2-D Haar IDWT of x2 by
% row-column method
for k = 1:scale % compute the 2-D IDWT scale by scale

for j = 1:lim
    X = x2(j,1:lim);  % copy each row of the coefficients
    x = pm_haar_inv(X,1); % Invoke pm_haar_inv.m to compute each row IDWTs
    x2(j,1:lim) = x;
end
%
for j = 1:lim
    X = x2(1:lim,j);  % copy each column of the results of row IDWTs
    x = pm_haar_inv(X',1); % Invoke pm_haar_inv.m to compute column IDWTs
    x2(1:lim,j) = x';
end
lim = lim * 2; % transform size increased for each scale
end
%
x2 % the given input after required stages of IDWT computation