A Comparative Simulation Study of Wavelet Based Denoising Algorithms M. C. E. Rosas-Orea, M. Hernandez-Diaz, V. Alarcon-Aquino, and L. G. Guerrero-Ojeda Universidad de las Américas, Puebla Departamento de Ingeniería Electrónica y Mecatrónica Sta. Catarina Mártir C.P. 72820, Cholula, Puebla. MEXICO Email: [email protected]

Abstract In this paper we present a comparative simulation study of three denoising algorithms using wavelets. The denoising algorithms (i.e., universal threshold, minimax threshold and rigorous SURE threshold) have been used to remove white Gaussian noise from synthetic and real signals. The analysis is done by applying soft and hard thresholds to signals with different sample sizes. The mean squared error (MSE) is used to evaluate the performance of these algorithms. The results show that the rigorous SURE algorithm with a hard threshold has a better performance than other algorithms in synthetic signals. On the other hand, the universal threshold algorithm with a soft threshold shows the best performance in real signals when using the Daubechies wavelet with 5 vanishing moments.

1. Introduction It is well known that due to traffic information signals present a noise problem that complicates the transmission and reception of data. The analysis of noise and its elimination requires special mathematical tools. Since any undesirable electrical noise falls within the passband of the signal [1], the Fourier transform is the most common tool to perform this analysis. The main source of noise is the thermal noise because it is random and continuous, and occurs at all frequencies. Its power spectral density is flat, hence, it is also called "white noise" [1]. In order to eliminate this noise different denoising algorithms have been proposed (see e.g., [2][9]). In this paper we present a comparative simulation study of three wavelet based denoising algorithms applied to synthetic and real signals. The analysis is done by applying soft and hard thresholds to signals with different sample sizes. The remainder of this paper is organized as follows. A review of wavelet theory is

presented in Section 2. Section 3 presents a description of three wavelet based denoising algorithms. Simulation results are reported in Section 4. Finally, Section 5 summarizes the conclusions drawn from previous sections.

2. A review of wavelet theory Wavelet transforms involve representing a general function in terms of simple, fixed building blocks at different scales and positions. These building blocks are generated from a single fixed function called mother wavelet by translation and dilation operations. The continuous wavelet transform considers a family [10]

"

> , (1) + where + − ‘ , , − ‘, with + Á !, and b œ È <Œ +

the admissibility condition. For discrete wavelets the scale (or dilation) and translation parameters in (1) are chosen such that at level 4 the wavelet +!4 <ˆ+!
4 >‰ is +!4 times the width of b œ +!
4Î# <ˆ+!
4 >

8,! ‰ß

so the discrete version of wavelet transform is .4ß8 œ ØBÐ>Ñß <4ß8 Ð>ÑÙ _

(2) (3)

BÐ>Ñ<ˆ+!
4 > 8,! ‰.>
_ Ø † ß † Ù denotes the P#-inner product. The mother wavelet function <Ð>Ñ, scaling +! and translation ,! parameters are specifically chosen such that <4ß8 Ð>Ñ constitute orthonormal bases for P# Ð‘Ñ [10], [11]. To œ +!
4Î# (

form orthonormal bases with good time-frequency localisation properties, the time-scale parameters Ð,ß +Ñ

Proceedings of the 15th International Conference on Electronics, Communications and Computers (CONIELECOMP 2005) 0-7695-2283-1/05 $20.00 © 2005 IEEE

are sampled on a so-called in the time-scale plane, namely, +! œ # and ,! œ ", [10]-[12]Þ Thus, from Eq. (2) substituting these values, we have a family of orthonormal bases, dyadic grid

<4ß8 Ð>Ñ œ # 4Î# <ˆ# 4 >

8‰ß

BÐ>Ñß <4ß8 Ð>Ñ¡ œ #

4Î# (


_ _

BÐ>Ñ<4ß8 ˆ# 4 > 8‰.>ß (5)

A formal approach to constructing orthonormal bases is provided by multiresolution analysis (MRA) [12]. The idea of MRA is to write a function BÐ>Ñ as a limit of successive approximations, each of which is a smoother version of BÐ>Ñ. The successive approximations thus correspond to different resolutions [12]. In other words, after choosing an initial resolution N , any signal BÐ>Ñ − P# Ð‘Ñ can be expressed as [10], [12]: _

BÐ>Ñ œ "-Nß8 9Nß8 Ð>Ñ
"".4ß8 <4ß8 Ð>Ñß 8−™

