Coding wavelets the easy way

% clear workspace
clear all
% Prepare data
load PieceRegSig;
a0.s = sig;
a0.d = 0;
% Prepare the filters
load daub4;
h.s = h;
h.d = 0; % degree is 3
g.s = g;
g.d = 0;
rh.s = rh;
rh.d = 3; % degree is -3
rg.s = rg;
rg.d = 3;
% compute the decomposition
% LTI filters first
a1 = myconv(a0,h);
d1 = myconv(a0,g);
% subsample
a1 = subsample(a1);
d1 = subsample(d1);
% we are now in the middle of the diagram and have computed the signals a1
% and d1
% Reconstrauction
% insert zeros first
a1 = zeroinsert(a1);
d1 = zeroinsert(d1);
% LTI filters now
a1 = myconv(a1,rh);
d1 = myconv(d1,rg);
% add the signals
a2 = addSignals(a1,d1);
% restrict the reconstruction to the support of a0
l = length(a0.s);
a3 = a2.s(a2.d + 1:end);
a3 = a3(1:l);
figure(1);
plot(sig-a3);
% now check the difference of the support of a2 (which includes the support
% of a0)
a4 = [ zeros(1,a2.d) sig];
extraLa4 = length(a2.s) - length(a4);
a4 = [a4 zeros(1,extraLa4)];
figure(2);
plot(a4-a2.s);