5. Filter Banks and the Discrete Wavelet Transform

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 $\phi_{j,k}(t)$ and $w_{j,k}(t)$ 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 $1/2^{J}$ 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.

5.1 Analysis: From Fine Scale to Coarser Scale

In the Discrete Wavelet Transform (DWT) we have $V_{j} = V_{j-1} + W_{j-1}$. That is, each signal $x(t)$ in $V_{j}$ can be expressed in two ways using the basis functions in each of the spaces.

$$\begin{eqnarray*} x(t) &=& \sum_{k} cA_{0}(k) \phi_{j,k}(t) \\ &=& \sum_{k} cA_{1}(k) \phi_{j-1,k}(t) + \sum_{k} cD_{1}(k) w_{j-1,k}(t) \end{eqnarray*}$$

We start with the coefficients $A_{0}(k)$ at scale index $j$ and produce the two sets of coefficients $A_{1}(k)$ and $D_{1}(k)$ at scale index $j-1$ (Analysis). Alternately, we can start with the two sets of coefficients $A_{1}(k)$ and $D_{1}(k)$ at scale index $j-1$ and produce the coefficients $A_{0}(k)$ at scale index $j$ (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 $cA_{1}(k)$ and $D_{1}(k)$ by the usual inner product formula:

$$\begin{eqnarray*} cA_{1}(k) &=& \left< x(t), \phi_{j-1,k}(t) \right>\\ &=& \l... ...\sum_{n} cA_{0}(n) \left< \phi_{j,n}(t), \phi_{j-1,k}(t) \right> \end{eqnarray*}$$

To complete this calculation we have to compute the inner product:

$$\begin{eqnarray*} \left< \phi_{j,n}(t), \phi_{j-1,k}(t) \right> &=& \int_{-\inf... ... ds \mbox{ (integral is 0 unless } m = n-2k) \\ &=& h_{0}(n-2k) \end{eqnarray*}$$

We can now complete the calculation previously started:

$$\begin{eqnarray*} cA_{1}(k) &=& \sum_{n} h_{0}(n-2k) cA_{0}(n) \end{eqnarray*}$$

The detail coefficients can be computed similarily.

$$\begin{eqnarray*} cD_{1}(k) &=& \left< x(t), w_{j-1,k}(t) \right>\\ &=& \left... ...=& \sum_{n} cA_{0}(n) \left< \phi_{j,n}(t), w_{j-1,k}(t) \right> \end{eqnarray*}$$

To complete this calculation we have to compute the inner product:

$$\begin{eqnarray*} \left< \phi_{j,n}(t), w_{j-1,k}(t) \right> &=& \int_{-\infty}... ... ds \mbox{ (integral is 0 unless } m = n-2k) \\ &=& h_{1}(n-2k) \end{eqnarray*}$$

Upon substitution of this formula into the previous calculation we get:

$$\begin{eqnarray*} cD_{1}(k) &=& \sum_{n} h_{1}(n-2k) cA_{0}(n) \end{eqnarray*}$$

5.1.1 Filtering and Downsampling

The two formulas for the approximation and detail coefficients look similar to convolution but there is a downsampling involved.

$$\begin{eqnarray*} cA_{1}(k) &=& \sum_{n} h_{0}(n-2k) cA_{0}(n)\\ cD_{1}(k) &=& \sum_{n} h_{1}(n-2k) cA_{0}(n) \end{eqnarray*}$$

Downsampling a discrete time signal $x(n)$ consists of omitting every other value. We can think of a system whose input is $x(n)$ and whose output is $y(n) = x(2n)$.

To understand the approximation and detail formulas it will help to define the time reversed filters $\tilde{h}_{0}(n) = h_{0}(-n)$ and $\tilde{h}_{1}(n) = h_{1}(-n)$. We temporarily use $m=2k$ to see the convolution.

$$\begin{eqnarray*} u(m) &=& \sum_{n} h_{0}(n-m) cA_{0}(n) \\ &=& \sum_{n} \tilde{h}_{0}(m-n) cA_{0}(n)\\ &=& \tilde{h}_{0}(m) * cA_{0}(m) \end{eqnarray*}$$

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 $h_{1}(-n)$ followed by downsampling produces the detail coefficients at the next level.

5.1.2 The One-Stage Analysis Filter Bank

We really should think of the two filtering operations followed by downsampling as a filter bank.

We are analyzing a function $x(t)$ in $V_{j}$ into a detail $D_{1}(t)$ in $W_{j-1}$ and an approximation, $A_{1}(t)$, in $V_{j-1}$, using a filter bank to calculate the coefficients $cD_{1}(k)$ and $cA_{1}(k)$.

$$\begin{eqnarray*} x(t) &=& \sum_{k} cA_{0}(k) \phi_{j,k}(t) \\ &=& \sum_{k} cA... ...) + \sum_{k} cD_{1}(k) w_{j-1,k}(t)\\ &=& A_{1}(t) + D_{1}(t) \end{eqnarray*}$$

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, $lx$ be the length of the input vector and assume that the filters both have length $lf$. The length of the convolution is $lx + lf -1$ so that the lengths of $cA$ and $cD$ are $\mbox{floor}((lx + lf-1)/2)$. The overall size of the data emerging from the filterbank is increased by the length of the filter minus 1.

5.1.3 The Analysis Filter Bank

We can further decompose $A_{1}(t)$ to get:

$$\begin{eqnarray*} x(t) &=& A_{1}(t) + D_{1}(t)\\ &=& \sum_{k} cA_{2}(k) \phi_... ...k} cD_{1}(k) w_{j-1,k}(t)\\ &=& A_{2}(t) + D_{2}(t) + D_{1}(t) \end{eqnarray*}$$

