Determining Decomposition Levels for Wavelet Denoising Using Sparsity Plot

We present a method to select decomposition levels for noise thresholding in wavelet denoising. It is essential to determine the accurate decomposition levels to avoid inadequate noise reduction and/or signal distortion by noise thresholding. We introduce the concept of sparsity plot that captures the abrupt transition from noisy to noise-free Detail component, readily revealing the cut-off for the maximum decomposition levels. The method uses the sparsity parameter to determine the noise presence in each detail component and measures the magnitude change in the sparsity values to distinguish between noisy and noise-free Detail components. The method is tested on both model and experimental signals, and proves effective for various signal lengths and types, as well as different Signal-to-Noise Ratios (SNRs). The method can be embedded with any wavelet denoising method to improve its performance. The code is available via GitHub and denoising.cornell.edu, as well as the corresponding author’s group website (http://signalsciencelab.com).


I. INTRODUCTION
Wavelet denoising [1]- [8] is a widely used approach to remove noise and is applied in a large number of fields such as magnetic resonance spectroscopy [9], speech recognition [10], time series analysis [11], traffic volume prediction [12], wireless ultraviolet communication [13], time-difference-of-arrival localization [14], fluorescence imaging [15], surface electromyography [16], computed tomography [17], high resolution surface scans [18], nanopore and ion-channel signals [19], and simulated chest wall motion [20]. Effective denoising relies heavily on selecting the decomposition levels to apply noise thresholding for retrieving the desired signal. Selecting fewer wavelet components can lead to inadequate denoising, whereas noise thresholding wavelet components with large signal presence can lead to signal distortion. Hence, it is essential to accurately identify decomposition level where wavelet components contain noise.
Currently, decomposition levels are either selected arbitrarily via trial and error or using a predefined cut-off [21]. The current wavelet denoising methods lack a general objective The associate editor coordinating the review of this manuscript and approving it for publication was Hengyong Yu . criterion for selecting decomposition levels that can be applied to a wide range of signals. In this paper, we address the problem by introducing ''Sparsity Plot'' to select the appropriate decomposition levels. The Sparsity Plot calculates the sparsity of each Detail component and plots its value across decomposition levels, revealing the level at which the sparsity value abruptly changes from low to high. The sparsity parameter reflects the presence of noise in a Detail component, with small value reflecting high noise presence and large value representing high signal presence. The abrupt switch from low to high value informs the decomposition level cut-off where noise containing Detail components separate from the signal containing Detail components.
We demonstrate through theorems that: 1) The sparsity of both noise-free and noisy Detail components gradually increases from decomposition level one to higher levels.
2) The sparsity of noise-free Detail component is greater than noisy Detail component at a given decomposition level.
3) The sparsity of noise-free Detail component at decomposition level j + 1 is substantially greater than noisy Detail component at j. Using model and experimental data, we show how to utilize the Sparsity plot and its change in values between decomposition levels to select the Detail components for noise thresholding. We also present an objective cut-off for automation of selecting the decomposition level. We provide codes to carry out the decomposition level selection procedure as well as the codes that we used to generate model data and generate the sparsity plot.
The method can be easily inserted into any wavelet denoising procedure to improve its performance.
The paper is organized as follows. In section II, we provide a mathematical formulation of the sparsity parameter and theorems that demonstrate the existence of a sparsity-based characteristic marker separating noisy and noise-free Detail components. In section III, we illustrate the procedure for selecting the decomposition level. In section IV, we describe the model and experimental data that is used to test our new decomposition level selection method. Section V discusses the results and Section VI summarizes the method and findings in the Conclusion.

