Channels ▼
RSS

Parallel

The Fast Wavelet Transform

Source Code Accompanies This Article. Download It Now.


APR92: THE FAST WAVELET TRANSFORM

THE FAST WAVELET TRANSFORM

Beyond Fourier transforms

This article contains the following executables: WAVELET.ARC

Mac A. Cody

Mac is a digital systems manager for Sunair Electronics Inc. in Ft. Lauderdale, Fla. He designs modems for shortwave radios using DSP and embedded C programming. He can be contacted on GENIE MAIL as M.CODY1.


The Fourier transform has been used as a reliable tool in signal analysis for many years. Invented in the early 1800s by its namesake, Jean-Baptiste-Joseph Fourier, the Fourier transform has become a cornerstone of modern signal analysis.

The Fourier transform has proven incredibly versatile in applications ranging from pattern recognition to image processing. Nevertheless, it suffers from certain limitations. Recently, a new kind of transform, the wavelet transform, has been shown to be as powerful and versatile as the Fourier transform, yet without some of its limitations.

The wavelet transform is the result of work by a number of researchers. Initially, a French geophysicist, Jean Morlet, came up with an ad hoc method to model the process of sound waves traveling through the earth's crust. Unlike Fourier analysis, he did not use sine and cosine curves, but simpler ones which he called "wavelets." Yves Meyer, a mathematician, recognized this work to be part of the field of harmonic analysis, and came up with a family of wavelets that he proved were most efficient for modeling complex phenomena. This work was improved upon by two American researchers, Stephane Mallat of New York University and Ingrid Daubechies of Bell Labs. Since 1988, there has been a small explosion of activity in this area, as engineers and researchers apply the wavelet transform to applications ranging from image compression to fingerprint analysis. The wavelet transform has even been implemented in silicon, in the form of a chip from Aware Inc.

But before discussing how the wavelet transform works, let's first examine its predecessor.

How the Fourier Transform Works

A function in the time domain is translated by the Fourier transform into a function in the frequency domain, where it can be analyzed for its frequency content. This translation occurs because the Fourier transform expands the original function in terms of orthonormal basis functions of sine and cosine waves of infinite duration. The Fourier coefficients of the transformed function represent the contribution of each sine and cosine wave at each frequency.

The Fourier transform works under the assumption that the original time-domain function is periodic in nature. As a result, the Fourier transform has difficulty with functions having transient components, that is, components localized in time. This is especially apparent when a signal has sharp transitions. Another problem is that the Fourier transform of a signal does not convey any information pertaining to translation of the signal in time. Applications that use the Fourier transform often work around the first problem by windowing the input data so that the sampled values converge to 0 at the endpoints. Attempts to solve the second problem, such as the development of the short-time Fourier transform, have met with marginal success.

In recent years, new families of orthonormal basis functions have been discovered that lead to transforms which overcome the problems of the Fourier transform. These basis functions are called "wavelets," and unlike the sine and cosine wave of the Fourier transform, they need not have infinite duration. They can be nonzero for only a small range of the wavelet function. This "compact support" allows the wavelet transform to translate a time-domain function into a representation that is localized not only in frequency (like the Fourier transform) but in time as well. This ability has brought forth new developments in the fields of signal analysis, image processing, and data compression.

This article provides a basic understanding of wavelets and their use in representing signals through the wavelet transform. A description and implementation of a computationally efficient fast wavelet transform is also presented. It does for the wavelet transform what the Fast Fourier Transform (FFT) does for the Discrete Fourier Transform (DFT).

Meeting the Wave Head On

Before implementing a fast wavelet transform, a brief explanation of wavelets and the wavelet transform is in order. There are two fundamental equations upon which wavelet calculations are based: the scaling function (also called the "basic dilation equation" or "fundamental recursion"), see Example 1(a); and the basic (or "primary wavelet") function, see Example 1(b), where Z is the set of integers and the ak are the wavelet coefficients. Both functions are two-scale difference equations. They are the prototypes of a class of orthonormal basis functions of the form shown in Example 1(c), where the parameter j controls the dilation or compression of the function in time scale and amplitude. The parameter k controls the translation of the function in time. The set of basis functions formed by (t) and (t) is a system of scaled and translated wavelets.

