home *** CD-ROM | disk | FTP | other *** search
/ Media Share 13 / mediashare_13.zip / mediashare_13 / ZIPPED / PROGRAM / APR94_1.ZIP / WAVE3.ASC < prev   
Text File  |  1994-02-27  |  11KB  |  226 lines
  1. _The Wavelet Packet Transform_
  2. by Mac A. Cody
  3.  
  4. Listing One
  5.  
  6. /* CONVOLVE.C */
  7. /* Copyright (C) 1993 Mac A. Cody - All rights reserved */
  8.  
  9. #include "convolve.h"
  10.  
  11. /* DualConvDec2 - Convolution of data array with decomposition filters
  12.                   followed by decimation by two.
  13.    Input(s): REAL_TYPE *src    - Pointer to source data sample array.
  14.              REAL_TYPE *htilda - Pointer to lowpass filter array.
  15.              REAL_TYPE *gtilda - Pointer to highpass filter array.
  16.              short srclen      - length of source data array.
  17.              short filtlen     - length of filter arrays.
  18.    Output(s): REAL_TYPE *adst  - Pointer to approximation data sample array.
  19.               REAL_TYPE *ddst  - Pointer to detail data sample array.
  20. */
  21. void DualConvDec2(REAL_TYPE *src, REAL_TYPE *adst, REAL_TYPE *ddst,
  22.            REAL_TYPE *htilda, REAL_TYPE *gtilda, short srclen, short filtlen)
  23. {
  24.   short i, j, brklen, endlen;
  25.   REAL_TYPE adot_p, ddot_p;
  26.   REAL_TYPE *head_src, *lp_fltr, *hp_fltr;
  27.   brklen = 1; /* initial break in dot product is after first element */
  28.   /* perform truncated dot products until filter no longer hangs off end of
  29.      array; decimation by two handled by two-element shift; break count
  30.      increases by two on every iteration */
  31.   for(j = 0; j < filtlen; j += 2, brklen += 2)
  32.   {
  33.     head_src = src + j; /* point to appropriate initial element at head */
  34.     lp_fltr = htilda;   /* set up pointer to lowpass filter */
  35.     hp_fltr = gtilda;   /* set up pointer to highpass filter */
  36.     adot_p = *head_src * *lp_fltr++;   /* initial lowpass product of head */
  37.     ddot_p = *head_src-- * *hp_fltr++; /* initial highpass product of head */
  38.     for(i = 1; i < brklen; i++) /* perform remaining products of head */
  39.     {
  40.       adot_p += *head_src * *lp_fltr++;
  41.       ddot_p += *head_src-- * *hp_fltr++;
  42.     }
  43.     *adst++ = adot_p; /* save the completed lowpass dot product */
  44.     *ddst++ = ddot_p; /* save the completed highpass dot product */
  45.   }
  46.   endlen = srclen + filtlen - 2; /* find total length of array */
  47.   /* perform convolution to the physical end of the array
  48.      with a simple dot product loop */
  49.   for(; j <= endlen; j += 2)
  50.   {
  51.     head_src = src + j; /* point to appropriate initial element at head */
  52.     lp_fltr = htilda;   /* set up pointer to lowpass filter */
  53.     hp_fltr = gtilda;   /* set up pointer to highpass filter */
  54.     adot_p = *head_src * *lp_fltr++;   /* initial lowpass product */
  55.     ddot_p = *head_src-- * *hp_fltr++; /* initial highpass product */
  56.     for(i = 1; i < filtlen; i++) /* perform remaining products */
  57.     {
  58.       adot_p += *head_src * *lp_fltr++;
  59.       ddot_p += *head_src-- * *hp_fltr++;
  60.     }
  61.     *adst++ = adot_p; /* save the completed lowpass dot product */
  62.     *ddst++ = ddot_p; /* save the completed highpass dot product */
  63.   }
  64.   /* perform convolution off the physical end of the array
  65.      with a partial dot product loop, like at the beginning */
  66.   for(brklen = filtlen - 2, j = 2; brklen > 0; brklen -= 2, j += 2)
  67.   {
  68.     head_src = src + endlen; /* point to last element in array */
  69.     lp_fltr = htilda + j;    /* set up pointer to lowpass filter offset */
  70.     hp_fltr = gtilda + j;    /* set up pointer to highpass filter offset*/
  71.     adot_p = *head_src * *lp_fltr++;   /* initial lowpass product */
  72.     ddot_p = *head_src-- * *hp_fltr++; /* initial highpass product */
  73.     for(i = 1; i < brklen; i++) /* perform remaining products */
  74.     {
  75.       adot_p += *head_src * *lp_fltr++;
  76.       ddot_p += *head_src-- * *hp_fltr++;
  77.     }
  78.     *adst++ = adot_p; /* save the completed lowpass dot product */
  79.     *ddst++ = ddot_p; /* save the completed highpass dot product */
  80.   }
  81. } /* End of DualConvDec2 */
  82. /* DualConvInt2Sum - Convolution of data array with reconstruction
  83.                      filters with interpolation by two and sum together.
  84.    Input(s): REAL_TYPE *asrc   - Pointer to approximation data sample array.
  85.              REAL_TYPE *dsrc   - Pointer to detail data sample array.
  86.              REAL_TYPE *h      - Pointer to lowpass filter array.
  87.              REAL_TYPE *g      - Pointer to highpass filter array.
  88.              short srclen      - length of source data array.
  89.              short filtlen     - length of filter arrays.
  90.    Output(s): REAL_TYPE *dst  - Pointer to output data sample array.
  91. */
  92. void DualConvInt2Sum(REAL_TYPE *asrc, REAL_TYPE *dsrc, REAL_TYPE *dst,
  93.                      REAL_TYPE *h, REAL_TYPE *g, short srclen, short filtlen)
  94. {
  95.   short i, j, endlen;
  96.   REAL_TYPE dot_pe, dot_po;
  97.   REAL_TYPE *head_asrc, *head_dsrc, *lp_fltr, *hp_fltr;
  98.   endlen = srclen + filtlen - 2; /* find total length of array */
  99.   filtlen /= 2;         /* adjust filter length value for interpolation */
  100.   j = filtlen - 1;      /* start with filter 'covering' end of array */
  101.   head_asrc = asrc + j; /* point to initial element at head */
  102.   head_dsrc = dsrc + j; /* point to initial element at head */
  103.   lp_fltr = h + 1;      /* set up pointer to lowpass filter */
  104.   hp_fltr = g + 1;      /* set up pointer to highpass filter */
  105.   /* initial lowpass and highpass odd product */
  106.   dot_po = *head_asrc-- * *lp_fltr + *head_dsrc-- * *hp_fltr;
  107.   lp_fltr += 2;  hp_fltr += 2; /* skip over even filter elements */
  108.   for(i = 1; i < filtlen; i += 2) /* perform remaining products */
  109.   {
  110.     dot_po += *head_asrc-- * *lp_fltr + *head_dsrc-- * *hp_fltr;
  111.     lp_fltr += 2;  hp_fltr += 2; /* skip over even filter elements */
  112.   }
  113.   /* save the completed lowpass and highpass odd dot product */
  114.   *dst++ = dot_po;
  115.   /* perform initial convolution with a simple dot product loop */
  116.   for(j++; j <= endlen; j++)
  117.   {
  118.     head_asrc = asrc + j; /* point to appropriate initial element at head */
  119.     head_dsrc = dsrc + j; /* point to appropriate initial element at head */
  120.     lp_fltr = h;          /* set up pointer to lowpass filter */
  121.     hp_fltr = g;          /* set up pointer to highpass filter */
  122.     /* initial lowpass and highpass even product */
  123.     dot_pe = *head_asrc * *lp_fltr++ + *head_dsrc * *hp_fltr++;
  124.     /* initial lowpass and highpass odd product */
  125.     dot_po = *head_asrc-- * *lp_fltr++ + *head_dsrc-- * *hp_fltr++;
  126.     for(i = 1; i < filtlen; i++) /* perform remaining products */
  127.     {
  128.       dot_pe += *head_asrc * *lp_fltr++ + *head_dsrc * *hp_fltr++;
  129.       dot_po += *head_asrc-- * *lp_fltr++ + *head_dsrc-- * *hp_fltr++;
  130.     }
  131.     /* save the completed lowpass and highpass even dot product */
  132.     *dst++ = dot_pe;
  133.     /* save the completed lowpass and highpass odd dot product */
  134.     *dst++ = dot_po;
  135.   }
  136. } /* End of DualConvInt2Sum */
  137.  
  138.  
  139. Listing Two
  140.  
  141. /* AWPT.C */
  142. /* Copyright (C) 1993 Mac A. Cody - All rights reserved */
  143.  
  144. #include "wp_types.h"
  145. #include "awpt.h"
  146. #include "convolve.h"
  147.  
  148. /* AWPT - Aperiodic Wavelet Packet Transform: Data is assumed to be 
  149.           non-periodic, so convolutions do not wrap around arrays.
  150.           Convolution data past end of data is generated and retained
  151.           to allow perfect reconstruction of original input.
  152.    Input(s): REAL_TYPE *indata - Pointer to input data sample array.
  153.              REAL_TYPE *htilda - Pointer to lowpass filter array.
  154.              REAL_TYPE *gtilda - Pointer to highpass filter array.
  155.              short filtlen     - Length of filter arrays.
  156.    Output(s): WPTstruct *out - Pointer to transform data structure.
  157.    Note: Structure pointed to by 'out' contains:
  158.          out->levels - Number of levels in transform (short).
  159.          out->length - Length of input data sample array (short).
  160.          out->data   - Pointer to pointer of arrays of data (REAL_TYPE ***).
  161. */
  162. void AWPT(REAL_TYPE *indata, WPTstruct *out,
  163.           REAL_TYPE *htilda, REAL_TYPE *gtilda, short filtlen)
  164. {
  165.   short i, j, nodes, datalen;
  166.   REAL_TYPE *src;
  167.   levels = out->levels;
  168.   datalen = out->length; /* start with length of input array */
  169.   /* loop for all levels, halving the data length on each lower level */
  170.   for (i = 0, nodes = 2; i < levels; i++, nodes <<= 1)
  171.   {
  172.     for(j = 0; j < nodes; j += 2)
  173.     {
  174.       if(out->data[i][j] == 0) continue;
  175.       if(i == 0)      /* ... source for highest level is input data */
  176.         src = indata;
  177.       else            /* ... source is corresponding node of higher level */
  178.         src = out->data[i-1][j >> 1];
  179.       DualConvDec2(src, out->data[i][j], out->data[i][j+1],
  180.                                     htilda, gtilda, datalen, filtlen);
  181.     }
  182.     datalen /= 2; /* input length for next level is half this level */
  183.   }
  184. } /* End of AWPT */
  185. /* IAWPT - Inverse Aperiodic Wavelet Packet Transform: Data is assumed to be
  186.            non-periodic, so convolutions do not wrap arround arrays.
  187.            Convolution data past end of data is used to generate perfect
  188.            reconstruction of original input.
  189.    Input(s): WPTstruct *in     - Pointer to transform data structure.
  190.              REAL_TYPE *htilda - Pointer to lowpass filter array.
  191.              REAL_TYPE *gtilda - Pointer to highpass filter array.
  192.              short filtlen     - Length of filter arrays.
  193.    Note: Structure pointed to by 'in' contains:
  194.          in->levels - Number of levels in transform (short).
  195.          in->length - Length of output data sample array (short).
  196.          in->data   - Pointer to pointer of arrays of data (REAL_TYPE ***).
  197.    Output(s): REAL_TYPE *indata - Pointer to input data sample array.
  198. */
  199. void IAWPT(WPTstuct *in, REAL_TYPE *outdata,
  200.            REAL_TYPE *htilda, REAL_TYPE *gtilda, short filtlen)
  201. {
  202.   short i, j, levels, nodes, datalen;
  203.   REAL_TYPE *dst;
  204.   levels = in->levels;
  205.   /* start with length of lowest level input array */
  206.   datalen = in->length >> levels;
  207.   /* loop for all levels, doubling the data length on each higher level;
  208.      destination of all but the highest branch of the reconstruction
  209.      is the next higher node */
  210.   for (i = levels - 1, nodes = 1 << levels; i >= 0; i--, nodes >>= 1)
  211.   {
  212.     for(j = 0; j < nodes; j += 2)
  213.     {
  214.       if(out->data[i][j] == 0) continue;
  215.       if(i == 0) /* ... destination for highest level is input data */
  216.         dst = outdata;
  217.       else       /* ... destination is corresponding node of higher level */
  218.         dst = in->data[i - 1][j >> 1];
  219.       DualConvInt2Sum(in->data[i][j], in->data[i][j+1], dst,
  220.                                       htilda, gtilda, datalen, filtlen);
  221.     }
  222.     datalen *= 2; /* input length for next level is half this level */
  223.   }
  224. }/* End of IAWPT */
  225.  
  226.