Digital Signal Processing Using MATLAB and Wavelets

Chapter 9.11 - Wavelet Transform Theory

The discrete wavelet transform convolves the input by the shifts (translation in time) and scales (dilations or contractions) of the wavelet. Below are variables commonly used in wavelet literature:

g represents the highpass (wavelet) filter

h represents the lowpass (scaling) filter

J is the total number of octaves

j is the current octave (used as an index. $1 <= j <= J$ )

N is the total number of inputs

n is the current input (used as an index. $1 <= n <= N$ )

L is the width of the filter (number of taps)

k is the current wavelet coefficient

$W_f(a, b)$ represents the continuous wavelet transform (CWT) of function $f$

$W_h[j, n]$ represents the discrete wavelet transform of function $f$

$W[j, n]$ represents the discrete scaling function of $f$, except:

$W[0, n]$ which is the input signal.

The continuous wavelet transform is represented by:


\begin{displaymath}W_f(a,b) = \int f(t) \psi(at+b) dt \end{displaymath}

where $f(t)$ is the function to analyze, $\psi$ is the wavelet, and $\psi( at + b)$ is the shifted and scaled version of the wavelet at time $b$ and scale $a$. An alternate form of the equation is:

\begin{displaymath}W_f(s,u) = \int_{-\infty}^{\infty} f(t) \sqrt{s} \psi(s(t-u)) dt \end{displaymath}

again where $\psi$ is the wavelet, while the wavelet family is shown above as $\sqrt{s}\psi(s(t-u))$, shifted by $u$ and scaled by $s$. We can rewrite the wavelet transform as an inner product:

\begin{displaymath}W_f(s,u) = \left< f(t), \sqrt{s} \psi(s(t-u)) \right>. \end{displaymath}

This inner product is essentially computed by the filters.

So far, this background focuses on how to get the wavelet transform given a wavelet, but how does one get the wavelet coefficients and implement the transform with filters? The relationship between wavelets and the filter banks that implement the wavelet transform is as follows. The scaling function, $\phi(t)$, is determined through recursively applying the filter coefficients, since multiresolution recursively convolutes the input vector after shifting and scaling. All the information about the scaling and wavelet functions is found by the coefficients of the scaling function and of the wavelet function, respectively. The scaling function is given as:


\begin{displaymath}\phi(t) = \sqrt{2} \sum_{k} h[k] \phi(2t-k) . \end{displaymath}

The wavelet function is given as:


\begin{displaymath}\psi(t) = \sqrt{2} \sum_{k} g[k] \phi(2t-k) . \end{displaymath}

There is a finite set of coefficients $h[k]$. Once these coefficients are found, allowing us to design the lowpass filter, then the highpass filter coefficients are easy to find. Daubechies chose the following coefficients, with many desirable (and required) properties, and solved the above equations with them.

Lowpass (scaling) coefficients:


\begin{displaymath}h[0], h[1], h[2], h[3] = \frac{1 + \sqrt{3}}{4 \sqrt{2}}, \fr...
...c{3 - \sqrt{3}}{4 \sqrt{2}}, \frac{1 - \sqrt{3}}{4 \sqrt{2}} . \end{displaymath}

These produce a space, $V_j$, that has invariance to shift and scale, which is needed for multiresolution analysis. Now that we have coefficients for the lowpass filter, the corresponding highpass (wavelet) filter coefficients can be calculated, as follows:


\begin{displaymath}g[0], g[1], g[2], g[3] = \frac{1 - \sqrt{3}}{4 \sqrt{2}}, \fr...
...{3 + \sqrt{3}}{4 \sqrt{2}}, \frac{-1 - \sqrt{3}}{4 \sqrt{2}} . \end{displaymath}

Notice that $g[0]=h[3]$, $g[1]=-h[2]$, $g[2]=h[1]$, and $g[3]=-h[0]$.

This pattern appears in the filter bank shown in Figure 9.4. The significance of the above coefficients is that they are used in a filter bank to generate the wavelet transform.

Writing the wavelet equation with regard to the discretely sampled sequence $x[n]$:


\begin{displaymath}W_f[a,b] = \frac{1}{\sqrt{a}} \sum_{n=b}^{aL+b-1} x[n] g\left(\frac{n - b}{a}\right) \end{displaymath}