Now, not just any set of ak can be used as wavelet coefficients. Certain wavelet conditions must be satisfied for a set of ak to qualify as wavelet coefficients. These conditions are shown in Example 1(d) where ak is the conjugate of ak. (ak can be complex.)

Wavelet systems can, as a result, be either real or complex valued. Most references have used real-valued wavelet systems. Wavelet systems may or may not have compact support. (As mentioned earlier, wavelets have compact support if, and only if, they have a finite number of nonzero coefficients.) It is compact support that allows wavelets to localize in both time and frequency, so we will deal with wavelets of this type only.

Several techniques have been used to create wavelet systems. These include cubic splines, complex exponentials, and parameter-space constructions. The wavelet coefficient generation routine in Listing Two uses formulas to generate a parameter space for the general, real-valued wavelet system with two to six coefficients. These parameter-space formulas were developed by David Pollen of Aware Inc. and are listed in Table 1. The parameter space modeled by the formulas includes the wavelet systems introduced by Ingrid Daubechies of Bell Laboratories and the original wavelet basis, the Haar basis.

Once a wavelet system is created, it can be used to expand a function g(t) in terms of the basis functions shown in Example 2(a), with the coefficients calculated by inner products, as seen in Example 2(b). If the wavelet system has compact support and an upper limit is placed upon the degree of dilation j, then the expansion equation becomes that shown in Example 2(c).

The expansion coefficients c(l) represent the approximation of the original signal g(t) with a resolution of one point per every 2J points of the original signal. The expansion coefficients d(j,k) represent details of the original signal at different levels of resolution. These coefficients completely and uniquely describe the original signal and can be used in a way similar to the Fourier transform. The wavelet transform, then, is the process of determining the values of c(l) and d(j,k) for a given g(t) and wavelet system.

Digital Filers, Trees, and Fast Transforms

The Fourier transform was made practical to implement on computers by the development of the Fast Fourier Transform by J.W. Cooley and J.W. Tukey in the mid-sixties. (Actually, it appears that Cooley and Tukey rediscovered the FFT. An unpublished treatise, written by Carl Friedrich Gauss in the early 1800s, describes an algorithm similar to the one developed by Cooley and Tukey. Gauss's algorithm seems to predate Fourier's 1807 work on harmonic analysis! You can draw your own conclusions.) The FFT eliminated redundancies that exist in the Discrete Fourier Transform. While a DFT of length N requires N2 multiplications and additions, the FFT requires only N log2 N such operations.

The expansion equation naturally leads to a recursive algorithm for the wavelet transform, given certain assumptions. First, the function g(t) is taken as a sequence of discrete points, YJ,p, sampled at 2J points per unit interval. These points can be viewed as the inner products of J,p and g(t). That is, the sample points are an approximation, or c(l) coefficients, of the continuous function g(t). This allows c(l) and d(j,k) terms to be calculated by direct convolution of g(t) samples with the coefficients ak.

Daubechies has discovered that the wavelet transform can be implemented with a specially designed pair of Finite Impulse Response (FIR) filters called a "Quadrature Mirror Filter" (QMF) pair. FIR filters and QMFs come from the field of digital signal processing. A FIR filter performs the dot product (or sum of products) between the filter coefficients and the discrete samples in the tapped delay line of the filter. The act of passing a set of discrete samples, representing a signal, through a FIR filter is a discrete convolution of the signal with the filter's coefficients; see Figure 1.

QMF filters are used in digital speech analysis and spectral decomposition. Briefly, QMFs are distinctive because the frequency responses of the two FIR filters separate the high-frequency and low-frequency components of the input signal. The dividing point is usually halfway between 0 Hertz (DC) and half the data Sampling rate (the Nyquist frequency).

