% 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