COMPRESSION, ESTIMATION, AND ANALYSIS OF ULTRASONIC SIGNALS

BY GUILHERME CARDOSO DE CARDOSO

Submitted in partial fulfillment of the requirements for the degree of Doctor in Philosophy in Electrical Engineering in the Graduate College of the Illinois Institute of Technology

Approved _________________________ Adviser

Chicago, Illinois May 2005

© Copyright by GUILHERME CARDOSO DE CARDOSO May 2005

ACKNOWLEDGEMENT “O importante é alegria no coração”

iii

TABLE OF CONTENTS Page ACKNOWLEDGEMENT

.......................................................................................

iii

LIST OF TABLES

...................................................................................................

vi

LIST OF FIGURES

.................................................................................................

vii

ABSTRACT

............................................................................................................. xiii

CHAPTER 1. INTRODUCTION

..............................................................................

1.1 Thesis Summary

1

....................................................................

3

...................................................................

5

2.1 Data Compression .................................................................. 2.2 Ultrasonic Imaging ................................................................ 2.3 Transform Coding, Subband Coding, and Data Compression .......................................................................... 2.4 Conclusions ............................................................................

5 7 10 18

3. THRESHOLDING TECHNIQUES FOR DENOISING AND COMPRESSING ULTRASONIC SIGNALS ...........................

19

2. LITERATURE REVIEW

3.1 3.2 3.3 3.4

Introduction ............................................................................ Noise Suppression .................................................................. Data Compression .................................................................. Conclusions ............................................................................

19 20 52 67

4. OPTIMAL DWT KERNEL FOR ULTRASONIC DATA REPRESENTATION ..............................................................

69

4.1 4.2 4.3 4.4 4.5

Introduction ............................................................................ Discrete Wavelet Transform .................................................. Optimal Design of the Echo Wavelet .................................... Performance Evaluation of Wavelet Kernels ......................... Conclusions ............................................................................

iv

69 70 72 75 86

5. PARAMETER ESTIMATION ALGORITHM FOR GAUSSIAN DECOMPOSITION OF ULTRASONIC SIGNALS ..........................

87

5.1 Introduction ............................................................................ 87 5.2 Successive Parameter Estimation Algorithm ......................... 89 5.3 Performance Evaluation with Simulated and Experimental Ultrasonic Signals .................................................................. 111 5.4 Conclusions ............................................................................ 120 6. PARAMETER ESTIMATION ALGORITHM FOR CHIRPLET DECOMPOSITION OF ULTRASONIC SIGNALS .......................... 122 6.1 Introduction ............................................................................ 6.2 Successive Parameter Estimation Algorithm ......................... 6.3 Performance Evaluation Using Simulated and Experimental Echoes .................................................................................... 6.4 Conclusions ............................................................................

122 123 127 133

7. RECONFIGURABLE ARCHITECTURE FOR ULTRASONIC SIGNAL COMPRESSION AND TARGET DETECTION ............... 134 7.1 7.2 7.3 7.4

Introduction ............................................................................ Ultrasonic Target Detection ................................................... An Unified Architecture Design ............................................ Conclusions ............................................................................

8. CONCLUSIONS BIBLIOGRAPHY

134 135 138 145

................................................................................. 147

.................................................................................................... 150

v

LIST OF TABLES Table 3.1 Coefficients Sorted by Position

Page ..................................................................

25

..............................................................

25

3.2 Coefficients Sorted by Amplitude

3.3 Original Parameters of Ten Interfering Simulated Echoes

........................

3.4 SNR of Broadband and Narrowband Signals Compressed by 95%

33

...........

61

3.5 Kernel Size and Energy Compaction

.........................................................

64

4.1 Daub4 and Echo4 Wavelet Kernels

...........................................................

74

4.2 Coif6 and Echo6 Wavelet Kernels

.............................................................

74

..........................................................

75

4.3 Symm8 and Echo8 Wavelet Kernels

4.4 Monte Carlo Simulation Result Summary 4.5 B-Scan Compression Ratio Comparison

.................................................

76

....................................................

82

5.1 Two Interfering Echoes, Original and Estimated Parameters

.................... 109

5.2 Original and Estimated Parameters of Multiple Interfering Simulated Echoes ....................................................................................... 113 6.1 Parameter Estimation of Two Interfering Echoes

vi

...................................... 128

LIST OF FIGURES Figure

Page

2.1

Ultrasonic Imaging System Block Diagram

2.2

Transform Coding System

2.3

Subband Coding System

2.4

Discrete Wavelet Transform

...................................................................

16

3.1

Thresholding Techniques, Before (Dashed Lines) and After Thresholding (Solid Lines). a) Hard Threshold, b) Soft Threshold, and c) Garrote Threshold .........................................................................

22

3.2

Flowchart of ATF Algorithm

..................................................................

26

3.3

a) 2048 Samples of Gaussian Process Realization, b) Sorted Samples of the Random Noise Realization in a) ...................................................

29

a) 2048 Samples of Uniform Process Realization, b) Sorted Samples of the Random Noise Realization in a) ...................................................

30

SNR Enhancement of Adaptive Thresholding Function (ATF), Hard Thresholding (HT), Soft Thresholding (ST), and Garrote Thresholding (GT) techniques with DCT Applied to a Single Ultrasonic Echo Corrupted with Uniform Noise .....................................

34

SNR Enhancement of Adaptive Thresholding Function (ATF), Hard Thresholding (HT), Soft Thresholding (ST), and Garrote Thresholding (GT) techniques with DWT Applied to a Single Ultrasonic Echo Corrupted with Uniform Noise .....................................

35

3.4 3.5

3.6

............................................

8

.......................................................................

11

.........................................................................

12

3.7

SNR Enhancement of Adaptive Thresholding Function (ATF), Hard Thresholding (HT), Soft Thresholding (ST), and Garrote Thresholding (GT) techniques with WHT Applied to a Single Ultrasonic Echo Corrupted with Uniform Noise ..................................... 36

3.8

SNR Enhancement of Adaptive Thresholding Function (ATF), Hard Thresholding (HT), Soft Thresholding (ST), and Garrote Thresholding (GT) techniques with DCT Applied to 10 Ultrasonic Echoes Corrupted with Uniform Noise ............................

3.9

SNR Enhancement of Adaptive Thresholding Function (ATF), Hard Thresholding (HT), Soft Thresholding (ST), and Garrote vii

37

Thresholding (GT) techniques with DWT Applied to 10 Ultrasonic Echoes Corrupted with Uniform Noise .................................

38

3.10 SNR Enhancement of Adaptive Thresholding Function (ATF), Hard Thresholding (HT), Soft Thresholding (ST), and Garrote Thresholding (GT) techniques with WHT Applied to 10 Ultrasonic Echoes Corrupted with Uniform Noise ...................................................

39

3.11 SNR Enhancement of Adaptive Thresholding Function (ATF), Hard Thresholding (HT), Soft Thresholding (ST), and Garrote Thresholding (GT) techniques with DCT Applied to a Single Ultrasonic Echo Corrupted with Gaussian Noise .....................................................

40

3.12 SNR Enhancement of Adaptive Thresholding Function (ATF), Hard Thresholding (HT), Soft Thresholding (ST), and Garrote Thresholding (GT) techniques with DWT Applied to a Single Ultrasonic Echo Corrupted with Gaussian Noise .....................................................

41

3.13 SNR Enhancement of Adaptive Thresholding Function (ATF), Hard Thresholding (HT), Soft Thresholding (ST), and Garrote Thresholding (GT) techniques with WHT Applied to a Single Ultrasonic Echo Corrupted with Gaussian Noise .....................................................

42

3.14 SNR Enhancement of Adaptive Thresholding Function (ATF), Hard Thresholding (HT), Soft Thresholding (ST), and Garrote Thresholding (GT) techniques with DCT Applied to 10 Ultrasonic Echoes Corrupted with Gaussian Noise ..................................................

43

3.15 SNR Enhancement of Adaptive Thresholding Function (ATF), Hard Thresholding (HT), Soft Thresholding (ST), and Garrote Thresholding (GT) techniques with DWT Applied to 10 Ultrasonic Echoes Corrupted with Gaussian Noise ..................................................

44

3.16 SNR Enhancement of Adaptive Thresholding Function (ATF), Hard Thresholding (HT), Soft Thresholding (ST), and Garrote Thresholding (GT) Techniques with WHT Applied to 10 Ultrasonic Echoes Corrupted with Gaussian Noise ..................................................

45

3.17 SNR Enhancement of Thresholding Techniques in an Ultrasonic Signal with 10 Interfering Echoes Corrupted by Uniform Noise with SNR = -2.5dB using DCT. a) Original Signal, b) Noisy Signal, Reconstructed Signal Using c) Adaptive Thresholding Function (ATF), d) Hard Thresholding (HT), e) Soft Thresholding (ST), and f) Garrote Thresholding (GT) ...................................................................................

46

3.18 SNR Enhancement of Thresholding Techniques in an Ultrasonic viii

Signal with 10 Interfering Echoes Corrupted by Uniform Noise with SNR = -2.5dB using DWT. a) Original Signal, b) Noisy Signal, Reconstructed Signal Using c) Adaptive Thresholding Function (ATF), d) Hard Thresholding (HT), e) Soft Thresholding (ST), and f) Garrote Thresholding (GT) ...................................................................................

47

3.19 SNR Enhancement of Thresholding Techniques in an Ultrasonic Signal with 10 Interfering Echoes Corrupted by Uniform Noise with SNR = -2.5dB using WHT. a) Original Signal, b) Noisy Signal, Reconstructed Signal Using c) Adaptive Thresholding Function (ATF), d) Hard Thresholding (HT), e) Soft Thresholding (ST), and f) Garrote Thresholding (GT) ...................................................................................

48

3.20 SNR Enhancement of Thresholding Techniques in an Ultrasonic Signal with 10 Interfering Echoes Corrupted by Gaussian Noise with SNR = 0dB using DCT. a) Original Signal, b) Noisy Signal, Reconstructed Signal Using c) Adaptive Thresholding Function (ATF), d) Hard Thresholding (HT), e) Soft Thresholding (ST), and f) Garrote Thresholding (GT) ...................................................................................

49

3.21 SNR Enhancement of Thresholding Techniques in an Ultrasonic Signal with 10 Interfering Echoes Corrupted by Gaussian Noise with SNR = 0dB using DWT. a) Original Signal, b) Noisy Signal, Reconstructed Signal Using c) Adaptive Thresholding Function (ATF), d) Hard Thresholding (HT), e) Soft Thresholding (ST), and f) Garrote Thresholding (GT) ...................................................................................

50

3.22 SNR Enhancement of Thresholding Techniques in an Ultrasonic Signal with 10 Interfering Echoes Corrupted by Gaussian Noise with SNR = 0dB using WHT. a) Original Signal, b) Noisy Signal, Reconstructed Signal Using c) Adaptive Thresholding Function (ATF), d) Hard Thresholding (HT), e) Soft Thresholding (ST), and f) Garrote Thresholding (GT) ...................................................................................

51

3.23 Ordered Transform Coefficients of Ultrasonic Echoes: a) Narrowband and b) Broadband ...........................................................

56

3.24 Narrowband Original Echo (Dashed Line) Superimposed to Most Energetic Coefficient of a) DCT, b) DWT Haar, c) WHT, and d) DWT Daub20 ...............................................................................

57

3.25 Broadband Original Echo (Dashed Line) Superimposed to Most Energetic Coefficient of a) DCT, b) DWT Haar, c) WHT, and d) DWT Daub20 ......................................................................................

58

3.26 Relation Between NBW and the Five Most Energetic Transform ix

Coefficients

.............................................................................................

59

3.27 a) Broadband signal, b) ultrasonic signal compressed 95% using DCT, c) DWT (Haar), d) DWT (Daubechies), e) WHT; f) Broadband signal, g) ultrasonic signal compressed 95% using DCT, h) DWT (Haar), i) DWT (Daubechies), j) WHT ................................................................

60

3.28 Wavelet Kernels: a) Haar, b) Daubechies 20, c) Beylkin, d) Coiflet, e) Symmlet, and f) Vaidyanathan ............................................................

62

3.29 Energy Accumulated Among the Five Most Dominant DWT Coefficients Using the Following Kernels: a) Haar, b) Daubechies, c) Beylkin, d) Coiflet, e) Symmlet, and f) Vaidyanathan ........................

63

3.30 a) Original Experimental Signal, and Reconstruction with 98% Compression Ratio: b) DCT, c) DWT Haar, d) DWT Daub, and e) WHT .............................................................................................

65

3.31 a) Original Experimental Signal, and Reconstruction with 92% Compression Ratio: b) DCT, c) DWT Haar, d) DWT Daub, and e) WHT .............................................................................................

66

4.1

Discrete Wavelet Transform

...................................................................

70

4.2

Objective Function

..................................................................................

72

4.3

Experimental A-Scan Data from a Thin Metal Sample

..........................

76

4.4

Fitting Results of MCS with Daub4 and Echo4. a) Echo4 with HT, b) Echo4 with ST, c) Daub4 with HT, and d) Daub4 with ST ................

78

Relation between Normalized MSE and Compression Ratio Using Daub4 and Echo4 Wavelet Kernels .........................................................

79

Relation between Normalized MSE and Compression Ratio Using Coif6 and Echo6 Wavelet Kernels ..........................................................

80

Relation between Normalized MSE and Compression Ratio Using Symm8 and Echo8 Wavelet Kernels .......................................................

81

4.8

Original B-Scan Image

83

4.9

Original Image Compressed 87% Using Echo4

......................................

84

4.10 Original Image Compressed 87% Using Daub4

.....................................

85

4.5 4.6 4.7

............................................................................

x

5.1

Effect of Uncertainty of the Values of Parameters on Reconstruction Error, Er ...................................................................................................

91

5.2

Successive Parameter Estimation Flowchart

........................................... 102

5.3

a) Original Signal (Dashed Line, Θ = [6, 1, 4, 0.6, 3.5] ), ˆ = [6, 1, 4, 0.6, 3.5] ); Estimated Signal (Solid Line, Θ b) Time×Frequency Representation of Signal in (a); c) Phase×Bandwidth Representation of Signal in (a) .............................. 103

5.4

Input and Output SNR

5.5

a) Echo 1 with Moderately Poor SNR, b) Estimated Echo 1 Superimposed with Actual Echo, c) TF Representation of Echo 1, d) TF Representation of Estimated Echo 1, e) Echo 2 with Severely Poor SNR, f) Estimated Echo 2 Superimposed with Actual Echo, g) TF Representation of Echo 2, and h) TF Representation of Estimated Echo 2 ................................................................................ 106

5.6

a) TF Representation of Two Interfering Echoes, b) Projection in the Frequency Domain, and c) Projection in the Time Domain .................... 108

5.7

Original (Solid Line) and Estimated (Dashed Line) Signals, a) Two Interfering Echoes, b) First Echo, and c) Second Echo ........................... 110

5.8

a) Original Signal, b) TF Representation of Signal in (a), c) Original Signal Corrupted with Noise, d) TF Representation of Signal in (c), e) Estimated Signal, and f) TF Representation of Signal in (e) .............. 112

5.9

Sequence of Single Echo Estimations from Highest Energy Echo Shown in (a) to the Lowest Energy Echo Shown in (j) Superimposed with the Original Echoes. The Solid Line Represents the Original Echoes, and the Dashed Line Represents the Estimated Echoes ............ 114

............................................................................. 105

5.10 Energy of Estimated Echoes from the Highest Energy (Echo “a”) to the Lowest Energy (Echo “k”) ............................................................ 115 5.11 SNR Enhancement of Successive Parameter Estimation Algorithm as Function of Compression Ratio ............................................................... 118 5.12 a) Experimental Signal, b) TF Representation of Signal in (a), c) Estimated Signal, and d) TF Representation of Signal in (c) 6.1

.............. 119

a) TF Representation of Two Interfering Echoes, b) Projection in the Frequency Domain, c) Projection in the Time Domain .......................... 129

xi

6.2

Actual (Solid Line) and Estimated (Dashed Line) Signals, a) Two Interfering Echoes, b) First Echo, and c) Second Echo ........................... 130

6.3

a) Experimental Signal, b) TF Representation of Signal in (a), c) Estimated Signal, d) TF Representation of Signal in (c) ..................... 132

7.1

Ultrasonic Target Detection System

7.2

Target Detection Results. a) Ultrasonic Data, b) Target-to-Clutter Enhancement Using DCT, c) Target-to-Clutter Enhancement Using DWT ............................................................................................. 137

7.3

Timing Requirements

.............................................................................. 138

7.4

System Components

................................................................................ 140

7.5

Processing Element (PE) Architecture

7.6

Forward DWT Implementation Using Multiple PE Arrays

7.7

Recursive Structures for DCT. a) Even Coefficient, b) Odd Coefficient Computation ......................................................................... 145

xii

........................................................ 135

.................................................... 142 .................... 143

ABSTRACT Ultrasonic imaging in medical and industrial applications often requires a large amount of data collection.

Consequently, it is desirable to use data compression

techniques to reduce data size and to facilitate the analysis and remote access of ultrasonic information. Hence, the locally obtained ultrasonic signals can be transferred efficiently through wireless or wired communication channels to the remotely located experts. In this research, we analyze different signal processing techniques to compress and denoise ultrasonic signals. Due to the increasing demand on hardware capabilities for ultrasonic applications, we also developed a reconfigurable hardware architecture implementation of an ultrasonic signal processor. The precise ultrasonic data representation is paramount to the accurate analysis of the shape, size, and orientation of ultrasonic reflectors, as well as to the determination of the properties of the propagation path. We address the ultrasonic signal representation by introducing a successive parameter estimation technique that is also able to compress and denoise ultrasonic signals. This technique uses a modified version of the continuous wavelet transform to decompose the ultrasonic signal in Gaussian shaped echoes. Furthermore, a chirplet transform is employed to decompose the ultrasonic signal in chirp-shaped echoes. This technique provides a high resolution and accurate estimation of the echo parameters, which allows the solution of the ultrasonic deconvolution problem. The algorithm also performs well in noisy environments where signal-to-noise ratio (SNR) enhancements beyond 60dB are feasible. Overall, the signal modeling and parameter estimation algorithm presented in this study not only offers data compression

xiii