The outputs of the QMF filter pair are decimated (or de-sampled) by a factor of two; that is, every other output sample of the filter is kept, and the others are discarded. The low-frequency (low-pass) filter outputs is fed into another identical QMF filter pair. This operation can be repeated recursively as a tree or pyramid algorithm (Figure 2), yielding a group of signals that divides the spectrum of the original signal into octave bands with successively coarser measurements in time as the width of each spectral band narrows and decreases in frequency; see Figure 3.

Mallat has shown that the tree or pyramid algorithm can be applied to the wavelet transform by using the wavelet coefficients as the filter coefficients of the QMF filter pairs. The same wavelet coefficients are used in both low-pass and high-pass (actually, band-pass) filters. The low-pass filter coefficients are associated with the ak of the scaling function . The output of each low-pass filter is the c(l), or approximation components, of the original signal for that level of the tree. The high-pass filter is associated with the ak of the wavelet function . (Note the alternating sign change.) The output of each high-pass filter is the d(j,k), or detail components, of the original signal at resolution 2j. The c(l) of the previous level are used to generate the new c(l) and d(j,k) for the next level of the tree. Decimation by two corresponds to the multiresolutional nature (the j parameter) of the scaling and wavelet functions; see Figure 4 .

The reverse fast wavelet transform essentially performs the operations associated with the forward fast wavelet transform in the opposite direction. The expansion coefficients are combined to reconstruct the original signal. The same ak coefficients are used as in the forward transform, but in reverse order. The process works down the branches of the tree combining the approximation and detail signals into approximation signals with higher levels of detail; see Figure 5. Instead of decimation, the signals are interpolated: 0s are placed between each approximation and detail sample, and the signals are then passed through the low-pass and high-pass filters. The intermediate 0 values are replaced by "estimates" derived from the convolutions. The filters' outputs are then summed to form the approximation coefficients for the next-higher level of resolution; see Figure 6. The final set of approximation coefficients at the tree's top level in the reverse transform is a reconstruction of the original signal data points.

The fast wavelet transform is actually more computationally efficient than the Fast Fourier Transform. As mentioned previously, an FFT of length N (where N is an integral power of 2) takes on the order of N log2 N operations. A fast wavelet transform of length N requires approximately N operations--the best possible. The FIR filters are linear processing elements. In the tree algorithm, each consecutive level performs only half as many operations as its predecessor.

Catch a Wave...

Wavelet filter coefficients are stored in arrays of type double, six elements in length. Unused array elements are set to 0. Storage for the sampled data signal is also an array of type double; see Figure 7(a). The only constraint on the length of the array is that the number of elements must be an integer multiple of the number of sample points per unit interval, 2J. (J, once again, is the number of levels in the wavelet transform to be calculated.) An additional five storage locations are allocated to allow for convolution overrun. Transform coefficient storage consists of a dynamically created array of arrays of type double; see Figure 7(b). The pairs of arrays store the approximation and detail coefficients for each given level. The length of each pair of arrays is set according to the length of the array containing the sampled data signal and the number of levels in the transform. If the length of the input data array is N+5, then the topmost pair of arrays has a length of N/2+5 data elements, the next pair has a length of N/4+5 data elements, and so on.

Listings One and Two, page 100, provide a set of C functions for implementing the forward and inverse fast wavelet transforms, the generation of wavelet filter coefficients, and data-structure management. Headers for the functions are provided in Listing Three; see page 101. The functions are written in Turbo C 2.0 and conform to standard C conventions. The code should be portable to most compilers and host systems.

The functions in Listing One provide utilities for managing the storage structure for wavelet transform coefficients. BuildTreeStorage uses calloc to create a dynamic storage structure (array of arrays) for wavelet coefficients and intermediate results. DestroyTreeStorage uses free to delete an existing wavelet-coefficient storage structure. TreeCopy copies wavelet coefficients from one storage structure to another. Only the c(l) (approximation) coefficients on the last level are copied along with all the d(j,k) (detail) coefficients. TreeZero sets all elements of the transform coefficient storage structure to 0. ZeroTreeDetail sets a given number of detail-coefficient levels to 0. The routine starts with the top level of detail coefficients and works down the required number of levels.