We can then decompose $A_{2}(t)$ to get:

$$\begin{eqnarray*} x(t) &=& A_{1}(t) + D_{1}(t)\\ &=& \sum_{k} cA_{2}(k) \phi_... ...sum_{k} cD_{2}(k) w_{j-2,k}(t) + \sum_{k} cD_{1}(k) w_{j-1,k}(t) \end{eqnarray*}$$

The coefficients, $cA_{m}{k}$ and $cD_{m}(k)$ form $m=1,2,3$ 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 $C$ is the set of details coefficients at each level together with the approximation coefficients at the final level.

$$C = \left[ cA_{N}, cD_{N}, cD_{N-1}, \cdots , cD_{2}, cD_{1} \right]$$

The vector $L$ is the vector of lengths of each of the entries in $C$ and the length of $x$ (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.

5.2 Synthesis: From Course Scale to Fine Scale

The decomposition of a signal into an approximation and a detail can be reversed. That is, we start with the two sets of coefficients $A_{1}(k)$ and $D_{1}(k)$ at scale index $j-1$ and produce the coefficients $A_{0}(k)$ at scale index $j$ (Synthesis). We have:

$$\begin{eqnarray*} x(t) &=& \sum_{k} cA_{0}(k) \phi_{j,k}(t) \\ &=& \sum_{k} cA... ...) + \sum_{k} cD_{1}(k) w_{j-1,k}(t)\\ &=& A_{1}(t) + D_{1}(t) \end{eqnarray*}$$

Using the fact that $\phi_{j,n}(t)$ is an orthogonal basis for $V_{j}$ we have:

$$\begin{eqnarray*} cA_{0}(n) &=& \left< x(t), \phi_{j,n}(t) \right>\\ &=& \lef... ... \sum_{k} cA_{1}(k) h_{0}(n-2k) + \sum_{k} cD_{1}(k) h_{1}(n-2k) \end{eqnarray*}$$

This synthesis formula can be understood in terms of upsampling and filtering.

5.2.1 Upsampling and Filtering

The expressions:

$$\begin{eqnarray*} u(n) &=& \sum_{k} cA_{1}(k) h_{1}(n-2k)\\ v(n) &=& \sum_{k} cD_{1}(k) h_{0}(n-2k) \end{eqnarray*}$$

look like convolutions but upsampling is involved. Upsampling of a discrete time signal $x(n)$ consists of inserting zeros between the values. We can think about a system with input $x(n)$ and output $y(n) = x(n/2)$ for even values of $n$ and $y(n) = 0$ for odd values of $n$.

The expressions for $u(n)$ and $v(n)$ consist of upsampling followed by filtering.

5.2.2 The One-Stage Synthesis Filter Bank

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');

5.2.3 Perfect Reconstruction Filter Bank

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, $N$ to each filter to make it causal before applying the input signal. This is the same as delaying the input signal by $N$ before applying it to the filter bank.

To make things simple, assume that the filters $h_{0}(n)$ and $h_{1}(n)$ both have length $N+1$. If we delay the time reversed filter $h_{i}(-n)$ by $N$ then we obtain the ``causal flip'' of the filter $h_{i}(n)$.

Denote the reversed filters by $rh_{i}(n) = h_{i}(N-n)$, $i=0,1$. Filtering a signal $x(n)$ by $rh_{i}(n)$ is the same a delaying the signal by $N$ then filtering by $h_{i}(-n)$.

Hence, we can have a filter bank consisting of causal filters which gives perfect reconstruction with an overall delay of $N$. The analysis filters are $rh_{1}(n)$ and $rh_{0}(n)$ while the synthesis filters are $h_{1}(n)$ and $h_{0}(n)$.

5.2.4 The Synthesis Filter Bank

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 $cD_{1}(k)$, $cD_{2}(k)$, $cD_{3}(k)$, and $cA_{3}(k)$. These are fed into the 3 level synthesis filter bank as shown:

5.2.5 Approximations and Details

We have seen that we can reconstruct the signal $x(t)$ in $V_{j}$ from the approximation and detail coefficients, $cA_{1}(k)$ and $cD_{1}(k)$ at level 1.