II. MATHEMATICAL FORMULATION
We formulate the sparsity parameter to calculate the noise presence in each Detail component. Using the vector of detail coefficients (i.e., Detail component), D j , from wavelet decomposition at each of the j = 1, . . . , N levels, we compute the sparsity component for each level, S j , as follows: where q j is the length of the Detail component (i.e. detail coefficient vector) at a given decomposition level j. We now develop three arguments: Theorem 1: Given S j and S j+1 being the sparsity values of the Detail component at decomposition level j and j + 1, respectively, S j < S j+1 : Proof: By equation 1, we rewrite this inequality as: For a given signal f [t m ], the Detail component (D j ) at decomposition level j can be computed as: where ψ is the wavelet function and p is the length of the input signal f [t m ]. Incorporating equation 4 into equation 3, we can rewrite the inequality as: As 2 j/2 and 2 j+1/2 are common in the numerator and denominator in their respective terms, they can be cancelled out and we can rewrite the inequality as: Since ψ(2 j t m − k) and ψ(2 j+1 t m − k)| on their respective side of the inequality represent index values, the maximum magnitude associated with the summation As a result, we can conclude that the numerators of each term are approximately equal and rewrite the inequality as: VOLUME 9, 2021 Suppose we let: Using the definition in equations 8 and 9, we can rewrite the inequality in equation 7 as, From equation 10, we now know that this inequality holds if and only if d j > d j+1 . Thus, we want to show that d j > d j+1 : We note that ψ(2 j t m − k) and ψ(2 j+1 t m − k) are same function with different index values sampled, and hence have same values for a given index. In ψ(2 j t m − k), twice the number of values are sampled from the function compared to is sampling every 4 th data point. This implies that for a given input signal length p, Replacing equation 12 in equation 11, we get As a result, the term on the RHS of our inequality sums roughly half the number of equal values as that on the LHS. Therefore, we can say that d j , is greater than d j+1 , and conclude that our claim S j < S j+1 is true.
Theorem 2: Given that S S j and S S+N j are the sparsity value of the noise-free signal and noisy signal at level j, respectively, S S j > S S+N j : Proof: Using equation 1, we rewrite this inequality as: where is n[k] is the noise value at k th index. For SNR 1, signal is dominant over noise. Thus, for the detail coefficients associated with maximum signal value: max|D j | max| ± n|, we can conclude that: Thus, equation 15 can be rewritten as, 1 Suppose we let: The inequality in equation 17 holds if and only if We can divide a noisy Detail component into three parts: 1) signal coefficients 2) signal coefficients + noise coefficients 3) signal coefficients − noise coefficients The length and contribution of each will vary depending upon the type of signal and the noise level. The LHS ( 20 can be expanded into the three parts and rewritten as: where π j is the length of the detail coefficient vector with positive noise at decomposition level j and ρ j − (1 + π j ) is the length of the detail coefficient vector with negative noise at decomposition level j.
Since we assume that the signal is sparse and noise is random Gaussian-centered around mean, as well as i.i.d, the summation of . We can rewrite the equation 21 as: Because the signal is sparse in the wavelet domain, the RHS in equation 20 can written as, Next, rewrite the RHS of our inequality as:

Replacing equations 22 and 23 in equation 20, we show that
Since the LHS has a noise term in addition to the summation of the detail coefficient vector, we know that it must be larger than the RHS. Therefore, d S+N > d S , and conclude that initial claim S S j > S S+N j is true.
By Theorem 2, we show that: As the noise-free Detail component is at the decomposition level j + 1 and the noisy Detail component is at the j th level, we combine equations 25 and 26 and show that: The gains associated with the two inequalities in equation 27 is additive, where increase in the sparsity value associated higher decomposition levels is coupled with the higher sparsity value of noise-free Detail component. Hence, the equation 27 can be rewritten as, Hence, equation 27 proves that Theorem 3 holds and reveals that an abrupt increase in the sparsity value between decomposition levels is a characteristic marker separating the noisy and noise-free Detail components.  components saturate. This suggests that the wavelet filter coefficients dominate the Detail component, instead of the signal and/or noise coefficients. It can be measured by finding the ratio of lengths of Detail component and the wavelet function filter as follows,

III. SELECTION OF DECOMPOSITION LEVEL
where |D j | is the length of the Detail component at decomposition level j and L filter is the length of the wavelet filter (i.e., high pass filter associated with Detail components).
To avoid selecting the Detail components dominated by the wavelet filter coefficients, we select the ratio cut-off of 1.5, where the Detail components above this value are only considered. In other words, maximum decomposition level N is defined as, As the wavelet transform is obtained by performing convolution, the cut-off value of 1.5 ensures that the signal and/or noise still has substantial contribution in the Detail component under study.  Fig. 1a) (B) Plot of Change in Sparsity values with respect to previous decomposition level for noise-free signal (cf. Fig. 1a). (C) Sparsity plots of noise-free (blue, 1a) and noisy signal (red, (cf. Fig. 1c)) at SNR = 10 (D) Plot of Change in Sparsity values from j to j + 1 decomposition level for j th noise-free signal (blue, cf. Fig. 1a) and j th noisy signal (black, cf. Fig. 1c) at SNR = 10 with respect to j + 1 th sparsity value of the noise-free signal.

B. STEP 2: SPARSITY PLOT
We calculate the sparsity of each Detail component up to the maximum decomposition level selection. This sparsity value enables us to separate noisy Detail components from noisefree ones. The Detail components with high noise presence have low sparsity values and vice versa [21]. The sparsity value is calculated as, where S j is the sparsity value and 1 ≤ j ≤ N . The sparsity plot displays the sparsity values for N levels. As described through Theorems 1 to 3, the sparsity values should increase from level one to level N , albeit at different rates and with occasional abrupt increases in values between the levels.