WaveletCoeffs, see Listing Two, implements the formulas to generate a parameter space for the general, real-valued wavelet system with two to six coefficients. The formulas are listed in Table 1. Varying the values of alpha and beta yields an infinite variety of wavelet coefficient sets (ak). MakeWaveletFilters takes the coefficients created by WaveletCoeffs and creates the low-pass and high-pass filters for either the forward or the inverse fast wavelet transform.

To make the implementation of the forward and inverse fast wavelet transforms easier to comprehend, the algorithms have been coded into several simple routines. These routines represent different levels of the algorithms functionality.

DotP performs the dot product for the forward fast wavelet transform. Each coefficient of the wavelet filter is used in calculating a dot product. ConvolveDec2 implements convolution and decimation for the transform. The dot products of the wavelet filter coefficients and the sample data in the input array are calculated via calls to DotP. As shown in Figure 4,there is a shift of two samples between each dot product, thus accomplishing the decimation. The dot products become the transform coefficients of the branch of the next level. Special test conditions are taken for data samples at the beginning and end of the input array.

DecomposeBranches decomposes a level of the forward fast wavelet transform. The tree's current approximation level is transformed into two branches containing the detail and approximation decimations for the next level. This is accomplished with two calls to ConvolveDec2. One call passes the low-pass wavelet filter coefficients to generate the approximation coefficients for the next-lower transform level. The second call passes the high-pass wavelet filter coefficients to generate the detail coefficients for the next-lower transform level.

WaveletDecomposition is the driver routine for the forward fast wavelet transform algorithm. The code performs the necessary number of decompositions for the desired number of levels via calls to DecomposeBranches. For the first pass of the loop, InData is the sampled input signal. For subsequent loops, the approximation coefficients of the previously calculated level are used.

DotpEven and DotpOdd use the even- and odd-numbered filter coefficients to calculate dot products for the inverse fast wavelet transform. ConvolveInt2 implements the convolution, interpolation, and summation operations for the inverse fast wavelet transform. The dot products of the wavelet filter coefficients and groups of samples in the input array are calculated via calls to DotpEven and DotpOdd. Alternating application of DotpEven and DotpOdd on the same subset of data points emulates interpolation of the data with 0s inserted between each data point, as shown in Figure 7. Special test conditions are taken for data samples at the beginning and end of the array.

ReconstructBranches performs reconstruction of an inverse wavelet transform level. ConvolveInt2 is called twice to convolve and sum the approximation and detail branches to form the next, more-detailed level of approximation coefficients. The first call of ConvolveInt2 convolves the low-pass wavelet filter with the level's approximation branch. The second call convolves the high-pass wavelet filter with the level's detail branch and sums the results of both convolutions. WaveletReconstruction is the driver routine for inverse fast wavelet transform. It performs the necessary number of reconstructions for the number of wavelet coefficient levels using ReconstructBranches.

Making Waves...

To experiment with the wavelet transform, use the WAVEDEMO program, which employs the functions in Listings One and Two. WAVEDEMO is written in Turbo C 2.0 and uses Borland's BGI graphics to provide support for CGA, EGA, VGA, and Hercules graphics adapters. The display shows the original signal and the expansion coefficients at decreasing resolutions. An alternate display shows the expansion coefficients and the reconstructed signal. The program can perform up to six levels of coefficient expansion. The screen acts as a window that can be scrolled over the displayed coefficients. Another display shows the signal data and expansion coefficients in numerical format. A menu interface allows the user to generate wavelet coefficients; decompose and reconstruct signals; view and modify the expansion coefficients; and store and retrieve signal data and expansion coefficients from disk files. Because of its size, the program is available electronically; see page 3. The electronic version includes the WAVEDEMO executable and source, sample input data and wavelet coefficient files, and documentation.

Wave Goodbye!

