Wavelet Analysis

Creating the Wavelet Decomposition

The following SAS code starts the wavelet analysis:

  
    %wavginit; 
    proc iml; 
       %wavinit; 
 

Notice that the previous code segment includes two SAS macro calls. You can use the IML wavelet functions without using the WAVGINIT and WAVINIT macros. The macros are called to initialize and load IML modules that you can use to produce several standard wavelet diagnostic plots. These macros have been provided as autocall macros that you can invoke directly in your SAS code.

The WAVGINIT macro must be called prior to invoking PROC IML. This macro defines several macro variables that are used to adjust the size, aspect ratio, and font size for the plots produced by the wavelet plot modules. This macro can also take several optional arguments that control the positioning and size of the wavelet diagnostic plots. See the section "Obtaining Help for the Wavelet Macros and Modules" for details about getting help about this macro call.

The WAVINIT macro must be called from within PROC IML. It loads the IML modules that you can use to produce wavelet diagnostic plots. This macro also defines symbolic macro variables that you can use to improve the readability of your code.

The following statements read the absorbance variable into an IML vector:

  
       use quartzInfraredSpectrum; 
       read all var{Absorbance} into absorbance;
 

You are now in a position to begin the wavelet analysis. The first step is to set up the options vector that specifies which wavelet and what boundary handling you want to use. You do this as follows:

  
       optn              =  &waveSpec;   /* optn=j(1,4,.); */ 
       optn[&family]     =  &daubechies; /* optn[3] = 1;   */ 
       optn[&member]     =  3;           /* optn[4] = 3;   */ 
       optn[&boundary]   =  &polynomial; /* optn[1] = 2;   */ 
       optn[&degree]     =  &linear;     /* optn[2] = 1;   */
 

These statements use macro variables that are defined in the WAVINIT macro. The equivalent code without using these macro variables is given in the adjacent comments. As indicated by the suggestive macro variable names, this options vector specifies that the wavelet to be used is the third member of the Daubechies wavelet family and that boundaries are to be handled by extending the signal as a linear polynomial at each endpoint.

The next step is to create the wavelet decomposition with the following call:

  
       call wavft(decomp,absorbance,optn);
 

This call computes the wavelet transform specified by the vector optn of the input vector absorbance. The specified transform is encapsulated in the vector decomp. This vector is not intended to be used directly. Rather you use this vector as an argument to other IML wavelet subroutines and plot modules. For example, you use the WAVPRINT subroutine to print the information encapsulated in a wavelet decomposition. The following code produces the output in Figure 16.2.

  
       call wavprint(decomp,&summary); 
       call wavprint(decomp,&detailCoeffs,1,4);
 

Decomposition Summary Decomposition Name DECOMP Wavelet Family Daubechies Extremal Phase Family Member 3 Boundary Treatment Recursive Linear Extension Number of Data Points 850 Start Level 0   Wavelet Detail Coefficients for DECOMP Translate Level 1 Level 2 Level 3 Level 4 0 -1.71343E-9 1.36819E-10 -6.6097E-12 5.23868E-11 1 1340085.30 -128245.70 191.084707 4501.36 2   62636.70 6160.27 -1358.23 3   -238445.36 -54836.56 -797.724143 4     39866.95 676.034389 5     -28836.85 -5166.59 6     223421.00 -6088.99 7       -5794.67 8       30144.74 9       -3903.53 10       638.063264 11       -10803.45 12       33616.35 13       -50790.30

Figure 16.2: Output of WAVPRINT CALLS

Usually such displayed output is of limited use. More frequently you want to represent the transformed data graphically or use the results in further computational routines. As an example, you can estimate the noise level of the data by using a robust measure of the standard deviation of the highest-level detail coefficients, as demonstrated in the following statements:

  
       call wavget(tLevel,decomp,&topLevel); 
       call wavget(noiseCoeffs,decomp,&detailCoeffs,tLevel-1); 
  
       noiseScale=mad(noiseCoeffs,"nmad"); 
       print "Noise scale = " noiseScale;
 

The result is shown in Figure 16.3.

NOISESCALE Noise scale = 169.18717

Figure 16.3: Scale of Noise in the Absorbance Data

The first WAVGET call is used to obtain the top level number in the wavelet decomposition decomp. The highest level of detail coefficients is defined at one level below the top level in the decomposition. The second WAVGET call returns these coefficients in the vector noiseCoeffs. Finally, the MAD function computes a robust estimate of the standard deviation of these coefficients.