**3. Discrete wavelet transform**

## **3.1 Background**

132 Advances in Wavelet Theory and Their Applications in Engineering, Physics and Technology

surplus smoothing may noticeably broaden the solution and conceal its important features (a parameter α that is too high leads to a solution that is not very sensitive to noise, but it is very far from the measured data). It is therefore necessary to find the optimum α and, hence, a smoothing factor at which the solution to problem (6) is wellstabilized but still close to the real distribution. Moreover, in the iterative algorithms, if the regularization parameter is not well calculated, the oscillations created at iteration n are amplified at iteration n + 1, which degenerates the final solution. The value of α, the regularization parameter associated with D, is very important in the regularization process, and its value determines the quality of the final solution. This can easily be understood if one analyzes the generalized matrix H+. As α increases, the matrix (H+)-1 tends toward a diagonal matrix, while the vector HTy, which is broadened in comparison to the initial data vector y due to the multiplication by the transposed distortion matrix, remains unchanged. As a result, with an increase of α, the shape of the solution tends to HTy, i.e., to the considerably broadened initial data. Figure 3 shows an example of the evolution of the spectrum of the matrix H+ for various values of the regularization

According to Fig. 3, the determination of the classical regularization parameter αc for the SIMS profile led to a value of 5.9290×10-4. For this value, the spectrum of the generalized matrix H+ is oscillatory (the matrix is not well-conditioned), leading to an unstable solution. In order to stabilize the system more, Mancina et al (Mancina et al, 2000) proposed multiplying the regularization parameter by a positive factor. The multiplication of this parameter by a factor K (K = 10, 102, 103) leads to increasingly regularized matrix H+ (Fig. 2). Nevertheless, this multiplication is arbitrary and it is not

Fig. 3. Evolution of the spectrum of the matrix *H*+ according to kαc (K =1, 10, 102, 103).

parameter α.

based on any physical support

Wavelet theory is widely used in many engineering disciplines (Rashed et al, 2007; Garcia-Talavera et al, 2003; Starck et al, 2003; Jammal et al, 2004; Rucka et al, 2006), and it provides a rich source of useful tools for applications in time-scale types of problems. The attention to study of wavelets becomes more attractive when Mallat (Mallat, 1989) established a connection between wavelets and signal processing. Discrete wavelet transform (DWT) is an extremely fast algorithm that transforms data into wavelet coefficients at discrete intervals of time and scale, instead of at all scales. It is based on dyadic scaling and translating, and it is possible if the scale parameter varies only along the dyadic sequence (dyadic scales and positions). It is basically a filtering procedure that separates high and low frequency components of signals with high-pass and low-pass filters by a multiresolution decomposition algorithm (Mallat, 1989). Hence, the DWT is represented by the following equation:

$$\mathcal{W}(j,k) = \sum\_{j} \sum\_{k} \mathcal{y}(k) 2^{-j} \mathcal{Y}(2^{-j}n - k), \tag{11}$$

where y is discretized heights of the original profile measurements, ψ is the discrete wavelet coefficients, and n is the sample number. The translation parameter k determines the location of the wavelet in the time domain, while the dilatation parameter j determines the location in the frequency domain as well as the scale or the extent of the space-frequency localization.

The DWT analysis can be performed using a fast, pyramidal algorithm by iteratively applying low-pass and high-pass filters, and subsequent down-sampling by 2 (Mallat, 1989). Each level of the decomposition algorithm then yields to low-frequency components of the signal (approximations) and high-frequency components (details). This is computed with the following equations:

$$\text{y}\_{\text{low}}[\text{k}] = \sum\_{n} y[n] f[\text{2}k - n]. \tag{12}$$

Multi-Scale Deconvolution of Mass Spectrometry Signals 135

Fig. 4. (a) *Sym4*: scaling function, wavelet function, and the associated filters. (b) *Coif3*:

Due to the compression and dilatation properties of the wavelets in representing a signal, wavelet-based filters can easily follow the sharp edges of the input signal. In other words, they restore any discontinuity in the input signal, or, in terms of the frequency domain interpretation, they pass high frequency components of the input signal. This is a very appealing feature of the wavelet-based methods in many applications, such as finding the location of discontinuities and abrupt changes in a signal. However, this feature may have adverse effects, especially when one wants to get rid of impulsive noise (outliers and gross

We notice that the lower level (high-frequency) wavelet components are similar to a random process, while the higher level (low-frequency) ones are not [Fig. 5(a)]. The noise in SIMS analysis is Gaussian, and one considers that, if there was no signal but the noise alone, the variance of the details would decrease by a factor of 2 at each resolution. Analysis at each level of detail (from small to large) separately on the same signal is shown in Fig. 5(b).

scaling function, wavelet function, and the associated filters.

errors).

$$\text{ything}[\mathbf{k}] = \sum\_{n} y[n] \text{lg}[2k - n],\tag{13}$$

where ylow[k] and yhigh[k] are the outputs of the low-pass (f) and high-pass (g) filters, respectively, after down sampling by 2. Due to down-sampling during decomposition, the number of resulting wavelet coefficients at each level is exactly the same as the number of input points for this level. It is sufficient to keep all detail coefficients and the final approximation coefficients (at the roughest level) in order to reconstruct the original data.

The approximation and details at the resolution 2-(j+1) are obtained from the approximation signal at resolution 2-j. In the matrix formalism, eqs. (12) and (13) can be written as

$$\mathbf{y}\_a^{(j+1)} = \mathbf{F} \ \mathbf{y}\_a^{(j)} \ \ \ \ \ \mathbf{y}\_d^{(j+1)} = \mathbf{G} \ \ \ \mathbf{y}\_a^{(j)} \ \tag{14}$$

where F and G are Toeplitz matrices constructed from the filters f and g, respectively.

The reconstruction algorithm involves up-sampling, i.e., inserting zeros between data points, and filtering with dual filters. By carefully choosing filters for the decomposition and reconstruction that are closely related, we can achieve perfect reconstruction of the original signal in the inverse orthogonal wavelet transform (Daubechies, 1990). The reconstructed signal is obtained from eq. (14) by

$$\tilde{y} = \tilde{F}y\_a^{(l)} + \tilde{G}y\_d^{(j)}\ ,\ \mathbf{j} = \mathbf{1}\,,\dots\ \mathbf{J}\,,\tag{15}$$

where *F* and *G* are Toeplitz matrices constructed from the reconstruction filters *f* and *g* , respectively. For a general introduction to discrete wavelet transform and filter banks, the reader is referred to refs. (Mallat, 1989; Daubechies, 1990).

The Mallat algorithm (Mallat, 1989) is a fast linear operation that operates on a data vector whose length is an integer power of two, transforming it into numerically different vectors of the same length. Many wavelet families are available. However, only orthogonal wavelets (such as Haar, Daubichies, Coiflet, and Symmlet wavelets) allow for perfect reconstruction of a signal by inverse discrete wavelet transform, i.e., the inverse transform is simply the transpose of the transform. Indeed, the selection of the most appropriate wavelet is based on the similarity between the derivatives of the signal and the number of wavelet vanishing moments. In practice, wavelets with a higher number of vanishing moments give higher coefficients and more stable performance. In this study, we restrict ourselves to Symlet and Coiflet families; after some experimentation, we chose "Sym4" wavelet with four vanishing moments and "Coif3" wavelet with three vanishing moments. Figure 4 shows the wavelet function, scaling function and the four filters of the wavelets "Sym4" and "Coif3". The decomposition on a wavelet basis (to the level 5) of a SIMS profile containing four delta-layers of boron in silicon to approximation and detail signals is illustrated in Fig. 5.

signal (approximations) and high-frequency components (details). This is computed with

*<sup>y</sup> nf k n* (12)

*<sup>y</sup> <sup>n</sup> <sup>g</sup> k n* (13)

*aa a <sup>d</sup> <sup>y</sup> Fy y Gy* , (14)

, j = 1,…, J, (15)

and *g* ,

*n*

*n*