The wavelet transform allows for simultaneous analysis of time and frequency in signals. An infinite variety of wavelets enables the user to tailor that analysis to extract this information of interest. An efficient implementation for computers has been made possible through the fast wavelet transform algorithm. The software in this article now makes this tool available to the masses.

Special thanks to Mr. Karl Jagler of Aware Inc. (Cambridge, Mass.), for giving me the insights and understanding into implementing the fast wavelet transform. His time and patience were at least as important in helping me understand wavelets as the reference articles and algorithm examples he sent me.

Bibliography

Elliot, Douglas F. ed. Handbook of Digital Signal Processing: Engineering Applications. San Diego, Calif.: Academic Press, 1987.

Mallat, Stephane G. "A Theory for Multiresolution Signal Decomposition: The Wavelet Representation." IEEE Transactions on Pattern Analysis and Machine Intelligence (July, 1989).

Resnikoff, H.L. "Foundations of Arithmeticum Analysis: Compactly Supported Wavelets and the Wavelet Group." Aware Report No. AD890507.1.

Resnikoff, H.L. and C.S. Burrus. "Relationships Between the Fourier Transform and the Wavelet Transform." Aware Report No. AD900609.

Strang, Gilbert. "Wavelets and Dilation Equations: A Brief Introduction." SIAM Review (December, 1989).

The Ultra Wave Explorer User's Manual. Cambridge, Mass.: Aware Inc., 1989.


_THE FAST WAVELET TRANSFORM_
by Mac A. Cody


[LISTING ONE]
<a name="00af_000c">

#define WAVE_MGT
#include <alloc.h>
#include "wave_mgt.h"
double **BuildTreeStorage(int inlength, int levels)
{
  double **tree;
  int i, j;
  /* create decomposition tree */
  tree = (double **) calloc(2 * levels, sizeof(double *));
  j = inlength;
  for (i = 0; i < levels; i++)
  {
    j /= 2;
    if (j == 0)
    {
      levels = i;
/*    printf("\nToo many levels requested for available data\n");
      printf("Number of transform levels now set to %d\n", levels); */
      continue;
    }
    tree[2 * i] = (double *) calloc((j + 5), sizeof(double));
    tree[2 * i + 1] = (double *) calloc((j + 5), sizeof(double));
  }
  return tree;
}
void DestroyTreeStorage(double **tree, int levels)
{
  char i;
  for (i = (2 * levels - 1); i >= 0; i--)
    free(tree[i]);
  free(tree);
}
void TreeCopy(double **TreeDest, double **TreeSrc, int siglen, int levels)
{
  int i, j;
  for (i = 0; i < levels; i++)
  {
    siglen /= 2;
    for (j = 0; j < siglen + 5; j++)
    {
      if ((i + 1) == levels)
        TreeDest[2 * i][j] = TreeSrc[2 * i][j];
      else
        TreeDest[2 * i][j] = 0.0;
      TreeDest[(2 * i) + 1][j] = TreeSrc[(2 * i) + 1][j];
    }
  }
}
void TreeZero(double **Tree, int siglen, int levels)
{
  int i, j;
  for (i = 0; i < levels; i++)
  {
    siglen /= 2;
    for (j = 0; j < siglen + 5; j++)
    {
      Tree[2 * i][j] = 0.0;
      Tree[(2 * i) + 1][j] = 0.0;
    }
  }
}
void ZeroTreeDetail(double **Tree, int siglen, int levels)
{
  int i, j;
  for (i = 0; i < levels; i++)
  {
    siglen /= 2;
    for (j = 0; j < siglen + 5; j++)
      Tree[(2 * i) + 1][j] = 0.0;
  }
}





<a name="00af_000d">
<a name="00af_000e">
[LISTING TWO]
<a name="00af_000e">

