In the previous section we found that the a Multiresolution Analysis allows us to decompose a signal into approximations and details.
On the theoretical level this is an Analysis-Synthesis situation. That is, we have bases and and we use these bases to decompose our signal.
On the practical level, we assume that our signal is represented by its approximation coefficients at some scale and we decompose it in terms of its coefficients at larger scale.
Both points of view are necessary for a real understanding of the subject.
In this section we will show that the approximation and detail coefficients can be computed using the filters previously mentioned. As we must compute these coefficients at many different scales we will need a filter bank.
In the Discrete Wavelet Transform (DWT) we have . That is, each signal in can be expressed in two ways using the basis functions in each of the spaces.
We start with the coefficients at scale index and produce the two sets of coefficients and at scale index (Analysis). Alternately, we can start with the two sets of coefficients and at scale index and produce the coefficients at scale index (Synthesis).
We can show that the two operations of Analysis and Synthesis are produced by certain filter banks.
As the wavelets and the scales at each index level are orthogonal we can compute the coefficients and by the usual inner product formula:
To complete this calculation we have to compute the inner product:
We can now complete the calculation previously started:
The detail coefficients can be computed similarily.
To complete this calculation we have to compute the inner product:
Upon substitution of this formula into the previous calculation we get:
The two formulas for the approximation and detail coefficients look similar to convolution but there is a downsampling involved.
Downsampling a discrete time signal consists of omitting every other value. We can think of a system whose input is and whose output is .
To understand the approximation and detail formulas it will help to define the time reversed filters and . We temporarily use to see the convolution.
If we follow this filter by the downsampler we get the approximation coefficients at the next level.
The same calculation holds for the detail coefficients. That is, convolution with the time reversed filter followed by downsampling produces the detail coefficients at the next level.
We really should think of the two filtering operations followed by downsampling as a filter bank.
We are analyzing a function in into a detail in and an approximation, , in , using a filter bank to calculate the coefficients and .
The dwt command in the Wavelet toolbox of Matlab carries out the calculations in this one stage filter bank. The syntax is:
[cA, cD] = dwt(x, Lo, Hi); = dwt(x, 'wname');
Note that the number of data values produced by the filter bank is about the same as the number of data values entering the system. To see this let, be the length of the input vector and assume that the filters both have length . The length of the convolution is so that the lengths of and are . The overall size of the data emerging from the filterbank is increased by the length of the filter minus 1.
We can further decompose to get:
We can then decompose to get:
The coefficients, and form can be calculated by iterataing or cascading the single stage filter bank to obtain a multiple stage filter bank.
The wavedec command in the Wavelet toolbox of Matlab carries out the calculations in the multi-stage filter bank. The syntax is:
[C, L] = wavedec(x, N, Lo, Hi); = wavedec(x, N, 'wname');
The vector is the set of details coefficients at each level together with the approximation coefficients at the final level.
The vector is the vector of lengths of each of the entries in and the length of (see the wavedec reference page).
Matlab provides functions, appcoef and detcoef compute any of the approximation or detail coefficients from the output of the wavedec command.
The decomposition of a signal into an approximation and a detail can be reversed. That is, we start with the two sets of coefficients and at scale index and produce the coefficients at scale index (Synthesis). We have:
Using the fact that is an orthogonal basis for we have:
This synthesis formula can be understood in terms of upsampling and filtering.
look like convolutions but upsampling is involved. Upsampling of a discrete time signal consists of inserting zeros between the values. We can think about a system with input and output for even values of and for odd values of .
The expressions for and consist of upsampling followed by filtering.
It follows that the synthesis formula consists of adding the outputs of the upsampled and filtered approximation and detail coefficients.
The idwt command in the Wavelet toolbox of Matlab carries out the calculations in this one stage filter bank. The syntax is:
x = idwt(cA, cD, Lo, Hi); = idwt(cA, cD, 'wname');
If we feed the output of the one-stage analysis filterbank to the input of the one-stage synthesis fitler bank then we get the original coefficients back. We say the we have a perfect reconstruction filter bank.
Note that the non-causal nature of the Analysis Bank filters is does not cause a problem in practice as we are dealing with FIR filters. This means that we can apply a fixed delay, to each filter to make it causal before applying the input signal. This is the same as delaying the input signal by before applying it to the filter bank.
To make things simple, assume that the filters and both have length . If we delay the time reversed filter by then we obtain the ``causal flip'' of the filter .
Denote the reversed filters by , . Filtering a signal by is the same a delaying the signal by then filtering by .
Hence, we can have a filter bank consisting of causal filters which gives perfect reconstruction with an overall delay of . The analysis filters are and while the synthesis filters are and .
The outputs of the multiple stage analysis filter bank can be fed into a multiple stage synthesis filter bank to reproduce the original coefficients. For example, a 3 level analysis bank produces outputs , , , and . These are fed into the 3 level synthesis filter bank as shown:
We have seen that we can reconstruct the signal in from the approximation and detail coefficients, and at level 1.
As and are in we can resolve them as:
Using the same reasoning as before:
That is, we obtain the approximation coefficients at level 0 by upsampling the approximation coefficients, at level 1 and then filtering with the low pass filter .
Similarily, we obtain the detail coefficients at level 0 by upsampling the detail coefficients, at level 1 and then filtering with the high pass filter .
We have previously analyzed the one-stage analysis filterbank to count
the number of coefficients produced.
We can also count the number of floating point operations involved. Assume that we are feeding in coefficients and the filter have length . Each convolution takes approximately operations. So, the system requires total floating point operations.
Now if we use two-stage filter bank then the input to the next bank has length (approximately) so that the second stage adds only operations to the first stage operations.
Continuing in this way, a multi-level bank requires:
That is, the computational complexity of the algorithm is linear in the size of the data.
Compare this to the FFT algorithm which has complexity . We really should not make too much of this apparent computational advantage as we know the CWT reflects a logarithmic division of frequency while the FFT uses equally spaced frequency divisions. The point is that for some applcations the logarithmic frquency divisions are sufficent for analyzing the situation so that the DWT has a computational advantage is these cases.
The m-file checkfb.m produced the following diary file.
echo on % m file to generate and test the analysis and synthesis banks % % | % +------+ +------+ | +------+ +------+ % | | y1 | | v1 | | | u1 | | w1 % +-----| rh1 |-----| dwn2 |----|----| up2 |-----| h1 |----+ % | | | | | | | | | | | % | +------+ +------+ | +------+ +------+ | %x | | | z %----+ | +---- % | | | % | +------+ +------+ | +------+ +------+ | % | | | y0 | | v0 | | | u0 | | w0 | % +-----| rh0 |-----| dwn2 |----|----| up2 |-----| h1 |----+ % | | | | | | | | | % +------+ +------+ | +------+ +------+ % | % Analysis Bank | Synthesis Bank % % % N = 2; % which db filters to use len = 2*N; % filter length p = 10; % size of input vectors % construct the string for dbN wave = ['db',int2str(N)]; % get some filters to work with [rh0 rh1 h0 h1] = wfilters(wave); % Have a look at the filters [rh0' rh1' h0' h1'] ans = -0.1294 -0.4830 0.4830 -0.1294 0.2241 0.8365 0.8365 -0.2241 0.8365 -0.2241 0.2241 0.8365 0.4830 -0.1294 -0.1294 -0.4830 x=randn(p,1); % white noise % first handle the analysis bank y0 = conv(rh0',x); y1 = conv(rh1',x); v0 = dyaddown(y0); v1 = dyaddown(y1); % print them out [y0 y1] ans = -0.0381 -0.1422 0.2389 0.8916 -0.1457 -1.5287 -1.0255 0.0748 0.4057 1.7050 1.4370 -1.4494 0.2355 0.0570 0.8709 1.7159 1.2926 -1.0294 -1.1244 -1.2863 -1.8471 1.0070 -0.2182 0.0585 0.2758 -0.0739 [v0 v1] ans = 0.2389 0.8916 -1.0255 0.0748 1.4370 -1.4494 0.8709 1.7159 -1.1244 -1.2863 -0.2182 0.0585 % now do the synthesis bank u0 = dyadup(v0); u1 = dyadup(v1); w0 = conv(h0,u0); w1 = conv(h1,u1); % print them out [u0 u1] ans = 0 0 0.2389 0.8916 0 0 -1.0255 0.0748 0 0 1.4370 -1.4494 0 0 0.8709 1.7159 0 0 -1.1244 -1.2863 0 0 -0.2182 0.0585 0 0 [w0 w1] ans = 0 0 0.1154 -0.1154 0.1998 -0.1998 -0.4417 0.7362 -0.8888 -0.4474 0.4642 0.2502 1.3348 0.2887 0.7427 -1.4345 0.5426 0.3154 -0.3478 1.6018 -1.0533 -0.5404 -0.3574 -1.0836 -0.0370 0.6081 -0.0489 0.0489 0.0282 -0.0282 0 0 % now add up the streams z = w0 + w1; % to compare x to z we pad the end with zeros cx = [x ; zeros(2*(len-1),1)]; [cx z] ans = 0.2944 0 -1.3362 0.0000 0.7143 -0.0000 1.6236 0.2944 -0.6918 -1.3362 0.8580 0.7143 1.2540 1.6236 -1.5937 -0.6918 -1.4410 0.8580 0.5711 1.2540 0 -1.5937 0 -1.4410 0 0.5711 0 -0.0000 0 0.0000 0 0 >> diary off
The m-file wavdemo.m produced the following plots.
% File: wavdemo.m % demonstrates various commands in the wavelet toolbox load sumsin x=sumsin; % Compute the wavelet decomposition at level 3 [C,L] = wavedec(x,3,'db2'); figure(1) subplot(2,1,1); plot(x); title('original signal');axis([0 1000 -5 5]); subplot(2,1,2); plot(C); title('wavelet decomposition');axis([0 1000 -5 5]); orient landscape; print wavdemofig1.eps % Now extract the coefficients cA3 = appcoef(C, L, 'db2', 3); cD3 = detcoef(C, L, 3); cD2 = detcoef(C, L, 2); cD1 = detcoef(C, L, 1); % plot these figure(2) subplot(4,2,1); plot(cA3); title('cA3');axis([0 1000 -5 5]); subplot(4,2,3); plot(cD3); title('cD3');axis([0 1000 -5 5]); subplot(4,2,5); plot(cD2); title('cD2');axis([0 1000 -5 5]); subplot(4,2,7); plot(cD1); title('cD1');axis([0 1000 -5 5]); % Now compute the reconstructed coefficients A3 = wrcoef('a', C, L, 'db2', 3); D3 = wrcoef('d', C, L, 'db2', 3); D2 = wrcoef('d', C, L, 'db2', 2); D1 = wrcoef('d', C, L, 'db2', 1); % plot these subplot(4,2,2); plot(A3); title('A3');axis([0 1000 -5 5]); subplot(4,2,4); plot(D3); title('D3');axis([0 1000 -5 5]); subplot(4,2,6); plot(D2); title('D2');axis([0 1000 -5 5]); subplot(4,2,8); plot(D1); title('D1');axis([0 1000 -5 5]); orient landscape; print wavdemofig2.eps
The function examine_recon.m can be useful in examining the reconstruction or synthesis of a signal.
It was used to produce the plots of approximations and details for the speech signal in the multiresolution section.
function [cA,A,cD,D] = examine_recon(x, N, wname) %[cA,A,cD,D] = examine_recon(x,N,wname) calculate and plot approx and details %up to level N for a signal x using wavelet type given by wname [m,n]=size(x); if m == 1 x=x'; end % perform wavelet decompostion at level N [C,L] = wavedec(x,N,wname); A = zeros(length(x),N); D = zeros(length(x),N); cA = zeros(length(x),N); cD = zeros(length(x),N); lena = zeros(1,N); lend = zeros(1,N); % compute the approximations at various levels for i=1:N A(:,i) = wrcoef('a',C,L,wname,i); end % compute the approximation coeff at various levels for i=1:N temp = appcoef(C,L,wname,i); lena(i)=length(temp); cA(1:lena(i),i) = temp; end % compute the details at various levels for i=1:N D(:,i) = wrcoef('d',C,L,wname,i); end % compute the detail coeff at various levels for i=1:N temp = detcoef(C,L,i); lend(i)=length(temp); cD(1:lend(i),i) = temp; end % make plots figure; subplot(N+1,1,1); plot(x); title(['Signal and approximations 1 to ', num2str(N)]); for i=1:N subplot(N+1,1,i+1); plot(A(:,i)); end figure; subplot(N+1,1,1); plot(x); title(['Signal and approximation coef 1 to ', num2str(N)]); for i=1:N subplot(N+1,1,i+1); plot(cA(1:lena(i),i)); end figure; subplot(N+1,1,1); plot(x); title(['Signal and details 1 to ', num2str(N)]); for i=1:N subplot(N+1,1,i+1); plot(D(:,i)); end figure; subplot(N+1,1,1); plot(x); title(['Signal and detail coef 1 to ', num2str(N)]); for i=1:N subplot(N+1,1,i+1); plot(cD(1:lend(i),i)); end
Given a signal , to apply the Discrete Wavelet Transform our starting point must be a signal in one of the spaces . Since the spaces, , are getting larger and larger as goes to infinity we can approximate our signal, , closely by choosing a large enough value of and projecting the signal into using the orthogonal basis , (all values of ). That is:
The projection coefficients, , in this summation can be obtained by sampling the output of a continuous time filter. That is,
This shows that we obtain the coordinates of the projection of onto the space by passing through the filter with impulse response and then sampling at times .
We will show that if is bandlimited then by choosing an appropriate sampling rate and a scale index, , we can recognize the projection coefficents in terms of a discrete time filtering process.
Assume that we have sampled a band limited signal at 20% above the Nyquist rate. After normalization we can assume that the sampling rate is 1Hz so that our signal contains only frequencies below 0.5/1.2 = 0.417. That is, we can represent our signal as
Where the Discrete Fourier Transform, , of is zero for .
We can now pass this signal into the above filter to find the projection onto . After rearranging the summation and itegration we find:
Where evaluated at . Hence the coordinates of the projection of our signal onto are just the samples, , filtered by .
Fact: For almost all orthogonal wavelets the filter is approximately magnitude 1 over frequencies and the phase response is linear. That is, the filter acts as a delay on our signal and so aside from that delay, we may use the samples of the signal as the scaling coefficients.
Example: The m-file initdwt.m can be used to investigate the filter for various wavelets.
% This m file investigates the filter which carries out the projection of % a band limited signal into V0 % % This filter is alpha(n) = sinc(r) * phi(-r), at r = n % % alpha(n) = integral phi(-s) sinc(n-s) ds % = integral phi(t) sinc(t+n) dt % % Since the support of phi(t) is 0 <= t <= N-1 the values of this % integral will become quite small outside the interval -2N <= n <= 2N-1. % So, we approximate alpha(n) by an FIR filter of length 4N. iter = 10; dt = 1/2^iter; wave = 'db7'; [f0 f1 h0 h1] = wfilters(wave); [phi, psi, t] = wavefun(wave,iter); N = length(h0); alpha = zeros(1,4*N); for n=(-2*N:(2*N-1)) alpha(n+2*N+1) = sinc(t + n) * phi' * dt; end subplot(3,1,1) ;stem(alpha); title(['Projection Filter for ' wave]); [Alpha,W]=freqz(alpha); subplot(3,1,2); plot(W/(2*pi), abs(Alpha)); axis([0 0.5 0 1.2]); title('Magnitude Response'); subplot(3,1,3); plot(W/(2*pi), unwrap(angle(Alpha))); title('Phase Response');