TECHNICAL REPORT SAS/IML Software: Changes And Enhancements, Release 8

8 Chapter 1. Wavelet AnalysisThe WAVGINIT macro must be called prior to invoking PROC IML. This macrodefines several macro variables that are used to adjust the size, aspect ratio, and fontsize for the plots produced by the wavelet plot modules. This macro can also takeseveral optional arguments that control the positioning and and size of the waveletdiagnostic plots. See “Obtaining Help for Wavelet Modules and Macros” on page 32for details on 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 definessymbolic 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 theoptions vector that specifies which wavelet and what boundary handling you want touse. You do this as ry]optn[°ree] /*/*optn j(1,4,.);optn[3] 1;optn[4] 3;optn[1] 3;optn[2] 1;*/*/*/*/*/These statements use macro variables that are defined in the WAVINIT macro. Theequivalent code without using these macro variables is given in the adjacent comments. As indicated by the suggestive macro variable names, this options vectorspecifies that the wavelet to be used is the third member of the Daubechies waveletfamily 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 inputvector 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 anargument to other IML wavelet subroutines and plot modules. For example, youuse the WAVPRINT subroutine to print the information encapsulated in a waveletdecomposition. The following code produces output in Figure 1.2.call wavprint(decomp,&summary);call wavprint(decomp,&detailCoeffs,1,4);

Creating the Wavelet DecompositionDecomposition SummaryDecomposition NameWavelet FamilyFamily MemberBoundary TreatmentNumber of Data PointsStart LevelDECOMPDaubechies Extremal Phase3Recursive Linear Extension8500Wavelet Detail Coefficients for DECOMPTranslateLevel 1Level 2Level 3Level 616.35-50790.30Figure 1.2.Output of WAVPRINT CALLSUsually such displayed output is of limited use. More frequently you will want torepresent the transformed data graphically or use the results in further computationalroutines. As an example, you can estimate the noise level of the data using a robustmeasure of the standard deviation of the highest level detail coefficients, as demonstrated in the following statements:call wavget(tLevel,decomp,&topLevel);call noiseScale mad(noiseCoeffs,"nmad");print "Noise scale " noiseScale;The result is shown in Figure 1.3;NOISESCALENoise scale Figure 1.3.169.18717Scale of Noise in the Absorbance DataThe first WAVGET call is used to obtain the top level number in the wavelet decomposition decomp. The highest level of detail coefficients are defined at one levelbelow the top level in the decomposition. The second WAVGET call returns thesecoefficients in the vector noiseCoeffs. Finally, the MAD function computes a robustestimate of the standard deviation of these coefficients. 9

10 Chapter 1. Wavelet AnalysisWavelet Coefficient PlotsDiagnostic plots greatly facilitate the interpretation of a wavelet decomposition. Onestandard plot is the detail coefficients arranged by level. Using a module included bythe WAVINIT macro call, you can produce the plot shown in Figure 1.5 as follows:call coefficientPlot(decomp, , , , ,"Quartz Spectrum");The first argument specifies the wavelet decomposition and is required. All otherarguments are optional and need not be specified. You can use the WAVHELP macroto obtain a description of the arguments of this and other wavelet plot modules. TheWAVHELP macro is defined in autocall the WAVINIT macro. For example, invokingthe WAVHELP macro as follows writes the calling information shown in Figure 1.4to the SAS log.%wavhelp(coefficientPlot);coefficientPlot ModuleFunction: Plots wavelet detail coefficientsUsage: call n - (required) valid wavelet decompostion producedby the IML subroutine WAVFTthreshopt- (optional) numeric vector of 4 elementsspecifying thresholding to be usedDefault: no thresholdingstartLevel- (optional) numeric scalar specifying the lowestlevel to be displayed in the plotDefault: start level of decompositionendLevel- (optional) numeric scalar specifying the highestlevel to be displayed in the plotDefault: end level of decompositionhowScaled- (optional) character: ’absolute’ or ’uniform’specifies coefficients are scaled uniformlyDefault: independent level scalingheader- (optional) character string specifying a headerDefault: no headerFigure 1.4.Log Output Produced by %wavhelp(coefficientPlot) Call

Wavelet Coefficient PlotsFigure 1.5.Detail Coefficients Scaled by LevelIn this plot the detail coefficients at each level are scaled independently. The oscillations present in the absorbance data are captured in the detail coefficients at levels7, 8, and 9. The following statement produces a coefficient plot of just these higherlevel detail coefficients and shows them scaled uniformly.call coefficientPlot(decomp, ,7, ,’uniform’,"Quartz Spectrum");The plot is shown in Figure 1.6. 11

