Wavelet Analysis
/****************************************************************/
/* S A S S A M P L E L I B R A R Y */
/* */
/* NAME: WAVELET */
/* TITLE: Wavelet Analysis */
/* PRODUCT: IML */
/* SYSTEM: ALL */
/* KEYS: WAVELETS SIGNAL PROCESSING */
/* PROCS: IML */
/* DATA: */
/* */
/* SUPPORT: FRWICK UPDATE: FEB2017 */
/* REF: */
/* MISC: */
/* */
/****************************************************************/
/* This example demonstrates wavelet analysis applied to an */
/* FT-IR spectrum of quartz (Sullivan, 2000). */
/* */
/* Data Source: http://www.che.utexas.edu/~dls/ir/ir_dir.htm */
/* */
/* The following DATA step creates a data set containing the */
/* spectrum,expressed as an absorbance value for each of */
/* 850 wave numbers. */
data quartzInfraredSpectrum;
WaveNumber=4000.6167786 - _N_ *4.00084378;
input Absorbance @@;
datalines;
4783 4426 4419 4652 4764 4764 4621 4475 4430 4618
4735 4735 4655 4538 4431 4714 4738 4707 4627 4523
4512 4708 4802 4811 4769 4506 4642 4799 4811 4732
4583 4676 4856 4868 4796 4849 4829 4677 4962 4994
4924 4673 4737 5078 5094 4987 4632 4636 5010 5166
5166 4864 4547 4682 5161 5291 5143 4684 4662 5221
5640 5640 5244 4791 4832 5629 5766 5723 5121 4690
5513 6023 6023 5503 4675 5031 6071 6426 6426 5723
5198 5943 6961 7135 6729 5828 6511 7500 7960 7960
7299 6484 7257 8180 8542 8537 7154 7255 8262 8898
8898 8263 7319 7638 8645 8991 8991 8292 7309 8005
9024 9024 8565 7520 7858 8652 8966 8966 8323 7513
8130 8744 8879 8516 7722 8099 8602 8729 8726 8238
7885 8350 8600 8603 8487 7995 8194 8613 8613 8408
7953 8236 8696 8696 8552 8102 7852 8570 8818 8818
8339 7682 8535 9038 9038 8503 7669 7794 8864 9163
9115 8221 7275 8012 9317 9317 8512 7295 7623 9021
9409 9338 8116 6860 7873 9282 9490 9191 7012 7392
9001 9483 9457 8107 6642 7695 9269 9532 9246 7641
6547 8886 9457 9457 8089 6535 7537 9092 9406 9178
7591 6470 7838 9156 9222 7974 6506 7360 8746 9057
8877 7455 6504 7605 8698 8794 8439 7057 7202 8240
8505 8392 7287 6634 7418 8186 8229 7944 6920 6829
7499 7949 7831 7057 6866 7262 7626 7626 7403 6791
7062 7289 7397 7397 7063 6985 7221 7221 7199 6977
7088 7380 7380 7195 6957 6847 7426 7570 7508 6952
6833 7489 7721 7718 7254 6855 7132 7914 8040 7880
7198 6864 7575 8270 8229 7545 7036 7637 8470 8570
8364 7591 7413 8195 8878 8878 8115 7681 8313 9102
9185 8981 8283 8197 8932 9511 9511 9101 8510 8670
9686 9709 9504 8944 8926 9504 9964 9964 9627 9212
9366 9889 10100 9939 9540 9512 9860 10121 10121 9828
9567 9513 9782 9890 9851 9510 9385 9339 9451 9451
9181 9076 9015 8960 9014 8957 8760 8760 8602 8584
8584 8459 8469 8373 8279 8327 8282 8341 8341 8155
8260 8260 8250 8350 8245 8358 8403 8355 8490 8490
8439 8689 8689 8621 8680 8661 8897 9028 8900 8873
8873 9187 9377 9377 9078 9002 9147 9635 9687 9535
9127 9242 9824 9928 9775 9200 9047 9572 10102 10102
9631 9024 9209 10020 10271 9830 9062 9234 10154 10483
10453 9582 9011 9713 10643 10701 10372 9368 9857 10865
10936 10572 9574 9691 10820 11452 11452 10623 9903 10787
11931 12094 11302 10604 11458 12608 12808 12589 11629 11795
12863 13575 13575 12968 12498 13268 14469 14469 13971 13727
14441 15334 15515 15410 14986 15458 16208 16722 16722 16618
17061 17661 18089 18089 18184 18617 19015 19467 19633 19830
20334 20655 20947 21347 21756 22350 22584 22736 22986 23412
24126 24498 24501 24598 24986 25729 26356 26356 26271 26754
27624 28162 28162 28028 28305 29223 30073 30219 30185 30308
31831 32699 32819 32793 33320 34466 35600 36038 36086 36518
37517 38765 39462 39681 40209 41243 42274 42772 42876 43172
43929 44842 45351 45395 45551 46035 46774 47353 47353 47362
47908 48539 48936 48978 49057 49497 50101 50670 50914 51134
51603 52276 53007 53399 53769 54281 54815 54914 55365 55874
56180 56272 56669 57076 57422 57458 57525 57681 57679 57318
57318 57181 57417 57409 57144 57047 56377 56551 56483 56098
56034 55598 55364 55364 55146 54904 54990 55501 55533 55362
54387 55340 55240 54748 53710 55346 55795 55795 55060 55945
55945 55753 56759 56859 57509 56741 56273 56961 58566 58566
58104 59275 59275 59051 59090 59461 60362 60560 61103 61272
61380 61878 62067 62237 62214 61182 61532 62173 62253 60473
61346 63143 63378 61519 61753 63078 63841 63841 62115 61227
63237 63237 61338 63951 63951 63604 63633 64625 65135 64976
63630 63494 63834 63338 63218 62324 64131 64234 65122 64551
64127 64415 64621 64621 63142 65344 65585 65476 65074 64714
63803 65085 65085 65646 65646 64851 65390 65390 64997 65541
65587 65682 65952 65952 65390 65702 65846 65734 65734 65628
65509 65571 65636 65636 65620 65487 65544 65547 65738 65758
65711 65360 65362 65362 65231 65333 65453 65473 65435 65302
65412 65412 65351 65242 65242 65170 65221 65297 65297 65202
65177 65183 65184 65179 65209 65209 65144 65134 65113 65009
64919 64945 64988 64988 64856 64686 64529 64370 64282 64233
64169 63869 63685 63480 63373 63349 63307 63131 63017 62885
62736 62736 62706 62666 62622 62671 62781 62853 62950 63106
63135 63141 63220 63263 63489 63807 63966 64132 64294 64612
64841 64985 65159 65204 65259 65540 65707 65749 65732 65719
65820 65895 65925 65925 65888 65937 66059 66109 66109 66078
66007 65897 65897 65747 65490 64947 64598 64363 64140 63801
63571 63395 63333 63442 63442 63339 63196 62911 62118 61795
61454 61456 61607 62025 62190 62190 62023 61780 61502 61482
61458 61320 61015 60852 60708 60684 60522 60488 60506 60640
60797 60995 61141 61141 61036 60664 60522 60017 59681 59129
58605 58035 57192 56137 54995 53586 52037 50283 48565 45419
43341 41111 36131 35377 34431 31679 29237 26898 24655 22417
19876 17244 15176 12575 10532 8180 6040 4059 2210 575
;
/*---plot of spectrum---*/
proc sgplot data=quartzInfraredSpectrum;
series x=WaveNumber y=Absorbance;
xaxis reverse min=0;
yaxis values=(0 to 70000 by 10000);
run;
ods graphics on; /* ODS graphics enabled */
proc iml;
%wavinit; /* define modules for ODS graphics */
use quartzInfraredSpectrum;
read all var "absorbance";
close;
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[°ree] = &linear; /* optn[2] = 1; */
/*---compute the discrete wavelet transform---*/
call wavft(decomp, absorbance, optn);
/*---demonstrate wavprint routine---*/
call wavprint(decomp, &summary);
call wavprint(decomp, &detailCoeffs, 1, 4);
/*---demonstrate wavget routine---*/
call wavget(tLevel, decomp, &topLevel);
call wavget(noiseCoeffs, decomp, &detailCoeffs, tLevel-1);
noiseScale = mad(noiseCoeffs, "nmad");
print noiseScale[label="Noise Scale"];
/*---demonstrate coefficientPlot module---*/
call coefficientPlot(decomp) header="Quartz Spectrum";
/*---show wavhelp macro usage---*/
%wavhelp(coefficientPlot);
call coefficientPlot(decomp, , 7, , 'uniform', "Quartz Spectrum");
call coefficientPlot(decomp, &SureShrink, 6, , , "Quartz Spectrum");
/*---demonstrate multiresolution approximation plot module---*/
call mraApprox(decomp, , 3, , "Quartz Spectrum");
call mraApprox(decomp, , 7, 7, "Quartz Spectrum");
call mraApprox(decomp, &SureShrink, 10, 10, "Quartz Spectrum");
*---demonstrate multiresolution decomposition plot module---*/
call mraDecomp(decomp, ,5, , , "Quartz Spectrum");
/*---demonstrate scalogram plot module---*/
call scalogram(decomp, &SureShrink, , , 0.25, 'log',"Quartz Spectrum");
call scalogram(decomp, &SureShrink, 6, , , , "Quartz Spectrum");
call scalogram(decomp, , 6, , , , "Quartz Spectrum");
/*---demonstrate inverse wavelet transformation---*/
call wavift(reconstructedAbsorbance, decomp);
errorSS=ssq(absorbance-reconstructedAbsorbance);
print errorSS[label="Reconstruction Error Sum of Squares"];
call wavift(smoothedAbsorbance, decomp, &Sureshrink);
create temp var "smoothedAbsorbance";
append;
close temp;
quit;
data Spectrum;
merge quartzInfraredSpectrum temp;
run;
proc sgplot data=Spectrum;
series x=WaveNumber y=smoothedAbsorbance;
xaxis reverse min=0;
yaxis values=(0 to 70000 by 10000);
run;
proc iml;
%wavinit;
use quartzInfraredSpectrum;
read all var "Absorbance";
close;
optn = j(1, 4, .);
optn[1] = 0;
optn[3] = 1;
optn[4] = 3;
optn = &waveSpec;
optn[&boundary] = &zeroExtension;
optn[&family] = &daubechies;
optn[&member] = 3;
call wavft(decomposition, absorbance, optn);
call wavget(n, decomposition, 1);
call wavget(fWavelet, decomposition, 8);
call wavget(n, decomposition, &numPoints);
call wavget(fWavelet, decomposition, &fatherWavelet);