signal at resolution 2-j. In the matrix formalism, eqs. (12) and (13) can be written as

where F and G are Toeplitz matrices constructed from the filters f and g, respectively.

( )*J* ( )*j <sup>a</sup> <sup>d</sup> y Fy Gy*

respectively. For a general introduction to discrete wavelet transform and filter banks, the

The Mallat algorithm (Mallat, 1989) is a fast linear operation that operates on a data vector whose length is an integer power of two, transforming it into numerically different vectors of the same length. Many wavelet families are available. However, only orthogonal wavelets (such as Haar, Daubichies, Coiflet, and Symmlet wavelets) allow for perfect reconstruction of a signal by inverse discrete wavelet transform, i.e., the inverse transform is simply the transpose of the transform. Indeed, the selection of the most appropriate wavelet is based on the similarity between the derivatives of the signal and the number of wavelet vanishing moments. In practice, wavelets with a higher number of vanishing moments give higher coefficients and more stable performance. In this study, we restrict ourselves to Symlet and Coiflet families; after some experimentation, we chose "Sym4" wavelet with four vanishing moments and "Coif3" wavelet with three vanishing moments. Figure 4 shows the wavelet function, scaling function and the four filters of the wavelets "Sym4" and "Coif3". The decomposition on a wavelet basis (to the level 5) of a SIMS profile containing four delta-layers of boron in silicon to approximation and detail

where *F* and *G* are Toeplitz matrices constructed from the reconstruction filters *f*

reader is referred to refs. (Mallat, 1989; Daubechies, 1990).

where ylow[k] and yhigh[k] are the outputs of the low-pass (f) and high-pass (g) filters, respectively, after down sampling by 2. Due to down-sampling during decomposition, the number of resulting wavelet coefficients at each level is exactly the same as the number of input points for this level. It is sufficient to keep all detail coefficients and the final approximation coefficients (at the roughest level) in order to reconstruct the original data. The approximation and details at the resolution 2-(j+1) are obtained from the approximation

( 1) ( ) ( 1) ( ) , *jj j <sup>j</sup>*

The reconstruction algorithm involves up-sampling, i.e., inserting zeros between data points, and filtering with dual filters. By carefully choosing filters for the decomposition and reconstruction that are closely related, we can achieve perfect reconstruction of the original signal in the inverse orthogonal wavelet transform (Daubechies, 1990). The reconstructed

ylow [k] = [ ] [2 ],

yhigh [k] = [ ] [2 ],

the following equations:

signal is obtained from eq. (14) by

signals is illustrated in Fig. 5.

Fig. 4. (a) *Sym4*: scaling function, wavelet function, and the associated filters. (b) *Coif3*: scaling function, wavelet function, and the associated filters.

Due to the compression and dilatation properties of the wavelets in representing a signal, wavelet-based filters can easily follow the sharp edges of the input signal. In other words, they restore any discontinuity in the input signal, or, in terms of the frequency domain interpretation, they pass high frequency components of the input signal. This is a very appealing feature of the wavelet-based methods in many applications, such as finding the location of discontinuities and abrupt changes in a signal. However, this feature may have adverse effects, especially when one wants to get rid of impulsive noise (outliers and gross errors).

We notice that the lower level (high-frequency) wavelet components are similar to a random process, while the higher level (low-frequency) ones are not [Fig. 5(a)]. The noise in SIMS analysis is Gaussian, and one considers that, if there was no signal but the noise alone, the variance of the details would decrease by a factor of 2 at each resolution. Analysis at each level of detail (from small to large) separately on the same signal is shown in Fig. 5(b).

Multi-Scale Deconvolution of Mass Spectrometry Signals 137

Wavelets have multiscale and local properties that make them very effective in analyzing the class of locally varying signals. Together the locality and multiscale properties enable the wavelet transform to efficiently match signals organized into levels or scales of localized variations. Thus DWT transforms the noisy signal in the wavelet domain, and by denoising we obtain a sparse representation with a few large dominating coefficients (Donoho et al, 1994, 1995). A large part of the wavelet coefficients does not carry significant information [see absolute wavelet coefficients for n = 1 to 5 in Fig. 5(b)]. We select the significant ones by a thresholding procedure, which is addressed in the