12 Chapter 1. Wavelet AnalysisFigure 1.6.Uniformly Scaled Detail CoefficientsAs noted earlier, noise in the data is captured in the detail coefficients, particularly inthe small coefficients at higher levels in the decomposition. By zeroing or shrinkingthese coefficients, you can get smoother reconstructions of the input data. This isdone by specifying a threshold value for each level of detail coefficients and thenzeroing or shrinking all the detail coefficients below this threshold value. The IMLwavelet functions and modules support several policies for how this thresholdingis performed as well as for selecting the thresholding value at each level. See the“WAVIFT Call” on page 26 for details.An options vector is used to specify the desired thresholding; several standard choicesare predefined as macro variables in the WAVINIT module. The following statementsproduce the detail coefficient plot with the “SureShrink” thresholding algorithm ofDonoho and Johnstone (1995).call coefficientPlot(decomp,&SureShrink,6,, ,"Quartz Spectrum");The plot is shown in Figure 1.7.

Multiresolution Approximation PlotsFigure 1.7.Thresholded Detail CoefficientsYou can see that “SureShrink” thresholding has zeroed some of the detail coefficientsat the higher levels but the larger coefficients that capture the oscillation in the data arestill present. Consequently, reconstructions of the the input signal using the thresholded detail coefficients will still capture the essential features of the data, but will besmoother as much of the very fine scale detail has been eliminated.Multiresolution Approximation PlotsOne way of presenting reconstructions is in a multiresolution approximation plot.In this plot reconstructions of the input data are shown by level. At any level thereconstruction at that level uses only the detail and scaling coefficients defined belowthat level.The following statement produces such a plot, starting at level 3:call mraApprox(decomp, ,3, ,"Quartz Spectrum");The results are shown in Figure 1.8. 13

14 Chapter 1. Wavelet AnalysisFigure 1.8.Multiresolution ApproximationYou can see that even at level 3, the basic form of the input signal has been captured.As noted earlier, the oscillation present in the absorbance data is captured in the detailcoefficients above level 7. Thus, the reconstructions at level 7 and below are largelyfree of these oscillation since they do not use any of the higher detail coefficients. Youcan confirm this observation by plotting just this level in the multiresolution analysisas follows:call mraApprox(decomp, ,7,7,"Quartz Spectrum");The results are shown in Figure 1.9.

Multiresolution Approximation PlotsFigure 1.9.Level 7 of the Multiresolution ApproximationYou can also plot the multiresolution approximations obtained with thresholded detailcoefficients. For example, the following statement plots the top level reconstructionobtained using the “SureShrink” threshold:call mraApprox(decomp,&SureShrink,10,10,"Quartz Spectrum");The results are shown in Figure 1.10. 15

16 Chapter 1. Wavelet AnalysisFigure 1.10.Top Level of Multiresolution Approximation with SureShrink Thresholding AppliedNote that the high frequency oscillation is still present in the reconstruction even with“SureShrink” thresholding applied.Multiresolution Decomposition PlotsA related plot is the multiresolution decomposition plot, which shows the detail coefficients at each level. For convenience, the starting level reconstruction at the lowestlevel of the plot and the reconstruction at the highest level the plot are also included.Adding suitably scaled versions of all the detail levels to the starting level reconstruction recovers the final reconstruction. The following statement produces such a plot,yielding the results shown in Figure 1.11.call mraDecomp(decomp, ,5, , ,"Quartz Spectrum");

Wavelet ScalogramsFigure 1.11.Multiresolution DecompositionWavelet ScalogramsWavelet scalograms communicate the time frequency localization property of the discrete wavelet transform. In this plot each detail coefficient is plotted as a filled rectangle whose color corresponds to the magnitude of the coefficient. The location andsize of the rectangle are related to the time interval and the frequency range for thiscoefficient. Ccoefficients at low levels are plotted as wide and short rectangles toindicate that they localize a wide time interval but a narrow range of frequencies inthe data. In contrast, rectangles for coefficients at high levels are plotted thin andtall to indicate that they localize small time ranges but large frequency ranges in thedata. The heights of the rectangles grow as a power of 2 as the level increases. Ifyou include all levels of coefficients in such a plot, the heights of the rectangles atthe lowest levels are so small that they will not be visible. You can use an option toplot the heights of the rectangles on a logarithmic scale. This results in rectanglesof uniform height but requires that you interpret the frequency localization of thecoefficients with care.The following statement produces a scalogram plot of all levels with “SureShrink”thresholding applied:call scalogram(decomp,&SureShrink, , ,0.25,’log’,"Quartz Spectrum");The sixth argument specifies that the rectangle heights are to be plotted on a logarithmic scale. The role of the fifth argument (0:25) is to amplify the magnitude of thesmall detail coefficients. This is necessary since the detail coefficients at the lower 17