4œN 8−™

(6)

where the detail or wavelet coefficients Ö.4ß8 × are defined in Eq. (5), while the approximation or scaling coefficients Ö-4ß8 × are defined by

-4ß8 œ #

4Î# (


_ _

BÐ>Ñ94ß8 ˆ# 4> 8‰.>Þ

(7)

where 94ß8 a † b denotes the scaling function. Equations (5) and (7) express that a signal BÐ>Ñ is decomposed in details Ö.4ß8 × and approximations Ö-4ß8 × to form a multiresolution analysis of the signal [12]. An attractive property of wavelets is that there exists a recursive relationship between scaling Ö-4ß8 × and wavelet coefficients Ö.4ß8 × at successive levels of resolution. That is, using Eq. (7) and the [10] yields dilation equation

-4ß8 œ "16-4

"ß#8 6

.4ß8 œ "26-4

"ß#8 6 Þ

6−™

6−™

-4

(4)

Using Eq. (3), the orthonormal wavelet transform is thus given by  

opposite direction, approximation coefficients at level of resolution 4 " are computed from scaling and wavelet coefficients at the coarser level of resolution 4 according to

Ð 8Ñ

where 16 represents the coefficients of a low-pass filter and 26 denotes the coefficients of a band-pass filter. Equation (8) denotes approximation and details coefficients respectively at level of resolution 4 and these are obtained from approximation coefficients at finer level of resolution 4 ". The sequences Ö-4ß8 × and Ö.4ß8 × are generated by down-sampling by a factor of two the output of the corresponding filters. In the

"ß8

œ "16 -4ß#8 6
"26.4ß#8 6Þ 6−™

6−™

(9)

The sequence Ö-4 "ß8 × is generated by up-sampling the output of the corresponding filters. This operation is achieved by inserting one zero every two samples. Equations (8)-(9) can be computed by a pyramid algorithm [12] called discrete wavelet transform. Equations (8) and (9) are used in the following section.

3. Denoising algorithms This section presents a description of three denoising algorithms as well as the general process for denoising signals using wavelet transforms. The algorithms are the Universal threshold, the Rigorous SURE threshold and the Minimax threshold [9]. A description of linear and non-linear methods for denoising signals is also presented. 3.1 General process

Three basic steps are required for denoising signals: decomposition, thresholding and reconstruction. The first and last steps are described in Section 2 (see Eqs. (8) and (9)). The thresholding step is divided into two steps. The first step defines an operator ?Ö † × to determine a threshold -. The second step defines another operator HÖ † × that works in the denoising process in order to obtain the modified ^Ð+ß ,Ñ coefficients, namely, - œ ?ÖWÐ+ß ,Ñ× (10) (11)

^Ð+ß ,Ñ œ HÖWÐ+ß ,Ñß -×

where +, ,, denotes the time-scale parameters (see Section 2) and WÐ+ß ,Ñ denotes the discrete wavelet coefficients computed by Eq. (8). The reconstruction step is computed by Eq. (9) using the modified ^Ð+ß ,Ñ coefficients. The computed thresholds need a noise level estimator which is expressed as [9], [13], [16] 5œ

medianÐÖlWÐN

"ß,Ñl À , œ !ß"ßÞÞÞß#N !Þ'(%&

"

"×Ñ(12)

where N denotes the number of decomposition levels. This estimator is used in the following sections.

Proceedings of the 15th International Conference on Electronics, Communications and Computers (CONIELECOMP 2005) 0-7695-2283-1/05 $20.00 © 2005 IEEE

3.2 Minimax threshold

3.4 Rigorous SURE threshold

The Minimax threshold denoising algorithm was proposed in [9]. It consists in an optimal threshold that is derived from minimising the constant term in an upper bound of the risk involved in the estimation of the function. The optimal threshold is defined as (13) -Q œ 5-‡8

The denoising algorithms described previously use global thresholds. That is, the computed threshold is applied to all wavelet coefficients. The Rigorous SURE threshold reported in [9], [13], [14], [16] describes a scheme that uses a threshold value -4 at each resolution level 4 of the wavelet coefficients. The Rigorous SURE threshold denoising algorithm is also known as SureShrink and uses the Stein's Unbiased Risk Estimate criterion to get an unbiased estimate. The threshold is computed as follows [9], [13]:

where -‡8 is defined as the value of - and it is obtained as follows:

V- Ð.Ñ 8 " V Ð.Ñ › (14) where V- Ð.Ñ œ IÐ$- Ð.Ñ
.Ñ# and V Ð.Ñ is the A‡8 À œ inf- sup. š

