Discrete Wavelet Transform via linear convolution matrix

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

Haar's ortogonal wavelets are used.

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
%        by means of a linear convolution matrix
%
% Author: David Valencia
% UPIITA IPN
%
% Posted at: http://www.dsprelated.com/showcode/11.php
%
% Description:
% 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 (linear convolution)
%
% 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 = 'db3';

% 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 = number of stages, it generates (2^n_stages) branches

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

fprintf('\nLinear convolution matrix processing\n');

%% STEP 3: Build analysis matrices

% -------------------------------------------------------------
% High-pass analysis matrix (single scale) [linear convolution]
% As the DWT has a simple high-pass branch, the matrix is easy 
% to build with the base filter.

% Append zeros to each side of the vector to use it as a sliding
% part to build the matrix
movil_1 = [zeros(1,N-1),fliplr(h1),zeros(1,N-1)];
lm = length(movil_1);
n_rows = ceil((L+N-1)/2); %Compute the number of rows
W1 = zeros(n_rows,N); % Allocate memory for the high pass matrix
dec_factor = 2; %Decimate factor for the high-pass matrix

% Build matrix (with sequential shifting)
for i = 0:(n_rows-1)
    W1(i+1,:) = movil_1(lm-(i*dec_factor)-N+1 : lm-(i*dec_factor)); %Rama pasaaltas
end

disp('High-pass analysis matrix: ');
W1

disp('High-pass branch results');
y1 = W1 * x

% -------------------------------------------------------------
% Low-pass several stage analysis matrix [linear convolution]

% Get an equivalent filter to get its length
hx = formfilter(n_stages, 1, h0, h1);
Lhx = length(hx);

% Compute the number of rows per branch
rowsperfilter = ceil((N+Lhx)/dec_fact);

% Allocate memory for the low-pass matrix
W0 = zeros((n_branches/2)*rowsperfilter,N);

tic; %Start cronometer

% Build low pass filter
for i = 0:(n_branches/2)-1
    % Get the equivalent filter depending on the number of stages
    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
    % Append zeros to the vector
    movil = [zeros(1,N-1),fliplr(hx),zeros(1,N-1)];
    lm = length(movil);
    % Shift and add the vector to the matrix
    for j = 0:rowsperfilter-1
        W0(i*rowsperfilter + j + 1,:) = movil(lm-(dec_fact*j)-N + 1 : lm-(dec_fact*j));
    end
end

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

disp('Low-pass filter bank results: ');
y0 = W0*x

%% STEP 4: Recover signal

disp('Recovered signal: ');
xn1 = (W1')*y1;
xn0 = (W0')*y0;
xn = xn1+xn0

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

fprintf('\nProcessing time = %d seg \n',toc);
fprintf('\nProgram end >>> \n');