capabilities, but also provides parameters that can be used for signal deconvolution, target detection, pattern recognition, and material characterization. Ultrasonic data is often embedded in noise. Hence, we introduce a technique to estimate an adaptive thresholding function that uses the statistical parameters of the noise embedded in the signal. Then, the statistical parameters are used to generate a thresholding function based on the probability distribution function of the noise. We analyze the performance of adaptive and classical thresholding techniques when applied to the discrete wavelet transform (DWT), discrete cosine transform (DCT), and WalshHadamard transform (WHT) coefficients. The results presented show that the adaptive thresholding technique is a very powerful method that allows the detection of low SNR ultrasonic backscattered echoes. This adaptive thresholding technique also achieves SNR improvements above 10dB over the classical thresholding techniques. We also develop subband and transform coding techniques to compress ultrasonic signals. In particular, the data compression performance of the DCT, the WHT, and the DWT are examined using simulated and experimental ultrasonic data. The results obtained show that the DWT is better in the representation of broadband signals, while the DCT and the WHT are more suitable in the representation of narrowband signals. For a narrowband ultrasonic signal the five most energetic coefficients of the DCT accumulate 95% of the total signal energy, while the DWT and the WHT accumulate 50% and 60% of the total signal energy respectively. On the contrary, for a broadband ultrasonic signal the 5 most dominant DCT coefficients accumulate 30% of the total signal energy, while the DWT and the WHT accumulate 97% and 15% of the total signal energy respectively.

xiv

The advances in the field programmable gate array (FPGA) area have brought this technology to a point where very complex architectures can be implemented at a relatively low cost. Thus, we also introduce the implementation of an ultrasonic data compressor that meets the stringent requirements of ultrasonic signal processing such as high speed, high data volume, and reconfigurability. This reconfigurable architecture for ultrasonic signal compression and target detection relies on a unified hardware implementation that is made possible since the DWT and DCT algorithms are designed to share subband decomposition logic and adaptable thresholding. This architecture is a flexible and efficient solution for real-time ultrasonic imaging systems where low-power and compactness are critical.

xv

1 1.

CHAPTER I

1 INTRODUCTION The field of ultrasonic engineering uses high frequency sound waves above the human hearing range (around 20kHz) [Fre65]. In particular, medical and industrial ultrasonic imaging applications use frequencies in the range between 500kHz and 20MHz. The goal of the ultrasonic application is to image the internal structure of materials by transmitting acoustic waves into the specimen using an ultrasound transducer. The acoustic waves travel inside the material under test, and the same transducer is used to measure the backscattered echoes. The analyses of the reflected echoes provide an important insight on the physical properties of the material. The properties of interest in the ultrasonic analysis include the microstructure, location, geometric shape, size, and orientation of the propagation path. Some of the echo parameters are fundamental for the evaluation of the specimen under test. The time-offlight is defined as the time that takes the echo to travel inside the material and bounce back. This parameter can be correlated to the position of the target. The amplitude of the backscattered echo shows the target orientation. Since the material is not necessarily linear, the frequency and bandwidth of the reflected echo shows how the material can modify the frequency of a sinusoidal wave in its propagation path. The phase of the reflected echo can give us indications on the angle of the target [San95]. This set of parameters can provide a good determination of the material properties. The ultrasonic measurement system can be modeled as the convolution of the input pulse with the propagation path. The propagation path in this model is similar to a filter,

2 and the output of such filter is the backscattered echo. The challenge of the ultrasonic signal echo estimation is to successfully deconvolve the single echoes that constitute the signal. Since the material medium is composed of several reflectors other than the target, the convolution of this multiple reflectors increase the level of complexity of our analysis. In this study, we address the parameter estimation, compression, and denoising of ultrasonic signals using different techniques. The transformation of the ultrasonic signal to different domains using the discrete cosine transform (DCT), discrete wavelet transform (DWT), and the Walsh-Hadamard transform (WHT) is analyzed. Once the ultrasonic signal is represented in these transform domains, we apply different thresholding techniques to compress and denoise the data. The effects of thresholding each of the transforms are analyzed, and the results with simulated and experimental data are presented. Furthermore, we analyze the different thresholding techniques and how they affect the compression and denoising performance. Applications where the ultrasonic signal is embedded in a high noise environment are not uncommon. In order to be able to recover the ultrasonic echoes in such applications, we developed a technique to estimate an adaptive thresholding function with denoising capabilities that outperform classical thresholding techniques. The DWT is a very important tool in signal representation. One important characteristic of the DWT transform is that we can select the correlation kernel (this is not the case for DCT, and WHT). The success of the DWT in data compression is closely tied to the correlation of the wavelet kernel to the ultrasonic signal. Therefore, we study the design of the wavelet kernel via optimization techniques for high data compression

3 ratios. The goal is to improve the similarity of the wavelet kernel to the ultrasonic signal, keeping some important properties of the wavelet kernel such as uniqueness and orthogonality. Given that these properties are very strict, we analyze the effects of relaxing some of these features in the data compression performance. For higher compression ratios, a model-based estimation algorithm is used to estimate parameters of the echo wavelets within an ultrasonic signal. This algorithm assumes that the ultrasonic signal is composed of a sequence of Gaussian or chip-shaped echoes with different parameters. This model approach has been proven to represent realistic ultrasonic signals with high level of accuracy [Dem01, Car04, Car05]. High compression ratios are achieved since a finite number of parameters are able to represent a complex ultrasonic signal. An ultrasonic data compressing and denoising implementation is also analyzed, and a reconfigurable architecture for such application is presented. This architecture is designed to address several problems of ultrasonic imaging applications, such as high data volume, high acquisition speeds, and flexibility requirements. 1.1

Thesis Summary

Chapter 2 presents a literature review on data compression and ultrasonic signal imaging. This literature review also includes the relationship between the different transform and subband coding techniques and their data compression performance. Chapter 3 studies thresholding techniques for denoising and compressing ultrasonic signals using adaptive and classical thresholding functions. In Chapter 4, we present a technique to design an optimal DWT kernel for ultrasonic data compression and denoising. For high fidelity ultrasonic signal representation, Chapter 5 introduces a

4 successive ultrasonic parameter estimation algorithm that relies on a Gaussian-shaped echoes decomposition of the ultrasonic signal. In Chapter 6, we expand the analysis in Chapter 5 for chirp-shaped ultrasonic echoes. In Chapter 7, we present a reconfigurable architecture implementation of an ultrasonic data compressor. Chapter 8 summarizes the thesis and the contributions of the research.

5 2.

CHAPTER II

2 LITERATURE REVIEW

2.1

Data Compression

Many of the modern technological accomplishments were enabled in part by improvements in data compression algorithms [Say00]. Among the applications that make extensive use of data compression are digital TV, mobile communications (audio and video), modems, telephony, satellite communications, medical diagnostic imaging, and file transfer over networks. Data compression allows data to flow over media where the information size exceeds the media bandwidth or storage capabilities. The subject of data compression has been in the literature for many years. The goal of data compression algorithms is to reduce the size required to store or transmit information. In order to accomplish this goal, the first step is to identify structures (e.g. statistical patterns) in the data. The correct determination of data structures is paramount for the success of data compression algorithms. Once the data is correctly identified, one can represent the “building blocks” of the data instead of the original information. On the other hand, not all data sources produce structured data. In this case a limited set of algorithms can be used to compress the data. These algorithms do not rely on data structure; hence limited data compression can be achieved. Data compression algorithms can also use previous knowledge of the user as a means to increase information compression. In this case, information that the user will not miss can be removed from the data without loss of quality.

6 Given that network speeds have increased substantially in the past years, one could argue that motivation for data compression could dim with time. Although this is a valid argument, reality has shown that Parkinson’s law [Par57] prevails. This law states that resource use expands to meet the resources available. In practice the need for transmission and storage has grown twice as fast as the transmission and storage capabilities improve [Say00]. Other than that, some transmission media, such as copper cables, fiber optics, and the atmosphere, have physical limitations to the maximum speed information can flow through them. In this scenario many application will always benefit from data compression algorithms. Although data compression algorithms have been the center of the scientific community in the past decades, data compression has been around since the beginning of human communication. Examples of data compression can be found most languages in the way information is represented in symbols or words. The words used more often are usually expressed with fewer letters, as less frequent words demand more letters. In Morse code letters used more often (e.g. “a” and “e”) are represented with shorter sequences of dots and traces. Braille code is another example of compressed data representation. There are 64 possible symbols in Braille coding, of which 26 are used to represent the letters of the alphabet, and the other 38 symbols represent words often used (e.g. “and” and “for”). The set of data compression algorithms where the original information can exactly be retrieved from the compressed data is called lossless algorithms. Since no information is lost in the compression process, these algorithms are widely used in application such as medical diagnosing, astronomical imaging, and text compression. One can realize that the

7 loss of information in any the mentioned situations could compromise the outcome of the data analysis. On the other hand, there are many applications that can handle some loss of information. In these cases lossy compression can be used, and higher compression ratios can be achieved. This algorithm has been implemented for the compression of speech for mobile and fixed telephony systems, images and video for portable digital cameras, and for many other applications. Furthermore, most lossy compression algorithms also allow the compression to be increased or decreased, depending on the user’s data reconstruction quality requirements. 2.2

Ultrasonic Imaging

Ultrasound are very useful in the identification of impedance discontinuities inside inhomogeneous materials. The ultrasonic waves can be aimed in one direction inside the material, and they can penetrate deep inside the specimen under test. The technique is also rapid and relatively low cost. Flaw detection by using ultrasound is possible because sound waves are reflected from a crack or other abrupt changes in the elastic properties of the material in which the waves are traveling. The ultrasonic waves excite the same transducer in their way back, which converts the acoustic waves back to electrical signals. The acquisition of a series of such echoes allows the composition of images of the specimen. The block diagram of a generic ultrasonic imaging system is shown in Figure 2.1. This figure also shows the ultrasonic echo transmitted inside the specimen under test. The transmitted echo, in this example, reflects on the material discontinuities and on the bottom of the specimen.

8 • Transducer: converts the electrical signal to acoustic pressure and vice-versa, used as

the measuring probe. • Electrical pulse generator: generates the electrical signal of ultrasonic echo, used to

excite the transducer. • Display unit: displays the transmitted and backscattered echo. The display unit is

usually part of the data analysis system, but can also be a standalone oscilloscope. • Data analysis system: processes the bounced echo and applies algorithms for the correct

identification of the echoes. Also used to concatenate the signals received to generate ultrasonic images. • Digitizer: converts the analog signal from the transducer to digital signal sent to the

data analysis system.

Figure 2.1. Ultrasonic Imaging System Block Diagram The medium through which the ultrasonic echo was propagated modifies the parameters of this echo, so that these modifications can be used to characterize specimen under test. Industrial applications use ultrasonic imaging to determine the different layers, defects, and/or grains of random size and distributions inside materials. Since the defects

9 and discontinuities inside the material reflect the acoustic wave, this method can also be used in medical diagnostic applications. Different tissue and organ composition represent differences in the acoustic impedance of the specimen; hence the backscattered echoes provide information about the composition of the features embedded in such tissues and organs, which are related to the differences in tissue and organs types and biological states [San81]. Ultrasonic imaging applications in general deal with reflected echoes that are usually complex due to heterogeneities of the material. These applications thus rely on sophisticated signal processing algorithms to discern the backscattered echoes. The ultrasonic wave is modified in general due to three different factors [Sun99]: absorption, scattering, and beam spreading. As the acoustic wave travels through the material, some of the ultrasonic energy is absorbed in the material generating heat. This energy cannot be recovered and cannot be used in the inspection of the material. The many discontinuities inside the specimen will scatter the ultrasonic wave, and that’s the way ultrasonic imaging is performed. The geometry of the transducer determines the ultrasonic beam spreading. The types of displays of such backscattered echoes were named after similar displays used in radar. Depending on the application, the ultrasonic imaging system can produce three different types of images. A-Scan images (amplitude scan, also known as radio frequency) show the transmitted and reflected echoes in the time domain. In the A-Scan image the transducer is stationary, so the reflected echoes represent one line through the material, from the top to the bottom of the specimen. The image formed by moving the transducer in one direction is called B-Scan (brightness scan). The B-Scan image can be

10 interpreted as a collection of A-Scan images in one direction in the face of the specimen. The transducer can also move back and forth in the x-y plane producing a series of BScan images called C-Scan (constant depth scan). 2.3

Transform Coding, Subband Coding, and Data Compression

In the past decades the use of transform and subband coding has been introduced to compress and analyze data from many different sources. These applications range from medical imaging [Muy95] to earthquake detecting [Kie95]. The use of transform has been shown successful in the decomposition of the signal into transform structures, converting spatial signal coefficients to transform coefficients. Since this is a linear process and no information is lost, the number of coefficients produced is identical to the number of signal coefficients transformed. The transforms analyzed in this chapter are unitary, so the energy of the signal is conserved in the transformed domain. Furthermore, transforms provide decorrelation of the signal, and also increases the energy concentration of the transform domain. Figure 2.2 shows an orthogonal transform coding system. Orthogonality in this system implies that matrix T is orthogonal, so that T-1=TT. The decimators remove one coefficient of the input signal x(n) every M coefficients, while the expander introduces a zero every M samples. Perfect reconstruction ( x(n) = xˆ (n) ) is achieved if no quantization is applied to the down sampled signal.

11

x(n)

z-1 z-1

↓M ↓M

x(Mn)

y0(n)

x(Mn-1)

y1(n)

Q0 Q1

q0(n) q1(n)

↑M

x(Mn-M+1)

Decimators

z-1

T

-1

↓M

z-1

T

T z

↑M

z-1 yM-1(n)

QM-1

qM-1(n)

Quantizers

↑M

xˆ (n)

Expanders

Figure 2.2. Transform Coding System The purpose of subband coding is to separate the original signal into different frequency bands using digital filters. Figure 2.3 shows a subband coding system based on filter banks, where Hk(z) and Fk(z) are the Z transforms [Opp75] of the analysis and synthesis filters respectively. Since the filter bank is usually implemented as FIR filters, Hk(z) and Fk(z) are also the impulse response of the filters. Once the input signal x(n)

has been divided into separate frequency bands, different coding strategies can be used on the different bands in order to improve data compression performance. Since different frequency components affect the reconstructed signal differently, each frequency band may have a different quantizer. Therefore, a smaller number of bits can be used to encode the frequency components that are perceptually less important. After the quantization stage the signal qx(n) is up-sampled, and the output of the expanders is input to the synthesis bank. The outputs of all filters in the synthesis bank are added together to form the reconstructed signal xˆ (n) .

12 x(n)

H0(z) H1(z)

HM-1(z) Analysis bank

s0(n) s1(n)

sM-1(n)

↓M ↓M

↓M Decimators

v0(n) v1(n)

vM-1(n)

Q0 Q1

QM-1 Quantizers

q0(n) q1(n)

qM-1(n)

↑M

F0(z)

↑M

F1(z)

↑M

FM-1(z)

Expanders

Synthesis bank

xˆ (n)

Figure 2.3. Subband Coding System It’s important to realize that the optimal transform coding system should be redesigned as the statistics of the input signal changes [Vai93], while the optimal subband coding system is independent of the input signal statistics. Nevertheless, the desired effect in both transform and subband coding is that the energy of the signal will be contained in a few large coded coefficients. If this is true, data compression is achieved by reducing the number of coefficients (and ultimately the number of bits) necessary to represent the desired information. Lossless coding algorithms can further code the larger coefficients. The smaller coefficients can be coarsely coded or even deleted (in this case lossy coding) without damaging the quality of the reproduced signal. Different coding techniques differ in the degree of concentration of energy in a few coefficients, the region of influence of each coefficient in the reconstruction of the signal, and the appearance and visibility of coding noise due to coarse quantization of the coefficients. In this research, we analyze the data compression performance of the discrete cosine transform (DCT), discrete wavelet transform (DWT), and the Walsh-Hadamard transform (WHT). The information of the original signal (in our case ultrasonic signals) is concentrated in the transform domain. Depending on the characteristics of the transform, this concentration is going to be higher or lower. The data compression quality and performance is closely related to the energy concentration in the transform domain.

13 2.3.1

Discrete Cosine Transform

The derivation of the discrete cosine transform (DCT) with the Fourier cosine transform (FCT), defined as X c ( w) =

2

π

∞

∫ x(t ) cos(wt )dt

(2.1)

0

and the inverse FCT [Rao90] is given by x(t ) =

2

π

∞

∫X

c

( w) cos( wt )dw, t ≥ 0

(2.2)

0

We see that the FCT has a kernel given by Kc(w,t)=cos(wt). This kernel in the discrete domain is given by Kc(wm,tn), where wm=2πmΔf and tn=nΔt (m and n are integers), Δf and Δt are the unit sample intervals for the frequency and time respectively. Hence Kc(wm,tn) can be written as Kc(wm,tn)= Kc(2πmΔf,nΔt)=cos(2πmnΔfΔt)= Kc(m,n). This equation represents the discretized FCT kernel. If we define ΔfΔt =

1 , where N is 2N

the size of the input data vector, the DCT kernel can be written as ⎛ πmn ⎞ K c (m, n) = cos⎜ ⎟ ⎝ N ⎠

(2.3)

We can also regard the kernel Kc(m,n) as elements in an N+1 by N+1 matrix (denoted by [M]). The DCT was first reported by [Ahm74], where the forward transform is defined as

X ( m) =

N 2 ⎛ mnπ ⎞ km ∑ k n x(n) cos⎜ ⎟, N n =0 ⎝ N ⎠

and the inverse DCT is given by

m = 0,1,..., N

(2.4)

14

x ( n) =

N 2 ⎛ mnπ ⎞ k n ∑ k m X (m) cos⎜ ⎟, N n =0 ⎝ N ⎠

n = 0,1,..., N

(2.5)

where if j ≠ 0 or N ⎧1, ⎪ 1 kj = ⎨ , if j = 0 or N ⎪⎩ 2

(2.6)

If we define the kernel matrix [CN+1]m,n, such that the elements of the matrix are given by

[C N +1 ]m,n =

2 N

⎡ ⎛ mnπ ⎢k m k n cos⎜ N ⎝ ⎣

⎞⎤ ⎟⎥, ⎠⎦

m, n = 0,1,..., N

(2.7)

The forward and inverse DCT are defined as X = [C N ]x

(2.8)

x = [C N ] X

(2.9)

−1

Since the matrix [CN] has to be inverted, unitarity is desirable. It has been shown in [Rao90] that this matrix is indeed unitary, hence

[C N ]−1 = [C N ]T = [C N ]

(2.10)

simplifying the calculation of the inverse DCT matrix kernel. The discrete sine transform (DST) can be derived in a similar fashion. The kernel of the DST is given by

[S N −1 ]m,n = where

2 ⎛ mnπ ⎞ sin ⎜ ⎟, N ⎝ N ⎠