where the wavelet is replaced by function $g$, which is obtained by sampling the continuous wavelet function. For the discrete case, we let $a=2k$ and require that the parameters ($a, b$) as well as $k$ be integers. To make this clearer, we replace $a$ and $b$ (which represent real numbers above) with $j$ and $n$ (which represent whole numbers). With the discrete case, redundancy is not required in the transformed signal. To reduce redundancy, the wavelet equation must satisfy certain conditions. First, we must introduce an orthogonal scaling function, which allows us to have an approximation at the last level of resolution (the initial level of resolution is simply the input signal). As we saw in previous sections, the coefficients for the scaling filter (lowpass) and the wavelet filter (highpass) are intimately related. The functions that produce these coefficients are also dependent on one another. Second, the representation of the signal at octave $j$ must have all the information of the signal at octave $j+1$, which makes sense: the input signal has more information than the first approximation of this signal. The function $x[n]$ is thus changed to $W[j-1,m]$, which is the scaling function's decomposition from the previous level of resolution, with $m$ as an index. Third, after $J$ levels of resolution, the result of the scaling function on the signal will be 0. After repeatedly viewing a signal at coarser and coarser approximations, eventually the scaling function will not produce any useful information. Fourth, the scaling function allows us to approximate any given signal with a variable amount of precision. The scaling function, $h$, gives us an approximation of the signal via the following equation. This is also known as the lowpass output:


\begin{displaymath}W[j,n] = \sum_{m=0}^{N-1} W[j-1,m] h[2n-m] . \end{displaymath}

The wavelet function gives us the detail signal, also called highpass output:


\begin{displaymath}W_{h}[j,n] = \sum_{m=0}^{N-1} W[j-1,m] g[2n-m] . \end{displaymath}

The $n$ term gives us the shift, the starting points for the wavelet calculations. The index $2n - m$ incorporates the scaling, resulting in half the outputs for octave $j$ compared to the previous octave $j-1$.



Example:
Suppose $g = \{ \frac{1}{\sqrt{2}}, \frac{-1}{\sqrt{2}} \}$, and $x = \{34, 28, 76, 51, 19, 12, 25, 43\}$. Find the detail signal produced by the wavelet transform on this signal.

Answer:

\begin{displaymath}W[0, n] = x[n] \end{displaymath}

\begin{displaymath}W_h[1, n] = W[0, 0] g[2n-0] + W[0, 1] g[2n-1] + W[0, 2] g[2n-2] + W[0, 3] g[2n-3] \end{displaymath}

\begin{displaymath}+ W[0, 4] g[2n-4] + W[0, 5] g[2n-5] + W[0, 6] g[2n-6] + W[0, 7] g[2n-7] \end{displaymath}


For this example, the value of $g$ with any index value outside of [0, 1] is zero. Also, the value of $W[j, n] = 0$ for any value of $n$ less than 0 or greater or equal to $N$. So we can find the individual values for $W_h$ as follows.

Looking at octave 1:

\begin{displaymath}W_h[1, n] = 34 g[2n-0] + 28 g[2n-1] + 76 g[2n-2] + 51 g[2n-3] + 19 g[2n-4] \end{displaymath}

\begin{displaymath}+ 12 g[2n-5] + 25 g[2n-6] + 43 g[2n-7] \end{displaymath}


\begin{displaymath}W_h[1, 0] = 34 g[0] + 28 g[-1] + 76 g[-2] + 51 g[-3] + 19 g[-4] + 12 g[-5] + 25 g[-6] + 43 g[-7] \end{displaymath}

\begin{displaymath}W_h[1, 0] = 34 \frac{1}{\sqrt{2}} \end{displaymath}


\begin{displaymath}W_h[1, 1] = 34 g[2] + 28 g[1] + 76 g[0] + 51 g[-1] + 19 g[-2] + 12 g[-3] + 25 g[-4] + 43 g[-5] \end{displaymath}

\begin{displaymath}W_h[1, 1] = 28 \frac{-1}{\sqrt{2}} + 76 \frac{1}{\sqrt{2}} \end{displaymath}


\begin{displaymath}W_h[1, 2] = 34 g[4] + 28 g[3] + 76 g[2] + 51 g[1] + 19 g[0] + 12 g[-1] + 25 g[-2] + 43 g[-3] \end{displaymath}

\begin{displaymath}W_h[1, 2] = 51 \frac{-1}{\sqrt{2}} + 19 \frac{1}{\sqrt{2}} \end{displaymath}


\begin{displaymath}W_h[1, 3] = 34 g[6] + 28 g[5] + 76 g[4] + 51 g[3] + 19 g[2] + 12 g[1] + 25 g[0] + 43 g[-1] \end{displaymath}

\begin{displaymath}W_h[1, 3] = 12 \frac{-1}{\sqrt{2}} + 25 \frac{1}{\sqrt{2}} \end{displaymath}


\begin{displaymath}W_h[1, 4] = 34 g[8] + 28 g[7] + 76 g[6] + 51 g[5] + 19 g[4] + 12 g[3] + 25 g[2] + 43 g[1] \end{displaymath}

\begin{displaymath}W_h[1, 4] = 43 \frac{-1}{\sqrt{2}} \end{displaymath}


