6. Properties of the Filters, and the Scale and Wavelet Functions

In this section we will examine properties of the filters, scale and wavelet using the two-scale equation:

$$\phi(t) = \sum_{k} h_{0}(k) \sqrt{2} \phi(2t-k)$$

The fact that the scale function and the low pass filter satisfy this equation has some remarkable consequences.

Keep in mind that we are making a number of assumptions:

1. Orthogonality for the Wavelet: The baby wavelets are orthgonal and to have unit energy.

   
   

   $$\int_{-\infty}^{\infty} w_{j,k}(t) w_{m,n}(t)dt = \delta(j-m) \delta(k-n)$$
   
   

2. Orthogonality for the Scale: The translates $\phi(t-k)$, $-\infty < k < \infty$, are orthogonal and have unit energy ( $<\phi(t), \phi(t)> = 1$).

   
   

   $$\int_{-\infty}^{\infty} \phi(t)\phi(t-n) dt = \delta(n)$$
   
   

3. Completness: The translates $\phi(t-k)$, $-\infty < k < \infty$, span the same space as the wavelets $w_{j,k}(t)$, $-\infty < j < 0$ and $-\infty < k < \infty$.

6.1 Double Shift Orthogonality of the Filters

The two-scale equation shows that coordinates of $\phi(t)$ with respect to the basis $\sqrt{2} \phi(2t - k)$ are the filter coefficients $h_{0}(k)$.