m, n = 1,2,..., N − 1

(2.11)

15 ⎧ 1, if n ≠ N ⎪ kn = ⎨ 1 , if n = N ⎪⎩ 2

(2.12)

The forward and inverse transforms for the DST family can be expressed as

X ( m) =

x ( n) =

N 2 ⎛ mnπ ⎞ k m ∑ k n x(n) sin⎜ ⎟, N n =0 ⎝ N ⎠

2 N ⎛ mnπ k n ∑ k m X (m) sin⎜ N n =0 ⎝ N

⎞ ⎟, ⎠

m = 0,1,..., N

(2.13)

n = 0,1,..., N

(2.14)

The DST kernel is also unitary, hence

[S N ]−1 = [S N ]T = [S N ] 2.3.2

(2.15)

Walsh-Hadamard Transform

The discrete Walsh-Hadamard transform (WHT) is one of the simplest transforms to implement in a digital system, since it’s composed only by 1’s and –1’s. The transform matrix is built from rearrangements of simpler transform matrices. The construction of WHT matrices with power of two dimensions is as follows ⎡H H 2N = ⎢ N ⎣H N

HN ⎤ − H N ⎥⎦

(2.16)

where HN is a N×N matrix designed in such a way that HHT=NI, I being a N×N identity matrix. The one and two dimensions WHT matrices are defined as H 1 = [1]

(2.17)

⎡1 1 ⎤ H2 = ⎢ ⎥ ⎣1 − 1⎦

(2.18)

16 The forward and inverse WHT are given by X = [H ]x

x=

2.3.3

(2.19)

1 [H ]X N

(2.20)

Discrete Wavelet Transform

The discrete wavelet transform (DWT) is implemented via filter banks. The filter bank representation of the DWT is shown in Figure 2.4. The low pass h0(n) and the high pass h1(n) filters are analysis filters, while g0(n) and g1(n) are low and high pass synthesis filters. The input signal x(n) is down-sampled by two (↓2) after the analysis filter, and up-sampled by two (↑2) before the synthesis filters. The inverse wavelet transform (IWT) reconstructs the signal from the up-sampled wavelet coefficients [Vai93, Dau90].

WT h0(n)

IWT 2

2

g0(n)

x(n)

+ h1(n)

2

2

xˆ (n)

g1(n)

Figure 2.4. Discrete Wavelet Transform

The DWT can be described in the Z domain as:

⎡ H ( z ) H 0 (− z )⎤ ⎡ X ( z ) ⎤ 1 Xˆ ( z ) = [G0 ( z ) G1 ( z )].⎢ 0 ⎥.⎢ ⎥ 2 ⎣ H 1 ( z ) H 1 (− z ) ⎦ ⎣ X (− z )⎦

(2.21)

17 X ( z) [G0 ( z ).H 0 ( z ) + G1 ( z ).H 1 ( z )] + X (− z ) [G0 ( z ).H 0 (− z ) + G1 ( z ).H 1 (− z )] Xˆ ( z ) = 2 2

(2.22)

Where G0(z), G1(z), H0(z), H1(z), X(z), and Xˆ ( z ) are the Z transforms of g0(z), g1(z), h0(z), h1(z), x(z), and xˆ (n) respectively. In order to have Xˆ ( z ) = c. X ( z ).z − n (such that 0

Xˆ ( z ) is a delayed and amplitude scaled version of X(z)) the following conditions for

perfect reconstruction have to be satisfied

⎧ G ( z ) = H 1 (− z ) G0 ( z ).H 0 (− z ) + G1 ( z ).H 1 (− z ) = 0 → ⎨ 0 ⎩G1 ( z ) = − H 0 (− z )

G0 ( z ).H 0 ( z ) + G1 ( z ).H1 ( z ) = c.z −n0 → H 0 ( z ) = H 1 (− z )

(2.23)

(2.24)

Many solutions for H0(z) and H1(z), satisfy the above conditions. One way to solve for H0(z) and H1(z) is to use quadrature mirror filters (QMF) [Vai93]. Perfect reconstruction using QMF is obtained if the filters are orthogonal to the even-shifted version of themselves M 2

∑ h ( n) h ( n + 2 k ) = δ ( k ) , k =0

i

i

i = 0,1

(2.25)

where M is the length of the QMF. The above equation sets the condition for perfect reconstruction. Although this is an important requirement in the solution of a wavelet kernel, there are other issues that affect the performance of the kernel in representing data. In [Car01] we investigated wavelet kernels that not only result in near perfect reconstruction, but that also improves signal data compression. There are many different wavelet kernels available in the literature [Say00]. The choice of which wavelet kernel to

18 use is application dependent, since each wavelet kernel has different characteristics that may improve or deteriorate signal analysis. 2.4

Conclusions

In this chapter, we reviewed fundamental concepts pertaining to the imaging, compression, and analysis of ultrasonic signals. Transform and subband coding techniques have been employed in the past for processing ultrasonic imaging signals. The use of such techniques has been shown successful in the decomposition of the signal into transform structures, converting spatial signal coefficients to transform coefficients. In the next chapter, we will introduce methods to improve data compression and denoising of ultrasonic signals using techniques to threshold the transform and subband coefficients.

19 3. CHAPTER III 3 THRESHOLDING TECHNIQUES FOR DENOISING AND COMPRESSING ULTRASONIC SIGNALS

3.1

Introduction

In many ultrasonic imaging applications, the signal acquired is embedded in noise, and situations with very small signal-to-noise ratio (SNR) are not uncommon. Thus, before any data analysis can be applied to the signal some level of noise removal is necessary. In this chapter we analyze the denoising and compression performances of the discrete wavelet transform (DWT), discrete cosine transform (DCT), and WalshHadamard transform (WHT). In section 3.2 we introduce a procedure for obtaining an adaptive thresholding function (ATF) that, when applied to DWT, DCT and WHT coefficients, improves the signal-to-noise (SNR) of ultrasonic signals embedded in noise. In particular, the ATF technique is very successful in denoising low SNR ultrasonic signals. Furthermore, the ATF approach outperforms the classical techniques when the ultrasonic signal has a low SNR (below 5 dB). The classical thresholding techniques analyzed (hard, soft, and Garrote thresholding) use a fixed threshold to separate noise from signal. On contrary, the ATF technique removes the noise from the signal by generating a thresholding function that is obtained from the statistical noise parameters. Hence, each transform coefficient is compared to its own “threshold”. Moreover, in section 3.3 the data compression performances of the DWT, DCT, and WHT are examined using simulated and experimental ultrasonic signals. Specifically, we

20 study the relation between the transforms’ kernel length, shape, energy compactness, and localization and their data compression performance. In this section we also study how the bandwidth of the ultrasonic echo (thus the shape of the echo) influences the energy packing of the transform coefficients. 3.2

Noise Suppression

Ultrasonic imaging signals are composed of ultrasonic echoes that are produced by the reflection of a transmitted pulse in scatterers inside the material under test. Often, the material under test has multiple structural scatterers. Hence, the ultrasonic signal is composed by echoes not only from the flaw in the material (in the case there is flaw in the material), but also by randomly distributed scatters. In addition to these multiple scatters, the received signal is also corrupted by thermal noise in the electronic components in the data acquisition system. Therefore, the goal of the denoising techniques presented in this chapter is improving the contribution of the echoes from flaws while worsening the noise contribution from multiple scatters and electronic systems. 3.2.1

Classical Thresholding Techniques

In this section we present three classical thresholding techniques: hard thresholding (HT), soft thresholding (ST) [Don92, Don93], and Garrote thresholding (GT) [Bre95, Gao98]. In the HT approach all transform coefficients smaller than a given threshold, T, are set to zero. All coefficients greater than T are kept at the same original value. In the ST approach, on the contrary, instead of introducing an abrupt change in the values of the transform coefficients, ST “smoothes” the transition from zero to non-zero coefficients. GT introduces a relation to the thresholding coefficients that emphasizes the larger

21 transform coefficients compared to the smaller transform coefficients. Hence, GT offers a compromise between HT and ST. The HT, ST, and GT can be described by the following equations

⎧⎪0, ⎪⎩ x,

ζ Hard ( x) = ⎨

x ≤T x >T

⎧0, ⎪ ζ Soft ( x) = ⎨ x − T , ⎪x + T , ⎩

ζ

Garrote

⎧0, ⎪ ( x) = ⎨ T 2 , ⎪x − x ⎩

(3.1)

x ≤T x >T x < −T

(3.2)

x ≤T x >T

(3.3)

where ζ represents the transform coefficients after thresholding, T is the threshold value, and x represents the original transform coefficients. Figure 3.1 shows the transform coefficients before and after thresholding.

22 0.5 a) Hard Threshold

0

−0.5 −0.5

−0.4

−0.3

−0.2

−0.1

0

0.1

0.2

0.3

0.4

0.5

−0.2

−0.1

0

0.1

0.2

0.3

0.4

0.5

−0.1

0

0.1

0.2

0.3

0.4

0.5

−T

0.5

T

b) Soft Threshold

0

−0.5 −0.5

−0.4

−0.3

0.5 c) Garrote Threshold

0

−0.5 −0.5

−0.4

−0.3

−0.2

Figure 3.1. Thresholding Techniques, Before (Dashed Lines) and After Thresholding (Solid Lines). a) Hard Threshold, b) Soft Threshold, and c) Garrote Threshold

23 In the case where the signal is corrupted by additive white Gaussian noise (AWGN) the observed signal, x, can be modeled as

x=s+n

(3.4)

where s is the original signal and n is the AWGN. To derive the best threshold T to denoise the signal x, the approach used by Donoho and Johnstone [Don92, Don93] is based on the minimization of the mean-square-error (MSE) between the original, s , and the reconstructed signal, sˆ , MSE ( sˆ, s ) =

1 N

N

∑ E (sˆ

i

i =1

− si )

2

(3.5)

The minimization of MSE in the previous equation leads to the “universal” threshold, T, given by

T=

σ N

2 ln( N − 1)

(3.6)

where N is the total number of transform coefficients and σ is the standard deviation of the AWGN. The statistical parameters of the AWGN can be estimated from the original signal using the unbiased estimation [Kay93a]

σ=

1 N (xi − μ )2 ∑ N − 1 i =1

(3.7)

where the mean, μ, is estimated by

1 μ= N

N

∑x i =1

i

(3.8)

24 3.2.2

Adaptive Thresholding Function

In this section, we present the procedure for estimating the adaptive thresholding function (ATF) in order to improve SNR of ultrasonic signals embedded in noise. Different from the classical thresholding techniques shown in the previous section, the ATF technique relies on a thresholding function, T(x), instead of a thresholding coefficient, T. The shape of this thresholding function is determined by the source of random noise in the signal. More specifically, T(x) has the shape of the probability distribution function of the noise source. Thus, instead of modifying each transform coefficient based on a single threshold coefficient, in the ATF technique a thresholding function is subtracted from the sorted transform coefficients, or

ζ sATF ( xs ) = xs − T ( xs )

(3.9)

where ζ sATF ( xs ) represents the transform coefficients after ATF thresholding, T ( xs ) is the ATF, and xs is the sequence of sorted coefficients of x. This result is then resorted and inverse transformed. The following sections show that the ATF technique leads to a reconstructed signal in the time domain with higher SNR than using the classical thresholding techniques. The performance of the ATF and classical thresholding techniques in denoising ultrasonic signals using different transforms (DWT, DCT, and WHT) is also the subject of the following section. 3.2.3

Description of the ATF Technique

The first step in the implantation of the ATF technique is to apply the appropriate transform to the original signal, x. The next step in the algorithm is to sort the transform coefficients to obtain xs, as shown in the flowchart in Figure 3.2. The sorting algorithm organizes the transform coefficients from smallest to greatest in amplitude while keeping

25 the position (or order) information of each of these coefficients. To better illustrate the sorting algorithm, Table 3.1 shows a set of 10 unsorted coefficients. Table 3.1. Coefficients Sorted by Position Coefficient Position

0

1

2

3

4

5

6

7

8

9

Coefficient Amplitude

9

2

4

3

2

7

1

8

0

-2

The result of running these 10 coefficients through the sorting algorithm is shown in Table 3.2. Table 3.2. Coefficients Sorted by Amplitude Coefficient Position

9

8

6

1

4

3

2

5

7

0

Coefficient Amplitude

-2

0

1

2

2

3

4

7

8

9

26

Ultrasonic Signal

Apply Transform DCT, DWT, or WHT to obtain x

Sort Transform Coefficients to obtain xs

Estimate Noise Statistical Parameters (µ and σ)

Generate Thresholding Function T(xs)

Subtract Thresholding Function from Original Transform Coefficients

ζ sATF ( xs ) = xs − T ( xs )

Resort Transform Coefficients

Apply Inverse Transform IDCT, IDWT, or IWHT

Estimated Ultrasonic Echoes Figure 3.2. Flowchart of ATF Algorithm

27 The recovery of the original sequence is done using the same sorting algorithm, but applied to the coefficient’s position instead of the coefficient’s amplitude. The result of the resorting process is shown in Table 3.2. In the ATF technique the transform coefficients are sorted such that the observed signal, xs, is given by

x s = s s + ns

(3.10)

where ss and ns are the sorted sequence of signal and noise transform coefficients, respectively. The shape of ss is unknown, as the number, position, and distribution of the ultrasonic echoes in s is by definition unknown. On the contrary, if the type of the underlying random process that generates ns is known, as we assume it is, the shape of ns is the probability distribution function of n [Dav81, Fuk90]. Hence, the ATF is generated from the noise parameters so that ns and T(xs) are similar to each other, and as a consequence

ζ sATF ( xs ) = ss + ns − T ( xs ) ζ sATF ( xs ) ≅ ss

(3.11)

or the thresholded transform coefficients are approximately equal to the transform coefficients of the original signal ss. Two examples of the outcome of sorting the realizations of random processes are shown in Figure 3.3 and Figure 3.4. Figure 3.3a shows 2048 samples of a Gaussian process with zero mean and unit variance (µ = 0 and σ = 1, i.e., N(0,1)). Figure 3.3b shows the sorted sequence of the samples in Figure 3.3a, where the coefficient amplitude is plotted in the horizontal axis and the coefficient number is plotted in the vertical axis. This figure shows that the sorted coefficients in

28 Figure 3.3a produce the probability distribution function (PDF) of the Gaussian random process, which has the shape of the error function. Similarly, Figure 3.4a shows 2048 samples of a uniform process with µ = 0 and σ2 = 0.1. Figure 3.4b shows the sorted sequence of the samples in Figure 3.4a, where the coefficient amplitude is plotted in the horizontal axis and the coefficient number is plotted in the vertical axis. Similarly, this figure shows that sorting the DWT coefficients in Figure 3.4b lead to the PDF of the realization in Figure 3.4a.

29

4

Coefficient Amplitude

a) 2

0

−2

−4

200

400

600

800 1000 1200 Coefficient Number

1400

1600

1800

2000

2000

Coefficient Number

b) 1500 1000 500

−4

−3

−2

−1

0 1 Coefficient Amplitude

2

3

4

Figure 3.3. a) 2048 Samples of Gaussian Process Realization, b) Sorted Samples of the Random Noise Realization in a)

30

Coefficient Amplitude

0.6

a)

0.4 0.2 0 −0.2 −0.4 −0.6 200

400

600

800 1000 1200 Coefficient Number

1400

1600

1800

2000

Coefficient Number

2000

b) 1500 1000 500

−0.6

−0.4

−0.2

0 0.2 Coefficient Amplitude

0.4

0.6

Figure 3.4. a) 2048 Samples of Uniform Process Realization, b) Sorted Samples of the Random Noise Realization in a)

31 The noise statistical parameters can be estimated from the sorted transform coefficients, xs, or the original transform coefficients, x. The noise estimated parameters (mean, µ, and standard deviation, σ) are used to generate the ATF, T(xs), given by T (x s )z =

⎡ (k − μ )2 ⎤ exp ∫ ⎢⎣− 2σ 2 ⎥⎦ dk 2π 0 1

z

(3.12)

in the case that the noise source is Gaussian (a common source of noise in experimental ultrasonic signals), and z represents the individual samples of T(xs) and it is defined in the range from the minimum coefficient of x (i.e., xs(1)) to the maximum coefficient of x (i.e., xs(N)) in

xs ( N ) − xs (1) steps. N

In the case that the noise source is uniform, the ATF is given by T ( x s )z =

where μ =

z − x s (1) 1 z−μ = + x s (N ) − x s (1) 2 2 3σ

(3.13)

x (N ) − x s (1) x s (1) + x s (N ) and σ = s . 2 2 3

The next step is to subtract T(xs) from the original signal, xs, to obtain the thresholded transform coefficients ζ sATF . The thresholded coefficients are then resorted using the original position of the transform coefficients, and the last step of the algorithm is to apply the inverse transform to the resorted transform coefficients. 3.2.4

Performance of ATF and Classical Techniques

In this section the performance of the ATF technique is compared to the performance of the classical techniques, i.e., hard thresholding (HT), soft thresholding (ST), and Garrote thresholding (GT) in denoising ultrasonic signals with single and multiple echoes

32 using the DCT, DWT, and WHT. In these simulations the pulse-echo ultrasonic testing the backscattered echo from a single reflector can be modeled as

[

]

f Θ (t ) = β exp − α (t − τ )2 cos(2πf c (t − τ ) + φ )

(3.14)

where Θ = [α , β , f c , φ , τ ] denotes the parameter vector and ultrasonic signals consisting of multiple interfering echoes can be modeled as s Θ (t ) =

M −1

∑ j =0

f Θ j (t ) =

M −1

∑β j =0

j

[

] (

exp − α j (t − τ j ) 2 cos 2πf c j (t − τ j ) + φ j

)

(3.15)

where M is the total number of echoes in the signal. A Monte Carlo analysis was done to evaluate the performance of the thresholding techniques when the applied to single and multiple interfering echoes in the DCT, DWT, and WHT domains. The SNR of the original (input) signal was varied from around -5dB to 10dB. The MC analysis generated 25 noise realizations for each input SNR, thus in these figures, each point in the plot shows the mean of the improvement obtained and the error bars show one standard deviation around this mean. Figure 3.5-Figure 3.10 show the performance of the ATF, HT, ST, and GT techniques when applied to a single and multiple ultrasonic echoes embedded in uniform random noise, while Figure 3.11-Figure 3.16 show the performance of these thresholding techniques with ultrasonic signals embedded in AWGN. The single echo used in these simulations has the following parameters: bandwidth factor