C. STEP 3: CHANGE IN SPARSITY
To capture the rate of increase between the levels, we measure the change in the sparsity parameter at the adjacent levels. The change in sparsity values SC j is given as, The plot of SC j values informs about the decomposition level of abrupt change that distinguishes between Detail components with and without noise presence. We evaluate this change for each decomposition level, starting from level two (at level one, we have SC 1 = 0 since there is no prior level to measure sparsity value). We plot this value for each level to render the Change in Sparsity Plot. These plots are effective visual tools allowing us to determine the cut-off between the noisy and noise-free component in signals and can readily reveal the cut-off for the maximum decomposition levels needed for denoising.

D. SPARSITY CUT-OFF
We can determine the decomposition level k, where D 1 to D k contain noise by using both subjective and object approaches. In the subjective approach, the change in the sparsity values SC j can be plotted and the abrupt change in values can be easily observed (see Figure 4 (B3-D3) and Figure 5 (A3-B3)). The decomposition level j is selected where the first abrupt change happens between j and j+1. In the objective approach, a threshold change of 5% with respect to 1 is selected as cut-off between the noisy and noise-free Detail components. The change of 5% is substantial to capture the abrupt increase in sparsity towards all signal components. The decomposition level j−1 is selected where SC j > 0.05 in the first occurrence.

A. MODEL DATA GENERATION
We first generated model data using the Lorenz in function (see attached code). This noise-free data contains 1024 values with identical peak heights of approximately 0.2753 with each peak spaced out at intervals of around 185 index values and negative and positive peaks separated by 13 index values. The values are symmetric with negative and positive values centered around zero. As seen in Figure 4, this data was used as our underlying signal, to which we added Gaussian noise with varying signal-to-noise ratios (SNRs). We added white Gaussian noise to our signal using the awgn() function in MATLAB and validated its SNR as: Thus, we generated noise with SNR values of ten, fifty, and one hundred (see Figure 4 (B1-D1)) and added these each to the model data (see Figure 4 (B2-D2)).

B. EXPERIMENTAL DATA COLLECTION
We use continuous wave-Electron Spin Resonance (cw-ESR) spectroscopy for collecting experimental data to test the results. Cw-ESR is used extensively to study the dynamics and structure of biomolecules, and is the most commonly used ESR technique [22]. The cw-ESR spectrum is acquired in the magnetic field (B 0 ) domain, which is swept, i.e., B 0 = B 0 (t) in a linear fashion. For Example 1, the ESR experiments were performed at 20 • C on a commercial spectrometer (BRUKER ELEXYS-II E500) at a standard microwave frequency of 9.4 GHz (X-band) corresponding to a dc magnetic field of 0.34 Tesla. In these experiments the microwave frequency is held constant, while the magnetic field is swept through resonance to obtain the ESR spectrum [22]. The sample consisted of 4 µL of a 100 µM aqueous solution of the commonly used spin-probe molecule Tempol (4-Hydroxy-2,2,6,6-Tetramethylpiperidine 1-oxyl) [22]- [25]. It was placed in a glass capillary of 0.8 mm ID, which was then introduced into the microwave cavity situated between the pole caps of the dc magnet. The magnetic field was then swept over a range of 60 G corresponding to the resonant spectral range which took 2 minutes, and a 82 ms time constant was used. The spectral data consisted of 4096 points along the magnetic-field sweep. In addition, small coils placed at the sides of the resonator provided a small magnetic field modulation of ±0.02 G at a frequency of 100 kHz. The 100 kHz modulated ESR signal was detected with a lock-in detector at this frequency, providing the first derivative of the absorption signal [22]. Low power (0.2 mW) microwave radiation was used to avoid saturating the ESR signal. Multi-scan experiments were performed with a delay of 4 s between scans. The results of these scans were then averaged. The ESR spectrum displayed in Figure 5 was obtained under very similar conditions, also for the Tempol spin-probe.
In Example 2 shown in Figure 5, an ESR spectrum was obtained under different conditions to provide a more complex spectrum for denoising. It was obtained on a home-built (ACERT) 95 GHz ESR spectrometer [25] with a dc magnetic field of 3.3 Tesla at 25 • C. The sample here contained ca. 5 µL of phospholipid vesicles doped with 0.5% of a lipid spin label: 16-PC (1-acyl-2-[16-(4,4-dimethyloxazolidine-N-oxyl)stearoyl]-sn-glycero-3-phosphocholine) in the fluid phase that had been suspended in water. It was placed in a disc-like sample holder utilized for millimeter-wave ESR methodology [25]. The acquisition parameters were: sweep width of 250 G, sweep time of 2 min with a time constant of 100 ms. The millimeter-wave power was 16 mW and the spectrum consists of 512 points. The field modulation parameters were: 6 G modulation amplitude and 100 kHz modulation frequency. The time between scans was 3 s.