oracle

-WYVI œ arg min !
-
- WYVIŒ-ß YRM

oracle

ideal risk achieved with the help of an oracle. Two oracles are considered. The diagonal linear projection (DLP) and the diagonal linear shrinker (DLS). The former states when to “keep” or “kill” each empirical wavelet coefficient, while the latter states how much shrinking is applied to each wavelet coefficient. The ideal risks for these oracles are given by (15) V HPT Ð.Ñ À œ minÐ. # ß "Ñ oracle

V HPW Ð.Ñ À œ . #. " #

oracle

(16)

where mina † b denotes the minimum value of . # . The minimax principle is used in statistics to design estimators. Since the denoised signal can be assimilated to the estimator of the unknown regression function, the minimax estimator is the option that realizes the minimum, over a given set of functions, of the maximum mean square error [15].

3.3 Universal threshold

UNI

œ 5 È#logÐ8Ñ

(17)

where 8 denotes the length of the signal and 5 is given by Eq. (12). The implementation in software of this threshold requires no costly development of lookup tables. Nevertheless, the universal threshold is substantially larger than the Minimax threshold (14) for any particular value of 8 [13].

5

WYVIÐ-à \Ñ œ 8 #♣Ö3 À l\3 l Ÿ -×
minÐl\3 lß -Ñ‘# and

♣

(18)

SUREÐÑ is defined as

denotes the cardinality of the set

(19)

Ö3 À l\3 l Ÿ

-×.

Note that for discrete wavelet coefficients the variable is changed to

WÐ+ß,Ñ

5

\

.

3.5 Soft and hard threshold

The denoising algorithms can be divided into linear and non-linear methods. The linear method is independent of the size of empirical wavelet coefficients, and therefore the coefficient size by itself is not taken into account. It assumes that signal noise can be found mainly in fine scale coefficients and not in coarse scales ones [5], namely,

. œ œ .! ßß 44
  -4ß5

The Universal threshold denoising algorithm was also proposed in [9], and is also known as VisuShrink. This algorithm is an alternative to the Minimax threshold, however it uses a fixed threshold form, namely,

-

where

WÐ+ß ,Ñ

4ß5

(20)

The non-linear method is based on the idea that the white noise can be found in every coefficient and is distributed over all scales. It can be applied in two ways: hard thresholding and soft thresholding. The former cuts off coefficients below a certain threshold -, while the latter reduces all coefficients by this threshold. The hard and soft thresholds are respectively given by [5]

ß lBl  =ÐBÑ œ œ =ÐBÑ ! ß lBl Ÿ -

Proceedings of the 15th International Conference on Electronics, Communications and Computers (CONIELECOMP 2005) 0-7695-2283-1/05 $20.00 © 2005 IEEE

(21)

=ÐBÑ œ œ

ÐBÑÐlBl -Ñß lBl
! ß lBl Ÿ -

Sign

(22)

where =ÐBÑ is the analysed signal and - is the chosen threshold.

4. Simulation results In this section we compare the Universal threshold, the Rigorous SURE threshold and the Minimax threshold [9] using synthetic and real signals.

4.1 Simulation parameters The three denoising algorithms are simulated in Matlab® [15]. We have used the following synthetic signals [9], [13]: Step, Waves, Blip, Blocks, Bumps, Angles, Time Shifted Sine and Corner functions. Mathematical expressions and general features of these signals can be found in [9], [13]. The sample sizes of the synthetic signals are 8 œ "#)ß &"#ß "!#%ß #!%), %!*'. The analysis is done using N œ & levels of decomposition. We use a white Gaussian noise with a Ð!ß "Ñ and a SNR=3dB. We have also compared the algorithms using a real electric signal called "lelecum" [15]. The following wavelets are assessed: the Haar wavelet, Daubechies, Coiflet and Symmlet wavelets. The Mean Squared Error (MSE) is used to evaluate the performance of the algorithms and is expressed as follows:

!Ð0 8

MSE œ

3œ"

0/ Ñ# 8

(23)

where 8 denotes the length of the signal, 0 represents the original signal and 0/ is the estimated signal obtained from the denoised wavelet coefficients.