If we delay $\phi(t)$ by n then the coordinates of $\phi(t-n)$ are the filter coefficients delayed by $2n$ (i.e. $h(k-2n)$. This can be easily seen by substituting $t-n$ for $t$ in the two-scale equation.

$$\begin{eqnarray*} \phi(t-n) &=& \sum_{k} h_{0}(k) \sqrt{2} \phi(2(t-n) - k)\\ ... ...(2t - (m))\\ &=& \sum_{k} h_{0}(k-2n) \sqrt{2} \phi(2t - k)\\ \end{eqnarray*}$$

The same is true of the wavelet function and the filter $h_{1}(k)$.

$$\begin{eqnarray*} w(t) &=& \sum_{k} h_{1}(k) \sqrt{2} \phi(2t - k)\\ & & \\ w(t-n) &=& \sum_{k} h_{1}(k-2n) \sqrt{2} \phi(2t - k)\\ \end{eqnarray*}$$

Since the integer time shifts of the wavelet and the scale are orthogonal we have:

$$\begin{array}{lcccl} \displaystyle \int_{-\infty}^{\infty} \... ...0 &=& \displaystyle \sum_{k} h_{0}(k) h_{1}(k-2n) \end{array}$$

These equations are called the double shift orthogonality relations of the filters. They lead to a number of other properties of the filters.

1. If we take $n=0$ in the equations above we get:

   
   

   $$\sum_{k} h_{0}(k)^{2} = 1 \mbox{ and } \sum_{k}h_{1}(k)^{2} = 1$$
   
   

2. If we integrate both sides of the two-scale equations we get:

   
   

   $$\sum_{k} h_{0}(k) = \sqrt{2} \mbox{ and } \sum_{k}h_{1}(k) = 0$$
   
   

3. The even and odd sums of both filters are:

   
   

   $$\sum_{k} h_{0}(2k) = \sum_{k} h_{0}(2k+1) = \frac{1}{\sqrt{2}}$$
   
   

   
   

   $$\sum_{k} h_{1}(2k) = -\sum_{k} h_{1}(2k+1) = \pm \frac{1}{\sqrt{2}}$$
   
   

   For a proof of this equation see Proof 3 in Appendix A, page 247 of the textbook (WWT).

4. The filter $h_{1}(k)$ is an alternating flip of the filter $h_{0}(k)$. That is, there is an odd integer $N$ such that:

   
   

   $$h_{1}(k) = (-1)^{k} h_{0}(N-k)$$
   
   

   For a proof of this fact see Proof 10 in Appendix A, page 250 of the textbook (WWT).

Example

To illustrate these equations we can use the m-file dshift.m as follows:

[h0 h1] = wfilters('db3','r')

h0 =

    0.3327    0.8069    0.4599   -0.1350   -0.0854    0.0352


h1 =

    0.0352    0.0854   -0.1350   -0.4599    0.8069   -0.3327

sum(h0)

ans =

    1.4142

sum(h1)

ans =

  -2.9981e-12

h0even = dyaddown(h0,1)

h0even =

    0.3327    0.4599   -0.0854

h0odd = dyaddown(h0,0)

h0odd =

    0.8069   -0.1350    0.0352

sum(h0even)

ans =

    0.7071

sum(h0odd)

ans =

    0.7071

h1even = dyaddown(h1,1)

h1even =

    0.0352   -0.1350    0.8069

h1odd = dyaddown(h1,0)

h1odd =

    0.0854   -0.4599   -0.3327

sum(h1even)

ans =

    0.7071

sum(h1odd)

ans =

   -0.7071

p0 = conv(h0,wrev(h0))

p0 =

  Columns 1 through 7 

    0.0117   -0.0000   -0.0977    0.0000    0.5859    1.0000    0.5859

  Columns 8 through 11 

    0.0000   -0.0977   -0.0000    0.0117

p1 = conv(h1, wrev(h1))

p1 =

  Columns 1 through 7 

   -0.0117   -0.0000    0.0977    0.0000   -0.5859    1.0000   -0.5859

  Columns 8 through 11 

    0.0000    0.0977   -0.0000   -0.0117

>> diary off

6.2 Frequency Domain Formulas

In the frequency domain we have a number of properties.

1. Frequency Domain Version of the Two-Scale Equation:

   
   

   $$\Phi(f) = \frac{1}{\sqrt{2}} H_{0}(f/2) \Phi(f/2)$$
   
   

   This is proven by taking the Fourier Transform of both sides of the two-scale equation.

   Note that this is a very unique equation because $\Phi(t)$ is the continuous time transform while $H(f)$ is the discrete time Transform.

2. If we iterate the previous equation we get:

   
   

   $$\Phi(f) = \Phi(0) \prod_{k=1}^{\infty} \left( \frac{1}{\sqrt{2}} H_{0}(f/2^{k}) \right)$$
   
   

3. $H_{0}(0) = \sqrt{2}$ and $H_{0}(1/2) = 0$. That is, $H_{0}(f)$ is a low pass filter since the frequency response is zero at the highest (normalized) frequency $f=1/2$.

4. $H_{1}(0) = 0$ and $H_{1}(1/2) = \sqrt{2}$. That is, $H_{1}(f)$ is a high pass filter since the frequency response is zero at the lowest frequency $f=0$.

5. $\vert H_{1}(f)\vert = \vert H_{0}(f+1/2)\vert$. That is, $\vert H_{1}(f)\vert$ is a mirror of $\vert H_{0}(f)\vert$ around the quadrature point $f=1/4$.

6. Frequency Domain Version of Double Shift Orthogonality

   
   

   $$\vert H_{0}(f)\vert^{2} + \vert H_{0}(f+1/2)\vert^{2} = 2$$
   
   

   To understand this equation we can define the filter

   
   

   $$p(n) = h_{0}(n) * h_{0}(-n) = \sum_{k} h_{0}(k) h_{0}(k-n)$$
   
   

   Notice that the double shift orthognality shows that $p(0) = 1$ and $p(2n) = 0$ for $n \ne 0$. That is, $p(n)$ is a centered, symmetric half-band filter.

   In the frequency domain we have $P(f) = \vert H_{0}(f)\vert^{2}$ so that we are trying to show that $P(f) + P(f+1/2) = 2$.

   The filter with frequency response $P(f+1/2)$ is the filter $p(n)$ with alternating signs: $(-1)^{n}p(n)$. We therefore have:

   
      Frequency Domain                         Time Domain
      
      P(f)     <----> [ . . . p(-3)   0  p(-1)  1  p(1)  0  p(3)  0  p(5) . . . ]
      P(f+1/2) <----> [ . . .-p(-3)   0 -p(-1)  1 -p(1)  0 -p(3)  0 -p(5) . . . ]
       sum     <----> [ . . .    0    0    0    2    0   0    0   0    0
   

   It follows that $P(f) + P(f+1/2) = 2$.

Example

To illustrate these equations we can use the m-file dshiftf.m as follows:

% File: dshiftf.m
% Illustrates the frequency domain properties of the filters
%echo on;
%diary dshiftf.diary;
[h0 h1] = wfilters('db3','r')

p0 = conv(h0,wrev(h0));

% alternate signs
ind = 1:2:length(p0);
p0shift = p0;
p0shift(ind) = -p0shift(ind);

[H0,W] = freqz(h0);
[H1,W] = freqz(h1);
[P0,W] = freqz(p0);
[P0shift,W] = freqz(p0shift);

subplot(2,2,1);
   plot(W/(2*pi), abs(H0));
   title('Low Pass Filter Response H_{0}(f)');
   xlabel('Frequency');
   axis([0 0.5 0 2.1]);

subplot(2,2,2);
   plot(W/(2*pi), abs(H1));
   title('High Pass Filter Response H_{1}(f)');
   xlabel('Frequency');
   axis([0 0.5 0 2.1]);

subplot(2,2,3);
   plot(W/(2*pi), abs(P0),W/(2*pi),abs(P0shift));
   title('P_{0}(f) and P_{0}(f + 1/2)');
   xlabel('Frequency');
   axis([0 0.5 0 2.1]);

subplot(2,2,4);
   plot(W/(2*pi), abs(P0)+abs(P0shift));
   title('P_{0}(f) + P_{0}(f + 1/2)');
   xlabel('Frequency');
   axis([0 0.5 0 2.1]);
%orient landscape;
%print -depsc 'dshiftf.eps';

6.3 Support of the Scale Function

The support of a function $\phi(t)$ is the range of values of $t$ where the function is non-zero.

The two-scale equation imposes a restriction on the support of the scale function.

Assume that the low pass filter $h(k)$ is causal with length $N+1$. Then we will show that the support of $\phi(t)$ is the interval $0 \le t \le N$.

The two-scale equation becomes:

$$\phi(t) = \sum_{k=0}^{N} h_{0}(k) \sqrt{2} \phi(2t-k)$$

If we assume that $\phi(t)$ has support on an interval $a \le t \le b$ then the first term in the sum, $\phi(2t)$, has support $a/2 \le t \le b/2$ and the last term in the sum, $\phi(2t - N)$, has support $(a+N)/2 \le t \le (b+N)/2$.

It follows that the largest possible support for the sum is $a/2 \le t \le (b+N)/2$. That is $a = a/2$ and $b = (b+N)/2$. Hence $a=0$ and $b = N$.

6.4 The Cascade Algorithm

The two-scale equation can be used to approximate the scale function using the low pass filter $h_{0}(k)$. We can think of the two-scale equation as an iteration. That is, starting with an approximation $\phi^{(0)}(t)$ we produce a sequence of approximations $\phi^{(n)}(t)$ by considering the two-scale equation as an update formula:

$$\phi^{(n+1)}(t) = \sum_{k} h_{0}(k) \sqrt{2} \phi^{(n)}(2t-k)$$

As a numerical algorithm we start with estimates $\phi^{(0)}(m)$ of $\phi(t)$ on the integer values $t=m$. Then, the first iteration produces an approximation to $\phi(t)$ at the half-integer points $t=m/2$ via:

$$\begin{eqnarray*} \phi^{(1)}(m/2) &=& \sum_{k} h_{0}(k) \sqrt{2} \phi^{(0)}(2(m/2)-k)\\ &=& \sum_{k} h_{0}(k) \sqrt{2} \phi^{(0)}(m-k) \end{eqnarray*}$$

The next iteration produces estimates of $\phi(t)$ at the quarter-integer points $t = m/4$ and so on.

If we carry out this algorithm for M steps then we get approximate values of $\phi(t)$ at points $t = \frac{k}{2^{M}}$ where $k$ ranges over $0 \le k \le N 2^{M} - 1$. That is, the times range from $t=0$ up to $t=N-\frac{1}{2^{M}}$.

Once we have an estimate for the scale function we can produce an estimate of the wavelet function from it and the high pass filter $h_{1}(k)$.

This algorithm can be implemented directly in Matlab but it is somewhat slow to carry out as Matlab is an interprted language.

This slowness is a general Matlab problem for any algorithm involving many nested for loops. If we can formulate our algorithm to use built in vector operations rather than for loops it will be much faster. There is a vector version of the algorithm which we can derive as follows.

We know that the subspaces $V_{0}$, $V_{1}$, $V_{2}$, ..., $V_{j}$ give a finer and finer decomposition of functions using the basis functions $\phi_{j,k}(t)$.

That is, any function $f(t)$ in $V_{1}$ can be written in terms of the basis in $V_{1}$ and also in terms of the basis in $V_{j}$:

$$\begin{eqnarray*} f(t) &=& \sum_{k} a_{1,k} \phi_{1,k}(t)\\ & & \\ f(t) &=& \sum_{k} a_{j,k} \phi_{j,k}(t) \end{eqnarray*}$$

We know that the coefficients $a_{j,k}$ are obtained by passing the coefficients $a_{1,k}$ through $(j-1)$ levels of the synthesis bank.

Now recall that the coordinates of the scale function, $\phi(t)$, in $V_{1}$ are exaclty the low pass filter coefficents, $h_{0}(n)$, and the coordinates of the wavelet, $w(t)$, in V1 are exactly the high pass filter coefficients, $h_{1}(n)$.

Hence we should upsample and low pass filter these two sets of coefficents $(j-1)$ times to get their coefficients in $V_{j}$.

At this point we must ``estimate'' $\phi_{j,k}(t)$. We know that

$$\phi_{j,0}(t) = 2^{j/2} \phi(2^{j} t)$$

So we can estimate this function by a box of height $2^{j/2}$ extending for $1/2^{j}$ units of time. We might as well center this box in the interval $0 \le t \le N)/2^{j}$, where $N+1$ = length of the filter $h_{0}(n)$.

Define $\mbox{offset} = (N-1)/2$. Then

$$\begin{array}{lclr} \phi_{j,k}(t) &=& 2^{j/2} & \mbox{ if } ... ...fset}+1)/2^{j}\\ &=& 0 & \mbox{ everywhere else} \end{array}$$