A. MODEL DATA
Using SNR values of ten, fifty, and one hundred, we applied our method and calculated the sparsity and change in sparsity for each decomposition level. The plots of these metrics are shown in Figure 4. We add a dotted vertical line in both plots to distinguish between noisy and noise-free components of the data and to reveal the decomposition level selection cut-off for each example. We then applied our method and calculated the sparsity for each decomposition level.  In Figure 4(B3), we show the sparsity plot for the model data plus a large amount of noise. We add a dotted vertical line to distinguish between noisy and noise-free components of the data and to reveal the decomposition level selection cut-off for each example.
In this first case, we see that there is a drastic change from all noise-containing wavelet coefficients to noise-free coefficients between levels 3 and 4. In Figure 4(C3) we show the sparsity plot for the model data plus some noise. We observe that there exists a similar change between levels 3 and 4, while the signal plus noise line approaches the signal line quicker than in the very noisy case. Finally, we show the sparsity plot for the model data plus minimal noise in Figure 4(D3). Here, we see that there exists a similar change between levels 2 and 3, since the signal plus noise line approaches the signal line much quicker than in the previous cases.
Next, we more directly quantified the features in the sparsity plots by plotting the changes in sparsity. These give us a more reliable way to precisely determine where to select the decomposition level cut-off. In Figure 4(B4), we see that the signal plus noise line reaches 0.05 between decomposition levels 3 and 4. This is important since it reveals a drastic change of five percent in the sparsity of the noisy signal, suggesting that we should choose our cut-off here. Likewise, in Figure 4(C4), we observe that the signal plus noise line reaches 0.05 between decomposition levels 3 and 4, informing our choose to select our cut-off between these two levels. Lastly, in the case of minimal noise addition, we see that the signal plus noise line reaches 0.05 between levels 2 and 3, suggesting that we should cut-off earlier.

B. EXPERIMENTAL DATA
In Figure 5, we apply our method to two different kinds of experimental data. In the Example 1, we generate three VOLUME 9, 2021 signals by with varying noise levels by averaging 1 scan, 4 scans and 500 scans, respectively. The sparsity plot for each of these noisy signals is shown in Figure 5(A2). It can be seen that the data undergoes the largest change in sparsity between decomposition levels 5 an 6 (see Figure 5(A3-B3)), revealing the cut-off for the decomposition level selection. We further confirm this is the case by plotting the change in sparsity for each of these data (see Figure 5(A3)). Here, we observe that there is quite a drastic change in sparsity between these two levels.
In the Example 2 (cf. Figure 5), we generate three signals by averaging 1 scan, 4 scans and 16 scans. The sparsity for each of these data in Figure 5(B2) and it is obvious that the data undergoes the largest change in sparsity between decomposition levels 4 and 5, and hence determining the decomposition level cut-off for our selection. This can be confirmed by plotting the change in sparsity for each of these data (see Figure 5(B3)), where there is a substantial change in sparsity, over ten percent, between these two levels.

VI. CONCLUSION
This paper presents a method for selecting decomposition levels for wavelet denoising. This new procedure facilitates the accurate determination of decomposition levels to ensure effective denoising and improved retrieval of the desired signal. It is demonstrated through rigorous mathematical justification that an objective sparsity-based criterion is an appropriate metric to utilize for decomposition level selection.
The presented method is advantageous to the standard practice of either arbitrary selection or trial and error to determine decomposition levels. While these latter strategies are prone to inadequate noise reduction and/or signal distortion, an objective sparsity-based criterion is shown to readily distinguish between noisy and noise-free Detail components. Signal in wavelet domain is sparse due to its coherent nature, whereas noise, being random and Gaussian in nature, is present throughout the noisy Detail component. This enables the sparsity value to separate the noisy and noise-free Detail components. The presented novel method can also be easily inserted into any wavelet denoising procedure and can be applied to a wide range of signals.
Our future work involves extending the decomposition level selection criterion for multidimensional signals.

CODE AVAILABILITY
We maintain a GitHub repository (https://github.com/ Signal-Science-Lab / Sparsity-Based-Decomposition-Level-Selection) containing MATLAB scripts to replicate our simulations and other analyses. The codes can also be found at https://denoising.cornell.edu/ and the Signal Science Lab website (https://signalsciencelab.com/). The code is also available by reaching out to the authors.