4.2 Results and comparisons All the analysed denoising algorithms show good signal estimation. Simulation results also show that when the sample size increases, the MSE decreases. The thresholds of the Universal threshold algorithm and the Minimax algorithm increase when the sample size increases, while the threshold of the Rigorous SURE algorithm decreases and oscillate near a certain value. The best estimations using the soft thresholding rule are presented in Table 1, while the best estimations using the hard thresholding rule are shown in Table 2. The results in Table 1 show that the soft thresholding rule is

better that the hard thresholding rule in terms of MSE (see also Table 2) for almost all synthetic signals. Table 2 shows that the best performance is achieved with the Coiflet wavelet when applied to different signals. Table 3 shows a comparison between the denoising algorithms using different wavelets. We note that these results are the averages of the MSE for all synthetic signals (see [17]). The results show that the best performance is achieved with the Coiflet 5. This is due to the fact that it has 10 vanishing moments. The worst wavelet is the Haar because it has one vanishing moment only. In the soft and hard thresholding cases we note that there is a better performance when using the hard thresholding rule. Further details on the simulation results reported in this paper can be found in [17].

Table 1. Estimations using a soft threshold.

Signal

Algorithm

Wavelet

Step UTST Coiflet 5 Waves MMST/RSST Coiflet 5 Blip UTST Symmlet 5 Blocks UTST Coiflet 5 Bumps UTST Coiflet 5 Angle MMST Coiflet 5 Sine RSST Coiflet 5 Corner RSST Daubechies 5 UTST: Universal Threshold & Soft Threshold. UTHT: Universal Threshold & Hard Threshold. MMST:Minimax Threshold & Soft Threshold. MMHT:Minimax Threshold & Hard Threshold. RSST: Rigorous Sure & Soft Threshold. RSHT: Rigorous Sure & Hard Threshold.

MSE

0.000843 0.000069 0.000317 0.000597 0.000882 0.000098 0.000134 0.000055

Table 2. Estimations using a hard threshold.

Signal

Algorithm

Wavelet

Step MMHT Coiflet 5 Waves MMHT Coiflet 5 Blip MMHT Symmlet 5 Blocks RSHT Coiflet 5 Bumps UTHT/MMHT Coiflet 5 Angle RSHT Daubechies 5 Sine RSHT Daubechies 5 Corner UTHT/RSHT Coiflet 5/Symmlet 5

MSE

0.000861 0.000077 0.000304 0.000603 0.000883 0.000106 0.000120 0.000054

We have also assessed the denoising algorithms using a real electrical signal called "lelecum" [15]. Figure 1 shows the original electrical signal. In this assessment we have tested only the algorithms with the best average in the MSE (see Table 3), that is, the UTST using the Daubechies wavelet with 5 vanishing moments and the MMHT using the Coiflet wavelet with 10 vanishing moments. In this experiment we have also used five levels of decomposition. The results are shown in Figure 2. It can be seen that although the Coiflet wavelet has 10 vanishing moments, the Daubechies

Proceedings of the 15th International Conference on Electronics, Communications and Computers (CONIELECOMP 2005) 0-7695-2283-1/05 $20.00 © 2005 IEEE

Table 3. Comparison between denoising algorithms using different wavelets. Wavelet/Algorithm Haar Daubechies 5 Coiflet 5 Symmlet 5

UTST 0.0047970 0.0037083 0.0037427 0.0038594

UTHT 0.0047660 0.0037931 0.0037294 0.0039042

MMST 0.00477408 0.00374358 0.00373075 0.00388270

MMHT 0.0047771 0.0038094 0.0037068 0.0038374

RSST 0.0047764 0.0037436 0.0037368 0.0038519

RSHT 0.0047320 0.0037920 0.0037201 0.0038268

MSE Ave. 0.0047704 0.0037650 0.0037278 0.0038604

MSE Average

0.0040269

0.0040482

0.0040328

0.0040327

0.0040272

0.0040177

0.0040309

wavelet with 5 vanishing moments performs better. This is due to the fact that the smoothness of the Daubechies wavelet increases as the order the filter increases. Thus, a better approximation of the signal may be obtained [10]. 550 500 450

Amplitude

400 350 300 250 200 150 100

0

200

400

600

800

1000 1200 Samples

1400

1600

1800

2000

Figure 1. Original electrical signal. 550 500 450

Amplitude

400

5. Conclusions