α = 4 (MHz )2 , arrival time τ = 1 μs , center frequency f c = 4.5 MHz , phase φ = 1 rad , and amplitude β = 1 . The multiple interfering echoes used have the parameters shown in Table 3.3.

33 Table 3.3. Original Parameters of Ten Interfering Simulated Echoes. Echo a b c d e f g h i j

f [MHz] 4.5 4 4 7 3 5 3.3 3 3 3

τ [μs]

α [MHz]2

β

φ [rad]

1 6 3 9 2 8.7 5 7 8 4

4 4.8 4 8 5 8 5 5 4 4

1 0.9 1 1 0.9 0.8 0.8 1 1 0.5

1 0.8 0.9 1.5 0.8 0.2 0 0.7 0.5 0.7

These plots show that the best transform to recover the signal embedded in low SNR is the DWT using the ATF. This result is a consequence of the time localization of the DWT kernel and its packing efficiency, as seen in more detail in the following section. Furthermore, these results show that for small SNR (-2dB) the ATF consistently outperforms the classical thresholding techniques, but the situation reverses as the input SNR is improved (15dB). As expected, the ATF is better suited for applications with poor SNR. To better appreciate the performance of the ATF and the classical thresholding techniques, Figure 3.17-Figure 3.19 show one realization of the multiple interfering echo signal embedded in uniform noise with SNR = -2.5dB, while Figure 3.20-Figure 3.22 show one realization of the multiple interfering echo signal embedded in AWGN with SNR = 0dB.

34

15

Output SNR [dB]

10

5

0

ATF HT ST GT

−5 −5

0 5 Input SNR [dB]

10

Figure 3.5. SNR Enhancement of Adaptive Thresholding Function (ATF), Hard Thresholding (HT), Soft Thresholding (ST), and Garrote Thresholding (GT) techniques with DCT Applied to a Single Ultrasonic Echo Corrupted with Uniform Noise

35

15

Output SNR [dB]

10

5

ATF HT ST GT

0

−5 −5

0 5 Input SNR [dB]

10

Figure 3.6. SNR Enhancement of Adaptive Thresholding Function (ATF), Hard Thresholding (HT), Soft Thresholding (ST), and Garrote Thresholding (GT) techniques with DWT Applied to a Single Ultrasonic Echo Corrupted with Uniform Noise

36

Output SNR [dB]

10

5 ATF HT ST GT

0

−5 −5

0 5 Input SNR [dB]

10

Figure 3.7. SNR Enhancement of Adaptive Thresholding Function (ATF), Hard Thresholding (HT), Soft Thresholding (ST), and Garrote Thresholding (GT) techniques with WHT Applied to a Single Ultrasonic Echo Corrupted with Uniform Noise

37

2

Output SNR [dB]

0 −2

ATF HT ST GT

−4 −6 −8 −5

0 5 Input SNR [dB]

10

Figure 3.8. SNR Enhancement of Adaptive Thresholding Function (ATF), Hard Thresholding (HT), Soft Thresholding (ST), and Garrote Thresholding (GT) techniques with DCT Applied to 10 Ultrasonic Echoes Corrupted with Uniform Noise

38

Output SNR [dB]

15

10

5 ATF HT ST GT

0

−2

0

2

4 6 8 Input SNR [dB]

10

12

Figure 3.9. SNR Enhancement of Adaptive Thresholding Function (ATF), Hard Thresholding (HT), Soft Thresholding (ST), and Garrote Thresholding (GT) techniques with DWT Applied to 10 Ultrasonic Echoes Corrupted with Uniform Noise

39

14 12

Output SNR [dB]

10 8 6 4

ATF HT ST GT

2 0 −2 −2

0

2

4 6 8 Input SNR [dB]

10

12

Figure 3.10. SNR Enhancement of Adaptive Thresholding Function (ATF), Hard Thresholding (HT), Soft Thresholding (ST), and Garrote Thresholding (GT) techniques with WHT Applied to 10 Ultrasonic Echoes Corrupted with Uniform Noise

40

8

Output SNR [dB]

6 4 2 0 ATF HT ST GT

−2 −4 −6 −8 −10 −10

−5

0 Input SNR [dB]

5

Figure 3.11. SNR Enhancement of Adaptive Thresholding Function (ATF), Hard Thresholding (HT), Soft Thresholding (ST), and Garrote Thresholding (GT) techniques with DCT Applied to a Single Ultrasonic Echo Corrupted with Gaussian Noise

41

Output SNR [dB]

10

5

0

ATF HT ST GT

−5

−10 −10

−5

0 Input SNR [dB]

5

Figure 3.12. SNR Enhancement of Adaptive Thresholding Function (ATF), Hard Thresholding (HT), Soft Thresholding (ST), and Garrote Thresholding (GT) techniques with DWT Applied to a Single Ultrasonic Echo Corrupted with Gaussian Noise

42

4

Output SNR [dB]

2 0 −2 −4 ATF HT ST GT

−6 −8 −10 −10

−5

0 Input SNR [dB]

5

Figure 3.13. SNR Enhancement of Adaptive Thresholding Function (ATF), Hard Thresholding (HT), Soft Thresholding (ST), and Garrote Thresholding (GT) techniques with WHT Applied to a Single Ultrasonic Echo Corrupted with Gaussian Noise

43

2

Output SNR [dB]

0 −2 −4

ATF HT ST GT

−6 −8 −5

0 5 Input SNR [dB]

10

Figure 3.14. SNR Enhancement of Adaptive Thresholding Function (ATF), Hard Thresholding (HT), Soft Thresholding (ST), and Garrote Thresholding (GT) techniques with DCT Applied to 10 Ultrasonic Echoes Corrupted with Gaussian Noise

44

15

Output SNR [dB]

10

5 ATF HT ST GT

0

−5 −5

0 5 Input SNR [dB]

10

Figure 3.15. SNR Enhancement of Adaptive Thresholding Function (ATF), Hard Thresholding (HT), Soft Thresholding (ST), and Garrote Thresholding (GT) techniques with DWT Applied to 10 Ultrasonic Echoes Corrupted with Gaussian Noise

45

6

Output SNR [dB]

4 2 0 −2

ATF HT ST GT

−4 −6 −6

−4

−2 0 Input SNR [dB]

2

4

6

Figure 3.16. SNR Enhancement of Adaptive Thresholding Function (ATF), Hard Thresholding (HT), Soft Thresholding (ST), and Garrote Thresholding (GT) Techniques with WHT Applied to 10 Ultrasonic Echoes Corrupted with Gaussian Noise

46

2 a) Original Signal 0 −2 0 2

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1 −5

x 10

b) Noisy Signal 0 −2 0 2

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1 −5

x 10

c) ATF 0 −2 0 2

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1 −5

x 10

d) HT 0 −2 0 2

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1 −5

x 10

e) ST 0 −2 0 2

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1 −5

x 10

f) GT 0 −2

0

0.1

0.2

0.3

0.4

0.5

Time [μs]

0.6

0.7

0.8

0.9

1 −5

x 10

Figure 3.17. SNR Enhancement of Thresholding Techniques in an Ultrasonic Signal with 10 Interfering Echoes Corrupted by Uniform Noise with SNR = -2.5dB using DCT. a) Original Signal, b) Noisy Signal, Reconstructed Signal Using c) Adaptive Thresholding Function (ATF), d) Hard Thresholding (HT), e) Soft Thresholding (ST), and f) Garrote Thresholding (GT)

47

2 a) Original Signal 0 −2 0 2

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1 −5

x 10

b) Noisy Signal 0 −2 0 2

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1 −5

x 10

c) ATF 0 −2 0 2

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1 −5

x 10

d) HT 0 −2 0 2

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1 −5

x 10

e) ST 0 −2 0 2

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1 −5

x 10

f) GT 0 −2

0

0.1

0.2

0.3

0.4

0.5

Time [μs]

0.6

0.7

0.8

0.9

1 −5

x 10

Figure 3.18. SNR Enhancement of Thresholding Techniques in an Ultrasonic Signal with 10 Interfering Echoes Corrupted by Uniform Noise with SNR = -2.5dB using DWT. a) Original Signal, b) Noisy Signal, Reconstructed Signal Using c) Adaptive Thresholding Function (ATF), d) Hard Thresholding (HT), e) Soft Thresholding (ST), and f) Garrote Thresholding (GT)

48

2 a) Original Signal 0 −2 0 2

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1 −5

x 10

b) Noisy Signal 0 −2 0 2

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1 −5

x 10

c) ATF 0 −2 0 2

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1 −5

x 10

d) HT 0 −2 0 2

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1 −5

x 10

e) ST 0 −2 0 2

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1 −5

x 10

f) GT 0 −2

0

0.1

0.2

0.3

0.4

0.5

Time [μs]

0.6

0.7

0.8

0.9

1 −5

x 10

Figure 3.19. SNR Enhancement of Thresholding Techniques in an Ultrasonic Signal with 10 Interfering Echoes Corrupted by Uniform Noise with SNR = -2.5dB using WHT. a) Original Signal, b) Noisy Signal, Reconstructed Signal Using c) Adaptive Thresholding Function (ATF), d) Hard Thresholding (HT), e) Soft Thresholding (ST), and f) Garrote Thresholding (GT)

49

2 a) Original Signal 0 −2 0 2

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1 −5

x 10

b) Noisy Signal 0 −2 0 2

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1 −5

x 10

c) ATF 0 −2 0 2

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1 −5

x 10

d) HT 0 −2 0 2

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1 −5

x 10

e) ST 0 −2 0 2

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1 −5

x 10

f) GT 0 −2

0

0.1

0.2

0.3

0.4

0.5

Time [μs]

0.6

0.7

0.8

0.9

1 −5

x 10

Figure 3.20. SNR Enhancement of Thresholding Techniques in an Ultrasonic Signal with 10 Interfering Echoes Corrupted by Gaussian Noise with SNR = 0dB using DCT. a) Original Signal, b) Noisy Signal, Reconstructed Signal Using c) Adaptive Thresholding Function (ATF), d) Hard Thresholding (HT), e) Soft Thresholding (ST), and f) Garrote Thresholding (GT)

50

2 a) Original Signal 0 −2 0 2

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1 −5

x 10

b) Noisy Signal 0 −2 0 2

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1 −5

x 10

c) ATF 0 −2 0 2

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1 −5

x 10

d) HT 0 −2 0 2

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1 −5

x 10

e) ST 0 −2 0 2

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1 −5

x 10

f) GT 0 −2

0

0.1

0.2

0.3

0.4

0.5

Time [μs]

0.6

0.7

0.8

0.9

1 −5

x 10

Figure 3.21. SNR Enhancement of Thresholding Techniques in an Ultrasonic Signal with 10 Interfering Echoes Corrupted by Gaussian Noise with SNR = 0dB using DWT. a) Original Signal, b) Noisy Signal, Reconstructed Signal Using c) Adaptive Thresholding Function (ATF), d) Hard Thresholding (HT), e) Soft Thresholding (ST), and f) Garrote Thresholding (GT)

51

2 a) Original Signal 0 −2 0 2

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1 −5

x 10

b) Noisy Signal 0 −2 0 2

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1 −5

x 10

c) ATF 0 −2 0 2

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1 −5

x 10

d) HT 0 −2 0 2

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1 −5

x 10

e) ST 0 −2 0 2

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1 −5

x 10

f) GT 0 −2

0

0.1

0.2

0.3

0.4

0.5

Time [μs]

0.6

0.7

0.8

0.9

1 −5

x 10

Figure 3.22. SNR Enhancement of Thresholding Techniques in an Ultrasonic Signal with 10 Interfering Echoes Corrupted by Gaussian Noise with SNR = 0dB using WHT. a) Original Signal, b) Noisy Signal, Reconstructed Signal Using c) Adaptive Thresholding Function (ATF), d) Hard Thresholding (HT), e) Soft Thresholding (ST), and f) Garrote Thresholding (GT)

52 3.3

Data Compression

Data compression is the process of obtaining a more efficient representation of a signal; consequently, it is desirable to use data compression techniques to reduce the ultrasonic data size while maintaining the signal integrity. Discrete transforms are used in data compression due to their data decorrelation, energy compaction, data independence, and computational speed properties. Several discrete transforms achieve these properties to various extents, such as the Karhunen-Loeve transform (KLT), the discrete cosine transform (DCT) [Ahm74], the Walsh-Hadamard transform (WHT) [Ahm71], and the discrete wavelet transform (DWT) [Dau92]. The KLT is optimal in relation to data decorrelation and energy compaction, which leads to optimal results in data compression. However, the data dependence of the KLT kernel deteriorates the performance of the KLT computation to the point of making it of limited practical interest. In contrast, the DCT, the DWT, and the WHT are data independent and can be implemented efficiently. The DCT has been used in the compression of images as part of the joint photographic experts group (JPEG) standard. More recently, JPEG has adopted the DWT as the basis of signal decomposition in the JPEG2000 standard. The DCT and the DWT have been applied to compress and decompose ultrasonic images [Muy95, Chi01, Cin01, Pun01, Oh03] and also to detect ultrasonic flaw echoes [Rob97, Car01]. In this section the compression properties of the DCT, the WHT, and the DWT for ultrasonic imaging are examined. In particular, we evaluate the relation between ultrasonic echo shape, bandwidth, and data compression performance.

53 3.3.1

Discrete Transforms

The data compression performance is closely related to how the energy of the transform is distributed among its coefficients (i.e., energy compaction); hence, the unitarity of the WHT, DCT, and DWT is an important property for determining their data compression performance. If the few most energetic (highest amplitude) coefficients of the transform concentrate most of the energy, a higher rate of data compression is achieved. In contrast, if the few most energetic coefficients of the transform do not have a significant portion of the energy, then a lower data compression performance is achieved. In this section, a hard thresholding technique is applied to the transform domain coefficients with the purpose of achieving data reduction. In the hard threshold technique all coefficients above the threshold are kept, while the coefficients below the threshold are set to zero. The DCT is widely used for signal compression and it is an excellent processing tool to perform signal analysis.

It has near optimal performance in terms of energy

compaction capability [Ahm71, Car01]. The Walsh-Hadamard transform (WHT) is one of the simplest transforms to be implemented, as the WHT is a unitary and orthogonal transform [Ahm71] composed by rectangular waveforms with values +1 and –1. Hence, no multiplications are necessary, and the transform is faster to implement than the DCT and DWT. Furthermore, the WHT can achieve both energy compaction and fast implementation, as it resembles the DCT kernel closely in terms of zero-crossing. The DWT is a unitary and orthonormal transform that is calculated by shifting and dilating the wavelet kernel and correlating it with the input signal. An important advantage of the DWT over the DCT and WHT is that the wavelet kernel is not restricted to a set of

54 functions (e.g., cosines in the DCT case). In principle, the number of different wavelet kernels is unlimited [Dau92]. 3.3.2

Data Compression Performance Analysis

In this section the data compression performance of the DWT (Daubechies with 20 coefficients, i.e., Daub20, and Haar), DCT, and WHT is analyzed. For a Gaussian echo, the 98% bandwidth (BW) that contains 98% of the signal energy is BW = 0.382 α . Normalizing this BW by the center frequency, fc, NBW = 0.382 α f c , allows the differentiation between narrowband and broadband echoes. Hence, a narrowband echo has a small NBW while a broadband echo has a large NBW. Moreover, the NBW is inversely proportional to the quality factor. The transform coefficients for a narrowband (NBW= 0.2) and a broadband (NBW= 0.5) ultrasonic echoes are shown in Figure 3.23. In this figure the transform coefficients have been ordered from largest to smallest energy, so that the first coefficient plotted is the most energetic one, and the last coefficient plotted is the least energetic. The echo energy is normalized to 1. Therefore, for the 2 narrowband echo the bandwidth factor α = 1.2 (MHz ) , time-of-arrival τ = 7 μs , center

frequency

f c = 4 MHz , phase φ = 0 rad , and amplitude β = 0.18 , and for the

broadband echo the bandwidth factor α = 27 (MHz ) , time-of-arrival τ = 7 μs , center 2

frequency fc = 4 MHz , phase φ = 0 rad , and amplitude β = 0.4 . Figure 3.23 shows that the DWT performs better for broadband echoes, while the WHT and DCT outperform the DWT for narrowband echoes. These results come from the fact that all transforms correlate their kernels to the ultrasonic echo. The DWT has a kernel localized in time; therefore it performs better when the echo has a short duration (broadband). On the

55 contrary, if the echo has a longer time duration (narrowband), the WHT and DCT kernels will show a higher correlation. The effect of this property is illustrated in Figure 3.24. In this figure the original echo is superimposed to the reconstructed signal using only the most energetic transform coefficient. The broadband echo can almost be completely reconstructed using one DWT coefficient (Figure 3.25b and Figure 3.25d). On the contrary, more DWT coefficients are necessary to represent a narrowband signal (Figure 3.24b and Figure 3.24d). The DCT and WHT kernels have higher amplitude (i.e., higher energy) if the signal is narrowband (Figure 3.24a and Figure 3.24c), and lower energy if the signal is broadband (Figure 3.25a and Figure 3.25c). Since all the analyzed transforms are unitary, a high-energy concentration of the transform coefficients means that a better data compression performance is achieved. The compression performance of the DCT, DWT, and WHT as a function of the ultrasonic echo bandwidth (i.e., NBW) is shown in Figure 3.26. This figure shows the total energy of the five most energetic transform coefficients. All signals are 512 16-bits samples long. For a broadband signal, the DWT Daub20 outperforms the DCT and WHT, as the DWT coefficients are able to recover over 90% of the signal energy. The DCT and the WHT outperform the DWT for narrowband signals. The time domain representation of reconstructed ultrasonic echoes is shown in Figure 3.27. This figure shows the reconstructed signals when the 25 most dominant coefficients (95% compression ratio) are used in the inverse transforms for a broadband (NBW=0.5) and a narrowband (NBW=0.1) signal. The signal-to-noise ratio (SNR) of the reconstructed signals is presented in Table 3.4.

Coefficient energy

56

0.6

a)

0.4 0.2 0 1

Coefficient energy

DWT Haar DWT Daub20 WHT DCT

0.6

2

3

4

5

6

7

8

9

10

DWT Haar DWT Daub20 WHT DCT

b)

0.4 0.2 0 1

2

3

4

5 6 7 Transform coefficient

8

9

10

Figure 3.23. Ordered Transform Coefficients of Ultrasonic Echoes: a) Narrowband and b) Broadband

57

0.2 a)

0.2 b)

0.1

0.1

0

0

−0.1

−0.1

−0.2

−0.2 6

7 Time [s]