$$\begin{eqnarray*} x(t) &=& \sum_{k} cA_{0}(k) \phi_{j,k}(t) \\ &=& \sum_{k} cA... ...) + \sum_{k} cD_{1}(k) w_{j-1,k}(t)\\ &=& A_{1}(t) + D_{1}(t) \end{eqnarray*}$$

As $A_{1}(t)$ and $D_{1}(t)$ are in $V_{j}$ we can resolve them as:

$$\begin{eqnarray*} A_{1}(t) &=& \sum_{k} cA_{1}(k) \phi_{j-1,k}(t)\\ &=& \sum_{... ...} cD_{1}(k) w_{j-1,k}(t)\\ &=& \sum_{k} D_{1}(k) \phi_{j,k}(t) \end{eqnarray*}$$

Using the same reasoning as before:

$$\begin{eqnarray*} A_{1}(n) &=& \sum_{k} cA_{1}(k) h_{0}(n-2k)\\ & &\\ D_{1}(n) &=& \sum_{k} cD_{1}(k) h_{1}(n-2k) \end{eqnarray*}$$

That is, we obtain the approximation coefficients at level 0 by upsampling the approximation coefficients, $cA_{1}(k)$ at level 1 and then filtering with the low pass filter $h_{0}(k)$.

Similarily, we obtain the detail coefficients at level 0 by upsampling the detail coefficients, $cD_{1}(k)$ at level 1 and then filtering with the high pass filter $h_{1}(k)$.

5.3 Numerical Complexity of the Discrete Wavelet Transform

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 $n$ coefficients and the filter have length $f$. Each convolution takes approximately $f n/2$ operations. So, the system requires $f n$ total floating point operations.

Now if we use two-stage filter bank then the input to the next bank has length (approximately) $n/2$ so that the second stage adds only $f n/2$ operations to the first stage operations.

Continuing in this way, a multi-level bank requires:

$$\mbox{Number of Operations } = f n + \frac{f n}{2} + \frac{f ... ... \left( 1 + \frac{1}{2} + \frac{1}{4} + \cdots \right) = 2 f n$$

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 $n \mbox{log}_{2}(n)$. 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.

5.4 Matlab Examples

5.4.1 One-Stage Perfect Reconstruction

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

5.4.2 Approximations and Details

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

5.4.3 A Useful Function

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

5.5 Initialization of the Discrete Wavelet Transform

Given a signal $x(t)$, to apply the Discrete Wavelet Transform our starting point must be a signal in one of the spaces $V_{j}$. Since the spaces, $V_{j}$, are getting larger and larger as $j$ goes to infinity we can approximate our signal, $x(t)$, closely by choosing a large enough value of $j$ and projecting the signal into $V_{j}$ using the orthogonal basis $\phi_{j,k}(t)$, (all values of $k$). That is:

$$x(t) \approx \sum_{k} a_{j,k} \phi_{j,k}(t) = \sum_{k} \left< x(t), \phi_{j,k}(t) \right> \phi_{j,k}(t)$$

The projection coefficients, $a_{j,k}$, in this summation can be obtained by sampling the output of a continuous time filter. That is,

$$\begin{eqnarray*} a_{j,k} &=& \left< x(t), \phi_{j,k}(t) \right> \\ &=& \int_... ...t) * \sqrt{2^{j}} \phi(-2^{j} t) \mbox{ evaluated at } t=2^{-j}k \end{eqnarray*}$$

This shows that we obtain the coordinates of the projection of $x(t)$ onto the space $V_{j}$ by passing $x(t)$ through the filter with impulse response $\sqrt{2^{j}} \phi(-2^{j}t)$ and then sampling at times $t=k/2^{j}$.

We will show that if $x(t)$ is bandlimited then by choosing an appropriate sampling rate and a scale index, $j$, 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

$$x(t) = \sum_{n} x(n) \mbox{sinc}(t-n)$$

Where the Discrete Fourier Transform, $X(f)$, of $x(n)$ is zero for $f > 0.417$.

We can now pass this signal into the above filter to find the projection onto $V_{0}$. After rearranging the summation and itegration we find:

$$\begin{eqnarray*} \left< x(t), \phi(t-k) \right> &=& \int_{-\infty}^{\infty} x(... ...ight)\\ &=& \sum_{n} x(n) \alpha(k-n)\\ &=& x(n) * \alpha(n) \end{eqnarray*}$$

Where $\alpha(n) = \mbox{sinc}(t) * \phi(-t)$ evaluated at $t=n$. Hence the coordinates of the projection of our signal onto $V_{0}$ are just the samples, $x(n)$, filtered by $\alpha(n)$.

Fact: For almost all orthogonal wavelets the filter $\alpha(n)$ is approximately magnitude 1 over frequencies $f < 0.4$ 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 $\alpha(n)$ 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');