/* WAVELET.C */
#include <math.h>
typedef enum { DECOMP, RECON } wavetype;
#include "wavelet.h"
void WaveletCoeffs(double alpha, double beta, double *wavecoeffs)
{
  double tcosa, tcosb, tsina, tsinb;
  char i;
  /* precalculate cosine of alpha and sine of beta to reduce */
  /* processing time */
  tcosa = cos(alpha);
  tcosb = cos(beta);
  tsina = sin(alpha);
  tsinb = sin(beta);
  /* calculate first two wavelet coefficients,  a = a(-2) and b = a(-1) */
  wavecoeffs[0] = ((1.0 + tcosa + tsina) * (1.0 - tcosb - tsinb)
                                      + 2.0 * tsinb * tcosa) / 4.0;
  wavecoeffs[1] = ((1.0 - tcosa + tsina) * (1.0 + tcosb - tsinb)
                                      - 2.0 * tsinb * tcosa) / 4.0;
  /* precalculate cosine and sine of alpha minus beta to reduce */
  /* processing time */
  tcosa = cos(alpha - beta);
  tsina = sin(alpha - beta);
  /* calculate last four wavelet coefficients  c = a(0), d = a(1), */
  /* e = a(2), and f = a(3) */
  wavecoeffs[2]  = (1.0 + tcosa + tsina) / 2.0;
  wavecoeffs[3]  = (1.0 + tcosa - tsina) / 2.0;
  wavecoeffs[4]  = 1 - wavecoeffs[0] - wavecoeffs[2];
  wavecoeffs[5]  = 1 - wavecoeffs[1] - wavecoeffs[3];
  /* zero out very small coefficient values caused by truncation error */
  for (i = 0; i < 6; i++)
  {
    if (fabs(wavecoeffs[i]) < 1.0e-15)
      wavecoeffs[i] = 0.0;
  }
}
char MakeWaveletFilters(double *wavecoeffs, double *Lfilter,
                               double *Hfilter, wavetype transform)
{
  char i, j, k, filterlength;
  /* find the first non-zero wavelet coefficient */
  i = 0;
  while(wavecoeffs[i] == 0.0)
    i++;
  /* find the last non-zero wavelet coefficient */
  j = 5;
  while(wavecoeffs[j] == 0.0)
    j--;
 /* form decomposition filters h~ and g~ or reconstruction filters h and g.
 Division by 2 in construction of decomposition filters is for normalization */
  filterlength = j - i + 1;
  for(k = 0; k < filterlength; k++)
  {
    if (transform == DECOMP)
    {
      Lfilter[k] = wavecoeffs[j--] / 2.0;
      Hfilter[k] = (double) (((i++ & 0x01) * 2) - 1) * wavecoeffs[i] / 2.0;
    }
    else
    {
      Lfilter[k] = wavecoeffs[i++];
      Hfilter[k] = (double) (((j-- & 0x01) * 2) - 1) * wavecoeffs[j];
    }
  }
  /* clear out the additional array locations, if any */
  while (k < 6)
  {
    Lfilter[k] = 0.0;
    Hfilter[k++] = 0.0;
  }
  return filterlength;
}
double DotP(double *data, double *filter, char filtlen)
{
  char i;
  double sum;
  sum = 0.0;
  for (i = 0; i < filtlen; i++)
    sum += *data-- * *filter++; /* decreasing data pointer is */
                                /* moving backward in time */
  return sum;
}
void ConvolveDec2(double *input_sequence, int inp_length,
          double *filter, char filtlen, double *output_sequence)
/* convolve the input sequence with the filter and decimate by two */
{
 int i;
 char shortlen, offset;
 for(i = 0; (i <= inp_length + 8) && ((i - filtlen) <= inp_length + 8); i += 2)
  {
    if (i < filtlen)
      *output_sequence++ = DotP(input_sequence + i, filter, i + 1);
    else if (i > (inp_length + 5))
    {
      shortlen = filtlen - (char) (i - inp_length - 4);
      offset = (char) (i - inp_length - 4);
      *output_sequence++ = DotP(input_sequence + inp_length + 4,
                                        filter + offset, shortlen);
    }
    else
      *output_sequence++ = DotP(input_sequence + i, filter, filtlen);
  }
}
int DecomposeBranches(double *In, int Inlen, double *Lfilter,
    double *Hfilter, char filtlen, double *OutL, double *OutH)
  /* Take input data and filters and form two branches of have the
     original length. Length of branches is returned */
{
  ConvolveDec2(In, Inlen, Lfilter, filtlen, OutL);
  ConvolveDec2(In, Inlen, Hfilter, filtlen, OutH);
  return (Inlen / 2);
}
void WaveletDecomposition(double *InData, int Inlength, double *Lfilter,
      double *Hfilter, char filtlen, char levels, double **OutData)