8

9

6

−6

x 10

0.2 c)

0.2 d)

0.1

0.1

0

0

−0.1

−0.1

−0.2

−0.2 6

7 Time [s]

8

9 −6

x 10

6

7 Time [s]

8

7 Time [s]

8

9 −6

x 10

9 −6

x 10

Figure 3.24. Narrowband Original Echo (Dashed Line) Superimposed to Most Energetic Coefficient of a) DCT, b) DWT Haar, c) WHT, and d) DWT Daub20

58

a)

b)

c)

d)

Figure 3.25. Broadband Original Echo (Dashed Line) Superimposed to Most Energetic Coefficient of a) DCT, b) DWT Haar, c) WHT, and d) DWT Daub20

59

1 0.9 0.8

Energy

0.7 WHT DCT DWT Daub20 DWT Haar

0.6 0.5 0.4 0.3 0.2 0.1

0.1

0.2

0.3 NBW

0.4

0.5

Figure 3.26. Relation Between NBW and the Five Most Energetic Transform Coefficients

60 0.5

a)

0 6 b)

7 Time [s]

6 c)

7 Time [s]

0 6 d)

7 Time [s]

0

0.2

g)

7 Time [s]

8

9 −6

x 10

DCT

0 8

9 −6

x 10 DWT Haar

6 0.2

h)

7 Time [s]

8

9 −6

DWT Haarx 10

0

8

9 −6

x 10 DWT Daub

6 0.2

i)

7 Time [s]

8

9 −6

x 10 DWT Daub

0 −0.2

−0.5

6 e)

7 Time [s]

8

9 −6

WHT

0 −0.5

x 10

6

−0.2

−0.5

0.5

9 −6

−0.2

−0.5

0.5

0

8 DCT

0

0.5

Narrowband

f)

−0.2

−0.5 0.5

0.2

Broadband

x 10

6 0.2

j)

7 Time [s]

8 WHT

9 −6

x 10

0 −0.2

6

7 Time [s]

8

9 −6

x 10

6

7 Time [s]

8

9 −6

x 10

Figure 3.27. a) Broadband signal, b) ultrasonic signal compressed 95% using DCT, c) DWT (Haar), d) DWT (Daubechies), e) WHT; f) Broadband signal, g) ultrasonic signal compressed 95% using DCT, h) DWT (Haar), i) DWT (Daubechies), j) WHT

61 Table 3.4. SNR of Broadband and Narrowband Signals Compressed by 95%.

DWT Daub20 DWT Haar DCT WHT

Broadband SNR [dB] 34 31 6 3

Narrowband SNR [dB] 22 8 50 5

Figure 3.28 shows the coefficients for six different wavelet kernels, Haar, Daubechies, Beylkin, Coiflet, Symmlet, and Vaidyanathan. The data compression performance of these six wavelet kernels is shown in Figure 3.29. This figure shows how much energy is concentrated in the five most dominant coefficients of the DWT as a function of the bandwidth of the ultrasonic signal. These results show that the Daub20 wavelet kernel has the best data compression performance, while the Haar wavelet kernel has the worst data compression performance. The size of each kernel and the energy compaction of the five most dominant coefficients of the different wavelet kernels are shown in Table 3.5. These values have been normalized by the best performing kernel; hence the Daub20 has energy compaction equal to 1.

62

2

3

4

5

Beylkin

−1 1

5

10

5

10 NBW

15

5

10

15

10

15

d)

0

1

e)

Symmlet

Daubechies

Coiflet

−1

15

0 −1

0

1

0

b)

−1

Energy

Energy

Haar

0 −1 1 1 c)

Energy

1

Energy

a)

Energy

Energy

1

5 f)

0 −1

Vaidyanathan 5

10 15 NBW

20

Figure 3.28. Wavelet Kernels: a) Haar, b) Daubechies 20, c) Beylkin, d) Coiflet, e) Symmlet, and f) Vaidyanathan

63

Figure 3.29. Energy Accumulated Among the Five Most Dominant DWT Coefficients Using the Following Kernels: a) Haar, b) Daubechies, c) Beylkin, d) Coiflet, e) Symmlet, and f) Vaidyanathan

64 Table 3.5. Kernel Size and Energy Compaction Wavelet kernel Haar Daubechies Beylkin Coiflet Symmlet Vaidyanathan

Kernel size 2 20 18 24 20 24

Energy compaction 0.8 1 0.95 0.97 0.98 0.94

The flaw detection performance of the DWT, WHT, and DCT is also evaluated using an ultrasonic experimental signal consisting of multiple interfering echoes. Figure 3.30a and Figure 3.31a show an experimental signal acquired from a steel sample block with a flat-bottom hole using a 5MHz transducer and sampling rate of 100MHz. The experimental signal has poor SNR and the flaw echo contains microstructure scattering and measurement noise. Figure 3.30b-e shows the reconstructed signal using 10 transform coefficients (i.e., 98% compression ratio), while Figure 3.31b-e shows the reconstructed signal when 41 transform coefficients are used. (i.e., 92% compression ratio). This figure shows that the dominant flaw echo has lower frequency content when compared to the backscattered grain echoes. It also shows that the DWT Haar is the best performing transform for target detection, followed by the DWT Daub20. For a compression ratio of 98% the flaw echo was not detected by the DCT or the WHT, while as the number of transform coefficients is increased to 41 both DCT and WHT are able to reconstruct the flaw echo. This result reveals that the DWT, due to its time-limited kernel, outperforms the DCT and WHT in detecting flaw echoes in the presence of microstructure scattering signals.

Amplitude

Amplitude

Amplitude

Amplitude

Amplitude

65

0.4

a) Original signal

0.2

Flaw

0 −0.2 0 0.4 0.2

0.5

1

1.5

2

Flaw

b) DCT

2.5 Time [s]

3

2.5 Time [s]

3

2.5 Time [s]

3

2.5 Time [s]

3

2.5 Time [s]

3

3.5

4

4.5 −6

x 10

0 −0.2 0 0.4 0.2

0.5

1

1.5

c) DWT Haar

2

Flaw

3.5

4

4.5 −6

x 10

0 −0.2 0 0.4 0.2

0.5

1

1.5

2

Flaw

d) DWT Daub

3.5

4

4.5 −6

x 10

0 −0.2 0 0.4 0.2

0.5

1

1.5

e) WHT

2

Flaw

3.5

4

4.5 −6

x 10

0 −0.2 0

0.5

1

1.5

2

3.5

4

4.5

Figure 3.30. a) Original Experimental Signal, and Reconstruction with 98% Compression Ratio: b) DCT, c) DWT Haar, d) DWT Daub, and e) WHT

−6

x 10

66

a)

b)

c)

d)

e)

Figure 3.31. a) Original Experimental Signal, and Reconstruction with 92% Compression Ratio: b) DCT, c) DWT Haar, d) DWT Daub, and e) WHT

67 3.4

Conclusions

In this chapter we have introduced an adaptive thresholding function technique that uses the statistical parameters of the noise embedded in the signal to generate a thresholding function based on the probability distribution function of the noise. This thresholding function is then subtracted from the sorted transform coefficients, leading to transform coefficients with a superior SNR. The results presented show that this is a very powerful technique that allows the detection of ultrasonic backscattered echoes embedded in low SNR environments. For signals with uniform noise added to DWT coefficients, the ATF technique achieves SNR improvements around 9dB over the classical thresholding techniques; these improvements are above 10dB for Gaussian noise. This chapter also compares the denoising performance of the ATF technique when implemented in combination with DCT, DWT, and WHT coefficients. The results show that the DWT has the best performance among the transforms in the group analyzed. The performance of the DWT is consistent with the results shown in the previous chapters. The energy packing capabilities of the DWT, along with the time localization of the wavelet kernel, allow that only a few DWT coefficients are necessary to correctly identify the location of the ultrasonic echo. Such localization is paramount for signal denoising and compression. In this chapter we have also analyzed ultrasonic signal compression using the DWT, the DCT, and the WHT. The results obtained show that the DWT is better in the representation of broadband signals, while the DCT and the WHT are more suitable in the representation of narrowband signals. For a narrowband ultrasonic signal the five

68 most energetic coefficients of the DCT accumulate 95% of the total signal energy, while the DWT and the WHT accumulate 50% and 60% of the total signal energy respectively. On the contrary, for a broadband ultrasonic signal the 5 most dominant DCT coefficients accumulate 30% of the total signal energy, while the DWT and the WHT accumulate 97% and 15% of the total signal energy respectively. Thus, the discrete transform methods analyzed in this chapter offer data compression and denoising capabilities for ultrasonic signals suitable for target detection, pattern recognition, and material characterization.

69 4.

CHAPTER IV

4 OPTIMAL DWT KERNEL FOR ULTRASONIC DATA REPRESENTATION

4.1

Introduction

The discrete wavelet transform (DWT) has many characteristics that make it a suitable tool for ultrasonic signal analysis. Among these features are good time and scale (frequency) localization, high signal similarity (one can tailor the wavelet structure to the signal being analyzed), constant relative bandwidth (this allows an estimation of frequency components with high resolution), and orthogonality (critical for matching a unique basis element to a portion of the signal). The filter bank structure of the DWT correlates the ultrasonic signal with the wavelet kernel in different channels. In this chapter, we present preliminary results on the determination of similarities between the wavelet kernel and the ultrasonic echoes to maximize data compression ratios. We show that the design of wavelet kernels via optimization methods improves data compression and denoising. Simulated data, as well as ultrasonic experimental signals, are used to verify the results of this study. High data compression ratios and echo detection in very low signal-to-noise ratio is achieved with this approach. Once the optimized wavelet kernel is designed, lossy data compression is achieved by applying thresholding techniques over the energy coefficients from each channel of the DWT filter bank. We also analyze the trade offs between data compression ratios and accuracy to the original signal. The DWT has also been shown to be an excellent technique for ultrasonic signal denoising once the wavelet kernel tracks the ultrasonic echoes and not the noise [Car01]. Since this algorithm allows perfect reconstruction, the

70 number of energy coefficients used in the inverse DWT can adjust the noise level of the reconstructed signal. This adjustment is achieved with thresholding in the wavelet domain. A set of different wavelet kernels including Coiflets, Daubechies, and Symmlets is used in the analysis and synthesis of ultrasonic signals. A selection of wavelet kernels from this set is examined for data compression of ultrasonic signals. 4.2

Discrete Wavelet Transform

The filter bank representation of the DWT [Vet92] is shown in Figure 4.1. The low pass h0(n) and the high pass h1(n) filters are analysis filters, while g0(n) and g1(n) are low and high pass synthesis filters. The input signal s(n) is down-sampled by two (↓2) before the thresholding block, and up-sampled by two (↑2) before the synthesis filters. Data compression and denoising takes place inside the thresholding block. The inverse discrete wavelet transform (IDWT) reconstructs the signal from the thresholded wavelet coefficients.

DWT 2

s(n) h1(n)

2

Thesholding Block

h0(n)

IDWT 2

+ 2

Figure 4.1. Discrete Wavelet Transform The DWT can be described in the Z domain as:

g0(n) g1(n)

sˆ(n)

71

⎡ H ( z ) H 0 (− z )⎤ ⎡ S ( z ) ⎤ 1 Sˆ ( z ) = [G0 ( z ) G1 ( z )].⎢ 0 ⎥.⎢ S (− z )⎥ ( ) ( ) − H z H z 2 ⎦ 1 ⎣ 1 ⎦⎣

(4.1)

S ( z) [G0 ( z ).H 0 ( z ) + G1 ( z ).H 1 ( z )] + Sˆ ( z ) = 2 S (− z ) [G0 ( z ).H 0 (− z ) + G1 ( z ).H 1 (− z )] 2

(4.2)

where G0 ( z ) , G1 ( z ) , H 0 ( z ) , H 1 ( z ) , S (z ) , and Sˆ ( z ) are the Z transforms of g 0 ( n ) , g 1 ( n ) , h0 (n) , h1 ( n ) , s (n) , and sˆ( n) respectively. In order to have

Sˆ ( z ) = c.S ( z ).z − n0 (such that Sˆ ( z ) is a delayed and amplitude scaled version of S (z) ) the following conditions for perfect reconstruction have to be satisfied [Say00]

G0 ( z ).H 0 (− z ) + G1 ( z ).H 1 (− z ) = 0 → ⎧ G0 ( z ) = H 1 (− z ) ⎨ ⎩G1 ( z ) = − H 0 (− z )

(4.3)

G0 ( z ).H 0 ( z ) + G1 ( z ).H 1 ( z ) = c.z

− n0

→

H 0 ( z ) = H 1 (− z ) Many solutions for H 0 ( z ) and H1(z) satisfy the above conditions. One way to solve for H 0 ( z) and H1(z) is to use quadrature mirror filters (QMF) [Vai93]. Perfect reconstruction in QMF is obtained if the filters are orthogonal to the even-shifted version of themselves: M 2

∑ h ( n) h ( n + 2k ) = δ ( k ) , k =0

i

i

i = 0,1

(4.4)

where M is the filter length (the filters in this chapter are always FIR filters unless stated otherwise). The above equation sets the condition for perfect reconstruction.

72 Although this is an important requirement in the solution of a wavelet kernel, there are other issues that affect the performance of the kernel in representing data. In this chapter, we investigate wavelet kernels that not only result in near perfect reconstruction, but also that improve signal data compression. The next section presents optimization techniques in estimating an optimal wavelet kernel, referred to as the “Echo” wavelet. The goal is to achieve a compromise between perfect reconstruction and data compression. 4.3

Optimal Design of the Echo Wavelet

In this chapter, a nonlinear optimization algorithm is used since both objective and constraint functions are nonlinear. In order to estimate the Echo wavelet kernel, an objective function is defined as the mean squared difference between the original signal s(n) and its compressed version ŝ(n) (the reconstructed signal with a given compression ratio). Depending on the desired compression ratio, the thresholding block keeps the correspondent number of wavelet coefficients. The smaller energy wavelet coefficients are discarded first.

Er =

N −1

∑ (sˆ(n) − s(n) )

2

(4.5)

n=0

The goal of the optimization is to generate a wavelet kernel as similar as possible to the ultrasonic backscattered signal. Figure 4.2 shows a block diagram of the objective function.

DWT

Threshold

IDWT

sˆ( n)

s(n)

Figure 4.2. Objective Function

∑(sˆ(n) − s(n))

2

Er

73 The solution for the Echo wavelet kernel minimizes the objective function in Equation 4.5 subject to the constraint described in Equation 4.4 [Col96]. This minimization problem can be described as

min

N −1

∑ (sˆ(n) − s(n) )

2

n =0

s.t.

(4.6)

M 2

∑ h ( n)h ( n + 2k ) = δ ( k ) , k =0

i

i

i = 0,1

Numerical methods using Lagrange multipliers were used to evaluate the optimization function. The Lagrange algorithm simplifies the calculation of nonlinear constrained optimization problems since it does not require the explicit solution of the conditions. The Lagrangian function can be written as

⎛ M2 ⎞ ⎜ ⎟ 2 L ( h, λ ) = (sˆ( n) − s ( n) ) − ∑ λi ⎜ ∑ hi ( n) hi ( n + 2k ) − δ ( k ) ⎟ i =1 ⎜ k =0 ⎟ ⎝ ⎠ m

(4.7)

where λi is the set of Lagrange multipliers. The necessary conditions so that h* is a local minimizer of the mean squared error function are stated as follows [Nas96]

∇ x L( h* , λ* ) = 0 Z ( h* )T ∇ xx L( h* , λ* ) Z ( h* ) is positive semi − definite 2

(4.8)

where Z(h*) is a null-space matrix for the Jacobian matrix of the constraint function. The necessary conditions to satisfy the minimization problem so that h* is a strict minimizer are

74

∇ x L( h* , λ* ) = 0 Z ( h* )T ∇ xx L( h* , λ* ) Z ( h* ) is positive definite 2

(4.9)

Conventional wavelet kernels are used as the starting point (seed) for the optimization routine. The wavelets used are the four-tap Daubechies (Daub4), six-tap Coiflet (Coif6) and eight-tap Symmlet (Symm8). The resultant wavelet kernels from the conventional wavelets are called Echo4, Echo6, Echo8, and are presented in the following tables. The Echo wavelets kernels shown are preliminary results of the optimization algorithm when the compression ratio is set to 78%. The tables also show the conventional wavelet kernel used as the seed for the optimization algorithm. 4.4

Performance Evaluation of Wavelet Kernels

In order to evaluate the performance of different wavelet kernels to noise supression, white Gaussian noise (WGN) was added to experimental ultrasonic data. The reconstructed signal is given by

sˆ(n) = IWT{Thresh [DWT (s(n) + g (n) )]} where g(n) ~ N(0,1). The signal s(n) is shown in Figure 4.3.

(4.10)

75

Table 4.1. Daub4 and Echo4 Wavelet Kernels h0(0) h0(1) h0(2) h0(3)

Daub4 0.48296292 0.83651631 0.22414387 -0.12940952

Echo4 0.57276109 0.78891378 0.19175563 -0.13921694

Table 4.2. Coif6 and Echo6 Wavelet Kernels h0(0) h0(1) h0(2) h0(3) h0(4) h0(5)

Coif6 0.03858078 -0.12696913 -0.07716156 0.60749164 0.74568756 0.22658427

Echo6 0.05387916 -0.17110494 -0.04582076 0.62389061 0.72505672 0.22831278

Table 4.3. Symm8 and Echo8 Wavelet Kernels h0(0) h0(1) h0(2) h0(3) h0(4) h0(5) h0(6) h0(7)

Symm8 -0.10714890 -0.04191097 0.70373907 1.13665824 0.42123453 -0.14031762 -0.01782470 0.04557034

Echo8 -0.07576571 -0.02963553 0.49761867 0.80373875 0.29785780 -0.09921954 -0.01260397 0.03222310

76

150

100

50

0

-50

-100

-150

0

50

100

150

200

250

Figure 4.3. Experimental A-Scan Data from a Thin Metal Sample

300

77 A Monte Carlo simulation (MCS) with 1,000 realizations was performed over this experimental data. A Gaussian fit of the set of realizations provides a comparative measure of the wavelet kernels’ performance. The fitting for the Daub4 and Echo4 is presented in Figure 4.4, where the vertical axis represents the number of events and the horizontal axis the mean squared reconstruction error. Table 4.4 summarizes the reconstruction error results for all wavelet kernels. These results show that the optimized wavelet kernels have a superior performance when compared to the conventional wavelet kernels. Table 4.4. Monte Carlo Simulation Result Summary