Substituting this for $\phi_{j,k}(t)$ we get that

$$f(t) = 2^{j/2}a_{j,k} \mbox{ for } (k+\mbox{offset})/2^{j} \le t < (k+\mbox{offset}+1)/2^{j}$$

The m-file ocascade.m carries out this algorithm.

function [phi,psi,x]=ocascade(h0,iter)
%[phi,psi,x]=ocascade(h0,iter)
%
%carrys out the cascade algorithm
%iter is the number of iterations used.

if nargin == 1 
   iter = 8;
end

N=length(h0)-1;

% compute the high pass filter (alternating flip)
h1 = fliplr(h0);
h1(2:2:(N+1)) = - h1(2:2:(N+1));


% do the first iteration
phi = h0;
psi = h1;

% now the rest of the iterations
for i=(2:iter)

   % upsample 
   top = 2*length(phi)-1;
   
   temp = phi;
   phi = zeros(1,top);
   phi(1:2:top) = temp;

   temp = psi;
   psi = zeros(1,top);
   psi(1:2:top) = temp;

   % filter
   phi = conv(h0,phi);
   psi = conv(h0,psi);
end


% scale the answers and pad with zeros on both ends to account for the offset
% note that the last x sample is exactly N - 1/2^iter. 

offset = (N-1)/2;
pad = zeros(1,offset);
%pad=[];
phi=2^(iter/2)*[pad phi pad];
psi=2^(iter/2)*[pad psi pad];

% compute the grid x
nf=length(phi); % same for all

%x = linspace(0, (nf-1)/2^iter,nf);
x = (0:(nf-1))/2^iter;

Example

The m-file ocascadedemo.m demonstrates the use the this m-file.