Noise is a phenomenon that affects all frequencies. Since the useful signal tends to dominate the low-frequency components, it is expected that the majority of high-frequency components above a certain level are due to noise. In the wavelet decomposition of signals, as has been described, the filter h is an averaging or smoothing filter, while its mirror g produces details. With the exclusion of the last remaining smoothed components, all wavelet coefficients in the final decomposition correspond to details. If the absolute value of a detail component is small (or set to zero), the general signal does not change much. Therefore, thresholding the wavelet coefficients is a good way of removing unimportant or undesired (insignificant) details from a signal. Thresholding techniques are successfully used in numerous data-processing domains, since in most cases a small number of wavelet coefficients with large amplitudes preserve most of the information about the original data

 Decomposition: Select the level N and type of wavelets and then determine the coefficients of SIMS signal by DWT. For wavelet denoising, we must decide from many possible selections, such as the type of mother wavelet, the decomposition levels, and the values of thresholds in the next step. In this study, decomposition at level 5 has been

 Thresholding: Estimating threshold values is based upon the analytical and empirical methods. For each level from 1 to 5, we use the estimated threshold values and set the detail coefficients below the threshold values to zero. Based on knowledge of the wavelet analysis in the data set, we use objective criteria to determine threshold values. Basically, the choice of mother wavelet appears not to matter much, while the values of thresholds do. Therefore, setting the values of the threshold is a crucial topic. According to the analysis described, we set threshold values based on the properties of SIMS data

 Reconstruction: We reconstruct the denoised signal using the original approximation coefficients of level N and the modified detail coefficients of levels from 1 to N by the

Wavelet denoising methods generally use two different approaches: hard thresholding and soft thresholding. The hard thresholding philosophy is simply to cut all the wavelet coefficients below a certain threshold to zero. However, soft thresholding reduces the value (referred to as "shrinkage") of wavelet coefficients towards zero if they are below a certain

A basic wavelet-based denoising procedure is described in the following:

following section.

set.

used.

sets.

inverse DWT.

**3.2 Denoising by wavelet coefficients thresholding** 

Fig. 5. Wavelet decomposition of SIMS profile measured at 8.5 keV/O2+, 38.1 rad. The wavelet used is *Sym4*; the decomposition level is 5.

(a) The original measured profile with the different approximation signals from level 1 to 5.

(b) Detail signals from level 1 to 5 with absolute wavelet coefficients

Fig. 5. Wavelet decomposition of SIMS profile measured at 8.5 keV/O2+, 38.1 rad.

(a) The original measured profile with the different approximation signals from level 1 to 5.

The wavelet used is *Sym4*; the decomposition level is 5.

(b) Detail signals from level 1 to 5 with absolute wavelet coefficients

Wavelets have multiscale and local properties that make them very effective in analyzing the class of locally varying signals. Together the locality and multiscale properties enable the wavelet transform to efficiently match signals organized into levels or scales of localized variations. Thus DWT transforms the noisy signal in the wavelet domain, and by denoising we obtain a sparse representation with a few large dominating coefficients (Donoho et al, 1994, 1995). A large part of the wavelet coefficients does not carry significant information [see absolute wavelet coefficients for n = 1 to 5 in Fig. 5(b)]. We select the significant ones by a thresholding procedure, which is addressed in the following section.

## **3.2 Denoising by wavelet coefficients thresholding**

Noise is a phenomenon that affects all frequencies. Since the useful signal tends to dominate the low-frequency components, it is expected that the majority of high-frequency components above a certain level are due to noise. In the wavelet decomposition of signals, as has been described, the filter h is an averaging or smoothing filter, while its mirror g produces details. With the exclusion of the last remaining smoothed components, all wavelet coefficients in the final decomposition correspond to details. If the absolute value of a detail component is small (or set to zero), the general signal does not change much. Therefore, thresholding the wavelet coefficients is a good way of removing unimportant or undesired (insignificant) details from a signal. Thresholding techniques are successfully used in numerous data-processing domains, since in most cases a small number of wavelet coefficients with large amplitudes preserve most of the information about the original data set.