18 Chapter 1. Wavelet Analysislevels are orders of magnitude larger than those at the higher levels. The amplification is done by first scaling the magnitudes of all detail coefficients to lie in theinterval [0; 1] and then raising these scaled magnitudes to the power 0:25. Note thatsmaller powers yield larger amplification of the small detail coefficient magnitudes.The default amplification is 1 3.The results are shown in Figure 1.12.Figure 1.12.Scalogram Showing All LevelsThe bar on the left-hand side of the scalogram plot indicates the overall energy ofeach level. This energy is defined as the sum of the squares of the detail coefficientsfor each level. These energies are amplified using the same algorithm for amplifyingthe detail coefficient magnitudes. The energy bar in Figure 1.12 shows that higherenergies occur at the lower levels whose coefficients capture the gross features ofthe data. In order to interpret the finer-scale details of the data it is helpful to focuson just the higher levels. The following statement produces a scalogram for levels 6and above without using a logarithmic scale for the rectangle heights, and using thedefault coefficient amplification.call scalogram(decomp,&SureShrink,6, , , ,"Quartz Spectrum");The result is shown in Figure 1.13.

Wavelet ScalogramsFigure 1.13.Scalogram of Levels 6 and Above Using SureShrink ThresholdingThe scalogram in Figure 1.13 reveals that most of the energy of the oscillation inthe data is captured in the detail coefficients at level 8. Also note that many of thecoefficients at the higher levels are set to zero by “SureShrink” thresholding. You canverify this by comparing Figure 1.13 with Figure 1.14, which shows the corresponding scalogram except that no thresholding is done. The following statement producesFigure 1.14:call scalogram(decomp, ,6, , , ,"Quartz Spectrum"); 19

20 Chapter 1. Wavelet AnalysisFigure 1.14.Scalogram of Levels 6 and Above Using No ThresholdingReconstructing the Signal from the Wavelet DecompositionYou can use the WAVIFT subroutine to invert a wavelet transformation computedusing the WAVFT subroutine. If no thresholding is specified, then up to numericalrounding error this inversion is exact. The following statements provide an illustration of this:call wavift(reconstructedAbsorbance,decomp);errorSS ssq(absorbance-reconstructedAbsorbance);print "The reconstruction error sum of squares " errorSS;The output is shown in Figure 1.15.ERRORSSThe reconstruction error sum of squares Figure 1.15.1.288E-16Exact Reconstruction Property of WAVIFTUsually you use the WAVIFT subroutine with thresholding specified. This yields asmoothed reconstruction of the input data. You can use the following statements tocreate a smoothed reconstruction of absorbance and add this variable to the QuartzInfraredSpectrum data set.call te temp from smoothedAbsorbance[colname ’smoothedAbsorbance’];append from smoothedAbsorbance;

Reconstructing the Signal from the Wavelet Decompositionclose temp;quit;data quartzInfraredSpectrum;set quartzInfraredSpectrum;set temp;run;The following statements produce the line plot of the smoothed absorbance datashown in Figure 1.16:symbol1 c black i join v none;proc gplot data quartzInfraredSpectrum;plot smoothedAbsorbance*WaveNumber/hminor 0vminor 0vaxis axis1hreverse frame;axis1 label ( r 0 a 90 );run;Figure 1.16.Smoothed FT-IR Spectrum of QuartzYou can see by comparing Figure 1.1 with Figure 1.16 that the wavelet smooth of theabsorbance data has preserved all the essential features of this data. 21

22 Chapter 1. Wavelet AnalysisSyntaxWavelet Analysis CallsWAVFT CallWAVGET CallWAVIFT CallWAVPRINT CallWAVTHRSH Callcomputes a specified wavelet transform of onedimensional datareturns requested information encapsulated in a wavelettransforminverts a wavelet transform after applying specified thresholding to the detail coefficientsdisplays requested information encapsulated in a wavelettransformapplies specified thresholding to the detail coefficients ofa wavelet transformWAVFT Callcomputes fast wavelet transformCALL WAVFT(decomp, data, opt , levels );The Fast Wavelet Transform (WAVFT) subroutine computes a specified discretewavelet transform of the input data, using the algorithm of Mallat (1989). This transform decomposes the input data into sets of detail and scaling coefficients defined ata number of scales or “levels.”The input data are used as scaling coefficients at the top level in the decomposition.The fast wavelet transform then recursively computes a set of detail and a set ofscaling coefficients at the next lower level by respectively applying “low pass” and“high pass” conjugate mirror filters to the scaling coefficients at the current level. Thenumber of coefficients in each of these new sets is approximately half the number ofscaling coefficients at the level above them. Depending on the filters being used, anumber of additional scaling coefficients, known as boundary coefficients, may beinvolved. These boundary coefficients are obtained by extending the sequence ofinterior scaling coefficients

with SAS Institute and the restrictions set forth in FAR 52.227-19 Commercial Computer Software-Restricted Rights (June 1987). SAS Institute Inc., SAS Campus Drive, Cary, North Carolina 27513. 1st printing, January 2001 SAS and all other SAS Institute Inc. product or service names are registered trademarks