Daub4 Echo4 Coif6 Echo6 Symm8 Echo8

Hard Threshold σ2 μ 311.16 28.53 288.44 33.82 287.42 31.38 287.76 31.89 253.11 40.71 252.34 42.37

Soft Threshold σ2 μ 362.11 60.83 346.70 63.24 350.81 62.57 351.63 63.23 316.82 67.48 317.39 68.23

Figure 4.5 shows that the optimal Echo wavelet achieves a higher compression ratio when compared to Daub4. Similar results are shown for Coif6, Echo6, Symm8, and Echo8 in Figure 4.6 and Figure 4.7. This figure also confirms that the reconstruction error is smaller using the optimal wavelet kernels than when using other conventional wavelet kernels.

78

Figure 4.4. Fitting Results of MCS with Daub4 and Echo4. a) Echo4 with HT, b) Echo4 with ST, c) Daub4 with HT, and d) Daub4 with ST.

79

Normalized Mean Squared Error

1

o - Daub4 . - Echo4

0.5

0 90

91

92

93

94

95

96

97

98

99

100

Compression Ratio [%] Figure 4.5. Relation between Normalized MSE and Compression Ratio Using Daub4 and Echo4 Wavelet Kernels

80

Normalized Mean Squared Error

1

o - Coif6 . - Echo6

0.5

0 90

92

94

96

98

100

Compression Ratio [%] Figure 4.6. Relation between Normalized MSE and Compression Ratio Using Coif6 and Echo6 Wavelet Kernels

81

Normalized Mean Squared Error

1

o - Symm8 . - Echo8

0.5

0 90

92

94

96

98

100

Compression Ratio [%] Figure 4.7. Relation between Normalized MSE and Compression Ratio Using Symm8 and Echo8 Wavelet Kernels

82 B-Scan images of ultrasonic signals are widely used in NDE applications. Figure 4.8 shows the B-Scan image from a coin with a rough surface. This image was compressed using Echo4 and Daub4, and the results are shown in Figure 4.9 and Figure 4.10. Table 4.5 presents a quantitative comparison of the wavelet kernels applied to the B-Scan image. This result demonstrates the superiority of the optimal Echo wavelet kernel when compared to conventional wavelet kernels. Table 4.5. B-Scan Compression Ratio Comparison Daub4 Echo4 Coif6 Echo6 Symm8 Echo8

Compression Ratio [%] 87 87 75 75 95 95

Reconstruction Error 1 0.84 1 0.89 1 0.91

83

Figure 4.8. Original B-Scan Image

84

Figure 4.9. Original Image Compressed 87% Using Echo4

85

Figure 4.10. Original Image Compressed 87% Using Daub4

86 4.5

Conclusions

In this chapter, we presented a method to design an optimal wavelet kernel based on the similarity between the wavelet kernel and the ultrasonic echoes. The optimal Echo wavelet kernel is obtained through a nonlinear constrained optimization method that allows near perfect reconstruction. This optimization algorithm is implemented numerically using Lagrange multipliers, where the optimization is achieved by minimizing the reconstruction error constrained to the even-shift orthogonality of the optimal Echo wavelet kernel. The similarity between the Echo wavelet kernel and the ultrasonic echo leads to the superior performance of these optimal wavelet kernels when compared to conventional kernels in ultrasonic signal compression and noise suppression. Experiments in simulated and experimental ultrasonic data indicate that the optimal Echo wavelet kernels outperform conventional wavelet kernels Daubechies, Symmlet, and Coiflet.

87 5. CHAPTER V 5 PARAMETER ESTIMATION ALGORITHM FOR GAUSSIAN DECOMPOSITION OF ULTRASONIC SIGNALS

5.1

Introduction

Signal modeling and parameter estimation for detecting and estimating multiple interfering echoes has been the subject of study in the field of ultrasonic imaging over the past two decades. In particular, Saniie [San81] and Saniie et al. [San88], [San89] have dealt with the ultrasonic signal analysis and modeling for non-destructive evaluation (NDE) applications when targets are reverberant and/or randomly distributed. In medical imaging, a large number of papers have dealt with modeling of ultrasonic signals and obtaining parameters for mean scatterer spacing and tissue characterization (e.g., [Abe96], [Cha98], [Che94], [Fel84], [Liu89], [Lan90], [Sim97], [Var95], [Wea93]). Modeling of superimposed signals and parameter estimation have also been studied by Feder and Weinstein [Fed88].

More recently, the modeling and estimation of the

ultrasonic signal parameters using the maximum likelihood estimation and the expectation maximization algorithm have been introduced by Demirli and Saniie [Dem01a] and [Dem01b]. In this chapter, we present a successive parameter estimation algorithm that relies on the assumption that any ultrasonic signal, no matter how complex it is, can be decomposed into the superposition of multiple single Gaussian echoes. The goal is to efficiently estimate the parameters of the individual echoes. Most importantly, with

88 ultrasonic signal parameters we can establish the analytical relationship between the signal model and the physical parameters of materials. The continuous wavelet transform (CWT) is an effective method to display the time×frequency (TF) information of signals, and it has been utilized for flaw detection in ultrasonic applications [Mal95] and [Abb97]. In particular, the Morlet wavelet [Car03] is used to successively estimate the echo parameters (amplitude, bandwidth, phase, time-ofarrival, and center frequency). In this chapter, it is shown analytically that the CWT (i.e. time×frequency representation) yields an exact solution for the time-of-arrival and a biased estimation of the center frequency. Consequently, a modified CWT (MCWT) based on the Gabor-Helstrom transform is introduced as a means to exactly estimate both time-of-arrival and center frequency of ultrasonic echoes. The parameter estimation method presented in this chapter uses the MCWT to perform the correlation of a mother wavelet with the ultrasonic signal [Gro84] and [Gou84]. Since this is a successive approach, the parameter estimation algorithm keeps searching until the estimation satisfies the error criteria. The error criteria can be generated based on the maximum number of echoes, the minimum echo energy, and/or the position of echoes. The parameters of the ultrasonic echo have different physical significance. Furthermore, the error in the estimation of each of the parameters affects the overall estimation accuracy differently. Thus, we have analyzed the sensitivity of the reconstruction error to the disparity of the estimated parameters. The successive parameter estimation algorithm presented in this chapter has several advantages over discrete time lossy compression techniques such as JPEG [Wal91] and SPIHT [Sai96]. Most data compression algorithms offer a compromise between

89 compression ratio and signal fidelity [Ger91], [Say00], [Car02]. Furthermore, discrete time techniques are not suitable for the estimation of the ultrasonic signal parameters (i.e., bandwidth, amplitude, center frequency, phase, and time-of-arrival). Therefore, in both medical and industrial imaging applications, a well-defined modeling of the ultrasonic echoes leading to the accurate parameter estimation is highly desirable for target detection, deconvolution, object classification, velocity measurement, and ranging system. 5.2

Successive Parameter Estimation Algorithm

The ultrasonic signal can be represented as a superposition of Gaussian echoes [Dem01a]. In pulse-echo ultrasonic testing the backscattered echo from a single reflector can be modeled as

[

]

f Θ (t ) = β exp − α (t − τ ) 2 cos(2πf c (t − τ ) + φ )

(5.1)

where Θ = [α , β , f c , φ , τ ] denotes the parameter vector. The parameters of this model are closely related to the physical properties of the ultrasonic signal propagating through the material. The time-of-arrival, τ, is related to the distance between the transducer and the reflector. The amplitude, β, is a function of the attenuation of the original signal and the size of the reflector relative to the beam field. The amplitude of the reflected signal also depends on the size of the reflector or scatterer compared with the wavelength. The parameters fc and α are the center frequency and bandwidth factor respectively. These parameters are governed by the transducer frequency characteristics and the propagation path. The phase of the signal, φ, accounts for the distance, impedance, size, and orientation of the reflector [San95].

90 The first step in building the parameter estimation algorithm is to identify the behavior of the signal parameters and the influence they have in the reconstruction error of the signal. Thus, we examined the reconstruction error as each of the parameters is altered. This simulates the situation where all the parameters, with the exception of one, are correctly estimated. The reconstruction error, Er, is calculated as

Er = f Θ (t ) − f Θˆ (t )

(5.2)

ˆ is the vector of estimated parameters. Figure 5.1 shows how Er behaves as where Θ

each of the estimated parameters deviates from -10% to 10% of its actual value. This figure reveals that the reconstruction error is more sensitive to time-of-arrival, τ, compared to other parameters. Hence τ is the most critical parameter to be estimated, followed by fc, β, φ, and α (i.e., successive parameter estimation algorithm). The successive parameter estimation algorithm is a recursive method that starts with a TF representation (MCWT) of the input signal. The successive estimation is achieved by applying a window to the TF representation in order to separate interfering echoes. The window process removes part of the echo’s energy that overlaps with neighboring echoes, therefore, it creates an incomplete echo. These windows are centered on each MCWT peak, and their size is determined by the proximity of the nearby echoes. Therefore, the parameter estimation algorithm relies on the detection of the ultrasonic echoes in the TF representation of the ultrasonic signal. The peaks of the TF representation provide information about the time-of-arrival and frequency of the multiple echoes embedded in the ultrasonic signal. Upon detection of a peak, an automatic windowing procedure is used to separate adjacent and interfering echoes.

91

←τ ↑f

0

c

r

Reconstruction Error, E [dB]

10

−10

←β

↓φ

←α

−20 −30 −40 −50 −10

−5

0 Θ variation [%]

5

10

Figure 5.1. Effect of Uncertainty of the Values of Parameters on Reconstruction Error, Er.

92 The window design strategy used to separate the ultrasonic echoes depends on the noise embedded on the signal. In the case where the SNR is high, the window procedure separates neighbor echoes using their respective projections in the time and frequency domains. If the SNR is low, it is desirable to constrain the window to a smaller region around the TF representation peak that represents the best SNR of the echo. The window must be small to minimize the introduction of noise in the estimation, but it must be large enough to contain sufficient information about the echo. 5.2.1

Echo Parameter Estimation Algorithms

The CWT of f Θ (t ) with respect to a wavelet kernel ψ(t) is defined as [Rao98] ∞

CWT ( a , b) =

∫ f (t )ψ (t )dt Θ

*

a ,b

(5.3)

t = −∞

The variable b represents time shifts in the wavelet kernel, while a is a positive variable and it is referred to as the scale of the dilation. The CWT maps f Θ (t ) into a 2dimensional TF representation. The inverse CWT is defined as ∞

∞

1 1 f Θ (t ) = ∫ ∫ 2 CWT ( a, b)ψ a ,b (t )dbda c a = −∞ b= −∞ a

(5.4)

where c is a scaling constant that satisfies the wavelet’s admissibility condition [Rao98] and

ψ a ,b (t ) =

1 ⎛t −b⎞ ψ⎜ ⎟ a ⎝ a ⎠

(5.5)

The distance between f Θ (t ) and ψa,b(t) is the norm of the difference of the two signals

93 f Θ (t ) −ψ a ,b (t ) = f Θ (t ) + ψ a ,b (t ) − 2 Re[CWT ( a, b)] 2

Minimizing

f Θ (t ) −ψ a ,b (t )

2

2

2

(5.6)

in terms of variables a and b will provide the best

similarity between f Θ (t ) and ψa,b(t). The quantities f Θ (t )

2

and ψ a ,b (t )

2

on the right-

hand side of the above equation are positive constants which are independent of a and b. Hence the minimization of f Θ (t ) −ψ a ,b (t ) Furthermore, f Θ (t ) + ψ a ,b (t )

2

2

implies the maximization of Re[CWT (a, b)] .

also measures the similarity between the echo, f Θ (t ) , and

the wavelet, -ψa,b(t). This similarity implies that the minimum of Re[CWT (a, b)] also represents the closeness between the ultrasonic echo and the wavelet. Therefore, the maximum of the absolute value of Re[CWT (a, b)] , abs{Re[CWT (a, b)]}, in terms of variables a and b represents the point of best similarity between the signal f Θ (t ) and the wavelet kernel ψa,b(t) or -ψa,b(t). The Morlet [Gou84] is the wavelet kernel of choice due to its similarity to the echoes and that brings many advantages in the decomposition of ultrasonic signals [Mal95], [Abb97], and [Car01]. The Morlet wavelet is defined as

ψ a ,b (t ) =

2 ⎡ 1 ⎛t −b⎞ ⎛ t − b ⎞⎤ exp ⎢ − γ 0 ⎜ ⎟ + iω0 ⎜ ⎟⎥ a ⎝ a ⎠ ⎝ a ⎠⎥⎦ ⎢⎣

(5.7)

where the variable a tracks the frequency and the variable b tracks the time-of-arrival of the echo, ω 0 and γ 0 are the center frequency and bandwidth factor of the Morlet wavelet kernel, respectively. The term

1 ensures that the energy of the wavelet kernel a

is the same for all a and b. The Morlet wavelet kernel is a complex function and it is one

94 sided in the frequency domain. Hence, the CWT of the ultrasonic echo shown in Equation 1

is

equivalent

to

the

CWT

of

the

ultrasonic

echo

represented

as

β exp[− α (t − τ ) 2 + i (2πf c (t − τ ) + φ )] . Then, the CWT (a, b) of a single ultrasonic echo becomes CWT ( a, b) = ⎡ 2⎛ γ0 ⎞ ⎛ 2bγ 0 iω0 ⎞ ⎛ 2 γ 0 b 2 ibω0 ⎞⎤ ⎜ ⎟⎟⎥dt exp α 2 ατ ω ατ ω τ φ i i − + − − − + + − + − + − t t i ⎜ ⎟ ⎜ ⎟ ⎢ c c 2 2 2 ∫ ⎜ a a a a a a t = −∞ ⎝ ⎠ ⎝ ⎠ ⎝ ⎠⎦ ⎣

β

∞

(5.8)

where ω c = 2πf c . The solution to the above equation is simplified to CWT ( a , b) = 2 ⎡ ⎛ ω0 ⎞ 4αγ 0 γ ⎞⎞ ⎤ ⎛ ⎛ αω0 γ 0ω c ⎞ ⎛ 2 + 2 ⎟(b − τ ) + 4φ ⎜ α + 02 ⎟ ⎟ ⎥ ⎢ − ⎜ ω c − ⎟ − 2 (b − τ ) + i ⎜ 4⎜ a ⎠ a a ⎠ a ⎠⎠ ⎥ aπ ⎝ ⎝ ⎝ ⎝ a β exp ⎢ ⎢ ⎥ γ0 ⎞ αa 2 + γ 0 ⎛ 4⎜ α + 2 ⎟ ⎢ ⎥ a ⎠ ⎝ ⎣⎢ ⎦⎥

(5.9)

The magnitude of the CWT (a, b) is given by 2 ⎡ ⎛ ⎤ ω0 ⎞ 4γ 0α 2 ⎟ − 2 (b − τ ) ⎥ ⎢ − ⎜ ωc − aπ a ⎠ a ⎥ exp ⎢ ⎝ CWT ( a , b) = β ⎢ ⎥ γ0 ⎞ αa 2 + γ 0 ⎛ 4⎜ α + 2 ⎟ ⎢ ⎥ a ⎝ ⎠ ⎣ ⎦

(5.10)

The maximum of CWT ( a, b) in terms of a and b is the same as the maximum of abs{Re[CWT ( a , b)]}, because abs{Re[CWT ( a , b)]} can be represented as abs{Re[CWT ( a , b)]} = CWT ( a , b) abs{cos(g (Θ))}

(5.11)

where g (Θ) is a function of the echo parameters. The term CWT ( a , b) in the right hand side of the above equation does not depend on the phase of the echo, φ . The cosine term, on the contrary, depends on φ . Hence, this extra degree of freedom can be used to

95 set the cosine term to unity, which yields that the maximum of CWT ( a, b) in terms of a and b is the same as the maximum of abs{Re[CWT ( a , b)]}. The maximization of CWT ( a , b) can be obtained by taking the partial derivatives of the above equation with

respect to variables b and a and setting the outcome to zero. 2 ⎡ ⎛ ⎤ ω0 ⎞ 4γ 0α 2 − − ω ⎜ c ⎟ − 2 (b − τ ) ⎥ ⎢ ∂ CWT ( a , b) aπ a ⎠ a ⎥⎛⎜ 2αγ 0 (b − τ ) ⎞⎟ = 0 = −β exp ⎢ ⎝ 2 ⎢ γ ⎞ ⎥⎜⎝ a 2α + γ 0 ⎟⎠ ∂b αa + γ 0 ⎛ 4⎜ α + 02 ⎟ ⎢ ⎥ a ⎠ ⎝ ⎣ ⎦

(5.12)

The solution to the above equation is b = τ , which proves that the CWT peak is the exact estimation of the time-of-arrival. Furthermore, ∂ CWT ( a , b) ∂a

πβ 2

(a α + γ )

2

2

0

2

=

2 ⎡ ⎛ ω0 ⎞ ⎤ ⎟ ⎥ ⎢ − ⎜ ωc − ⎞ a ⎠ ⎥⎛ a (− aγ 0ω c2 − a 2αωc ω0 + ω02 aα + γ 0ω cω0 ) ⎜ exp ⎢ ⎝ + (γ 0 − a 2α )⎟⎟ = 0 ⎢ ⎛ γ 0 ⎞ ⎥⎜⎝ a 2α + γ 0 ⎠ ⎢ 2⎜ α + a 2 ⎟ ⎥ ⎠ ⎝ ⎦ ⎣

(5.13)

which implies that ⎛ a (− aγ 0ωc2 − a 2αωcω0 + ω02aα + γ 0ωcω0 ) ⎞ ⎜⎜ + (γ 0 − a 2α )⎟⎟ = 0 2 a α +γ0 ⎝ ⎠

(5.14)

The solution of the above equation does not yield an exact estimation of the center frequency which is ωc =

ω0 a

. Hence, there is a bias in the estimation of the center

frequency when using the CWT. To calculate the bias, it is considered that

ωc =

ω0 a

5.14 leads to

+ δ , where δ is the estimation bias. Substituting this result into Equation

96

δ 2 (a 2γ 0 )+ δ (ω0 aγ 0 + ω0 a 3α ) + a 4α 2 − γ 02 = 0

(5.15)

The explicit solution for the bias is − ω0 (a 2α + γ 0 ) ± ω02 (a 4α 2 + 2a 2αγ 0 + γ 02 ) − 4γ 0 (a 4α 2 − γ 02 ) δ+ , δ− = 2aγ 0

(5.16)

There are two possible solutions to Equation 5.15, δ + and δ − . In the case where