A basic wavelet-based denoising procedure is described in the following:


Wavelet denoising methods generally use two different approaches: hard thresholding and soft thresholding. The hard thresholding philosophy is simply to cut all the wavelet coefficients below a certain threshold to zero. However, soft thresholding reduces the value (referred to as "shrinkage") of wavelet coefficients towards zero if they are below a certain value. For a certain wavelet coefficient k on scale j, the thresholded details coefficients are given by

$$\hat{y}\_d(k) = \text{sign}\left\|y(k)\right\| - \mathcal{X}\big|\,. \tag{16}$$

Multi-Scale Deconvolution of Mass Spectrometry Signals 139

Fig. 6. (a) Original SIMS profile superimposed on the denoised profile. (b) Thresholded

denoising by wavelet coefficient thresholding are as follows.

asymptotically closer towards the true signal.


**4. Multiscale deconvolution (MD): The proposed algorithms** 

deconvolution scheme is constructed by the following steps:


that it will lead to aberrant results.

**signal as model of solution** 

Because the few largest wavelet coefficients preserve almost the entire energy of the signal, shrinkage reduces noise without distorting the signal features. The main results after


Finally, we note that the result obtained [Fig. 6(a)] is of both good smoothness and regularity. Thus, we may exploit this advantage in deconvolution procedure without fearing

**4.1 First algorithm: Tikhonov-Miller regularization with a denoisy and deconvoluted** 

We have seen in § 2.2 eq. 10 that Mancina (Mancina, 2000) proposed to reiterate the algorithm of Barakat (Barakat et al, 1997) and to use a pre-deconvoluted signal as model of the solution with sufficient regularization. The accuracy of the solution is referred to the accuracy of the model, which suggests a reasonable formulation. It is obvious that a significant lack of precision in the a priori model leads to an error restoration more important than the usual one without the model. Moreover, if the pre-deconvoluted signal is a noisy signal (which is the case for SIMS signals) or contains aberrations, the iterative process worsens these aberrations and the result is an oscillatory signal. For this reason, it is important to remove noise components from the signal (the model of solution). The idea is to introduce a denoisy and deconvoluted signal as model of solution in Barakat's approach, which constitutes our first contribution in this field (Boulakroune, 2008). The first proposed

wavelet coefficients.

where the function "sign" returns the sign of the wavelet coefficient, and λ is the threshold value. In the case of Gaussian white noise (which is the kind of noise in SIMS analysis), Donoho and Johnstone (Donoho et al, 1994, 1995) modeled this threshold by

$$
\lambda = \sigma \sqrt{2 \log(N)}\,\,\,\,\,\tag{17}
$$

where N is the number of the observed data points, and σ is the standard deviation of noise. This standard deviation, in the case of white and Gaussian noise, is estimated by

$$
\hat{\sigma} = \text{median}\left( \left| cd^{(1)}(k) \right| \Big/ 0.6754 \right) \tag{18}
$$

where median(cd(1)(k)) is the median value of detail coefficients at the first level of decomposition which can be attributed to noise. After thresholding, the reconstructed signal of eq. (15) becomes:

$$\tilde{y} = \tilde{F}y\_a^{(I)} + \tilde{G}\hat{y}\_d^{(j)}\ ,\ \mathbf{j} = \mathbf{1}\ ,\ \dots \ \mathbf{J}\ .\tag{19}$$

By using this process, high-frequency components above a certain threshold can be removed. A raw SIMS profile and corresponding denoised profile are shown in Fig. 6(a). In particular, the figure shows that low-frequency components, which usually represent the main structure of the signal, are separated from high-frequency components. These preliminary results demonstrate the superior capabilities of the wavelet approach to SIMS profiles analysis over traditional techniques.

