Discrete Wavelet Transform via circular convolution matrix

The following code computes the Discrete Wavelet Transform via a circular convolution matrix, that is built from equivalent filters generated from basic filters by the function "filterform.m"

Required for this program are:

http://www.dsprelated.com/showcode/10.php

http://www.dsprelated.com/showcode/12.php

For more details on this code and its dependencies, please visit these blog posts:
http://www.dsprelated.com/showarticle/115.php
http://www.dsprelated.com/showarticle/116.php

% ----------------------------------------------------
% Title: Discrete Wavelet Transform
% Author: David Valencia
% UPIITA IPN 2010
%
% Computes the Discrete Wavelet Transform of
% n levels (and its branches) from a n sample input.
% Base filters can be of any lenght.
%
% Generalized DWT filter bank computed via the
% convolution matrix method (circular convolution)
%
% Published in: http://www.dsprelated.com/showcode/23.php
%
% Dependencies:
%
% formfilter.m
% http://www.dsprelated.com/showcode/12.php
% upsample2.m
% http://www.dsprelated.com/showcode/10.php
%
% Revisions:
%       v1.0a Commented and translated in English
% ----------------------------------------------------

close all; clear all; clc;

disp('== GENERALIZED MATRICIAL DWT ==')

%% STEP 1 - Base Filters
% Select type of base filters
typeofbasis = 'o';
typbior = 'bior2.2';
typor = 'db2';

% Obtain base filters
if(typeofbasis == 'b')
    [Rf,Df] = biorwavf(typbior);
    [h0,h1,g0,g1] = biorfilt(Df,Rf);
elseif (typeofbasis == 'o')
    [h0,h1,g0,g1] = wfilters(typor);
end;

%% STEP 2 - Parameter configuration
% One can declare or recover an input vector from another
% program, but here an example vector is built for testing.

N = 16; %Number of input samples
x = (1:N)';
L = length(h0); %Base filter lenght
n_stages = 2; %Of the low-pass stage
n_branches = 2^n_stages; %Number of branches
dec_fact = n_branches; %Decimating factor

% L = Basic filters lenght (h0 ó h1)
% N = Input vector length
% n_stages = stage quantity, it generates (2^n_stages) branches

fprintf('\nCircular convolution processing\n');

% ###############################################
% Analysis High-Pass 1-branch Circular convolution matrix

hx = formfilter(n_stages, 1, h0, h1);
Lhx = length(hx);

if Lhx>N
    beep; beep; beep;
    disp('ERROR: Input vector must be longer, increase n');
    disp('End of program ### ');
    return;
end

movil_1 = [fliplr(h1),zeros(1,N-L),fliplr(h1),zeros(1,N-L)];
lm = length(movil_1);
dec_fact = 2;
rowsperfilter = N/dec_fact; % Square matrix gets decimated
W1 = zeros(rowsperfilter,N); %MALLOC

for i = 0:(N/2)-1
    i;
    W1(i+1,:) = movil_1(lm-(N)-(i*dec_fact)+1:lm-(i*dec_fact)); %Rama pasaaltas
end
disp('High-pass matrix: ');
W1

y1 = W1*x;

% ###############################################
% Low-Pass n-branch Analyisis Circular convolution matrix
dec_fact = 2^n_stages;
rowsperfilter = N/dec_fact; % Square matrix gets decimated
% Space for a 3-dimensional array is reserved
% Each branch is assigned to a page of that array
W0 = zeros(rowsperfilter, N, n_branches); %MALLOC
Y0 = zeros(rowsperfilter, n_branches); %MALLOC

% Matrix generation
for i = 0:(n_branches/2)-1
    if n_stages < 3
        hx = formfilter(n_stages, (n_branches/2) - i, h0, h1);
    else
        hx = formfilter(n_stages, (n_branches/2) - i - 1, h0, h1);
    end
    Lhx = length(hx);
    movil = [fliplr(hx),zeros(1,N-Lhx),fliplr(hx),zeros(1,N-Lhx)];
    lm = length(movil);
    for j = 0:rowsperfilter-1
        W0( j+1 , : , i+1) = movil(lm-(dec_fact*j)-N + 1 : lm-(dec_fact*j));
    end
    Y0(:,i+1) = W0(:,:,i+1)*x;
end

disp('Low pass several-stage analysis matrix');
W0

%% Low-pass part recovering code
xnm = zeros(N, n_branches/2); %MALLOC
for i=0:(n_branches/2)-1
    xnm(:,i+1) = ((W0(:,:,i+1))')*Y0(:,i+1);
end

[n,m] = size(xnm);

xn1 = (W1')*y1; %High-pass part recovering

xn = zeros(n,1); %MALLOC
% Sum between the low-pass and high-pass recovered parts
for i=1:n
    xn(i,1) = sum(xnm(i,:)) + xn1(i);
end

disp('Original signal');
disp(x)
disp('Recovered signal');
disp(xn)

% Compute error
err = (x - xn).^2;
disp('Error :');
err = (sum(err)/N);
disp(err)