\begin{displaymath}W_h[1, 5] = 34 g[10] + 28 g[9] + 76 g[8] + 51 g[7] + 19 g[6] + 12 g[5] + 25 g[4] + 43 g[3] \end{displaymath}

\begin{displaymath}W_h[1, 5] = 0 . \end{displaymath}

It should be obvious that any value $W_h[1, n] = 0$ when $n \ge 5$. Comparing this to the filter bank approach, we look at Figure 9.25, which shows a channel from a filter bank. Keeping in mind that $g[0] = a$, and $g[1] = b$, we can trace the signal through as follows.


\begin{displaymath}y = \frac{34 - 0}{\sqrt{2}}, \frac{28 - 34}{\sqrt{2}}, \frac{...
...}{\sqrt{2}}, \frac{43 - 25}{\sqrt{2}}, \frac{0 - 43}{\sqrt{2}} \end{displaymath}

The output is a down-sampled version of $y$, so when we eliminate every other value, we get:

\begin{displaymath}output = \frac{34}{\sqrt{2}}, \frac{76 - 28}{\sqrt{2}}, \frac...
...1}{\sqrt{2}}, \frac{25 - 12}{\sqrt{2}}, \frac{-43}{\sqrt{2}} . \end{displaymath}

We end up with the same result, as expected. Of course, the scaling equation is the same basic equation, with different coefficients.

Figure 9.25: One channel of the analysis side of a filter bank.

Before we can look at octave 2, we need to know the $W[1, n]$ values, which we present as the following. Realize that these can be found in the same fashion as above, except that both filter coefficients are positive.


\begin{displaymath}W[1, 0] = \frac{34}{\sqrt{2}} \end{displaymath}

\begin{displaymath}W[1, 1] = \frac{28 + 76}{\sqrt{2}} = \frac{104}{\sqrt{2}}\end{displaymath}

\begin{displaymath}W[1, 2] = \frac{51 + 19}{\sqrt{2}} = \frac{70}{\sqrt{2}}\end{displaymath}

\begin{displaymath}W[1, 3] = \frac{12 + 25}{\sqrt{2}} = \frac{37}{\sqrt{2}}\end{displaymath}

\begin{displaymath}W[1, 4] = \frac{43}{\sqrt{2}} \end{displaymath}

Looking at octave 2:

\begin{displaymath}W_h[2, n] = W[1, 0] g[2n-0] + W[1, 1] g[2n-1] + W[1, 2] g[2n-2] + W[1, 3] g[2n-3] \end{displaymath}

\begin{displaymath}+ W[1, 4] g[2n-4] + W[1, 5] g[2n-5] + W[1, 6] g[2n-6] + W[1, 7] g[2n-7] \end{displaymath}

\begin{displaymath}W_h[2, n] = \frac{34}{\sqrt{2}} g[2n-0] + \frac{104}{\sqrt{2}...
...1] + \frac{70}{\sqrt{2}} g[2n-2] + \frac{37}{\sqrt{2}} g[2n-3] \end{displaymath}

\begin{displaymath}+ \frac{43}{\sqrt{2}} g[2n-4] + 0 g[2n-5] + 0 g[2n-6] + 0 g[2n-7] \end{displaymath}


\begin{displaymath}W_h[2, 0] = \frac{34}{\sqrt{2}} g[0] + \frac{104}{\sqrt{2}} g[-1] + \frac{70}{\sqrt{2}} g[-2] + \frac{37}{\sqrt{2}} g[-3] \end{displaymath}

\begin{displaymath}+ \frac{43}{\sqrt{2}} g[-4] + 0 g[-5] + 0 g[-6] + 0 g[-7] \end{displaymath}

\begin{displaymath}W_h[2, 0] = \left( \frac{34}{\sqrt{2}} \right) \left( \frac{1}{\sqrt{2}}\right) \end{displaymath}

\begin{displaymath}W_h[2, 0] = \frac{34}{2} = 17 \end{displaymath}


\begin{displaymath}W_h[2, 1] = \frac{34}{\sqrt{2}} g[2] + \frac{104}{\sqrt{2}} g[1] + \frac{70}{\sqrt{2}} g[0] + \frac{37}{\sqrt{2}} g[-1] \end{displaymath}

\begin{displaymath}+ \frac{43}{\sqrt{2}} g[-2] + 0 g[-3] + 0 g[-4] + 0 g[-5] \end{displaymath}

