% 1-D signal analysis

% biorwavf generates symmetric biorthogonal wavelet filters.
% The argument has the form biorNr.Nd, where
%    Nr = number of zeros at pi in the synthesis lowpass filter, s[n].
%    Nd = number of zeros at pi in the analysis lowpass filter, a[n].
% We find the famous Daubechies 9/7 pair, which have Nr = Nd = 4.  
% Note: In earlier versions of the Matlab Wavelet Toolbox (v2.1 and 
% below,) the vectors s and a are zero-padded to make their lengths equal:
%    a[-4] a[-3] a[-2] a[-1] a[0] a[1] a[2] a[3] a[4]
%      0   s[-3] s[-2] s[-1] s[0] s[1] s[2] s[3]   0
[s,a]= biorwavf('bior4.4');  

% Find the zeros and plot them.
close all
clf
fprintf(1,'Zeros of H0(z)')
roots(a)
subplot(1,2,1)
zplane(a)
title('Zeros of H0(z)')

fprintf(1,'Zeros of F0(z)')
roots(s)
subplot(1,2,2)
zplane(s)          % Note: there are actually 4 zeros clustered at z = -1.
title('Zeros of F0(z)')
pause


% Determine the complete set of filters, with proper alignment.
% Note: Matlab uses the convention that a[n] is the flip of h0[n].
%    h0[n] = flip of a[n], with the sum normalized to sqrt(2).
%    f0[n] = s[n], with the sum normalized to sqrt(2).
%    h1[n] = f0[n], with alternating signs reversed (starting with the first.)
%    f1[n] = h0[n], with alternating signs reversed (starting with the second.)
[h0,h1,f0,f1] = biorfilt(a, s);

clf
subplot(2,2,1)
stem(0:8,h0(2:10))
ylabel('h0[n]')
xlabel('n')
subplot(2,2,2)
stem(0:6,f0(2:8))
ylabel('f0[n]')
xlabel('n')
v = axis; axis([v(1) 8 v(3) v(4)])
subplot(2,2,3)
stem(0:6,h1(2:8))
ylabel('h1[n]')
xlabel('n')
v = axis; axis([v(1) 8 v(3) v(4)])
subplot(2,2,4)
stem(0:8,f1(2:10))
ylabel('f1[n]')
xlabel('n')
pause


% Examine the Frequency response of the filters.
N = 512;
W = 2/N*(-N/2:N/2-1);
H0 = fftshift(fft(h0,N));
H1 = fftshift(fft(h1,N));
F0 = fftshift(fft(f0,N));
F1 = fftshift(fft(f1,N));
clf
plot(W, abs(H0), '-', W, abs(H1), '--', W, abs(F0), '-.', W, abs(F1), ':')
title('Frequency responses of Daubechies 9/7 filters')
xlabel('Angular frequency (normalized by pi)')
ylabel('Frequency response magnitude')
legend('H0', 'H1', 'F0', 'F1', 0)
pause


% Load a test signal. 
load noisdopp
x = noisdopp;
L = length(x);
clear noisdopp

% Compute the lowpass and highpass coefficients using convolution and
% downsampling.
y0 = dyaddown(conv(x,h0));
y1 = dyaddown(conv(x,h1));

% The function dwt provides a direct way to get the same result.
[yy0,yy1] = dwt(x,'bior4.4');

% Now, reconstruct the signal using upsamping and convolution.  We only
% keep the middle L coefficients of the reconstructed signal i.e. the ones
% that correspond to the original signal.
xhat = conv(dyadup(y0),f0) + conv(dyadup(y1),f1);
xhat = wkeep(xhat,L);

% The function idwt provides a direct way to get the same result.
xxhat = idwt(y0,y1,'bior4.4');

% Plot the results.
subplot(4,1,1);
plot(x)
axis([0 1024 -12 12])
title('Single stage wavelet decomposition')
ylabel('x')
subplot(4,1,2);
plot(y0)
axis([0 1024 -12 12])
ylabel('y0')
subplot(4,1,3);
plot(y1)
axis([0 1024 -12 12])
ylabel('y1')
subplot(4,1,4);
plot(xhat)
axis([0 1024 -12 12])
ylabel('xhat')
pause

% Next, we perform a three level decomposition.  The following
% code draws the structure of the iterated analysis filter bank.
clf
t = wtree(x,3,'bior4.4');
plot(t)
pause
close(2)

% For a multilevel decomposition, we use wavedec instead of dwt.
% Here we do 3 levels.  wc is the vector of wavelet transform
% coefficients.  l is a vector of lengths that describes the
% structure of wc.
[wc,l] = wavedec(x,3,'bior4.4');

% We now need to extract the lowpass coefficients and the various
% highpass coefficients from wc.
a3 = appcoef(wc,l,'bior4.4',3);
d3 = detcoef(wc,l,3);
d2 = detcoef(wc,l,2);
d1 = detcoef(wc,l,1);

clf
subplot(5,1,1)
plot(x)
axis([0 1024 -22 22])
ylabel('x')
title('Three stage wavelet decomposition')
subplot(5,1,2)
plot(a3)
axis([0 1024 -22 22])
ylabel('a3')
subplot(5,1,3)
plot(d3)
axis([0 1024 -22 22])
ylabel('d3')
subplot(5,1,4)
plot(d2)
axis([0 1024 -22 22])
ylabel('d2')
subplot(5,1,5)
plot(d1)
axis([0 1024 -22 22])
ylabel('d1')
pause

% We can reconstruct each branch of the tree separately from the individual
% vectors of transform coefficients using upcoef.
ra3 = upcoef('a',a3,'bior4.4',3,1024);
rd3 = upcoef('d',d3,'bior4.4',3,1024);
rd2 = upcoef('d',d2,'bior4.4',2,1024);
rd1 = upcoef('d',d1,'bior4.4',1,1024);

% The sum of these reconstructed branches gives the full recontructed
% signal.
xhat = ra3 + rd3 + rd2 + rd1;

clf
subplot(5,1,1)
plot(x)
axis([0 1024 -10 10])
ylabel('x')
title('Individually reconstructed branches')
subplot(5,1,2)
plot(ra3)
axis([0 1024 -10 10])
ylabel('ra3')
subplot(5,1,3)
plot(rd3)
axis([0 1024 -10 10])
ylabel('rd3')
subplot(5,1,4)
plot(rd2)
axis([0 1024 -10 10])
ylabel('rd2')
subplot(5,1,5)
plot(rd1)
axis([0 1024 -10 10])
ylabel('rd1')
pause

clf
plot(xhat-x)
title('Reconstruction error (using upcoef)')
axis tight
pause

% We can also reconstruct individual branches from the full vector of
% transform coefficients, wc.
rra3 = wrcoef('a',wc,l,'bior4.4',3);
rrd3 = wrcoef('d',wc,l,'bior4.4',3);
rrd2 = wrcoef('d',wc,l,'bior4.4',2);
rrd1 = wrcoef('d',wc,l,'bior4.4',1);
xxhat = rra3 + rrd3 + rrd2 + rrd1;

clf
plot(xxhat-x)
title('Reconstruction error (using wrcoef)')
axis tight
pause

% To reconstruct all branches at once, use waverec.
xxxhat = waverec(wc,l,'bior4.4');

clf
plot(xxxhat-x)
axis tight
title('Reconstruction error (using waverec)')
pause