/* Assumes input data has 2 ^ (levels) data points/unit interval. First InData
is input signal; all others are intermediate approximation coefficients */
{
  char i;
  for (i = 0; i < levels; i++)
  {
    Inlength = DecomposeBranches(InData, Inlength, Lfilter, Hfilter,
                        filtlen, OutData[2 * i], OutData[(2 * i) + 1]);
    InData = OutData[2 * i];
  }
}
double DotpEven(double *data, double *filter, char filtlen)
{
  char i;
  double sum;
  sum = 0.0;
  for (i = 0; i < filtlen; i += 2)
    sum += *data-- * filter[i]; /* decreasing data pointer is moving */
                                /* backward in time */
  return sum;
}
double DotpOdd(double *data, double *filter, char filtlen)
{
  char i;
  double sum;
  sum = 0.0;
  for (i = 1; i < filtlen; i += 2)

    sum += *data-- * filter[i];  /* decreasing data pointer is moving */
                                 /* backward in time */
  return sum;
}
void ConvolveInt2(double *input_sequence, int inp_length, double *filter,
          char filtlen, char sum_output, double *output_sequence)
  /* insert zeros between each element of the input sequence and
     convolve with the filter to interpolate the data */
{
  int i;
  if (sum_output) /* summation with previous convolution if true */
  {
      /* every other dot product interpolates the data */
    for(i = (filtlen / 2) - 1; i < inp_length + filtlen - 2; i++)
    {
      *output_sequence++ += DotpOdd(input_sequence + i, filter, filtlen);
      *output_sequence++ += DotpEven(input_sequence + i + 1, filter, filtlen);
    }
    *output_sequence++ += DotpOdd(input_sequence + i, filter, filtlen);
  }
  else /* first convolution of pair if false */
  {
      /* every other dot product interpolates the data */
    for(i = (filtlen / 2) - 1; i < inp_length + filtlen - 2; i++)
    {
      *output_sequence++ = DotpOdd(input_sequence + i, filter, filtlen);
      *output_sequence++ = DotpEven(input_sequence + i + 1, filter, filtlen);
    }
    *output_sequence++ = DotpOdd(input_sequence + i, filter, filtlen);
  }
}
int ReconstructBranches(double *InL, double *InH, int Inlen,
        double *Lfilter, double *Hfilter, char filtlen, double *Out)
  /* Take input data and filters and form two branches of have
     original length. length of branches is returned */
{
  ConvolveInt2(InL, Inlen, Lfilter, filtlen, 0, Out);
  ConvolveInt2(InH, Inlen, Hfilter, filtlen, 1, Out);
  return Inlen * 2;
}
void WaveletReconstruction(double **InData, int Inlength, double *Lfilter,
              double *Hfilter, char filtlen, char levels, double *OutData)
  /* assumes that input data has 2 ^ (levels) data points per unit interval */
{
  double *Output;
  char i;
  Inlength = Inlength >> levels;
  /* Destination of all but last branch reconstruction is the next
     higher intermediate approximation */
  for (i = levels - 1; i > 0; i--)
  {
    Output = InData[2 * (i - 1)];
    Inlength = ReconstructBranches(InData[2 * i], InData[(2 * i) + 1],
                          Inlength, Lfilter, Hfilter, filtlen, Output);
  }
  /* Destination of the last branch reconstruction is the output data */
  ReconstructBranches(InData[0], InData[1], Inlength, Lfilter, Hfilter,
                                                      filtlen, OutData);
}
double CalculateMSE(double *DataSet1, double *DataSet2, int length)
{
  /* calculate mean squared error of output of reconstruction as
     compared to the original input data */
  int i;
  double pointdiff, topsum, botsum;
  topsum = botsum = 0.0;
  for (i = 0; i < length; i++)
  {
    pointdiff = DataSet1[i] - DataSet2[i];
    topsum += pointdiff * pointdiff;
    botsum += DataSet1[i] * DataSet1[i];
  }
  return topsum / botsum;
}