\begin{displaymath}W_h[2, 1] = \left( \frac{104}{\sqrt{2}} \right) \left( \frac{...
... \frac{70}{\sqrt{2}} \right) \left( \frac{1}{\sqrt{2}} \right) \end{displaymath}

\begin{displaymath}W_h[2, 1] = \frac{70 - 104}{2} = \frac{-34}{2} = -17 \end{displaymath}


\begin{displaymath}W_h[2, 2] = \frac{34}{\sqrt{2}} g[4] + \frac{104}{\sqrt{2}} g[3] + \frac{70}{\sqrt{2}} g[2] + \frac{37}{\sqrt{2}} g[1] \end{displaymath}

\begin{displaymath}+ \frac{43}{\sqrt{2}} g[0] + 0 g[-1] + 0 g[-2] + 0 g[-3] \end{displaymath}

\begin{displaymath}W_h[2, 2] = \left( \frac{37}{\sqrt{2}} \right) \left( \frac{-...
... \frac{43}{\sqrt{2}} \right) \left( \frac{1}{\sqrt{2}} \right) \end{displaymath}

\begin{displaymath}W_h[2, 2] = \frac{43 - 37}{2} = \frac{6}{2} = 3 \end{displaymath}


\begin{displaymath}W_h[2, 3] = \frac{34}{\sqrt{2}} g[6] + \frac{104}{\sqrt{2}} g[5] + \frac{70}{\sqrt{2}} g[4] + \frac{37}{\sqrt{2}} g[3] \end{displaymath}

\begin{displaymath}+ \frac{43}{\sqrt{2}} g[2] + 0 g[1] + 0 g[0] + 0 g[-1] \end{displaymath}

\begin{displaymath}W_h[2, 3] = 0 . \end{displaymath}

Again, we can verify this with the filter bank approach. Notice how the $g[2n-m]$ term takes care of down-sampling for us. We can continue on, with the third octave. First, we need the lowpass outputs of octave 2, found in the same way as above, only with different coefficients.


\begin{displaymath}W[2, 0] = \frac{34}{2} = 17 \end{displaymath}

\begin{displaymath}W[2, 1] = \frac{174}{2} = 87 \end{displaymath}

\begin{displaymath}W[2, 2] = \frac{80}{2} = 40 \end{displaymath}

\begin{displaymath}W[2, 3] = 0 \end{displaymath}

Looking at octave 3:

\begin{displaymath}W_h[3, n] = W[2, 0] g[2n-0] + W[2, 1] g[2n-1] + W[2, 2] g[2n-2] + W[2, 3] g[2n-3] \end{displaymath}

\begin{displaymath}+ W[2, 4] g[2n-4] + W[2, 5] g[2n-5] + W[2, 6] g[2n-6] + W[2, 7] g[2n-7] \end{displaymath}


\begin{displaymath}W_h[3, n] = 17 g[2n-0] + 87 g[2n-1] + 40 g[2n-2] \end{displaymath}

\begin{displaymath}W_h[3, 0] = 17 g[0] + 87 g[-1] + 40 g[-2] = \frac{17}{\sqrt{2}} \end{displaymath}

\begin{displaymath}W_h[3, 1] = 17 g[2] + 87 g[1] + 40 g[0] = \frac{40 - 87}{\sqrt{2}} = \frac{-47}{\sqrt{2}} \end{displaymath}

\begin{displaymath}W_h[3, 2] = 17 g[4] + 87 g[3] + 40 g[2] = 0 . \end{displaymath}

Finding the lowpass outputs for this octave:


\begin{displaymath}W[3, 0] = 17 h[0] + 87 h[-1] + 40 h[-2] = \frac{17}{\sqrt{2}} \end{displaymath}

\begin{displaymath}W[3, 1] = 17 h[2] + 87 h[1] + 40 h[0] = \frac{40 - 87}{\sqrt{2}} = \frac{-47}{\sqrt{2}} \end{displaymath}

\begin{displaymath}W[3, 2] = 17 h[4] + 87 h[3] + 40 h[2] = 0 . \end{displaymath}

We stop here, because we get no more useful information after this point. Why? Because our original signal has only 8 samples, and each octave uses only half the data of the previous octave. In other words, since we have $2^3$ values, we should not perform more than 3 octaves of resolution on the signal.



We end this chapter with a program to show the frequency magnitude response of a given wavelet.




Using the default db2 wavelet, we get the graphs shown in Figure 9.26 and Figure 9.27. In order to see the frequency magnitude responses of other wavelets, the wavelets toolbox must be installed. These figures show what we should expect: the lowpass filter (and inverse lowpass filter) keep low frequencies about the same, but attenuate high frequencies.

Figure 9.26: The first graph of the "show_wavelet' program.

Figure 9.27: The second graph of the "show_wavelet' program.
.

UNLIMITED FREE
ACCESS
TO THE WORLD'S BEST IDEAS

SUBMIT
Already a GlobalSpec user? Log in.

This is embarrasing...

An error occurred while processing the form. Please try again in a few minutes.

Customize Your GlobalSpec Experience

Category: DSP Boards
Finish!
Privacy Policy

This is embarrasing...

An error occurred while processing the form. Please try again in a few minutes.