% Program to compute FDWT
clear
f= [9;1;2;0];

J=2;
a(:,J+1)=f;

%Wavelet basis matrix
fprintf('The wavelet basis vectors are:\n');
T=[1 1 1 1; 1 1 -1 -1; 1 -1 0 0;0 0 1 -1]'

%Lowpass and Highpass filters
L=[.5 .5 0 0; 0 0 .5 .5];
H=[.5 -.5 0 0; 0 0 .5 -.5];


for k=J+1:-1:2,
	a(:,k-1)=[L*a(:,k);0;0]
	b(:,k-1)=[H*a(:,k);0;0]
%blur=[1 0; 1 0; 0 1; 0 1]*(L*a(:,k))
end

fprintf('The wavelet coefficients are:\n')
B= [a(1,1) b(1,1) b(1:2,2)']


figure(1)
subplot(421)
stem(f)
title('original signal f=a2')



subplot(424)
stem(T(:,3)*B(3) + T(:,4)*B(4))
title('detail lost in a2 to a1')

subplot(426)
stem(T(:,2)*B(2))
title('detail lost in a1 to a0')

subplot(423)
stem(T(:,2)*B(2)+T(:,1)*B(1))
title('blurred signal a1')

subplot(425)
stem(T(:,1)*B(1))
title('blurred signal a0')