In this paper we have presented a comparative simulation study of three denoising algorithms using wavelets. The universal threshold, the minimax threshold and the rigorous SURE threshold have been used to remove white Gaussian noise from synthetic and real signals. The analysis is done by applying soft and hard thresholds to signals with different sample sizes. For the algorithms settings and set of test signals carried out in this paper the evidence suggest that the Rigorous SURE algorithm with a hard threshold has a better performance than other algorithms in synthetic signals. On the other hand, the Universal threshold algorithm with a soft threshold shows the best performance in real signals when using the Daubechies wavelet with 5 vanishing moments. This is due to the fact the smoothness of the Daubechies wavelet increases as the order the filter increases. Thus, a better approximation of the signal may be obtained.

350

6. References

300 250 200 150 100

0

200

400

600

800

1000 1200 Samples

1400

1600

1800 2000

(a)

550

[2] H. W. Ott, Noise Reduction Techniques in Electronic Systems, AT&T Laboratories, Wiley-Interscience Publications. USA. 1988.

500 450 400

Amplitude

[1] W. Tomasi, Electronic Communications Systems: Fundamentals Through Advanced. Prentice Hall, 4th. ed. 2001. New Jersey, USA.

[3] M. Wachowiak, G. Rash, P. Quesada, A. Desoky, "Comparison of Wavelet-Based and Traditional Noise Removal Techniques", Gait and Biomechanics Laboratory, University of Louisville. August, 1998.

350 300 250 200 150 100

0

200

400

600

800

1000 1200 Samples

(b)

1400

1600

1800 2000

Figure 2. Reconstructed electrical signal using (a) the UTST with Daubechies 5 and (b) the MMHT with Coiflet 5.

[4] G. Cristobal, M. Chagoyen, B. Escalante, and J. Lopez, "Wavelets Based Denoising Methods: A Comparative Study with Applications in Microscopy", Universidad Autónoma de Madrid- UNAM, 1996. [5] F. Hess, M. Kraft, M. Richter, and H. Bockhorn, "Comparison and Assessment of Various Wavelet and

Proceedings of the 15th International Conference on Electronics, Communications and Computers (CONIELECOMP 2005) 0-7695-2283-1/05 $20.00 © 2005 IEEE

Wavelet Packet based Denoising Algorithms for Noisy Data", Universität Kaiserslautern, 1997.

Thesis, Department of Electrical and Electronic Engineering, Universidad de las Américas, Puebla, MEXICO 2003.

[6] R. Öktem, L. Yaroslavsky, K. Egiazarian, and J. Astola, "Transform Based Denoising Algorithms: Comparative Study", Tampere University of Technology. Finland. 1999. [7] R. Sathish, V. Anand, "Wavelet Denoising for Plane Wave DOA Estimation by MUSIC", IEEE Conference Tencon 2003. Bangalore, INDIA, October, 2003. [8] C. Taswell, "The What, How and Why of Wavelet Shrinkage Denoising", Computing in Science & Engineering, June 2000. [9] D. Donoho, I. M. Johnstone, "Ideal Spatial Adaptation by Wavelet Shrinkage", Department of Statistics, Stanford University, USA. April 1993. [10] I. Daubechies, Ten Lectures on Wavelets, New York: SIAM, 1992. [11] S. G. Mallat, A Wavelet Tour of Signal Processing, Academic Press, 2nd Edition, 1999. [12] S. G. Mallat S. G., A Theory of Multiresolution Signal Decomposition: The Wavelet Representation, IEEE Transactions on Pattern Analysis and Machine Intelligence, July 1989, 11(7), pp. 674-693. [13] B. Antoniadis and Sapatinas, "Wavelet Estimators in Non-parametric Regression: A Comparative Simulation Study", Laboratoire IMAG-LMC, University Joseph Fourier, France, 2001. [14] M. Jansen, A. Bultheel, "Asymptotic Behavior of The Minimum Mean Squared Error Threshold for Noisy Wavelet Coefficients of Piece-wise Smooth Signals", Katholieke Universiteit Leuven, Department of Computer Science, Belgium, October 4, 1999. [15] The MathWorks, Inc. "Wavelet Toolbox. Matlab: The Language of the Technical Computing". Version 6.0.0.88, Release 12, September 22, 2000. [16] S. Sardy, A. Antoniadis, and P. Tseng, "Automatic Smoothing with Wavelets for a Wide Class of Distributions", Journal of Computational and Graphical Statistics, January 2003. [17] M. Hernandez-Diaz, Analisis Comparativo de Algoritmos para Reducción de Ruido en Señales Utilizando Wavelets, BSc

Proceedings of the 15th International Conference on Electronics, Communications and Computers (CONIELECOMP 2005) 0-7695-2283-1/05 $20.00 © 2005 IEEE