Home

Wavelet Toolware - EBM English Books

1. Figure 3 8 DWT coefficient display after thresholding 3 3 4 Wavelet reconstruction We note that once a data set is decomposed its selected decomposition parameters are locked into the tool to prevent erroneous meaningless reconstruction When clicking on the M button a window for se lection of reconstruction parameters will pop up Toolware can perform reconstruction between different levels and their results will be stored 52 WAVELET TOOLWARE in its internal buffers with the possibility of overwriting their exist ing contents For example if an image is decomposed into four levels its DWT coefficients at levels 1 2 3 and 4 are kept in the internal buffers If a level 3 to level 0 reconstruction request is selected as in Fig 3 9 DWT coefficients stored in level 3 buffers are retrieved and then reconstructed back to level 2 subbands The new level 2 values will be saved into the level 2 buffers Then level 2 coefficients will be reconstructed into level 1 coefficients and new level 1 coefficients will be reconstructed back to level 0 If any of level 3 coefficients have been altered before the reconstruction request it is likely that old level 2 level 1 and level 0 coefficients will be overwritten in the reconstruction process Figure 3 9 Selection of DWT reconstruction As usual the user can view the new DWT values after reconstruc tion If the user exits Toolware at this point the original file on the
2. 1 6 2 Decomposition relation The relations 1 36 and 1 37 are often referred to as the reconstruc tive relations The decomposition relation is for separating a signal into components at different scales The subspace relation in an MRA Vier Vj W indicates that there exists a relation between the basis functions in these subspaces Since 2t and 2 1 are two distinct basis functions in 14 WAVELET TOOLWARE Two Scale Relation 25 20 0 500 1000 1500 2000 2500 3000 3500 4000 4500 Figure 1 7 Cubic B spline Viu there exist two sequences ap and bp in 2 such that P 2t azot k box ep t k k B 2t 1 So ax 10 t k boxap t k 1 38 k Writing these relations in their generalized form we have p 2t L X azx ep t k boxe t k 1 39 k for all Z For an arbitrary resolution level 7 this relation becomes Pej Q X 1 40 J ax ep Vt k boe_eh 2 t k 1 41 k 5 Maori Dak 1 42 k There exists a definite relationship between the two scale sequences pk qz and the decomposition sequences ag bk We will con sider these relations when different types of wavelets are discussed 1 6 MULTIRESOLUTION ANALYSIS 15 1 6 3 Development of the decomposition and reconstruction algorithms The decomposition analysis and reconstruction synthesis algorithms are the most often used algorithms in wavelet signa
3. functions are not resolvable If the time domain window is narrow enough to resolve the delta functions the frequency resolution is low We will show later that the CWT can resolve the delta functions and the frequencies simultaneously in only one pass by varying the scale parameter 2 3 ALGORITHM FOR COMPUTING THE CWT 31 STFT Toolware delta pts ATT TT all a ee eee Figure 2 1 STFT outputs for the window size equal 8 2 3 Algorithm for Computing the CWT The definition of the CWT given in part A is written here 1 se t b Wy ECD a zla z t y dt f tlt h b t dt x t h t 2 83 i where h t y i is used in 2 83 Using the convolution integral expression 2 83 can be computed by using FFT as in the case of STFT The steps are outlined as follows 1 Select a wavelet of your choice to be used for this computation 2 For a given scale a sample the wavelet at the same rate as the input data 3 Treat the data at the edge the same way as indicated in the STFT 4 Compute the FFT of the input data and that of the wavelet at scale a 5 Multiply the two FFTs and take the IFFT to obtain the CWT of the function at the scale a 32 WAVELET TOOLWARE STFT Toolware delta pts al ao ines 2048 il Figure 2 2 STFT outputs for the window size equal 64 6 Repeat steps 2 through 5 for as many values of a as the user wishes 7 Display data using a 3 D pl
4. you will see the Harr wavelet at iteration 1 Each time you click on the yellow rightward arrow the wavelet order increases by 1 That means that you will see the linear cubic and other splines and wavelets by clicking on the same button Fig 3 23 displays the cubic spline and its wavelet The Daubechies orthonormal wavelets are generated in a similar manner Some precomputed coefficients are already saved in the system for fast retrieval and export The values of these wavelets functions can be saved into text files for further processing For exam ple you could save the numerical values of wavelet functions into data files and then use the STFT tool to examine their spectra STFT Toolware delta pts Sars z uo 10 MLN LAE ar f il Figure 3 22 STFT coefficients for window size 64 Toolware Wavelet Builder Document Scaling Function Spline 0 666667 T 0 666667 f f i WaveletiSpline 4 Cubic 0 351287 4 f f 4 i ll 1 j i E S f t t L ra l j EEE wan z H Se i SN SS H 0 351287 Figure 3 23 Cubic spline and its wavelet 65 This page intentionally left blank Bibliography Haar A Zur Theorie der Orthogonalen Funktionensysteme Math Annal 69 1910 331 371 2 Chui C K Wavelet A Mathematical Tool for Signal Analysis SIAM Publ 1997 3 Chui C K An Introduction to Wavelets Academic Press 1992 i 4 Vetterli M and Kopvacevic Wavelet
5. 8 wavelet component 24 27 wavelet decomposition 23 wavelet packet 24 27 wavelet packet algorithm 24 wavelet packet series 25 wavelet packet subspace 24 wavelet packet transforms 23 wavelet packets 23 24 wavelet series 18 24 wavelet signal processing 6 18 wavelet space 21 38 wavelet subspace 10 11 19 wavelets 35 window area 9 window function 2 3 9 30 window measures 9 INDEX window width 9 window widths 9 10 Windowed Fourier Transform 2 windowing effect 3 z transform 21 73
6. PGM a II ni Figure 3 4 File selection window for 2 D images A loaded image is automatically displayed on the screen for DWT and processing see Fig 3 5 for a sculpture image file name art pgm taken from the main campus of Texas A amp M University The left most corner of this window lists the title of the tool and the cho sen file name The two buttons at the upper right corner as maximize minimize and close the current window The second row lists the available functions in the 2 D DWT tool where each of these headings is a pull down menu Although they could be activated from the pull down menu key operations of the 2 D tool are represented by thumbnail icons at the third row of its window Sse BAS ea 2c Sw The first group of four icons MEME enable opening saving and print ing a displayed image The second group of three icons start the decomposition reconstruction and processing of the selected im age Up to five levels of DWT coefficients can be stored in the internal buffer for display but only levels 0 to 3 can be processed by processing functions through the P button By clicking on one of the following icons OEI ZEI the user can display any particular decomposi 48 WAVELET TOOLWARE 2D DWT Toolware art pgm arag lalele 1AA 0 Ay ie pu p R Figure 3 5 A 2 D image in the 2 D tool 3 3 2 D DWT TOOL 49 tion level of wavelet coefficients All DWT coefficients at differe
7. Position the center of the window sample data or the center of the window at n 0 and multiply the window sample with the data set Take an N point FFT to obtain the spectrum of the function at n 0 5 Shift the window data to the right by one point and repeat step 4 to obtain the spectrum of the function at n 1 6 Repeat step 5 until the center of the window has reached the end point of the data set at n M 1 The resulting output data should be a rectangular array of M x N points 7 Display the data through either a 3 D plot or a gray scale plot Example We use an example similar to the one given in Dauchies book for illustration Let the function to be processed be f t sin 1000nt sin 2000t a 6 t t1 t ta 2 82 where a is a constant chosen to be 3 and t and tz are 0 192 and 0 196 msec respectively The function is sampled at 8000 points per second The frequencies are chosen so that when the time window width is small enough to resolve the delta functions the corresponding frequency win dow width is too wide to resolve the two frequencies accurately We use a rectangular window of different window widths 8 and 64 points respectively The STFTs of the functions for these windows are com puted using the Toolware and the results are shown in Fig 2 1 and Fig 2 2 One can easily see from these graphs that if the window in the time domain is wide enough to resolve the low frequencies the delta
8. and universal thresholding meth ods are quite simple One simply subtracts the threshold value from the magnitude of each coefficient If the difference is negative the co efficient is set to zero If the difference is positive no change is applied to the coefficient To implement the soft thresholding by using a linear function of unit slope the thresholding rule is y r ocifr o gt 0 0ifz o0 lt 0 2 99 2 10 Two Dimensional Extension of Wavelet Algo rithms We have mentioned in previous sections that the 2 D wavelets are tensor products of the 1 D scaling function and the wavelet Corresponding to the scaling function and the wavelet in one dimension are three 2 D wavelets and one 2 D scaling function at each level of resolution As a result the 2 D extension of the wavelet algorithms is the 1 D algorithm applied to both the z and y directions of the 2 D signal Let us consider a 2 D signal as a rectangular matrix of signal values In the case where the 2 D signal is an image we call these signal values Pixel values corresponding to the intensity of the optical reflection Consider the input signal m n as an N x N square matrix We may process the signal along the z direction first That is we decompose the signal row wise for every row using the 1 D decomposition algorithm Because of the downsampling operation the two resultant matrices are rectangular and of size N x x These matrices are then transposed and pr
9. bases functions have infinite time duration oco lt t lt 00 it is awkward to use any of these bases to represent finite time domain transient signals Haar s work demonstrated that a continuous function can be approximated arbitrarily closely by an orthonormal basis with local support the Haar basis Since then mathematicians have constructed various new bases to represent a variety of analog functions During the past three decades signal processing has grown to be come a major discipline in engineering and computing science Since most finite energy signals either natural or otherwise are transient or nonstationary in nature it is most natural and effective to represent these signals by localized finite energy bases Because wavelets belong to this class of bases the development of wavelet theory and its applica tions to signal processing in various engineering disciplines has gained tremendous popularity in recent years Fourier analysis has been the major mathematical tool for signal representation and processing for at least the past 50 years The dis crete version of this approach called Fourier series breaks down a given signal into sinusoidal functions with different harmonics of the funda mental frequency Since sinusoidal signals are periodic signals Fourier analysis is an excellent tool for analyzing this class of signals How ever it is inefficient for representing transient signals Recall that the definiti
10. hard disk is not affected but if the file save button is invoked then the latest level 0 image will replace the original image if the same file name is chosen Coefficient processing In most signal processing applications DWT coefficients are altered deleted or clustered to achieve certain effects The 2 D tool allows the user to add new functions for processing of DWT coefficients in its buffers Processing of DWT coefficients is triggered by the button and a list of processing functions will pop up 3 3 2 D DWT TOOL 53 Processing Functions Level 2 lt I uj m LI ka Figure 3 10 Selection of DWT processing functions In the example shown in Fig 3 10 the user has chosen to apply a simple thresholding routine to eight subbands of the decomposed art pgm im age When the OK button is clicked on this routine will be applied to each of the subbands one after another until all the selected subbands are processed Post processed DWT coefficients will be displayed on the screen automatically see Fig 3 11 Figure 3 11 Post processed DWT coefficients In addition to its own internal routines the 2 D tool will accept user defined processing functions provided that they conform to the following dll interface structure More detail on the integration of 54 WAVELET TOOLWARE processing functions to the Toolware will also be discussed in the section on the 1 D tool define PROC2D Tell the Toolwa
11. in the spectral domain while the time domain resolution is coarse For processing of a signal at low frequency one should use wavelets with a large scale Although there are striking similarities between the STFT and the CWT their results in signal processing are quite different One must compute the STFT of a broadband signal on the t f plane for several sizes of the same window to capture the high and low frequency char acteristics of the signal However one only needs to compute the CWT once to detect both the high and low frequency events in the signal 1 6 MULTIRESOLUTION ANALYSIS 9 by aot b a t b a t Figure 1 4 Window widths for wavelets The original signal is uniquely recovered by the double integral 1 f f dbda z Wb O 018 Cy 00 J 00 a where Cy is a finite constant given by the integral 2 b w 1 19 Cy a 2 du lt 00 1 19 1 6 Multiresolution Analysis We have shown that the CWT has a unique advantage because its win dow widths can be controlled by the scale parameter a However we also see that the computation load of the CWT is quite heavy in order to capture all the characteristics of the signal To alleviate this compu tational burden mathematicians have developed the Discrete Wavelet Transform DWT to minimize the redundancies existing in CWT Al though the algorithm of DWT is identical to that of the two channel filter bank analysis the underly
12. reads in a 1 D signal and gener ates the STFT coefficients for display The wavelet generation tool is a simple utility to generate B spline wavelets and Daubechies orthonor mal wavelets The 1 D and 2 D DWT tools have a similar operating structure They assist the user to decompose a signal using scaled wavelets at different levels process DWT coefficients and then reconstruct the signal in a very flexible manner A decomposed signal can be recon structed from any level to any other lower level with only seleced levels affected In addition to some simple processing functions already in cluded in Toolware the user can also add new functions to Toolware through a defined dll program interface Toolware creates new files only at the explicit request of the user Thus accidental crash or termination of Toolware execution should not affect existing files Toolware runs on the Microsoft Win95 or NT 3 1 3 5 3 51 4 operating system For performance and stability reasons we recommend systems with a Pentium 133Mhz or equivalent processor and 32MB of main memory The memory size should proportionally increase for processing of large files to keep intermediate results 3 1 1 Installation of Toolware Toolware installation is very simple The user just follows these steps and the installation program does the rest 1 Insert the CD ROM into the CD drive E which is the most typical CD drive label If the label of the CD drive on your machi
13. saved into level 3 60 WAVELET TOOLWARE Tem A ret seal at GA TA 24 28 gaj DB1 2B Low lt 0 534 0 7 DB2 2B High lt 0 062 1 092 gt Figure 3 16 A magnified view of small DWT coefficients Processing Function Figure 3 17 The window for 1 D DWT processing functions 3 4 1 D DWT TOOL 61 econstiuction Functions Figure 3 18 Selection of the reconstruction structure buffers Level 2 buffers are updated in a similar fashion This process repeats until all levels between the chosen reconstruction levels are up dated As usual you can view the modified signal after reconstruction Custom processing routines Like the 2 D tool the 1 D tool also accepts custom routines for pro cessing wavelet coefficients The following is an example routine for smoothing of the chosen signal define PROC1D Tell the Toolware this is a 1 D routine define NAME Simple Smoothing Name of the routine include extfunci h Toolware symbols This routine accepts floating point numbers placed in the array input process the array and then place the results at the array output The total number of data samples is defined by len TOOLBOX_FUNCTION process float huge input float huge output int len int i for i 1 i lt len 1 i output liJ input i 1 input i input i 1 3 0 62 WAVELET TOOLWARE 3 5 Continuous Wavelet Transform Tool The continous wavelet transfor
14. the CWT coefficients through a self explanatory sequence After the signal is processed its CWT coefficients will be displayed on the window In this example displayed in Fig 3 20 the data set file delta pts represents two compounded sinusoidal signals in addition to two im pulses The frequencies of the two sinusoidal signals are 500 and 1000 Hz respectively The magnitudes of the two impulses are located at t 0 192 and 0 196 respectively and their magnitude is 6 The sam pling rate for this signal is set at 8 kHz You can see from the image that the CWT coefficients of the two sinusoidal signals are clearly dis played horizontally and those of the two impulses are displayed verti cally at their time instants 3 6 SHORT TIME FOURIER TRANSFORM TOOL 63 CWT Toolware delta pts EEE Zick AO 290066 nnnng5 910353 anil E e Se e DER A Figure 3 20 CWT coefficients 2048 3 6 Short Time Fourier Transform Tool Similar to the CWT tool the STFT tool generates the short time Fourier transform coefficients of a one dimensional signal for a given window size Except for the selection of the window size which deter mines the size of a signal segment for FFT computations the STFT tool has an identical set of control buttons except for selection of the short time window and they serve the same purposes as those in the CWT tool By clicking on the file button you can choose a one dimension
15. the Toolware it is graphically displayed on the screen in a three row format see Fig 3 13 and the data is ready for DWT decomposition processing and reconstruction The three rows have three different buffers denoted as DB DBo and DBs which will be also used for processing and reconstruction of DWT co efficients In this example both DB and DB contain an unprocessed signal and DB3 is empty Details on management of the three display buffers will be explained later 3 4 1 1 D DWT decomposition The 1 D tool supports both pyramidal decomposition and wavelet packet decomposition Similar to the decomposition procedure of the 2D DWT Tool you need to first specify a decomposition level and a wavelet type for decomposition of a one dimensional signal The 1 D DWT decomposition is activated by clicking on MEM A list of wavelets will pop up on the screen together with decomposition struc tures see Fig 3 14 At the user s choice one needs to select the wavelet type the decomposition level and the decomposition structure i e pyramidal or wavelet packet and then click on the OK button to start decomposition In this example we have chosen a level 0 to level 4 3 4 1 D DWT TOOL 57 Error DB3 DB1 lt 129 000 255 000 gt DB3 0A lt 0 000 0 000 gt Figure 3 13 The three row display of an acoustic signal wavelet packet decomposition of the signal by the cubic spline wavelet 3 4 2 Viewing and processing the data
16. time and frequency domains simultaneously engineers have used the Short Time Fourier Transform STFT to accomplish their goal even if only approximately The STFT is formally defined by the integral transform ST f t w a i fu De dr 1 2 where the overbar indicates complex conjugation The function w t is called the window function to be chosen by the user This is the reason that the STFT is also sometimes called the Windowed Fourier Transform WFT The variables w are the transform variables from the single time domain variable 7 Hence the STFT transforms a one variable function into a two variable function It provides time frequency information of the function f t on the time frequency t f plane The complex spectrum of ST f t w gives approximated spec tral properties of the function near the time location t It is only an approximation since the expression in 1 2 is the Fourier transform of the product of the functions f T and w 7 t In other words 1 2 gives the spectrum of the product function f 7r w 7 t instead of the function f 7 alone Even if the window function w t is a rectangu lar window with unit amplitude and f r is a pure sinusoidal function cos wot the transform of the product f 7 w z7 t is not just a delta function in the spectrum Rather it is the convolution of the delta spectrum with the sinc function spectrum of the rectangular function This broadening of the delta fun
17. value of zero at w 0 In engineering terms the wavelet has no d c offset This spectral domain behavior makes the wavelet a bandpass filter The reader can use the Toolware to generate the graph of any wavelet listed in the Part III by following the instructions given later in this manual It is easy to see from any wavelet and its spectrum that the energy of the wavelet is concentrated in a certain region of both the t and w axes This localization property is an important feature of the wavelets If a wavelet is more localized that is the energy of the wavelet is concentrated in a smaller region it produces a better representation of the signal in the time frequency or time scale plane In our case better means higher resolution and requires less coefficients in the 6 WAVELET TOOLWARE D2 scaling function scale 1 2 T T 4 a Wg see ma ue eee 2 a a a e ai D2 scaling function scale 1 2 DO a ee S D2 scaling function scale 1 4 2m 2 0 1 2 3 Figure 1 2 The scaling functions of the Daubechies 2 wavelet at three different scales representation 1 4 1 Wavelets at different resolutions scales For a given wavelet y t a scaled and translated version is designated a afte tral Fe 1 12 The parameter a corresponds to the scale while b is the translation parameter The wavelet po t w t is called the basic wavelet or mother wavelet The Daubechies D2 scal
18. 1 regularity 35 Reisz basis 6 resolution 3 RMS 3 sample discretized values 17 scale 10 scaling 6 scaling function 6 10 16 18 20 34 35 scaling function component 24 scaling function representation 17 scaling function series 18 24 72 scaling functions 25 scaling parameter 10 semi orthogonal 19 20 Short Time Fourier Transform 2 signal compression 30 signal distortion 22 signal recognition 39 similarities 9 sinusoidal function 1 9 small wave 5 soft thresholding 39 speckle noise 41 spectral domain 3 5 9 27 spectral energy 3 spectral resolution 23 spectral domain windowing 3 spectrum 6 spline space 21 splitting trick 23 STFT 9 31 subspace relation 15 symmetric 22 tensor product 40 tensor product 25 thresholding 30 time domain 2 time frequency information 2 time dependent spectrum 5 time domain 9 time domain windowing 3 time frequency 6 Toolware 5 transient signals 1 5 translate 6 20 INDEX translate scaling function 36 truncation 3 Two dimensional wavelet pack ets 27 two dimensional wavelets 25 two scale relation 13 20 two scale relations 15 16 24 two scale sequences 15 22 two variable function 2 uncertainty principle 3 unique transformation 3 universal thresholding 39 up sampling 17 vector space 19 wavelet 5 6 11 16 wavelet algorithm 24 wavelet basis 6 wavelet coefficient
19. 3 4 1 D DWT Tool After the 1 D DWT tool is activated it does not display any signal before any data file is opened The 1 D DWT tool has the same set of control buttons as the 2 D DWT tool but their display and processing mechanisms are quite different The list of files that can be read into the 1 D DWT tool by clicking on the BEES button is displayed on the screen see Fig 3 12 Data file format The 1 D DWT tool accepts numerical integer or floating point files in the ASCII format Appending pts to a file name i e file name pts makes it visible to the 1 D DWT tool and the file name will appear on the file list of Toolware The folowing are three examples that have an acceptable data format for the Toolware 000000 000111 000444 000998 001773 002765 003973 005395 007027 008866 oooo co 0 00 0 0 56 WAVELET TOOLWARE 0 010906 0 013145 0 015577 0 000000 0 000111 0 000444 0 000998 0 001773 0 002765 0 003973 0 005395 0 007027 0 008866 0 010906 0 013145 0 015577 37 40 204 80 88 163 186 131 112 157 129 120 82 69 116 74 43 73 114 73 90 123 147 149 165 193 185 168 180 203 172 165 181 182 183 177 191 210 197 219 224 148 96 255 100 1 144 88 20 7 55 45 27 74 43 58 193 74 93 179 185 113 114 172 114 10 0 87 70 100 64 38 75 109 63 85 139 135 138 161 172 188 186 173 185 206 170 163 1 93 182 171 178 196 208 194 214 223 191 74 226 181 0 123 109 33 5 42 56 26 49 After a data file is loaded into
20. 635881557 1 115087052458 0 533728236886 0 266864118443 0 557543526229 0 591271763114 1 205898036472 0 078223266529 0 295635881557 0 057543526228 0 533728236886 0 016864118443 0 028771763114 0 091271763114 0 156446533058 0 026748757411 0 045635881557 0 033728236886 0 053497514822 1 8 Wavelet Packets The hierarchical wavelet decomposition produces signal components whose spectra reside in consecutive octave bands In certain applica tions the spectral resolution produced by this decomposition may not be fine enough One may want to use the CWT to obtain the nec essary finer resolutions by changing the scale parameter a by a small increment However this increases the compuation load by orders of magnitude In order to avoid this problem wavelet packets are general izations of wavelets in that each octave frequency band is itself further subdivided into finer frequency bands by wavelet packet transforms In other words the development of wavelet packets is a refinement of wavelets in the freqency domain and is based on a mathematical the orem proven by I Daubechies the splitting trick The theorem is stated as follows If f k kez forms an orthonormal basis and Fi z J mfz k 1 71 F z dlaf c k 1 72 then F 2k Fo 2k k Z is an orthonormal basis for E 1 8 WAVELET PACKETS 23 span f n n e Z The proof of this theorm is given in her book Ten Lecture
21. For any scaling function with its cor responding wavelet 4 there are three different 2 D wavelets and one 2 D scaling function using the tensor product approach We can write 1 9 TWO DIMENSIONAL WAVELETS 25 Bree 00000000 SHooee o Figure 1 10 Block diagram for the WP decomposition of a signal 26 WAVELET TOOLWARE the 2 D wavelets as W y glz ily j 1 78 wPl e y ple ijely j 1 79 TEle y ple ily j 1 80 and the 2 D scaling function as Slzy O a i b y J 1 81 In the spectral domain each of the wavelets and the scaling function occupies a different portion of the 2 D spectral plane When we analyze a 2 D signal the decomposition algorithm generates four components from the input signal With respect to the spectral domain the compo nents are labeled low high LH high low HL and high high HH corresponding to the wavelets vP y M 1 2 3 The component that corresponds to the scaling function signal is called the low low LL component The terms low and high refer to whether the pro cessing filter is lowpass or highpass The decomposition of a 2 D signal results in what is known as a hierarchacal pyramid Because of the downsampling operation each image is decomposed into four subim ages The size of each subimage is only a quarter of that of the original image Readers can refer to the last section of this manual to try to decom
22. The 1 D DWT tool can decompose a signal up to four levels As in the 2 D tool it also uses different buffers to store DWT coefficients For example if a signal is decomposed for three levels its DWT coefficients in any of the three levels can be retrieved for processing reconstruc tion and display When a particular subband at a level is chosen by clicking on one of the following level buttons the DWT coeffcients will be loaded into the three display buffers DB DB2 and DB for Figure 3 14 Selection of decomposition parameters 58 WAVELET TOOLWARE visualization based on the following rules When the user clicks on 0A the original signal is displayed on DB and the reconstructed signal if any is displayed on DB3 The recon struction difference between DB and DB is shown at DB If the system does not have any reconstructed data yet then DB and thus the third display row is empty and both DB and DB display the original data because DB keeps the difference between DB and DB3 For other decomposition levels DB displays the parent subband of the two chosen subbands displayed in DB and DB That is at level 1A which represents the two subbands decomposed from the original signal at 0A its low pass and bandpass subbands are displayed in DB and DB respectively The parent of 1A which corresponds to the original signal is displayed on DBs For following levels letters A and B respectively denote the l
23. Wavelet toolware Wavelet Toolware Software for Wavelet Training This page intentionally left blank Wavelet Toolware Software for Wavelet Training Andrew K Chan Department of Electrical Engineering Texas A amp M University Steve J Liu Department of Computer Science Texas A amp M University San Diego London Boston New York Sydney Tokyo Toronto This book is printed on acid free paper Copyright 1998 by Academic Press All rights reserved No part of this publication may be reproduced or transmitted in any form or by any means electronic or mechanical including photocopy recording or any information storage and retrieval system without permission in writing from the publisher ACADEMIC PRESS 525 B Street Suite 1900 San Diego CA 92101 4495 USA 1300 Boylston Street Chestnut Hill MA 02167 USA http www apnet com ACADEMIC PRESS LIMITED 24 28 Oval Road London NW1 7DX UK http www hbuk co uk ap Printed in the United States of America 98 99 00 01 02 IP 987654321 Disclaimer This eBook does not include the ancillary media that was packaged with the original printed version of the book Contents Part I 1 1 1 2 1 3 1 4 1 5 1 6 1 7 1 8 1 9 Part II 2 1 2 2 2 3 2 4 ii From Fourier Analysis to Wavelet Analysis 1 Short Time Fourier Transform 2 Window Measures 02 2 0004 ee 3 What Is a Wavelet e 0 4 4 4 4 4 8 da We ea
24. al data file in either the pts format or the standard wav format to be processed by the STFT tool Similar to the CWT tool the STFT coefficients of the signal are organized into a two dimensional matrix and then graphically displayed where the z axis of the display is the time line and the y axis is the frequency axis of the signal It is interesting to compare the results obtained from the CWT to those from the STFT The following two figures represent the results when the window sizes are set to be w 8 Fig 3 21 and w 64 Fig 3 22 respectively As expected for a given window size the STFT has either good time resolution or good frequency resolution but not both Using the CWT we can easily determine that the signal delta pts has two dominating frequencies and two impulses from the CWT coefficients On the other hand one must use different window sizes to test this signal in order to reach the same or a similar conclusion 64 STFT Toolware delta pts aang ACI o 5 tt i ii Tt i 2048 Figure 3 21 STFT coefficients for window size 8 3 7 Wavelet Builder Tool The wavelet builder tool constructs two different types of wavelets and their scaling functions In addition to their graphical display the nu merical values of those functions can also be exported to data files The B spline functions are generated using their closed form expres sions By choosing the B spline family from the file selection button
25. ble Wavelet subspaces Let the subspace V constitute an orthogonal sum for which we use the symbol of mutually orthogonal subspaces W of L as Va O72 Wj 1 22 j 00 Then Va 1 Bj coW iW D Wn Va D Wn 1 23 The subspace W the wavelet subspace is said to be the orthogo nal complementary subspace of Vp in Vasa This relationship is best 1 6 MULTIRESOLUTION ANALYSIS 11 Va Wa Figure 1 5 Approximation subspaces and wavelet subspaces described by Fig 1 5 We should take note of the following two obser vations 1 Subspaces V are nested 2 Subspaces W are mutually orthogonal Mathematically they are expressed by ViNVe V l gt j 1 24 W NW 0 7 2 1 25 ViNW 0 j lt 1 26 The subspaces W are generated by some function 7 t called a wavelet in the same way as V are generated by t In other words any function f t V can be written as 00 i oO fit gt C5 p0 2 t k 5 Cj Pikes 1 27 k 00 k and any function g t W can be written as g t d pp 2t k 5D di Pik 1 28 k 00 k 00 12 WAVELET TOOLWARE for some coefficient sets fc uz and d e N in Here we have adopted the notation Prj t k deg p t k We may restate here that a function f t can be decomposed into many components such that M m l where the function fm t is the approximation of f t from the
26. ction 6 w wo to a sinc w wo func tion represents the effect of the truncation due to the window function We call that a windowing effect on the function f r We will see in 1 3 WINDOW MEASURES 3 later sections that the accuracy of the information on the t f plane is controlled by the widths of the window in time and frequency domains In addition the uncertainty principle governing the area of the window in the t f plane places a limit on the resolution that one can achieve using STFT The spectral domain equivalent of 1 2 comes from the Paseval iden tity which gives ST flt w l fina He dr z erit E FOET 13 Expression 1 3 indicates that the time domain windowing process is also a spectral domain windowing process It provides information on the spectral energy of the signal near the center of the window The original signal f t can be recovered from the STFT uniquely Hence the STFT is a unique transformation The recovery formula for the STFT is given by f id w E ST 7 w w r the dwdr 1 4 where w is the norm of the function that will be defined in the next section 1 3 Window Measures In order to make an intelligent choice of a window function for the STFT some measures of goodness must be established to evaluate the suitability of the window function Engineering measurements such as the main lobe zero crossing and the 3db points do not gauge the concentration of energy of t
27. display on a small floating window see Fig 3 16 The zoom in window disappears at the release of the button 1D DWT Toolware hu pts Bl SARRE BE ex 22 za aa ec so aja acl ap ar ac asa 1250 DB1 3D Low lt 1 644 557 740 gt 55 ce DB2 3D High lt 1 447 778 035 gt DB3 2B High lt 20 161 669 170 gt Figure 3 15 The three row display of DWT coefficients at level 3 For a given display level the user can process data stored in the DB to DBs by clicking on the P processing button MEEI and then the list of processing functions will pop up see Fig 3 17 By clicking on a processing function together with one or more buffer buttons up permost signal DB middle signal DB bottom signal DB3 the user can process data stored in selected buffers by the chosen function Both the internal and display buffers are refreshed after the processing function is executed At any time the user can reconstruct all or some DWT coefficients Similar to the 2D Tool once a signal is decomposed its decomposition structure is registered into Toolware and the user can only use the se lected decomposition wavelet and structure for reconstruction In the example shown in Fig 3 18 only the reconstruction levels can be chosen by the user and all other parameters will be set by Toolware Once the reconstruction command is executed new level 3 DWT coefficients are reconstructed from level 4 DWT coefficients and
28. e ad 5 1 4 1 Wavelets at different resolutions scales 6 Continuous Wavelet Transform 7 Multiresolution Analysis 205 9 1 6 1 Two scale relations 12 1 6 2 Decomposition relation 04 13 1 6 3 Development of the decomposition and reconstruc tion algorithms 4 6044 oN ae wee we 15 1 6 4 Mapping of functions into the approximation space 16 Types of Wavelets 0 02 0 02 0000 0 08 18 1 7 1 Orthonormal wavelet bases 18 1 7 2 Semi orthogonal wavelet bases 19 1 7 3 Biorthogonal wavelet bases 21 Wavelet Packets n e e Se Sr e ate p 2 va cotit a a ae d 22 1 8 1 Wavelet packet algorithms 23 Two Dimensional Wavelets 004 24 1 9 1 Two dimensional wavelet packets 26 27 Wavelet Algorithms Overview 29 Algorithm for Computing the STFT 29 Algorithm for Computing the CWT 31 Computation and Display of Scaling Functions 32 v vi WAVELET TOOLWARE 2 5 Computation and Graphical Display of Wavelets 34 2 6 Mapping between B Splines and Their Duals 35 2 7 Decomposition ofa 1 D Signal 36 2 8 Reconstruction of 1 D Signals 37 2 9 Thresholding lt p Ho ce ate ok tai gi eve Sh Se Be ke SSS 37 2 9 1 Hard thresholding aaa 6 amp wae oe eS 38 2 9 2 Soft thresholding 4 sw lak a 3 9 Gals ose Be 38 2 9 3 Quantile th
29. ea of the STFT in the t f plane is given by u t t Aw Aw x w w Ad wtw Ae 1 11 The center of the window is located at t t w w on the t f plane where t is the time parameter being considered Once we have chosen the window function w t the window widths 2Aw and 2A are fixed everywhere in the t f plane since they are independent of t and w This effect is demonstrated on the t f plane as given in Fig 1 1 This characteristic presents a distinct disadvantage of the STFT since it creates difficulty in accurately processing signals with a time dependent spectrum such as a chirp signal In order to achieve a high degree of accuracy the STFT must be repeatedly applied to the signal with a varying window width each time 1 4 WHAT IS A WAVELET 5 Figure 1 1 STFT windows 1 4 What Is a Wavelet Wavelets are finite energy functions with localization properties that can be used very efficiently to represent transient signals Efficiency here means only a small finite number of coefficients are needed to rep resent a complicated signal In contrast with the sinusoidal functions of infinite extent big waves wavelet implies a a small wave Strictly speaking it also means mathematically that the area under the graph of the wavelet y t is zero i e f y t dt 0 In the spectral do main this property is equivalent to 0 0 i e the spectrum of the wavelet has a
30. efficients are usually very sparse For much better visual effects we display wavelet coefficients except for the coefficients in the lowest subband in reverse brightness That is except for DWT co efficients in the lowest subband i e LL LLLL LLLLLL etc larger 3 3 2 D DWT TOOL 51 DWT coefficients are represented by darker display pixels For conve nience of operation a set of buttons LEE a 2 is available for quick preview of thresholding effects without changing their internal values The two leftmost buttons decrease increase the smallest absolute values to be displayed The next two buttons de crease increase the largest absolute values to be displayed Next to them the ruler like icon indicates the selected lower upper bounds with respect to the range of true values Wavelet coefficients that fall into the lower upper bound setting are scaled to the display range of 0 to 255 The far right button resets all the threshold settings and returns the image display back to its original form We remark here that un like the processing functions the display thresholding does not have permanent effects on the internal values of the displayed DWT coeffi cients The user should use built in or custom processing functions to alter DWT coefficients for actual signal processing purposes Fig 3 8 shows an example in which some DWT coefficients are thresholded in the display buffer so that they become less visible on the screen
31. el Repeating the same procedure with coefficient set b yields the wavelet coeffi cients These procedures are repeated to yield the coefficients at lower resolution levels 2 8 Reconstruction of 1 D Signals The two scale relations for the approximation space and the wavelet space constitute the reconstruction algorithm The scaling function co efficients at a higher resolution level are computed by using the formula Che So De 2kCj k qe 2kd k k Each summation can be interpreted as a convolution process after up sampling This process is depicted in Fig 1 8 It can be repeated for coefficient sequences cj_px and dj pk p M M 1 0 The reader should use the Toolware to practice this algorithm for 1 D sig nals 2 9 Thresholding Thresholding is one of the most often used processing tools in wavelet signal processing It is used in noise reduction in signal and image compression and sometimes in signal recognition The four types of thresholding we use are l hard thresholding 2 soft thresholding 3 quantile thresholding and 4 universal thresholding The choice of thresholding methods depends on the application We discuss each type here briefly 38 WAVELET TOOLWARE 2 9 1 Hard thresholding Hard thresholding sometimes is called gating If a signal or a coeffi cient value is below a preset value it is set to zero That is y z forz gt o 2 96 O0 forrt lt o where is the threshold value
32. general learning and experimentation Typ ical numerical errors in the decomposition and reconstruction of a 512 x 512 gray scale image are listed in the following table The signal error ratio is defined as the absolute ratio of the reconstruction errors divided by a signal expressed in decibels Our test results show that the difference in numerical errors on different machines is very small 3 2 USING TOOLWARE 45 Daub7 2 37440e 02 dB Daub8 4 24760e 19 2 31849e 02 dB Daub9 7 74388e 20 2 39241e 02 dB Daubl0 8 43506e 19 2 28870e 02 dB BiDaub1_1 oo dB BiDaub1_3 oo dB BiDaub2_1 oo dB BiDaub2_2 0 00000e 00 oo dB Biaub2 3 0 00000e 00 co dB BiDaub2_4 8 46053e 13 1 68857e 02 dB BiDaub3_l oo dB BiDaub3_2 0 00000e 00 co dB BiDaub3_3 0 00000e 00 co dB BiDaub3_4 5 09022e 01 dB BiDaub3_5 1 68857e 02 dB BiDaub4_1 2 41220e 02 dB BiDaub4_2 2 57555e 02 dB BiDaub5 1 4 92118e 20 Linear BSpline 7 94940e 08 1 19127e 02 dB Cubic BSpline 3 08383e 01 5 32399e 01 dB Meyerzero 6 31288e 02 6 01285e 01 dB Meyerfirst 4 62319e 02 6 14814e 01 dB 3 2 Using Toolware Toolware is activated by double clicking on the Toolware group icon Fig 3 2 which is placed at the designed directory during installation to bring out the main menu Fig 3 3 After the main menu appears on the screen the user can activate any of the five tools simultaneously The five tools include one di
33. he function In the wavelet literature the RMS width of the function is used to measure energy concentration in the time and frequency domains These measurements are explicitly expressed respectively as Aw dal p t PA a 1 5 t Tal AN jw t P ae 1 6 4 WAVELET TOOLWARE and i ees ial fw apte ate da 1 7 gee ie ECOLA 1 8 Here the usual definition of the energy norms T wta 1 9 iol f 0a 1 10 applies to these formulas The measurement Aw is considered as half of the time domain window width which measures the concentration of energy in the window A smaller window width means more energy is concentrated in a smaller area One may use a small time window width to deduce the fine structure of a signal On the other hand Ai is the measurement of energy concentration in the spectral domain In addition t and w are the centers of gravity of the window in the time and spectral domains respectively If w t is real and even and the average value of the function is zero then w is also symmetric with respect to the origin and 0 0 In this case the integration in the frequency domain expression extends only from 0 to oo and by the definition of 1 8 w is located on the w axis other than the origin In general one must compute the window widths numerically since most of the window functions we are considering in this manual may not be described analytically We can easily see that the window ar
34. infinite sequences with exponential decay Explicit expressions for these sequences are II ak m k fim 1 yom SS e SO Nam 2h Daw e 0 2h k l m41 k 1 i m b Ma 2 D 2k 2m E T 2h m 1 7 TYPES OF WAVELETS 21 where the z transform of gn is the reciprocal of the z transform of the Euler Frobenius Laurent Polynomial Ey z namely 1 G z gt ghz F EG et Nom m h zh h m l Numerical values for these sequences are listed in Table 5 5 of the book Wavelets A Mathematical Tool for Signal Analysis SIAM Publ 1997 for m 2 linear splines and m 4 cubic splines 1 7 3 Biorthogonal wavelet bases It is known in the signal processing literature that finite impulse re sponse FIR filters with symmetric or antisymmetric coefficients have linear phase or psuedo linear phase characteristics Linear phase prop erties are important for minimizing signal distortion during processing The processing sequences or the FIR filters p and qx of finitely supported orthogonal wavelets are generally not symmetric In order to take advantage of the symmetric property one has to give up or thogonality Filter designs in a two channel filter banks are examples of a biorthogonal basis In this case orthogonality between filters is not required The only requirement for filter bank design is that the output of the filter bank be a delayed version of the input That is the output signal must be a pe
35. ing function at three different scales is shown in Fig 1 2 It is important to note that the shape of the wavelet remains the same under translation and scaling We will show below that with the proper choice of b a the set of wavelets Yra t constitutes the basis functions of a Reisz or stable basis in the Z or finite energy space Keeping this concept of basis in mind we will show that the idea of 1 5 CONTINUOUS WAVELET TRANSFORM 7 wavelet signal processing is not much different from that of Fourier pro cessing at least when orthonormal wavelet are considered Instead of decomposing a signal into sinusoidal functions of different frequen cies the Fourier basis wavelet signal processing seeks to decompose a transient signal into a linear combination of a scaled and translated version of the basic wavelet the wavelet basis 1 5 Continuous Wavelet Transform The continuous integral wavelet transform CWT IWT of a signal z t is a linear transform defined by the integral f i a t Wo t dt Jefe x t Pa t 1 13 The last expression of 1 13 represents the inner product of two functions defined in the L space by 6 9 FOD 1 14 Wyx 0 a Referring to 1 13 the CWT computes via the inner product formula the wavelet coefficient of x t associated with the wavelet Ysa t This coefficient indicates the correlation between the function x t to the wavelet w t Higher correlation produce
36. ing meanings of these algorithms are different We will point out these differences in later sections The most important feature of multiresolution analysis MRA MRA is the ability to separate a signal into many components at different 10 WAVELET TOOLWARE scales or resolutions For a specific choice of the scaling parameter such as a 2 j Z the decomposition algorithm is equivalent to putting signal components into successive frequency octaves Similar to multiband signal decomposition the goal here is to apply the divide and conquer strategy on the signal so that individual components may be processed by different algorithms We present here the essence of the MRA by considering the properties of the approximation subspaces and the wavelet subspaces Approximation subspaces A function t 12 called a scaling function or an approximate function generates a nested sequence of subspaces approximation subspaces V of L such that 10 C Vea CVichcUuchc 9L 1 20 and satisfies a dilation refinement equation p t X polat k 1 21 k for some a gt 0 a 1 and some coefficient sequence p 2 For the DWT to be discussed in some detail later we will consider a 2 that corresponds to a dyadic octave scale in the spectral domain More precisely Vo is generated by k k e Z In general V is generated by 2 k k Z The symbol stands for the dummy varia
37. l processing As mentioned earlier the idea consists of simply dividing the signal into components so that each component may be processed with a differ ent algorithm The important issue of this algorithm is to be able to reconstruct the signal perfectly if we apply all pass filters on all signal components These algorithms are based on the two scale relations and the decomposition relation as developed in the previous sections We rewrite these relations here for convenience Let filt Viu and fjalt LCi kj i klt fit e Vj and f t Dental g t W and g t dana t 1 43 Since the MRA requires that View Vi W 1 44 we have frit fi t glt 1 45 gt CpiwPjt1e t 2 Ci Pj a t gt dj Wy a t 1 46 We substitute the decomposition relation 1 40 in 1 46 and inter change the order of summations Comparing the coefficients of t and p g t we obtain Cie J aon ecjsie 1 47 dp gt dox ecj410 1 48 z where the right sides of 1 47 and 1 48 correspond to a down sample after convolution These formulas relate the coefficients of the scaling function from a given scale to the scaling function and wavelet coeffi cients at twice the given scale 16 WAVELET TOOLWARE EO OL 3 ds n Gk Figure 1 8 Decomposition and reconstruction The reconstruction algorithm is based on the two scale relations When we sum the functions at the jth resolution we have F
38. latter is used in filter bank processing digital processing in the frequency domain 1 7 Types of Wavelets In the preceding discussions of the approximate subspace and wavelet subspace the basis functions that span the subspace can have orthonor mal semi orthogonal or biorthogonal bases There are conditions to be satisfied by each type of basis We will state these conditions separately as follows 1 7 1 Orthonormal wavelet bases Orthogonality of functions depends on the definition of the inner prod uct in that vector space We have used the inner product defined in the L space Two different functions are orthogonal to each other if the inner product of these functions is zero That is lt f es e PROTOL 0 1 61 If we take g x f x and oF e Boa si 1 62 we say that the function f x is normalized Hence if the set of basis functions j Z spans an approximation space V and satisfies the condition lt Bibi gt gt de DAE Ido iy 163 it is a normalized orthogonal set or for short an orthonormal set This definition applies to the wavelet subspace as well A wavelet is said to be orthonormal if the set of basis functions 4 24 y 2 r k satisfies lt Vie Vik gt a Vielr Win z dx di je 1 64 1 7 TYPES OF WAVELETS 19 Since the wavelet subspace and the approximation subspace are orthog onal we have lt Qoi 0 3 gt ie boi Wo j x dx 0 y i j EZ 1 65 No
39. m tool CWT tool generates the con tinuous wavelet coefficients of a one dimensional signal It has eight different control buttons see Fig 3 19 Figure 3 19 CWT tool control buttons Starting from left to right the first button is for selection of a data file for processing the second and third are for saving the displayed CWT coefficients into an ASCII data file or a pgm image file respec tively The image of the CWT coefficient matrix can be printed to your printer by clicking on the printer button The relative display intensity of CWT coefficients can be adjusted by the next four buttons Clicking on the fifth button which has a straight line means that the display intensity directly reflects the magnitude of CWT coefficients The next two buttons are for de emphasis of small and large coefficients The eighth button allows you to display either only the magnitude which is the default mode or signed amplitudes with different colors assigned to positive and negative values The CWT coefficients of the signal which are organized into a two dimensional matrix are graphically shown on the display window of the CWT tool The z axis of the display is the time line and the y axis is the scale of the signal After you invoke the CWT tool from the main menu it will display an empty window By clicking on the file selection button you can select a pts or wav file as the signal input Then the CWT tool will produce
40. mensional discrete wavelet transform and processing BINE two dimensional discrete wavelet transform and processing 2w f continuous wavelet transform E short time Fourier transform STFT a and the generation of 46 WAVELET TOOLWARE amp D to olvware 98012 24 E a ea A Extfune h Extfuncl h flower2 pgm fountanpgm hydrant pom m MI o 9 hi koldus pgm lake pgm leaves pgm lupt manhole pgm mo ao n 3 i mecbowl2 pgm mscwall pgm o amp m pgm oul pts output data e i li le fa patk2 pgm resistor pgm ruddfoun pgm s1 pts s2pts Figure 3 2 The start up icon for Toolware programs Toolware Figure 3 3 Toolware main menu wavelets N These tools operate independently and they use different data formats 3 3 2 D DWT Tool After you open the 2 D DWT tool you need to choose a data file before proceeding with other processing functions This is done by clicking on the E button 3 3 2 D DWT TOOL 47 3 3 1 Data file selection Toolware accepts pgm image files meaning that any files appended with pgm can be read into Toolware for processing After the user clicks on the File button or its corresponding graphical icon the file selection window Fig 3 4 will pop out for the user to open an image file for processing From the file window one can identify any image file in PGM format in any directory and then load the file by clicking on the mE button image Format
41. nageable otherwise the data storage required is substantial The following alternative approach can be used 34 WAVELET TOOLWARE Iterative method We modify the two scale equation to fit an iterative technique by con sidering nyila SO pebn 2e k n 0 1 2 2 85 k The index n indicates the number of iterative loops For initial ization the user may use either first or second order B splines The iteration should converge after 5 to 10 iterations if the regularity of the function is high However for some sequences designed by the per fect reconstruction filter bank approach the algorithm may not even converge at all In other words there are filter banks that are not as sociated with scaling functions and wavelets The final data may be seen by using the usual 1 D graphic display The code to generate the scaling function in this Toolware is based on this algorithm We outline the procedure as follows 1 Initialize the program by setting all data files to zero 2 Set desired number of iteration say 10 and set the iteration index to 1 3 Input the initial trial function o x One may use an impulse function a rectangular pulse i e first order B spline a triangu lar pulse i e second order B spline etc 4 Carry out 3 13 by convolving the function with the pz sequence 5 Upsample the resulting sequence by inserting zeros in between every other data point This sequence is z 6 Increa
42. ne is not E use the proper CD drive label in subsequent steps 2 Click on the run icon on the desktop 3 Type in F setup command and the installation will start 44 WAVELET TOOLWARE 4 Interact with the installation software to load Toolware program files and data files into the directory that you want to host Tool ware 3 1 2 Wavelets Toolware supports several different wavelets including Daubechies or thonormal wavelets Daubl to Daub10 linear B spline wavelet and cubic B spline wavelet Battle Lemarie spline wavelets zeroth and first order Meyer wavelets and biorthogonal Daubechies wavelets 1 1 o 1 3 2 1 to 2 4 3 1 to 3 5 4 1 4 2 and 5 1 As shown in Fig 3 1 an image is decomposed using the tensor product of the 1 D transform of an image along its vertical and horizontal directions respectively Boundary data for both 1 D and 2 D data are treated by a simple wrap around or data reflection method to achieve perfect or nearly perfect reconstruction Row Filters gt k Column Filters gt LL Subband UC LH Subband eRe Original Image Lowpass Down or Previous Filter Sampling LL Subband ar Rua CORC E Bandpass Down CH 240 HH Subband Filter Sampling S Figure 3 1 The one lelvel DWT decomposition structure of a 2 D image in Toolware Toolware is not designed for solving high precision computational problems but is meant for
43. nt levels are buffered in Toolware and they are dynamically packed together to form one display image as needed Level 0 represents the original data or the image reconstructed from its DWT coefficients Obviously the reconstructed image can be different from the original one if certain DWT coefficients are manip ulated by processing functions in the tool The user needs to take precaution in preserving the original image while saving new results into the hard disk 3 3 2 Wavelet decomposition After an image is loaded into the 2 D tool it can be decomposed by clicking on the MESE button Then a list of wavelet decomposition se quences will pop up along with a choice of the decomposition structure Fig 3 6 shows an example in which a five level wavelet packet struc ture is selected for decomposition of an image using the Daubechies biorthogonal wavelet The decomposition process is started and termi nated automatically by clicking on the OK button Decomposition Functions Figure 3 6 Selection of the decomposition parameters After Toolware has performed the decomposition based on selected parameters the decomposed image will pop up on the window so that the user can view the absolute value of the wavelet coefficients see Fig 3 7 The displayed image is generated by using the truncated absolute 50 WAVELET TOOLWARE values of wavelet coefficients but raw DWT coefficients are retained in an internal buffer for actual proce
44. o that the reader can easily follow them to build codes using his her preferred computer languages The following algorithms are included in this sec tion Computation of STFT for 1 D signal Computation of CWT for 1 D signal Algorithm for generation of the graph of scaling functions Algorithm for generation of the graph of wavelets Mapping between B splines and their duals Decomposition of 1 D signals Reconstruction of 1 D signals Thresholding O wan mo A A WO N e Two dimensional extension of the wavelet algorithms O Other algorithms 2 2 Algorithm for Computing the STFT Let us rewrite the definition of STFT here for easy reference CO X STe f t w I f r w r tje dr 00 A straightforward approach to computing the STFT is to make use of the Fast Fourier Transform FFT algorithm that can be found in many numerical analysis books The procedure for computing the STFT is given here 1 Choose a symmetric window function and sample it at the same rate as the input data 2 Count the total number of samples of the window function say N 30 WAVELET TOOLWARE 3 To take care of the edge or boundary effect of the algorithm we may reflect if N is even or 231 if N is odd data points to the outside the edges of the 1 D data set at both ends of the set That is if the data set has M data points the modified data set will have a total of M N points from toM 4
45. ocessed row wise again to obtain four x x x square matrices namely cl m n di m n di m n and di 1 m n The subscripts of the d matrices correspond to the three different wavelets This procedure can be repeated an arbitrary number of times to the c m n matrix or the LL component and the total number of coefficients after the decomposition is always equal to the initial input coefficient N The reader may use the Toolware to practice the 2D decomposition and reconstruction algorithms If the coefficients are not processed the original data can be recov ered exactly through the reconstruction algorithm The procedure is simply the reverse of the decomposition except that the sequences are 40 WAVELET TOOLWARE Dx x instead of a bx Care should be taken to remember upsam pling first before convolution with the processing sequences 2 11 Other Algorithms Toolware also includes following codes for wavelet processing e Noise removal using hard thresholding e Image compression using thresholding e 1 D wavelet packet decomposition and reconstruction e 2 D wavelet packet decomposition and reconstruction Implementation of noise removal is straightforward and is not elab orated here The user can make use of the Toolware to see the effect of speckle noise being removed from an image using hard thresholding Setting the wavelet or wavelet packet coefficients at small ampli tudes to zero by thresholding compresses
46. olware contains among other algorithms and codes the one dimensional 1 D and two dimensional 2 D Discrete Wavelet Transforms DWTs and their inverse transforms as well as computations of the Continuous Wavelet Transform CWT and the Short Time Fourier Transform STFT Simple signal processing appli cations are also included for reader to practice with the software This manual for Wavelet Toolware leads the reader through the ba sics of wavelet theory and shows how it can be applied to a number of engineering problems It is organized into three major divisions Part I concerns the basics of wavelet theory and readers should get a very fundamental understanding of the development of the theory and gain appreciation of the time frequency or time scale representation of sig nals Part II contains several major algorithms that correspond to the topics discussed in Part I In most cases the algorithms are outlined in procedural steps Finally Part III provides hands on practice with the Wavelet Toolware It contains many practice sections ranging from installation of the software to denoising of images This manual is designed to be self contained for the reader who is unfamiliar with wavelet theory but wants to gain a basic understanding and the usage of wavelet theory After thoroughly studying this man ual and practicing with the Toolware readers should be able to apply the basic aspects of wavelet analysis to problems of their respec
47. on of the sequence is necessary with a reasonable filter length to ensure adequate accuracy For processing of large images e g mammograms at 4000 x 5000 x 12 bits per image or 3 D images e g multislice CAT scans the decomposition process is inefficient and the compuational load is large Since the B splines N 2 t k and its duals N 2 t k span the same subspace V it is convenient and efficient to map the B spline coefficients cj to the dual spline coefficients cpu After that we may process the signal using sequences p and g which are finite with relatively short filter lengths Let s t 5D Crab 2 87 E fibralt 2 88 k be a signal approximated by a B spline series t N 2t k and the dual B spline series Ge t Nm 2t k respectively cx are the spline coefficients that have been determined through the initial 36 WAVELET TOOLWARE mapping from the sampled values To obtain the dual B spline coef ficients one simply makes use of the orthogonality condition between the B spline and its dual Nalt Nm 25t k dx 2 89 Let j 0 by making the sampling interval unity We compute the dual coefficients cz by a Stati mt 9 2 90 k 2 CkNm t k Nin t o 2 91 k Since B splines have compact support the autocorrelation sequence 00 ip J Net k Nin t Oat 2 92 00 is a finite sequence where Tk Nom m k k m 1 m4 1 so
48. on of the Fourier transform of an analog function f t is defined by fw f j fedt 1 1 It is easy to show that if f t is replaced by a delta function t to at t to then the Fourier transform or spectrum becomes a sinusoidal function e 7 t cos tow j sin tow Hence it takes an infinite number of frequencies to represent a signal that exists at only one point in the time domain On the other hand if the function to be analyzed is a sinusoidal function of a single frequency the spectrum is a delta function A careful study of this consideration leads to discrete Fourier analysis Discrete Fourier analysis is very efficient for studying global periodic function However for transient signals such as music image speech acoustic noise seismic signals thunder and lightning it is more efficient to represent these signals by localized finite energy bases 2 WAVELET TOOLWARE Other drawbacks of the Fourier analysis include the following 1 The Fourier transform can be computed for only one frequency at a time 2 Exact representations cannot be computed in real time 3 The Fourier transform provides information only in the frequency domain but none in the time domain 1 2 Short Time Fourier Transform Fourier analysis being unable to provide localized time and frequency information simultaneously becomes the most serious drawback for pro cessing of transient signals In order to gain information from both
49. or the gate value 2 9 2 Soft thresholding Soft thresholding is defined as y f x o forr gt a O0 org 2 97 The function f x generally is a linear function a straight line with slope to be chosen However spline curves of third or fourth orders may be used to effectively weight the value greater than In some applications such as signal compression using a quadratic spline curve may increase the compression ratio by a certain amount 2 9 3 Quantile thresholding In certain applications such as image compression where a bit quota has been assigned to the compressed file it is more advantageous to set a certain percentage of wavelet coefficients to zero to satisfy the quota requirement In this case the setting of the threshold value o is based on the histogram and total number of coefficients The thresholding rule is the same as hard thresholding 2 9 4 Universal thresholding In some noise removal applications in which the noise statistics is known it may be more effective to set the threshold value based on the noise statistics For example Donoho and Johnstone 8 set the threshold value to be 2 log gy osl 2 98 where v is the standard deviation of the noise and is the cardinality of the data set This threshold value can be used in either hard or soft thresholding as shown earlier 2 10 TWO DIMENSIONAL EXTENSION OF WAVELET ALGORITHMS 39 Implementation Implementations of the hard quantile
50. ot or gray level display Example The same function used to demonstrate the STFT is also used here to demonstrate the characteristics of the CWT We choose the Morlet wavelet and compute the CWT via FFT The result shown in Fig 2 3 clearly shows the advantage of a flexible window of the CWT It requires only one computation to resolve both the frequencies and the delta functions in 2 82 One can compare this graph of the CWT to that of the STFT 2 4 Computation and Display of Scaling Functions We give two algorithms for computing the scaling functions Spectral domain method Let us recall the two scale relations for the scaling function olt E ped 2t k 2 4 COMPUTATION AND DISPLAY OF SCALING FUNCTIONS 33 CWT Toolware delta pts saaa ad 10 290066 nnangs 91 0353 w 8 ann NI 204 Figure 2 3 CWT outputs and consider the spectral domain expression of this relation namely TIP ei C Rar ea 2 84 where P x is the z transform of the sequence pp Hence the scaling function t can be obtained by finding the IFFT of the infinite prod uct given by 2 84 In general the infinite product must be truncated so that there is only a finite number of points in the IFFT program This approach is not easy to implement because the complexity of the product of the discrete Fourier transform of the sequence px If the number of elements in this sequence is small the problem is still ma
51. ow pass and bandpass children of a set of DWT coefficients When level 2A is chosen the low pass and bandpass children of the level 1A low pass subband which is displayed on DBs are respectively displayed on DB and DB Similarly at level 2B the low pass and bandpass children of the level 1A bandpass subband which is displayed on DBs are respectively displayed in DB and DB gt Rules for use of DB and DB at other levels are listed as follows Level DB and DB display 3A children from the low pass band of 2A 3C children from the low pass band of 2B children from the low pass band of 3B children from the bandpass band of 3B children from the low pass band of 3C 4F children from the bandpass band of 3C 4G children from the low pass band of 3D 4H children from the bandpass band of 3D Fig 3 15 shows DWT coefficients of an signal In its normal display mode the maximum and minimum values of all DWT coefficients in 3 4 1 D DWT TOOL 59 cluding the original data are used to set the display scales This usually leads to loss of visual details of small coefficients in some subbands as expected for the sake of display consistency To remedy this problem the 1 D tool is equipped with a sliding magnifier through which one can better view details of small coefficient values By simply pressing on the left mouse button a magnified version of data surrounding the cur sor will
52. pectrum 2 compression 39 continuous wavelet transform 6 convolution 2 16 38 correlation 8 cubic spline 13 CWT 9 d c offset 5 Daubechies 6 13 22 Daubechies orthonormal wavelet 20 70 decomposition 17 27 decomposition algorithm 17 decomposition relation 15 16 decomposition sequences 15 20 24 delta function 1 2 detection 37 diagonal bands 18 digital processing 18 dilation 11 discrete Fourier analysis 1 discrete Fourier transform 35 discrete wavelet transform 10 divide and conquer 10 downsample 16 27 38 dual B spline 21 dual B spline coefficients 37 dual B spline series 37 dual scaling function 20 21 dual spline 21 dual spline wavelet 21 dual wavelet 20 dual wavelet spaces 21 dual wavelets 21 dual spline 36 37 duals 36 energy 6 24 energy concentration 4 entropy 24 Euler Frobenius Laurent Polyno mial 21 exponential decay 21 Fast Fourier Transform FFT 30 Feauvears 22 INDEX FFT 31 filter bank 10 18 22 35 finite impulse response 22 finite energy functions 5 finitely supported 22 FIR filters 22 first order B spline 17 Fourier analysis 1 Fourier processing 6 Fourier series 1 Fourier transform 1 frequency bands 23 frequency domain 2 18 frequency octaves 10 fundamental frequency 1 gating 39 global periodic function 1 Haar 1 Haar wavelet 17 Haar wavelets 20 hard thresholding 39 hie
53. pose a 2 D image signal and see the resultant component display 1 9 1 Two dimensional wavelet packets Two dimensional wavelet packets are refinements of 2 D wavelets just as in the 1 D case Not only is the LL component the scaling function component decomposed to obtain further details of the image but the other wavelet components LH HL HH can also be further decom posed For example starting with an original image with size 128x128 a 2 D wavelet decomposition of this image will result in four subimages of size 64x64 Continuing the decomposition one gets 16 2 D wavelet packet subimages of size 32x32 The computational algorithm for 2 D wavelet packets is no different than that for the 2 D wavelets It requires an orderly bookkeeping for keeping track of the directions x or y the filters used highpass or lowpass and the resolutions The reader can make use of this wavelet Toolware to practice the wavelet packet decomposition and reconstruction of 2 D signals Part II Wavelet Algorithms This page intentionally left blank 2 1 WAVELET ALGORITHMS OVERVIEW 29 2 1 Wavelet Algorithms Overview This part discusses the approach we take to numerically implement the mathematical statements given in Part I In addition we also in clude several algorithms that have been used in many wavelet signal processing applications such as thresholding noise removal and signal compression The algorithms are written in procedural form s
54. rarchacal pyramid 27 high high HH 27 high low HL 27 highpass 27 histogram 39 infinite product 34 36 infinite sequences 21 inner product 8 19 interpolatory 17 large scale 9 linear B spline 17 linear combination 6 INDEX linear phase 22 linear transform 6 LL component 27 localization property 5 6 localized finite energy basis 2 low frequency 9 low high LH 27 low low LL 27 lowpass 27 main lobe zero crossing 3 Meyer 20 Meyer wavelet 20 Morlet wavelet 33 mother wavelet 6 MRA 10 15 multiband signal decomposition 10 neural network 37 noise reduction 39 noise removal 30 40 41 noise statistics 40 normalized 19 normalized orthogonal set 19 octave bands 23 one variable function 2 orthogonal 19 orthogonal basis 21 orthogonal complementary sub space 11 orthogonal projection 17 18 orthogonal subspace 11 orthogonal sum 11 orthogonality 19 20 22 37 orthonormal 19 71 orthonormal basis 1 18 20 orthonormal set 19 orthonormal wavelet basis 20 orthonormal wavelets 6 Paseval identity 3 perfect reconstruction 20 periodic signals 1 processing filter 27 processing sequences 22 psuedo linear phase 22 quadratic spline 39 quantile thresholding 39 recognition 37 reconstruction 17 reconstruction algorithm 16 reconstruction sequences 20 reconstructive relations 15 rectangular function 2 rectangular window 2 refinement 1
55. re that this is a 2 D processing routine define NAME Simple Smoothing The name of this routine include extfunci h Some Toolware symbols This routine accepts floating point numbers placed in the 2 D array input process the array and then place the results at the 2 D array output The dimension of the matrix is specified by Nyt and My amp TOOLBOX_FUNCTION process float huge input float huge output int x int y int i j for i 1 i lt x 1 i for j 1 j lt y 1 j output i j input i 1 j 1 input i j 1 input i 1 j 1 input i 1 j input i j input i 1 j input i 1 j 1 input i j 1 input i 1 j 1 9 0 One needs to use a compiler to compile a new routine into a Windows dll dynamic link library file and place it into the Toolware directory Otherwise the 1 D tool will not be able to find your custom routine for execution In using Borland C 5 0 one loads the ide project in C into the compiler and makes necessary functional changes Then simply build the dll libary to be saved into the Toolware directory If you use other compilers to generate codes in order to be compatible the dll used must have names exported explicitly see dil map If these files do not match Toolware will assume that the dil is not an external tool designed for it and will ignore it 3 4 1 D DWT TOOL 55 Figure 3 12 1 D data file selection
56. resholding 38 2 9 4 Universal thresholding 38 2 10 Two Dimensional Extension of Wavelet Algorithms 39 2 11 Other Algorithms me se e Gwe lee CaS RES YS GE 40 Part III 41 3 1 The Wavelet Toolware Overview 43 3 1 1 Installation of Toolware 43 312 Wavelets a O55 22 areal raceste e Gia ete Bs al Pe Ee a as 44 3 2 Using Toolware 6 8 4 c asi as eis oe ot ee 45 3 3 2 D DWT T ol oo a erat we ee eho Dek a Gen es 46 3 3 1 Data file selection a0 n cs 4G e we ce So 47 3 3 2 Wavelet decomposition 49 3 3 3 Display thresholding 50 3 3 4 Wavelet reconstruction 51 SA IBD DWT Tool sore uinea one a Pele Aa Toi A Cetea 55 3 4 1 1 D DWT decomposition 56 3 4 2 Viewing and processing the data 57 3 5 Continuous Wavelet Transform Tool 62 3 6 Short Time Fourier Transform Tool 63 3 7 Wavelet Builder Tool n tar ar 6 e act aa de Oe RR aod 64 Preface Wavelet Toolware is a companion software package to the book An Introduction to Wavelets by Charles K Chui It is designed for the reader to gain some hands on practice in the subject of wavelets The objective is to provide basic signal analysis and synthesis tools that are flexible enough for the reader to easily increase the capability of the Toolware by adding processing algorithms to tailor to specific areas of application Wavelet To
57. rfect copy perfect reconstruction condition of the input signal Biorthogonal wavelets of Cohen Daubechies and Feauvears are examples of this class of bases In addition there may not exist a scaling function or wavelet associated with filters designed for filter bank processing In other words there may not be a nested V spaces for arbitrary FIR filters as their two scale sequences For this reason the two channel filter bank is an excellent alternative for processing digital signals Decomposition and reconstruction sequences The two scale sequences of two biorthogonal wavelets of Daubechies et al are listed here as a reference For more details of biorthogonal wavelets readers are referred to Section 5 6 particularly Example 5 41 of the accompanying text An Introduction to Wavelets 3 the text Wavelets A Mathematical Tool for Signal Analysis Chapter 5 by 22 WAVELET TOOLWARE Chui 2 and the engineering texts Wavelets and Subband Coding by Vetterli and Kovacevic Prentice Hall 1995 4 and Wavelets and Filter Banks by Strang and Nguyen Wellesley Cambridge Press 1996 5 Coefficients of the B 97 Daubechies biorthogonal wavelet are given in the following table BiDaub9_7 qk 0 026748757411 0 016864118443 0 091271763114 0 053497514822 0 078223266529 0 045635881557 0 057543526228 0 033728236886 0 266864118443 0 028771763114 0 591271763114 0 156446533058 0 602949018236 0 295
58. s a larger coefficient An example of the CWT is shown in Fig 1 3 For the sake of comparison we write the STFT in the inner product form as ST 3 x t w t bye at x t w t bei 1 15 One immediately sees from 1 13 and 1 15 that the difference be tween CWT and STFT is that the kernels for these two transforms are different CWT and STFT are interchangeable if we switch w t bje o Yalt 1 16 8 WAVELET TOOLWARE CWT Toolware voice pt agag A ur e 8 433100 8 433100 TUNUL az ah NUI i WY NUD AU Y I po MINNI il 0 010000 Figure 1 3 The graphical display of the CWT coefficients of a voice signal We remark here that the basis function w t bei is a sinusoidal function at a single frequency modulated by the window function w t The window widths are fixed independent of b and w Based on the window measures discussed earlier in this manual the window area in the t f plane for the wavelet at a an scale a is a by t a t Ap t a t Ay x w Ay Lu A 1 17 We see clearly that the window area of the wavelet remains the same namely 4AwAw However the time and spectral window widths are dependent on the scale a See Fig 1 4 If a becomes large i e at low frequency the window width in the time domain is large while the window width in the spectral domain becomes small This is the case in processing signals in fine resolution
59. s and Subband Coding Pren tice Hall Inc 1995 5 Strang G and T Ngyuen Wavelets and Filter Banks Wellesley Cambridge Press 1996 6 Daubechies I Ten Lectures on Wavelets SIAM Publ 1992 326 7 Coifman R R and M V Vickerhauser Entropy Based Algorithms for Best Bases Selection IEEE Trans on Information Theory 38 1992 713 718 8 Donohoe D and Johnstone Minimaz estimation via wavelet shrinkage Technical Report Department of Statistics Stanford University 1992 67 Index Nm 20 pr 13 ae 13 ot 13 w t 13 z transform 34 2 D signal 27 2 D wavelet decomposition 27 2 D wavelet packet 27 2 D wavelets 25 27 3db points 3 all pass filters 15 analog processing 18 antisymmetric 22 approximate function 10 approximation space 19 24 38 approximation subspace 10 19 autocorrection sequence 18 autocorrelation matrix 18 autocorrelation sequence 37 B spline 20 35 37 B spline functions 21 B spline series 21 37 B wavelet 20 B wavelet x 20 bandpass filter 5 basic wavelet 6 69 basis function 9 13 basis functions 6 15 19 Battle Lamarie wavelets 20 Battle Lemarie 20 Best bases 24 best level 24 biorthogonal basis 19 22 biorthogonal wavelet 22 biorthogonality 20 boundary effect 31 broadband signal 10 cardinal B spline function 20 chirp signal 5 coefficient sequence 18 Cohen 22 compactly supported 18 20 complex s
60. s on Wavelets SIAM Publ 1992 p 326 6 This theorem is obviously true when f is the scaling function 4 since the two scale relations for and the wavelet y give V 3 90230 Z po k W gt v t Tai k If we apply this theorem to the W spaces we generate the wavelet packet subspaces The general recursive formulas for wavelet packet generation are Halt pebe 2t k 1 73 k Y gnue 2i k kez 1 74 k i Loe41 t where fo and pu w are the scaling function and the wavelet respectively For 1 we have the wavelet packets u2 and pus gen erated by the wavelet u y This process is repeated so that many wavelet packets can be generated from the two scale relations Because of the many components available any given signal can be represented by a choice of wavelet packets at different levels of resolution An op timized algorithm can be constructed to maximize a certain measure such as energy or entropy Best bases and best level are two popular algorithms for signal representations The reader can find these algo rithms in 7 The Toolware is designed so that the decomposition and reconstruction algorithms can use either wavelets or wavelet packets 1 8 1 Wavelet packet algorithms The wavelet packet algorithm follows the wavelet algorithm with the extension that the wavelet coefficients d are also being processed in the same way as the scaling function coefficients c Hence for e
61. se the iteration index by 1 and repeat steps 4 to 6 until 10 iteration cycles have been completed The reader may use the Toolware to familiarize himself herself with scaling functions included in the software 2 5 Computation and Graphical Display of Wavelets To generate a wavelet that satisfies the two scale relation b t Yo a 2t k 2 86 k we can use either the spectral or the time domain approach For the spectral domain approach we simply multiply the infinite product in 2 6 MAPPING BETWEEN B SPLINES AND THEIR DUALS 35 2 84 except that the index k runs from 2 to oo by one half of the z transform of q namely ae aG 5 gt ne HZP ei For the time domain approach we take the linear combination of the translate scaling function 2t k using the coefficients q We use the time domain approach in this Toolware to generate the graphs of different wavelets The reader should use the Toolware to view different wavelets as given in Part C Il p w 2 6 Mapping between B Splines and Their Duals When B splines are chosen to represent a signal s t we have shown in Part I that it is necessary to map the sample values s t to an approximation of the signal s t Ep chj x j t This procedure maps the signal into the B spline or the approximation space at the jth resolution Since the decomposition sequences ap and b for B splines and B spline wavelets are infinitely long a truncati
62. space Vm For example if a signal is sampled at the rate of 256 samples per second we may represent the signal by fa t 5 cs kol 2t k 1 30 k 00 This signal can also be represented by D w _m t fs t 1 31 1 6 1 Two scale relations The equation in 1 21 expresses the relation of a basis function between two different scales The two scale relation formalizes this relationship mathematically Since oft e wc 1 32 Yt Wc 1 33 we can write o t and y t in terms of the bases that generate Vi In other words there exist two sequences p and g4 L such that ot 9 2 t Lied t k Sopedis t 1 34 k Sa 2t k pielt 1 35 k k e x Il In general for any j Z the relationships between V W and Vj 1 are governed by the two scale relations 1 6 MULTIRESOLUTION ANALYSIS 13 Two Scale Relation 2 T T T t T 1 L l L 8 L 500 1000 1500 2000 2500 3000 3500 Figure 1 6 Daubechies D2 scaling function ll 2t Y t X pro 2 ttt k 1 36 X gott k 1 37 k II These formulas of the basis function sets 2 t and 4p 271 are built up by a linear combination of 2 1t For a given px sequence the relations 1 36 and 1 37 can be graphically describable We use the scaling functions of Daubechies D2 and the cubic spline as exam ples The two scale relation is graphically displayed in Fig 1 6 and Fig 1 7
63. ssing As needed any other decom position level within the decomposition range can be seen by clicking on one of the icons Ty When the level is clicked on all DWT coefficients from level 1 to level are displayed in the same win dow In a composed image of DWT coefficients the upper left block is the LL low low band the lower left block is the LH low high band the upper right block is the HL high low band and the lower right block is the HH high high band Although the image is decomposed for four levels the example in Fig 3 7 shows that we can still retrieve level 2 wavelet packet coefficients for processing and display since all coefficients at every level are saved into an internal buffer The screen remains unchanged if an invalid level of an image is chosen This can happen for example if the user tries to view the second level DWT coefficients of an image not yet decomposed wal TR z Hag nisi Sa DE A AL eb ay gt n gt pa ho aF 4 m Ag e A D a a is 7 Priza sub i ts oS m ete gt 3 Phy Li N TAS st a yy a gt i i B 1a abs oes QA SEZ amd nn LHT OY WUL ae a X 5 i lt I n amp rar we ee Taek s ie ae 4 mcr hae led narada enant it Fe Se CA en AI m R we 3 sd Se Y Se wen ee Bet ee a 2 Bi A nd aN Om i Sa yi Figure 3 7 DWT coefficient display 3 3 3 Display thresholding Wavelet co
64. t g t gt Ci ey a t 2 dj aVjx t fist 1 49 We substitute the two scale relations into 1 49 to yield gt Cj k do pebissalt SURF gt djk gt qebirup t 1 50 gt Ciu kBi 0 Again we compare the coefficents of 41 t on both sides of 1 50 to yield Cj 1e So Pe 2k6 k T Qe ir 1 51 k where the right sides of the equations correspond to an up sampling before convolution We use a block diagram Fig 1 8 to illustrate the procedure of these algorithms 1 6 4 Mapping of functions into the approximation space Before we go on any further on the extensions of the decomposition and reconstruction algorithms we first need to show how to obtain the scaling function representation of an input data That is for a given input signal s t we wish to find the representation s t s t D Cjk j rlt 5 Ci eh 2 t k 1 52 k k where cj are the input scaling function coefficients to be processed Many users of wavelets have used the sample discretized values of s t as the input We want to stress the fact that in general s t F Cik 1 53 1 6 MULTIRESOLUTION ANALYSIS 17 If we require s s then the two sides of 1 53 are equal to each other only if the scaling function is interpolatory such as first order B spline that corresponds to the Haar wavelet and the linear B spline It is therefore advisable to obtain c from the sample data s t E before applying the decomposi
65. that e cutii 2 93 k can be computed very efficiently If the final objective of the processing requires the reconstruction of the signal as in image compression one needs to eventually convert the processed coefficients back to B spline coefficients for image quality evaluation and display However if the processing objective is for detection and recognition the dual spline coefficients can be used for neural network training purposes 2 7 Decomposition of a 1 D Signal The decomposition algorithm for scaling function coefficients is de scribed by Cik 2 A2k Cj 1 We can rewrite the expression as Cik gt D Qm e m 2k Cit 2 94 2 8 RECONSTRUCTION OF 1 D SIGNALS 37 which is interpreted as convolution before downsampling by 2 In terms of digital signal processing symbolism it is given in Fig 1 8 The computation of the wavelet coeffieicnts at one level of resolution is carried out in a similar manner namely djk SS Om 2 maak Ci ri a 2 95 Comine these two steps to form the decomposition block as shown later This decomposition block can be repeated and sequentially applied to the scaling function coefficients c _p to yield the wavelet coefficient sequences dj p 1 x P 0 1 2 M Implementation of 2 94 is simple Taking a scaling coefficient set to convolve with the coefficient set ap and downsampling yields the scaling function coefficients at one lower resolution lev
66. the image information because it takes many fewer bits to represent clusters of zero coefficients The reconstructed image will lose fidelity because some high frequency com ponents of the image have been removed by thresholding The algorithmic logic of the 1 D and 2 D wavelet packet algorithms has been discussed in Part I The basic algorithm is no different from the wavelet algorithms except the user has to keep track of the de composition structure A procedural discussion of these algorithms is very lengthy and will not be presented here We encourage the user to practice this algorithm in the Toolware using the images provided or images from other sources Part III Toolware User s Manual This page intentionally left blank 3 1 THE WAVELET TOOLWARE OVERVIEW 43 3 1 The Wavelet Toolware Overview Wavelet Toolware is a 32 bit stand alone software program with all necessary library routines bundled together in a package It has five tools the one and two dimensional 1 D 2 D discrete wavelet trans form DWT tools the continuous wavelet transform CWT tool the short time Fourier transform STFT tool and the wavelet generation tool 1 D 2 D DWT tools respectively support DWT decomposition DWT coefficient processing and DWT reconstruction of one and two dimensional signals The CWT tool reads in a 1 D signal data file computes its CWT coefficients and then displays the coefficients on the screen Similarly the STFT tool
67. tice that the relationships in 1 63 to 1 65 are all referred to the same level of resolution The relationships between wavelets at different resolution levels are also orthogonal That is lt Weis Peg gt di j k 1 66 Examples of the orthonormal wavelet basis include the Haar wavelets the Daubechies orthonormal wavelet bases of all orders the Meyer wavelets of all orders and Battle Lemarie orthonormal spline wavelets not compactly supported We note here that both Meyer and Battle Lamarie wavelets are not compactly supported finite time duration Decomposition and resconstruction sequences The relations between the decomposition sequences ap bx and the reconstruction sequences px qx are the simplest for an orthonormal basis For perfect reconstruction of the signal we may use the following relations ge pia ak Pk bk 4k Once the p sequence has been found through the construction of the two scale relation of the scaling function the other sequences are de termined 1 7 2 Semi orthogonal wavelet bases In the case of cardinal B spline functions of orders higher than 1 and their corresponding compactly supported spline wavelets the integer translates of the scaling functions i e B splines as well as the B wavelets are not orthogonal That is ie Nn t DNa e jar 7 bi 1 67 ji Vn T 2 Py p z jide F i j 1 68 20 WAVELET TOOLWARE However the orthogonality bet
68. tion algorithm In the following we consider the orthogonal projection of s t onto the V space Assuming V is a subspace of L and s t L we wish to find s t V as an approximation to s t Write s t do cinbialt 1 54 k and consider s t to be the orthogonal projection of s t onto the V subspace Then s t s t is orthogonal to V so that lt s t 5 t fie gt 0 1 55 lt s t Pje gt a So cinbialt bj u t gt 1 56 k In matrix form this equation becomes ay ao ay Cik _ gj lt S Pj k gt 1 57 ap ai i ao where sa i f opt m dt 1 58 is the autocorrection sequence of the scaling function If the scaling function is compactly supported the corresponding autocorrelation ma trix in 1 57 is banded with a finite number of diagonal bands Also if the scaling function forms an orthonormal basis then m bm 0 1 59 If we have only the sample values of the signal s gt s t we have to approximate the integral in 1 56 by a sum Ce lt s t pje gt a ga Oat x s t Su 1 60 k 18 WAVELET TOOLWARE This demonstrates the difference between the coefficients of the scal ing function series and the sample values of the signal The former is used in wavelet signal processing analog processing in the time do main since the coefficient sequence of the scaling function series and the wavelet series are to be processed while the
69. tive disciplines Andrew K Chan College Station Texas Steve J Liu February 13 1998 This page intentionally left blank Acknowledgments The developmental work of this software for training in wavelet algo rithms is based on our experience in teaching short courses for the Texas Engineering Extension Services over the last few years Our gratitude goes to Professor Charles K Chui for his organization of the short courses help in the review of this manuscript and many valuable suggestions Some of the computer code evaluations and graphical illus trations were carried out by our former and current students Howard Choe Hung Ju Lee Nai wen Lin Dongkyoo Shin Hsien Hsun Wu Michael Yu Xigiang Zhi We thank them for their careful and detailed work Special thanks are due to Brandon Scott Wadsworth who spent countless hours perfecting the GUI of the software We are indebted to Margaret Chui for her help in the editorial work of this manuscript This page intentionally left blank To Jesus Christ our Lord and To our wives Sophia and Grace This page intentionally left blank Part I Wavelet Theory This page intentionally left blank 1 1 FROM FOURIER ANALYSIS TO WAVELET ANALYSIS 1 1 1 From Fourier Analysis to Wavelet Analysis Before the work of Haar 1 a continuous function or an analog signal f t was represented either by an entire series polynomial basis or by a Fourier series sinusoidal basis Since both
70. very decomposition process the number of components is double that of the previous resolution Suppose an input signal s n has been mapped into the approximation space V That is s n lt gt cj Using the de composition sequences a and b we obtain the scaling function series coefficients c _ and the wavelet series coefficients d _1 We can rep resent the scaling function component and the wavelet component at 24 WAVELET TOOLWARE s n level 0 PR natia SERI AN AN level 1 Soo So Sio N ala ZX VA S 100 So no m level 3 Figure 1 9 Wavelet packet decomposition tree resolution j 1 solk lt 7 Clik 1 75 s k os dig 1 76 The coefficient sequences c _1 and d _a are decomposed through an and b into four sequences dijo Vj 2k 1 77 w Pama RE ee aci DN cae i z ee eee Se in which the last two sequences uj_2 and v 2y are wavelet packet series coefficients This process can be repeated for the four sequences in 1 77 to produce eight coefficient sequences corresponding to eight wavelet packet components of the original signal s n The decomposi tion tree structure is given in Fig 1 9 The algorithm for WP decom position and reconstruction of a one dimensional signal is shown in Fig 1 10 1 9 Two Dimensional Wavelets When the input signal is two dimensional 2 D it is necessary to rep resent the signal components by two dimensional wavelets and two dimensional scaling functions
71. ween the scaling function Nm and the B wavelet 4 as well as translates of wavelets at different resolutions is still preserved In this case the compactly supported spline wavelet is called a semi orthogonal wavelet In order to facilitate efficient com putations of the expansion coefficients in the approximation spaces or the wavelet spaces a dual scaling function and a dual wavelet 4 can be found to satisfy the biorthogonality relation lt dni dej gt di ja 1 69 lt es Pez gt ijk b 1 70 where Nm and Y wy have been used in these two equations We remark here that both the dual scaling function the dual B spline and the dual spline wavelet are both functions in the spline space of the same order Hence these duals can be expanded into B spline series of the same spline order We will show in later sections how to use these relations to map coefficients between wavelet spaces and dual wavelet spaces Decomposition and reconstruction sequences The sequences for the two scale relations in the m order spline space are 1 p P 0skam 0 otherwise 1 F m qn S T Nom k 1 2 amp 0 3m 2 Qm 1 ari The two scale relations for the dual spline and dual wavelets are slightly more complicated Since the B spline functions of order m gt 2 are not orthogonal basis the sequences ap bp comes from the dual spline and dual wavelets which are not compactly supported but they are both

Wavelet Toolware - EBM English Books

Contents

Download Pdf Manuals

Related Search

Related Contents