γ 0 = α and a = 1 (i.e., ωc = ω0 ) the above equation simplifies to

δ+ =

− ω0 (a 2α + γ 0 ) + ω02 (a 4α 2 + 2a 2αγ 0 + γ 02 ) − 4γ 0 (a 4α 2 − γ 02 ) =0 2aγ 0

(5.17)

and − ω0 (a 2α + γ 0 ) − ω02 (a 4α 2 + 2a 2αγ 0 + γ 02 ) − 4γ 0 (a 4α 2 − γ 02 ) δ− = = −2ω0 2aγ 0

(5.18)

Thus the correct solution for the center frequency bias is given by δ+. Furthermore, the bias is unknown and cannot be used to correct the frequency estimation since the bandwidth factor and the phase are not known a priori. To circumvent the biasness of the CWT a modified version of the CWT (MCWT) has been developed. The MCWT is introduced as a means to exactly estimate all parameters of the ultrasonic echo. The estimation of the center frequency, time-of-arrival, bandwidth factor, phase, and amplitude is executed with an overcomplete Morlet wavelet kernel,

ψ Θˆ (t ) , that spans in γ (bandwidth factor) and θ (phase) space ψ Θˆ (t ) =

⎡ ⎤ ⎛t −b⎞ 2 exp ⎢ − γ (t − b ) + iω0 ⎜ ⎟ + iθ ⎥ ε ⎝ a ⎠ ⎣ ⎦

1

(5.19)

97

where the term

1

ε

normalizes the energy of the modified wavelet kernel, and

π . This overcomplete Morlet wavelet differs from the kernel used in the Gabor2γ

ε=

Helstrom transform ([Hel66] and [Mal95]) by including two additional parameters: phase ˆ = ⎡γ , βˆ , ω0 , θ , b⎤ represents the vector of estimated and bandwidth. The Θ ⎢⎣ ⎥⎦ 2πa

()

ˆ of a single echo is given by parameters. The MCWT Θ

( ) ∫f

ˆ = MCWT Θ

∞

Θ

(t )ψ Θ*ˆ (t )dt

(5.20)

t = −∞

=

β ε

∞

iω ⎞ ⎛ ibω0 ⎡ ⎛ ⎞⎤ exp ⎢ − t 2 (α + γ ) − t ⎜ − 2τα − iω c − 2bγ + 0 ⎟ − ⎜ ατ 2 + iω cτ − iφ + γb 2 − + iθ ⎟⎥ dt a a ⎝ ⎠ ⎝ ⎠⎦ ⎣ t = −∞

∫

The solution to the above equation results in

()

ˆ = MCWT Θ

β ε

2 ⎡ ⎛ ω0 ⎞ ⎡ ⎛ αω ⎤⎤ ⎞ 2 − − ω ⎟ − 4αγ (b − τ ) + i ⎢4⎜ 0 + γωc ⎟(b − τ ) + 4(α + γ )(φ − θ )⎥ ⎥ (5.21) ⎢ ⎜ c a ⎠ π ⎝ ⎠ ⎣ ⎝ a ⎦⎥ exp ⎢ ⎢ ⎥ α +γ 4(α + γ ) ⎢ ⎥ ⎣⎢ ⎦⎥

It is the objective of the parameter estimation algorithm to find the peaks of the TF representation of the ultrasonic echo to estimate the signal’s center frequency and time-

()

ˆ is used for the TF of-arrival. To accomplish this goal, the magnitude of the MCWT Θ representation of the signal, which is given by

()

ˆ = β MCWT Θ

ε

2 ⎤ ⎡ ⎛ ω0 ⎞ 2 ⎢ − ⎜ ωc − ⎟ − 4αγ (b − τ ) ⎥ π a ⎠ ⎥ exp ⎢ ⎝ ⎥ ⎢ α +γ 4(α + γ ) ⎥ ⎢ ⎦ ⎣

(5.22)

98 The maximum of the above equation can be obtained by taking partial derivatives as a function of a and b, as follows.

()

ˆ ∂ MCWT Θ ∂a

()

ˆ ∂ MCWT Θ ∂b

=

=

β ε

−β

ε

2 ⎤ ⎡ ⎛ ω0 ⎞ ⎛ω ⎞⎞ 2 ⎛ ω − − − 4αγ (b − τ ) ⎥⎜ ω0 ⎜ 0 − ω c ⎟ ⎟ ⎜ ⎟ ⎢ c π a ⎠ ⎠⎟ = 0 ⎥⎜ ⎝ a exp ⎢ ⎝ ⎥⎜ 2a 2 (α + γ ) ⎟ ⎢ α +γ 4(α + γ ) ⎟ ⎥⎜ ⎢ ⎠ ⎦⎥⎝ ⎣⎢

2 ⎤ ⎡ ⎛ ω0 ⎞ 2 ω − − − 4αγ (b − τ ) ⎥ ⎜ ⎟ ⎢ c a ⎠ π ⎥⎛⎜ 2αγ (b − τ ) ⎞⎟ = 0 exp ⎢ ⎝ ⎥⎜⎝ (α + γ ) ⎟⎠ ⎢ 4(α + γ ) α +γ ⎥ ⎢ ⎦⎥ ⎣⎢

The maximum of Equation 5.22 is reached when b = τ and

ω0 a

(5.23)

(5.24)

= ωc . The solutions to

Equations 5.23 and 5.24 show that the peak (maximum) of the

()

ˆ MCWT Θ

representation exactly (with no bias) estimates the time-of-arrival and center frequency of the ultrasonic echo. It is important to point out that these estimates are not a function of the phase and the bandwidth of the kernel, which is a desirable property. Consequently, the TF representation based on Equation 5.22 can be obtained by using γ = 1 and θ = 0 .

()

ˆ is proportional to the amplitude of the actual Furthermore, the peak value of MCWT Θ echo, and lead to the estimation of β. Similarly, the estimation of the phase and bandwidth factor of the ultrasonic echo is

{

( ) } as a function of θ

ˆ determined by taking partial derivatives of Re MCWT Θ

{

and γ ,

( ) } ) is used in this step because

ˆ respectively. The real part of the MCWT ( Re MCWT Θ

the phase information is not contained in the magnitude representation of the transformation.

99

{

( )}

ˆ = Re MCWT Θ 2 ⎡ ⎛ ω0 ⎞ ⎡ ⎛ αω ⎤ ⎞ 2⎤ − − ω ⎟ − 4αγ (b − τ ) ⎥ ⎢ ⎜ 0 + γωc ⎟(b − τ ) + (α + γ )(φ − θ )⎥ ⎢ ⎜ c π a a ⎠ ⎠ ⎥ cos⎢ ⎝ exp ⎢ ⎝ ⎥ (α + γ ) α +γ 4(α + γ ) ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥⎦ ⎣ ⎦ ⎣

β ε

(5.25)

The center frequency and time-of-arrival of the ultrasonic echo have already been estimated in the previous step of the parameter estimation algorithm. Hence, Equation 5.25 can be simplified to

{

( )}

ˆ Re MCWT Θ

b =τ

ω0 a

=ω c

β = ε

1

1

π

⎛ 2γ ⎞ 4 ⎛ π cos(φ − θ ) = β ⎜ ⎟ ⎜⎜ α +γ ⎝ π ⎠ ⎝α +γ

⎞ 2 ⎟⎟ cos(φ − θ ) ⎠

(5.26)

Therefore, the partial derivatives as a function of θ and γ result in

{

( )}

ˆ ∂ Re MCWT Θ

b =τ

ω0 a

∂θ

{

( )}

ˆ ∂ Re MCWT Θ ∂γ

=ωc

b =τ

ω0 a

=ωc

1

⎛ 2γ ⎞ 4 ⎛ π = β ⎜ ⎟ ⎜⎜ ⎝ π ⎠ ⎝α + γ

1

⎞ 2 ⎟⎟ sin (φ − θ ) = 0 ⎠

⎛ γ 2 = cos(φ − θ )(2π ) 4 ⎜⎜ 2 ⎝α + γ

β

1

1

⎞ ⎟ ⎟ ⎠

− 12

(5.27)

⎡ (α − γ )γ − 2 ⎤ =0 ⎢ 2 ⎥ ⎣ 2(α + γ ) ⎦ 1

(5.28)

The solution to Equation 5.27 leads to maximum of Equation 5.26 when

θ = φ ± 2πk , k = 0,1, 2... The solution to Equation 5.28 results in the maximization of

()

ˆ can not only estimate Equation 5.26 when γ = α . These results show that the MCWT Θ the time-of-arrival and center frequency with no bias, but it can also allow the exact estimation of the phase and the bandwidth factor of the ultrasonic echo. 5.2.2

Description of the Algorithm and Results

The successive parameter estimation technique can be applied to ultrasonic signals consisting of multiple echoes

100 N −1

N −1

j =0

j =0

[

] (

sΘ (t ) = ∑ f Θ j (t ) = ∑ β j exp − α j (t − τ j ) 2 cos 2πf c j (t − τ j ) + φ j

)

(5.29)

To search for an optimal result, the estimation method is iterated one echo at a time until the reconstruction error, Er, is below an acceptable value Emin. The value of Emin is application specific, since it varies based on the requirements of the reconstruction quality of the signal. The noise level of the input signal also influences in the determination of Emin, as the algorithm starts reconstructing noise after a certain number of iterations. If the error is not acceptable, the estimated echoes are subtracted from the original signal, and the estimation process is repeated for additional echoes until the error is within the acceptance level. If the total number of echoes N in the ultrasonic signal is unknown, then the reconstruction error is caused by two components. The first component is the error due to the incorrect estimation of the parameters for the first M echoes (M128). They use simple IIR filtering for obtaining DCT coefficients. Figure 7.7 shows the IIR filters used for even and odd coefficients. For an input signal x[n], DCT coefficient y[k], when k is even, is given by [Che04]

y[k ] =

2 Ek ⋅ (−1) k / 2 ⋅ g N / 2 −1 (k ) N

(7.1)

where Ek = 1 / 2 for k=0 and E k = 1 for k ≠ 1 , and j

1 g j [k ] = ∑ wk [ j − n] cos(n + )θ k 2 n=0

(7.2)

wk [n] = x[n] + (−1) k x[ N − 1 − n]

(7.3)

and

When k is odd, y[k] is given as

y[ k ] = where

2 E k ⋅ ( −1) k −1/ 2 ⋅ hN / 2−1 (k ) N

(7.4)

145 j 1 h j [k ] = ∑ wk [ j − n] sin( n + )θ k 2 n =0

a) W k [j]

+

1 Z

+

-1

+

θk

cos

gj[k]

b) W k [j]

+ Z-1

-1

+ Z

1

2

2cosθk

(7.5)

+ θ k sin

hj[k]

2

-1

2cosθk Z-1

-1

-1

-1

Figure 7.7. Recursive Structures for DCT. a) Even Coefficient, b) Odd Coefficient Computation Although sequential operation is not as fast as parallel operation, the hardware requirements and power consumption are significantly reduced. The computation of DCT coefficients is independent from each other. Therefore, the system throughput can be improved by introducing more IIR units into the system. Two programmable PEs are used to implement each IIR structure in Figure 7.7. For the inverse transform channels, all the available PE resources can be allocated to one channel or these resources can be distributed to each channel for parallel execution. Therefore, the configurability of the architecture plays an important role in enhancing the throughput of DCT realization. 7.4

Conclusions

The advances in the field programmable gate array (FPGA) area have brought this technology to a point where very complex architectures can be implemented at a relatively low cost. These improvements came with applications that demand everincreasing performance out of such devices. Ultrasonic applications are known for being

146 savvy for data volume and speed. In this chapter, we present a reconfigurable architecture for ultrasonic signal compression and target detection. A unified hardware implementation is made possible since both of the algorithms are designed to share subband decomposition logic and adaptable thresholding. This architecture is a flexible and efficient solution for real-time ultrasonic imaging systems where low-power and compactness are critical.

147 8.

CHAPTER VIII 8 CONCLUSIONS

In this research, we analyzed different signal processing techniques to compress and denoise ultrasonic signals. The design of an optimal wavelet kernel for representing ultrasonic signals is summarized in [Car01], while the data compression and denoising performance of the discrete wavelet transform (DWT) is summarized in [Car02, Car03]. The performance evaluation of the DWT, discrete cosine transform (DCT), and WalshHadamard transform (WHT) is presented in [Car04]. The successive parameter estimation algorithm is summarized in [Car05]. Moreover, the implementation of a reconfigurable architecture for ultrasonic signal compression and target detection is summarized in [Oru05]. In this study, we introduced a technique to estimate an adaptive thresholding function that uses the statistical parameters of the noise embedded in the signal to generate a thresholding function based on the probability distribution function of the noise. This thresholding function is then subtracted from the sorted transform coefficients, leading to transform coefficients with a superior SNR. The results presented show that this is a very powerful technique that allows the detection of ultrasonic backscattered echoes embedded in low SNR environments. We have also analyzed ultrasonic signal compression using the DWT, the DCT, and the WHT. The results obtained show that the DWT is better in the representation of broadband signals, while the DCT and the WHT are more suitable in the representation of narrowband signals.

148 Furthermore, we introduced a method to design an optimal wavelet kernel based on the similarity between the wavelet kernel and the ultrasonic echoes. The optimal Echo wavelet kernel is obtained through a nonlinear constrained optimization method that allows near perfect reconstruction. This optimization algorithm is implemented numerically using Lagrange multipliers, where the optimization is achieved by minimizing the reconstruction error constrained to the even-shift orthogonality of the optimal Echo wavelet kernel. The similarity between the Echo wavelet kernel and the ultrasonic echo leads to the superior performance of these optimal wavelet kernels when compared to conventional kernels in ultrasonic signal compression and noise suppression. We have also analyzed a signal modeling and successive parameter estimation technique to compress and denoise ultrasonic signals. It has been shown through computer simulations and analytical derivations that both the modified continuous wavelet transform (MCWT) algorithm and the chirplet transform (CT) algorithm can efficiently estimate all the echo parameters (Gaussian-shaped echo in MCWT, chirpshaped echo in CT). The performance of this data compression algorithm has been evaluated using both simulated and experimental ultrasonic signals. The algorithm is able to compress ultrasonic data by estimating the echo parameters with high accuracy, which leads to high fidelity signal reconstruction capabilities. Overall, the signal modeling and parameter estimation algorithm not only offers data compression capabilities, but also provides parameters that can be used for signal deconvolution, target detection, pattern recognition, and material characterization. Finally, the advances in the field programmable gate array (FPGA) area have brought this technology to a point where very complex architectures can be implemented at a

149 relatively low cost. These improvements came with applications that demand everincreasing performance out of such devices. Ultrasonic applications are known for being savvy for data volume and speed. Hence, we present a reconfigurable architecture for ultrasonic signal compression and target detection. A unified hardware implementation is made possible since both of the algorithms are designed to share subband decomposition logic and adaptable thresholding. This architecture is a flexible and efficient solution for real-time ultrasonic imaging systems where low-power and compactness are critical.

150 9 BIBLIOGRAPHY

[Abb97]

Abbate, A., Koay, J., Frankel, J., Schroeder, S.C., and Das, P., “Signal Detection and Noise Suppression Using a Wavelet Transform Signal Processor: Application to Ultrasonic Flaw Detection,” IEEE Transactions on Ultrasonics, Ferroelectrics and Frequency Control, Vol. 44, No. 1, pp. 14-26, January 1997.

[Abe96]

Abeyratne, U.R., Petropulu, A.P., and Reid, J.M., “On Modeling the Tissue Response from Ultrasonic B-scan Images,” IEEE Transactions on Medical Imaging, Vol. 15, No. 4, pp. 479-490, August 1996.

[Ahm71] Ahmed, N., Rao, K.R., Abdussattar, A.L., “BIFORE or Hadamard Transform,” IEEE Transactions on Audio and Electroacoustics, Vol. 19, No. 3, pp. 225-234, September 1971. [Ahm74] Ahmed, N., Natarajan, T., Rao, K.R., “Discrete Cosine Transform,” IEEE Transactions on Computers, Vol. C-23, pp. 90-93, January 1974. [And74]

Anderson, D., Sweeney, D., and Williams, T., Linear Programming for Decision Making, West Publishing, St. Paul, 1974.

[Bre95]

Breiman, L., “Better Subset Regression Using the Nonnegative Garrote,” Technometrics, No. 37, pp. 373-384, 1995.

[Car01]

Cardoso G., and Saniie, J., “Optimal Wavelet Estimation for Data Compression and Noise Suppression of Ultrasonic NDE Signals,” IEEE Ultrasonics, Ferroelectrics and Frequency Control Symposium, Vol. 1, pp. 675-678, Atlanta, October 2001.

[Car02]

Cardoso G., and Saniie, J., “Compression of Ultrasonic Data Using Transform Thresholding and Parameter Estimation Techniques,” IEEE Ultrasonics, Ferroelectrics and Frequency Control Symposium, Vol. 1, pp. 837-840, Munich, October 2002.

[Car03]

Cardoso G., and Saniie, J., “Data Compression and Noise Suppression of Ultrasonic NDE Signals Using Wavelets,” IEEE Ultrasonics, Ferroelectrics and Frequency Control Symposium, Vol. 1, pp. 250-253, Honolulu, October 2003.

[Car04]

Cardoso G., and Saniie, J., “Performance Evaluation of DWT, DCT, and WHT for Compression of Ultrasonic Signals,” IEEE Ultrasonics, Ferroelectrics and Frequency Control Symposium, Montreal, August 2004.

151 [Car05]

Cardoso G., and Saniie, J., “Ultrasonic Data Compression via Parameter Estimation,” To be published in the IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control, February 2005.

[Cha98]

Chaturvedi P., and Insana, M.F., “Errors in Biased Estimators for Parametric Ultrasonic Imaging,” IEEE Transactions on Medical Imaging, Vol. 17, No. 1, pp. 53-61, February 1998.

[Che94]

Chen, E.J., Jenkins, W.K., and O'Brien Jr, W.D., “The Impact of Various Imaging Parameters on Ultrasonic Displacement and Velocity Estimates,” IEEE Transactions on Ultrasonics, Ferroelectrics and Frequency Control, volume 41, no 3, pp. 293-301, May 1994.

[Che99]

Chen, C.T., Linear System Theory and Design, Oxford University Press, New York, 1999.

[Che04]

Chen, C., Liu, B., Yang, J., and Wang, J., “Efficient Recursive Structures for Forward and Inverse Discrete Cosine Transform,” IEEE Transactions on Signal Processing, Vol. 52, No. 9, pp. 2665-2669, September 2004.