<a name="00af_000f">
<a name="00af_0010">
[LISTING THREE]
<a name="00af_0010">

/* WAVE_MGT.H */
double **BuildTreeStorage(int inlength, int levels);
void DestroyTreeStorage(double **tree, int levels);
void TreeCopy(double **TreeDest, double **TreeSrc, int siglen, int levels);
void TreeZero(double **Tree, int siglen, int levels);
void ZeroTreeDetail(double **Tree, int siglen, int levels);
/* WAVELET.H */
void WaveletCoeffs(double alpha, double beta, double *wavecoeffs);
char MakeWaveletFilters(double *wavecoeffs, double *Lfilter,
                double *Hfilter, wavetype transform);
double DotP(double *data, double *filter, char filtlength);
void ConvolveDec2(double *input_sequence, int inp_length,
        double *filter, char filtlen, double *output_sequence);
int DecomposeBranches(double *In, int Inlen, double *Lfilter,
      double *Hfilter, char filtlen, double *OutL, double *OutH);
void WaveletDecomposition(double *InData, int Inlength, double *Lfilter,
      double *Hfilter, char filtlen, char levels, double **OutData);
double DotpEven(double *data, double *filter, char filtlength);
double DotpOdd(double *data, double *filter, char filtlength);
void ConvolveInt2(double *input_sequence, int inp_length, double *filter,
        char filtlen, char sum_output, double *output_sequence);
int ReconstructBranches(double *InL, double *InH, int Inlen,
      double *Lfilter, double *Hfilter, char filtlen, double *Out);
void WaveletReconstruction(double **InData, int Inlength, double *Lfilter,
      double *Hfilter, char filtlen, char levels, double *OutData);
double CalculateMSE(double *DataSet1, double *DataSet2, int length);


Copyright © 1992, Dr. Dobb's Journal


Related Reading


More Insights






Currently we allow the following HTML tags in comments:

Single tags

These tags can be used alone and don't need an ending tag.

<br> Defines a single line break

<hr> Defines a horizontal line

Matching tags

These require an ending tag - e.g. <i>italic text</i>

<a> Defines an anchor

<b> Defines bold text

<big> Defines big text

<blockquote> Defines a long quotation

<caption> Defines a table caption

<cite> Defines a citation

<code> Defines computer code text

<em> Defines emphasized text

<fieldset> Defines a border around elements in a form

<h1> This is heading 1

<h2> This is heading 2

<h3> This is heading 3

<h4> This is heading 4

<h5> This is heading 5

<h6> This is heading 6

<i> Defines italic text

<p> Defines a paragraph

<pre> Defines preformatted text

<q> Defines a short quotation

<samp> Defines sample computer code text

<small> Defines small text

<span> Defines a section in a document

<s> Defines strikethrough text

<strike> Defines strikethrough text

<strong> Defines strong text

<sub> Defines subscripted text

<sup> Defines superscripted text

<u> Defines underlined text

Dr. Dobb's encourages readers to engage in spirited, healthy debate, including taking us to task. However, Dr. Dobb's moderates all comments posted to our site, and reserves the right to modify or remove any content that it determines to be derogatory, offensive, inflammatory, vulgar, irrelevant/off-topic, racist or obvious marketing or spam. Dr. Dobb's further reserves the right to disable the profile of any commenter participating in said activities.

 
Disqus Tips To upload an avatar photo, first complete your Disqus profile. | View the list of supported HTML tags you can use to style comments. | Please read our commenting policy.