In the analysis of SIMS data, we find that most wavelet coefficients at high-frequency levels from 1 to 4 [see Fig. 5(b)], can be mostly ignored. However, we must be very cautious when manipulating the low-frequency components to keep as many true coefficients as possible after thresholding. According to the exploratory data analysis in the beginning of this section, we select a threshold value large enough to ignore most of the wavelet coefficients at levels 1-4, which represent the noise signals, especially in the beginning and the end of the profile. The denoising results show the good performance of wavelet application and exploratory data analysis.

The remaining wavelet coefficients after shrinkage are less than one tenth of those of the original SIMS. These thresholded wavelet coefficients [those "stuttered" 2n times on level n are concentrated in the zone where the signal is too noisy, see Fig. 6(b)] give us an idea of the remaining details in the approximation (denoised signal) of the original signal, which are higher than the determined threshold (significant details). For example, the estimated threshold of the previous SIMS signal, obtained using soft universal shrinkage [eq. (17)], is λ = 55.7831 cps. The estimated level of noise, using eq. (18), gives a signal-to-noise ratio (SNR) = 40.9212 dB.

value. For a certain wavelet coefficient k on scale j, the thresholded details coefficients are

ˆ ( ) *<sup>d</sup> y k* sign *y k*( )

Donoho and Johnstone (Donoho et al, 1994, 1995) modeled this threshold by

 

This standard deviation, in the case of white and Gaussian noise, is estimated by

where the function "sign" returns the sign of the wavelet coefficient, and λ is the threshold value. In the case of Gaussian white noise (which is the kind of noise in SIMS analysis),

where N is the number of the observed data points, and σ is the standard deviation of noise.

where median(cd(1)(k)) is the median value of detail coefficients at the first level of decomposition which can be attributed to noise. After thresholding, the reconstructed signal

By using this process, high-frequency components above a certain threshold can be removed. A raw SIMS profile and corresponding denoised profile are shown in Fig. 6(a). In particular, the figure shows that low-frequency components, which usually represent the main structure of the signal, are separated from high-frequency components. These preliminary results demonstrate the superior capabilities of the wavelet approach to SIMS

In the analysis of SIMS data, we find that most wavelet coefficients at high-frequency levels from 1 to 4 [see Fig. 5(b)], can be mostly ignored. However, we must be very cautious when manipulating the low-frequency components to keep as many true coefficients as possible after thresholding. According to the exploratory data analysis in the beginning of this section, we select a threshold value large enough to ignore most of the wavelet coefficients at levels 1-4, which represent the noise signals, especially in the beginning and the end of the profile. The denoising results show the good performance of wavelet application and

The remaining wavelet coefficients after shrinkage are less than one tenth of those of the original SIMS. These thresholded wavelet coefficients [those "stuttered" 2n times on level n are concentrated in the zone where the signal is too noisy, see Fig. 6(b)] give us an idea of the remaining details in the approximation (denoised signal) of the original signal, which are higher than the determined threshold (significant details). For example, the estimated threshold of the previous SIMS signal, obtained using soft universal shrinkage [eq. (17)], is λ = 55.7831 cps. The estimated level of noise, using eq. (18), gives a signal-to-noise ratio

( ) ( ) ˆ *<sup>J</sup> <sup>j</sup>*

, (16)

2log( ) *N* , (17)

ˆ median (1) ( ( ) 0.6754 *cd k* , (18)

*<sup>a</sup> <sup>d</sup> <sup>y</sup> Fy Gy* , j = 1, …, J. (19)

given by

of eq. (15) becomes:

exploratory data analysis.

(SNR) = 40.9212 dB.

profiles analysis over traditional techniques.

Fig. 6. (a) Original SIMS profile superimposed on the denoised profile. (b) Thresholded wavelet coefficients.

Because the few largest wavelet coefficients preserve almost the entire energy of the signal, shrinkage reduces noise without distorting the signal features. The main results after denoising by wavelet coefficient thresholding are as follows.


Finally, we note that the result obtained [Fig. 6(a)] is of both good smoothness and regularity. Thus, we may exploit this advantage in deconvolution procedure without fearing that it will lead to aberrant results.