[Chi01]

Chiu, E., Vaisey, J., Atkins, M.S., “Wavelet-Based Space-Frequency Compression of Ultrasound Images,” IEEE Transactions on Information Technology in Biomedicine, Vol. 5, No. 4, pp. 300-310, December 2001.

[Cin01]

Cincotti, G., Loi, G., Pappalardo, M., “Frequency Decomposition and Compounding of Ultrasound Medical Images with Wavelet Packets,” IEEE Transactions on Medical Imaging, Vol. 20, No. 8, pp. 764-771, August 2001.

[Cho01]

Chong, E., Zak, S., An Introduction to Optimization, John Wiley and Sons, New York, 2001.

[Col96]

Coleman, T.F. and Y. Li, “An Interior, Trust Region Approach for Nonlinear Minimization Subject to Bounds,” SIAM Journal on Optimization, Vol. 6, pp. 418-445, 1996.

[Dau90]

Daubechies, I., “The Wavelet Transform, Time-Frequency Localization and Signal Analysis,” IEEE Transactions on Information Theory, Vol. 36, No. 5, September 1990.

[Dau92]

Daubechies, I., Ten Lectures on Wavelets, SIAM, Philadelphia, 1992.

[Dav72]

Davisson, L.D., “Rate-Distortion Theory and Application,” Proceedings of the IEEE, Vol. 60, No. 7, pp. 800-808, July 1972.

[Dav81]

David, H.A., Order Statistics, John Wiley and Sons, New York, 1981.

152 [Dem97] Demirli, R., and Saniie, J., “High Resolution Parameter Estimation for NDE Applications,” IEEE Ultrasonics, Ferroelectrics and Frequency Control Symposium, pp. 661-664, Toronto, October 1997. [Dem98] Demirli, R., and Saniie, J., “Parameter Estimation of Multiple Interfering Echoes Using the SAGE Algorithm,” IEEE Ultrasonics, Ferroelectrics and Frequency Control Symposium, pp. 831-834, Sendai, Japan, October 1998. [Dem01a] Demirli, R., and Saniie, J., “Model-Based Estimation of Ultrasonic Echoes Part I: Analysis and Algorithms,” IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control, Vol. 48, No. 3, pp. 787-802, May 2001. [Dem01b] Demirli, R., and Saniie, J., “Model-Based Estimation of Ultrasonic Echoes Part II: Nondestructive Evaluation Applications,” IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control, Vol. 48, No. 3, pp. 803-811, May 2001. [Dem01c] Demirli, R., Model Based Estimation of Ultrasonic Echoes: Analysis, Algorithms, and Applications, Ph.D. Thesis, December 2001, Illinois Institute of Technology. [Don92]

Donoho, D.L., “Denoising by Soft Thresholding”, Technical Report, Department of Statistics, Stanford University, April 1992.

[Don93]

Donoho, D.L., Johnstone, I.M., “Ideal Spatial Adaptation by Wavelet Shrinkage,” Technical Report, Department of Statistics, Stanford University, April 1993.

[Don98]

Donoho, D.L., Vetterli, M., DeVore, R.A., Daubechies, I., “Data Compression and Harmonic Analysis,” IEEE Transactions on Information Theory, Vol. 44, No. 6, pp. 2435-2476, October 1998.

[Fan02]

Fan, W., Zou, H., Sun, Y., Li, Z., and Shi, R., “Decomposition of Seismic Signal via Chirplet Transform,” IEEE International Conference in Signal Processing, Vol. 2, pp. 1778-1782, August 2002.

[Fed88]

Feder M., and Weinstein, E., “Parameter Estimation of Superimposed Signals Using the EM Algorithm,” IEEE Transactions on Acoustics, Speech, and Signal Processing, Vol. 36, No. 4, pp. 477-489, April 1988.

[Fel84]

Fellingham L.L., and Sommer, F.G., “Ultrasonic Characterization of Tissue Structure in the In Vivo Human Liver and Spleen,” IEEE Transactions on Sonics and Ultrasonics, Vol. SU-31, No. 4, pp. 418-428, July 1984.

153 [Fen01]

Feng, A., Wu, X., and Yin, Q., “Fast Algorithm of Adaptive Chirplet-Based Real Signal Decomposition,” IEEE International Symposium on Circuits and Systems, Vol. 2, pp. 6-9, May 2001.

[Fre65]

Frederick, J., Ultrasonic Engineering, John Wiley and Sons, New York, 1965.

[Fuk90]

Fukunaga, K., Introduction to Statistical Pattern Recognition, Morgan Kaufmann, San Francisco, 1990.

[Gao98]

Gao, H.Y., “Wavelet Shrinkage Denoising Using the Non-Negative Garrote,” Journal of Computational and Graphical Statistics, Vol. 7, No. 4, pp. 469-488, 1998.

[Ger91]

Gersho A., and Gray, R.M., Vector Quantization and Signal Compression, Kluwer Academic Publishers, Boston, 1991.

[Gou84]

Goupillaud, P., Grossmann, A., and Morlet, J., “Cycle-Octave and Related Transforms in Seismic Signal Analysis,” American Geoexploration, Vol. 23, pp. 85-102, 1984.

[Gre94]

Grennberg, A., Sandell, M., “Estimation of Subsample Time Delay Differences in Narrowband Ultrasonic Echoes Using the Hilbert Transform,” IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control, Vol. 41, No. 5, pp. 588-595, September 1994.

[Gro84]

Grossmann A., and Morlet, J., “Decomposition of Hardy Functions into Square Integrable Wavelets of Constant Shape,” SIAM Journal on Mathematical Analysis, Vol. 15, pp 723-736, 1984.

[Hel66]

Helstrom, C.W., “An Expansion of a Signal in Gaussian Elementary Signals,” IEEE Transactions on Information Theory, Vol. 12, No. 1, pp. 81-82, January 1966.

[Kay93a] Kay, S.M., Fundamentals of Statistical Signal Processing Volume I: Estimation Theory, Prentice-Hall, Upper Saddle River, New Jersey, 1993. [Kay93b] Kay, S.M., Fundamentals of Statistical Signal Processing Volume II: Detection Theory, Prentice-Hall, Upper Saddle River, New Jersey, 1993. [Kie95]

Kiely, A.B., and Pollara, F., “A Seismic Data Compression System Using Subband Coding,” The Telecommunications and Data Acquisition Progress Report 42-121, January-March 1995, Jet Propulsion Laboratory, Pasadena, California, pp. 242-251, March 1995.

[Kit76]

Kitajima, H., “Energy Packing Efficiency of the Hadamard Transform,” IEEE Transactions on Communications, Vol. 24, No. 11, November 1976.

154 [Lan90]

Landini L., and Verrazzani, L., “Spectral Characterization of Tissues Microstructure by Ultrasounds: a Stochastic Approach,” IEEE Transactions on Ultrasonics, Ferroelectrics and Frequency Control, Vol. 37, No. 5, pp. 448-456, September 1990.

[Law00]

Law, A., and Kelton, W., Simulation Modeling and Analysis, McGraw-Hill, New York, 2000.

[Lia04]

Liao, H., Mandal, M.K., and Cockburn, B.F., “Efficient Architectures for 1-D and 2-D Lifting Based Wavelet Transforms,” IEEE Transactions on Signal Processing, Vol.52, No.5, pp. 1315-1326, May 2004.

[Liu89]

Liu, D.L., and Saito, M., “A New Method for Estimating the Acoustic Attenuation Coefficient of Tissue from Reflected Ultrasound Signals,” IEEE Transactions on Medical Imaging, Vol. 8, No. 1, pp. 107-110, March 1989.

[Mal95]

Malik, M.A., Unified Time-Frequency Analysis of Ultrasonic Signals, Ph.D. Thesis, Illinois Institute of Technology, July 1995.

[Mal96]

Malik, M.A., Saniie, J., “Performance Comparison of Time-Frequency Distributions for Ultrasonic Nondestructive Testing,” IEEE Ultrasonics, Ferroelectrics and Frequency Control Symposium, pp. 701-704, San Antonio, Texas, November 1996.

[Mal01]

Mallat, S., A Wavelet Tour of Signal Processing, Academic Press, San Diego, 2001.

[Man92]

Mann, S., and Haykin, S., “Time-Frequency Perspectives: The Chirplet Transform,” IEEE International Conference on Acoustics, Speech, and Signal Processing, Vol. 1, pp. 417-420, March 1992.

[Man95]

Mann, S., and Haykin, S., “The Chirplet Transform: Physical Considerations,” IEEE Transactions on Signal Processing, Vol. 43, No. 11, pp. 27452761, November 1995.

[Mic94]

Micheli, G.D., Synthesis and Optimization of Digital Circuits, McGraw-Hill, New York, 1994.

[Muy95] Muyshondt, R.A., and Mitra, S., “Visual Fidelity of Reconstructed Radiographic Images Using Wavelet Transform Coding and JPEG,” IEEE Symposium on Computer-Based Medical Systems, pp. 196-204, Lubbock, Texas, June 1995. [Nas96]

Nash, S., and Soffer, A., Linear and Nonlinear Programming, McGraw-Hill, New York, 1996.

155 [Oh03]

Oh, T.H., and Besar, R., “JPEG2000 and JPEG: Image Quality Measures of Compressed Medical Images,” Proceedings of the 4th National Conference on Telecommunication Technology, pp. 31-35, January 2003.

[Oia98]

Qian, S., Chen, D., and Yin, Q., “Adaptive Chirplet Based Signal Approximation,” IEEE International Conference on Acoustics, Speech, and Signal Processing, Vol. 1, pp. 1781-1784, May 1998.

[Opp75]

Oppenheim, A.V., and Schafer, R., Digital Signal Processing, Prentice-Hall, Upper Saddle River, New Jersey, 1974.

[Oru03]

Oruklu, E., and Saniie, J., “High Speed Design and Performance Evaluation of Frequency-Diverse Orthogonal Transforms for Ultrasonic Imaging Applications,” IEEE Ultrasonics, Ferroelectrics and Frequency Control Symposium, Vol. 1, pp. 258-261, October 2003.

[Oru04a] Oruklu, E., and Saniie, J., “Ultrasonic Flaw Detection Using Discrete Wavelet Transform for NDE Applications,” IEEE Ultrasonics, Ferroelectrics and Frequency Control Symposium, Montreal, August 2004. [Oru04b] Oruklu, E., Martinez, F., and Saniie, J., “A Reconfigurable Architecture for Target Detection in High Density Clutter Environments Using Subband Decomposition Algorithms,” GSPx Technical Conference, Santa Clara, California, October 2004. [Oru04c] Oruklu, E., and Saniie, J., “Performance Evaluation of Multiplierless Architectures for Frequency Diverse Target Detection Algorithms,” IEEE International Conference on Circuits, Signals and Systems, Clearwater Beach, Florida, November 2004. [Oru05]

Oruklu, E., Cardoso, G., and Saniie, J., “Reconfigurable Architecture for Ultrasonic Signal Compression and Target Detection,” to be published in the Proceedings of the IEEE International Conference on Acoustics, Speech, and Signal Processing, Philadelphia, March 2005.

[Pap02]

Papoulis, A., and Pillai, S.U., Probability, Random Variables and Stochastic Processes, McGraw-Hill, Boston, 2002.

[Par57]

Parkinson, C.N., Parkinson’s Law and Other Studies in Administration, Ballantine Books, New York, 1957.

[Per00]

Percival, D.B., and Walden, A.T., Wavelet Methods for Time Series Analysis, Cambridge University Press, Cambridge, 2000.

[Pro96]

Proakis, J.G., and Manolakis, D.G., Digital Signal Processing: Principles, Algorithms, and Applications, Prentice-Hall, Upper Saddle River, New Jersey, 1996.

156 [Pun01]

Puniene, J., Punys, V., and Punys, J., “Ultrasound and Angiographic Image Compression by Cosine and Wavelet Transforms and its Approval in Clinical Environment,” Proceedings of the 2nd International Symposium on Image and Signal Processing and Analysis, pp. 74-79, June 2001.

[Raj02]

Rajoub, B.A., “An Efficient Coding Algorithm for the Compression of ECG Signals Using the Wavelet Transform,” IEEE Transactions on Biomedical Engineering, Vol. 49, No. 4, pp. 355-362, April 2002.

[Rao90]

Rao, K.R., and Yip, P., Discrete Cosine Transform: Algorithms, Advantages, Applications, Academic Press, Boston, 1990.

[Rao96]

Rao, K.R., and Hwang, J.J., Techniques and Standards for Image, Video and Audio Coding, Prentice-Hall, Upper Saddle River, New Jersey, 1996.

[Rao98]

Rao, R.M., and Bopardikar, A.S., Wavelet Transforms, Addison-Wesley, Reading, Massachusetts, 1998.

[Ris03]

Risueno, G.L., Grajal, J., and Ojeda, O.A.Y., “CFAR Detection of Chirplets in Coloured Gaussian Noise,” IEEE International Conference on Acoustics, Speech, and Signal Processing, Vol. 1, pp. 741-744, April 2003.

[Rob97]

Robini, M.C., Magnin, I.E., Benoit-Cattin, H., and Baskurt, A., “TwoDimensional Ultrasonic Flaw Detection Based on the Wavelet Packet Transform,” IEEE Ultrasonics, Ferroelectrics and Frequency Control Symposium, Vol. 44, pp. 1384-1392, November 1997.

[Ros00]

Ross, S.M., Introduction to Probability Models, Harcourt Academic Press, San Diego, 2000.

[Sai96]

Said, A., and Pearlman, W.A., “A New Fast and Efficient Coder Based on Set Partitioning in Hierarchical Trees,” IEEE Transactions on Circuits and Systems for Video Technologies, pp. 243-250, June 1996.

[San81]

Saniie, J., Ultrasonic Signal Processing: System Identification and Parameter Estimation of Reverberant and Inhomogeneous Targets, Ph.D. Thesis, Purdue University, August 1981.

[San88]

Saniie, J., Wang, T., and Bilgutay, N.M., “Statistical Evaluation of Backscattered Ultrasonic Grain Signals,” Journal of the Acoustical Society of America, Vol. 84, No. 1, pp. 400-408, July 1988.

[San89]

Saniie, J., and Nagle, D.T., “Pattern Recognition in the Ultrasonic Imaging of Reverberant Multilayered Structures,” IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control, Vol. 36, No. 1, pp. 80-92, January 1989

157 [San90]

Saniie, J., Donohue, K., Bilgutay, N., “Order Statistic Filters as Postdetection Processors,” IEEE Transactions on Acoustics, Speech, and Signal Processing, Vol. 38, No. 10, pp. 1722-1732, October 1990.

[San92]

Saniie, J., and Nagle, D.T., “Analysis of Order Statistic CFAR Threshold Estimators for Improved Ultrasonic Flaw Detection,” IEEE Transactions on Ultrasonics, Ferroelectrics and Frequency Control, Vol. 39, No. 5, pp. 618-630, September 1992.

[San95]

Sandell, M., and Grennberg, A., “Estimation of the Spatial Impulse Response of an Ultrasonic Transducer Using a Tomographic Approach,” Journal of Acoustical Society of America, Vol. 98, No. 4, pp. 2094-2103, October 1995.

[Say00]

Sayood, K., Introduction to Data Compression, Morgan Kaufmann, San Francisco, 2000.

[Sch97]

Schick, I.C., and Krim, H., “Robust Wavelet Thresholding for Noise Suppression,” IEEE International Conference on Acoustics, Speech, and Signal Processing, Vol. 5, pp. 21-24, April 1997.

[Sim97]

Simon, C., Shen, J., Seip, R., and Ebbini, E.S., “A Robust and Computationally Efficient Algorithm for Mean Scatterer Spacing Estimation,” IEEE Transactions on Ultrasonics, Ferroelectrics and Frequency Control, Vol. 44, No. 4, pp. 882-894, July 1997.

[Ska96]

Skahill, K., VHDL for Programmable Logic, Addison-Wesley, New York, 1996.

[Sta97]

Staunstrup, J., and Wolf, W., Hardware/Software Co-Design: Principles and Practice, Kluwer Academic Publishers, Boston, 1997.

[Sun99]

Sun, H.C., Frequency-Diverse Ultrasonic Target Detection Using Bayes, Fuzzy Discriminant, and Neural Network, Ph.D. Thesis, May 1999, Illinois Institute of Technology.

[Swe96]

Sweldens, W., "The Lifting Scheme: A Custom-Design Construction of Biorthogonal Wavelets," Applied Computing and Harmonic Analysis, pp. 186-200, 1996.

[Vai93]

Vaidyanathan, P.P., Multirate Systems and Filter Banks, Prentice-Hall, Upper Saddle River, New Jersey, 1993.

[Var95]

Varghese, T., and Donohue, K.D., “Estimating Mean Scatterer Spacing with the Frequency-Smoothed Spectral Autocorrelation Function,” IEEE Transactions on Ultrasonics, Ferroelectrics and Frequency Control, Vol. 42, No. 3, pp. 451-463, May 1995.

158 [Vet92]

Vetterli, M., Herley, C., “Wavelets and Filter Banks: Theory and Design,” IEEE Transactions on Signal Processing, Vol. 40, No 9, pp. 2207-2231, September 1992.

[Wal91]

Wallace, G.K., “The JPEG Still Picture Compression Standard,” Communications of the ACM, Vol. 34, pp. 31-44, April 1991.

[Wea93]

Wear, K.A., Wagner, R.F., Insana, M.F., and Hall, T.J., “Application of Autoregressive Spectral Analysis to Cepstral Estimation of Mean Scatterer Spacing,” IEEE Transactions on Ultrasonics, Ferroelectrics and Frequency Control, Vol. 40, No. 1, pp. 50-58, January 1993.

[Yan99]

Yang, J., and Fan, C., “Recursive Cosine Transforms with Selectable FixedCoefficient Filters,” IEEE Transactions on Circuits and Systems-II, Vol. 46, No. 2, pp. 211-216, February 1999.

[Yip78]

Yip, P., Rao, K.R., “Energy Packing for the Generalized Discrete Transforms,” IEEE Transactions on Communications, Vol. 26, No. 8, August 1978.

[Zei00]

Zeigler, B.P., Praehofer, H., Kim, T.G., Theory of Modeling and Simulation, Academic Press, San Diego, 2000.

[Zha01]

Zhang, W., Zhao, X. H., “Wavelet Thresholding Using Higher-Order Statistics for Signal Denoising,” Proceedings of the 2001 International Conference on Info-tech and Info-net, Vol. 1, pp. 363-368, Beijing, November 2001.

10