Properties of the SURE Estimates When Using Continuous Thresholding Functions for Wavelet Shrinkage

Abstract

Wavelet analysis algorithms in combination with thresholding procedures are widely used in nonparametric regression problems when estimating a signal function from noisy data. The advantages of these methods lie in their computational efficiency and the ability to adapt to the local features of the estimated function. It is usually assumed that the signal function belongs to some special class. For example, it can be piecewise continuous or piecewise differentiable and have a compact support. These assumptions, as a rule, allow the signal function to be economically represented on some specially selected basis in such a way that the useful signal is concentrated in a relatively small number of large absolute value expansion coefficients. Then, thresholding is performed to remove the noise coefficients. Typically, the noise distribution is assumed to be additive and Gaussian. This model is well studied in the literature, and various types of thresholding and parameter selection strategies adapted for specific applications have been proposed. The risk analysis of thresholding methods is an important practical task, since it makes it possible to assess the quality of both the methods themselves and the equipment used for processing. Most of the studies in this area investigate the asymptotic order of the theoretical risk. In practical situations, the theoretical risk cannot be calculated because it depends explicitly on the unobserved, noise-free signal. However, a statistical risk estimate constructed on the basis of the observed data can also be used to assess the quality of noise reduction methods. In this paper, a model of a signal contaminated with additive Gaussian noise is considered, and the general formulation of the thresholding problem with threshold functions belonging to a special class is discussed. Lower bounds are obtained for the threshold values that minimize the unbiased risk estimate. Conditions are also given under which this risk estimate is asymptotically normal and strongly consistent. The results of these studies can provide the basis for further research in the field of constructing confidence intervals and obtaining estimates of the convergence rate, which, in turn, will make it possible to obtain specific values of errors in signal processing for a wide range of thresholding methods.

1. Introduction

Wavelet methods are widely used for various signal and image processing tasks. They are applied for the frequency–time analysis of signal spectra, the identification of singular points and data compression. Moreover, during the process of receiving, processing and transmitting, signals are often contaminated with noise due to system errors, environmental interference and the internal noise of receivers. To solve this problem, various noise reduction methods have been proposed, such as linear and median filtering. Recently, wavelet decomposition has become a popular tool in solving noise reduction problems. Compared with traditional linear and median filters, wavelet analysis methods offer advantages such as computational speed and adaptive spectral analysis capabilities: a flexible time–frequency window is automatically narrowed when considering high-frequency phenomena and expanded when studying low-frequency oscillations. This makes wavelets an effective tool for the analysis of non-stationary signals with varying regularity in their domain. Telecommunication channels and computer networks may contain additive Gaussian noise caused by thermal sources, solar radiation and other natural causes. Sensor noise, electronic component noise and poor lighting are also sources of Gaussian noise in digital images. The use of standard filtering methods may lead to the excessive blurring of small-scale image details. The use of wavelets in this situation also offers certain advantages, since applying the wavelet transform retains the Gaussian noise structure, and the useful signal is economically represented as a small number of wavelet coefficients. Thus, in the space of wavelet coefficients, it is possible to effectively separate the noise from the useful signal.

Many methods based on wavelet theory have been proposed [1]. Among these methods, threshold processing is the most widely used due to its simplicity and efficiency. Donoho and Johnstone [2] suggested that only a small number of wavelet coefficients are needed to restore the signal, and they proposed the classical method of hard threshold processing. When using the hard threshold processing function, coefficients with absolute values lower than the threshold are removed, while large coefficients remain unchanged. This function is discontinuous, and additional artifacts, such as an analog of the Gibbs effect, may appear in the restored signal. To make the threshold function continuous, Donoho proposed a soft threshold processing function that reduces large coefficients by the threshold value. However, this introduces an additional bias that can lead to distortions in the reconstructed signal (for more information, see Ref. [3]). Many authors have proposed various versions of thresholding functions that improve the classical methods to some extent [4,5,6,7,8,9,10,11,12,13,14,15]. Such functions, like the soft thresholding function, are continuous and allow one to construct an unbiased estimate of the mean square risk using Stein’s method [16]. In this case, the threshold value is usually proposed to be calculated by solving the problem of minimizing this estimate. Due to the appearance of such a variety of threshold processing methods, Johnstone [3] proposed to consider a fairly general class of thresholding functions that satisfy some natural requirements and allow the construction of an unbiased risk estimate. The study of risk estimate properties is also an important practical task, since this estimate allows one to quantitatively evaluate the error of the noise suppression method using only the observed data. However, while the properties of theoretical risk are well studied, the properties of its estimates have received less attention so far. The purpose of this paper is to fill this gap and study the statistical properties of thresholding risk estimates. This paper proves that, under some additional requirements for the threshold function, the unbiased estimate of the mean square risk is strongly consistent and asymptotically normal. These properties serve as the basis for the use of the risk estimate as a quality criterion for the thresholding method.

2. Wavelet Decomposition

For a function $f\in L^2\left(\mathbf{R}\right)$, the wavelet decomposition is a series of the form

$$f=\sum _{j,k\in \mathbf{Z}}\langle f,{\psi }_{j,k}\rangle {\psi }_{j,k},$$

where ${\psi }_{j,k}\left(x\right)=2^{j/2}\psi (2^{j}x-k)$ are shifts and dilations of the chosen wavelet function $\psi$. The index j in this expansion is called the scale, and the index k is called the shift. Under certain conditions, the sequence $\left\{{\psi }_{j,k}\right\}$ forms an orthonormal basis. The function $\psi$ can be chosen so that it is sufficiently smooth and has other useful properties [1].

In what follows, we will consider functions f with support on some finite interval $[a,b]$, uniformly Lipschitz regular with the exponent $\gamma >0$. If a wavelet function is M times continuously differentiable ($M⩾\gamma$), has M vanishing moments, i.e.,

$$\underset{-\infty }{\overset{\infty }{\int }}{x}^k\psi \left(x\right)dx=0,k=0,\dots ,M-1,$$

and rapidly decays at infinity together with its derivatives, i.e., for all $0\le k\le M$ and any $m\in N$, there is a constant $C_m$ such that, for all $x\in \mathbf{R}$,

$$|{\psi }^{\left(k\right)}\left(x\right)|\le {\displaystyle \frac{C_m}{1+{\left|x\right|}^m}},$$

then

$$|\langle f,{\psi }_{j,k}\rangle |\le {\displaystyle \frac{C_f}{2^{j\left(\gamma +1/2\right)}}},$$

where $C_f$ is some positive constant depending on f (a wavelet function with the given properties always exists—for example, a Daubechies wavelet of the corresponding order). This inequality shows that, in the space of wavelet coefficients, Lipschitz regular functions have a sparse representation.

In practice, functions are defined by their values on some grid of points. Suppose that this grid is uniform and the number of points is $N=2^J$ for some $J>0$. In the case where the number of observations differs from a power of two, some extension of the array to the nearest power of 2 is applied [17]. The discrete wavelet transform is the multiplication of the vector of function values by the orthogonal matrix W defined by the wavelet function $\psi$. This results in a set of discrete wavelet coefficients ${\mu }_{j,k}$, for which, for a sufficiently large N, the relation ${\mu }_{j,k}\approx \sqrt{N}\langle f,{\psi }_{j,k}\rangle$ is valid [1]. Thus, for ${\mu }_{j,k}$, the following inequality holds:

$$|{\mu }_{j,k}|\le {\displaystyle \frac{C_f{2}^{J/2}}{2^{j\left(\gamma +1/2\right)}}}.$$

(1)

3. Thresholding of Wavelet Coefficients

Real observations are usually contaminated by noise. In many cases, it is assumed that the noise is additive, white and Gaussian. Thus, we consider the following model:

$$X_i=f_i+z_i,i=1,\dots ,N,$$

where $f_i$ are values of the function f, and the random variables $z_i$ have the distribution $N(0,{\sigma }^2)$ and are independent. Since the matrix W is orthogonal, the discrete wavelet coefficients have the form

$$Y_{j,k}={\mu }_{j,k}+{ϵ}_{j,k},j=0,\dots ,J-1,k=0,\dots ,2^j-1,$$

where ${ϵ}_{j,k}$ have distribution $N(0,{\sigma }^2)$ and are independent. Classical methods of signal processing assume the additivity and Gaussianity of the noise. However, there are approaches (see, for example, [18,19]) that assume non-Gaussianity. The consideration of such a problem could be the subject of a separate study.

To suppress noise, a thresholding function is applied to the wavelet coefficients, the purpose of which is to remove sufficiently small coefficients that are considered noise. The most common are the hard thresholding function

$${\rho }_H(y,T)=\left\{\begin{array}{cc}y\hfill & \mathrm{for} \left|y\right|>T;\hfill \\ 0\hfill & \mathrm{for} \left|y\right|\le T\hfill \end{array}\right.$$

and the soft thresholding function

$${\rho }_S(y,T)=\left\{\begin{array}{cc}y-T\hfill & \mathrm{for} y>T;\hfill \\ y+T\hfill & \mathrm{for} y<-T;\hfill \\ 0\hfill & \mathrm{for} \left|y\right|\le T\hfill \end{array}\right.$$

with some threshold T. However, each of them has its own drawbacks. The function ${\rho }_H$ is discontinuous, which leads to a lack of stability and the appearance of additional artifacts, and the function ${\rho }_S$ leads to the appearance of an additional bias in the estimate of the function f. A number of studies have proposed some alternative thresholding functions $\rho (y,T)$, which, in fact, are a compromise between hard and soft thresholding [4,5,6,7,8,9,10,11,12,13,14,15]. These functions are continuous, like the function ${\rho }_S(y,T)$, but $\rho (y,T)\to y$ as $\left|y\right|\to \infty$, i.e., for large absolute values of the coefficients, they resemble the hard thresholding function.

Due to the emergence of various types of thresholding functions, it makes sense to consider some general classes of such functions. Let the function $h(y,T)$ have the following properties:

1. $h(-y,T)=-h(y,T)$ (oddness);

2. $0\le h(y,T)\le T$ for $y⩾0$ (boundedness);

3. $h(T,T)=T$ (continuity).

Then, the thresholding function is defined by the expression

$${\rho }_h(y,T)=\left\{\begin{array}{cc}y-h(y,T)\hfill & \mathrm{for} \left|y\right|>T;\hfill \\ 0\hfill & \mathrm{for} \left|y\right|\le T.\hfill \end{array}\right.$$

This form has, for example, the soft thresholding function, the hybrid thresholding function [5,15], the thresholding function based on the hyperbolic tangent [10], the sigmoid function [12] and some others.

Thus, the signal processing scheme is as follows. First, a discrete wavelet transform is performed using the orthogonal matrix W. Then, the threshold processing of the resulting wavelet coefficients is performed, which results in the suppression of most of the noise, since the useful signal has an economical representation in the space of wavelet coefficients, and the noise is uniformly distributed over all coefficients. Finally, the inverse wavelet transform is performed using the matrix $W^T$ to obtain the processed signal ($W^{-1}=W^T$, since W is orthogonal). Moreover, the task of the economical presentation of data is solved along the way for storage or transmission over communication channels.

The error (or risk) of threshold processing is defined as follows:

$$R_J\left(T\right)=\sum _{j=0}^{J-1}\sum _{k=0}^{2^j-1}\mathsf{E}{\left({\rho }_h(Y_{j,k},T)-{\mu }_{j,k}\right)}^2.$$

(2)

One of the key issues is the choice of the threshold value T (note that a “reasonable” threshold value should increase with J [1]). One of the first to appear was the so-called universal threshold $T_U=\sigma \sqrt{2ln{2}^J}$. The following statement demonstrates that this threshold is, in a certain sense, the maximum among the “reasonable” ones [1].

Lemma 1.

Let $Z_1,\dots ,Z_N$ be independent random variables with distribution $N(0,{\sigma }^2)$. Then,

$$\mathsf{P}\left(\underset{1\le i\le N}{max}\left|Z_i\right|>\sigma \sqrt{2lnN}\right)\to 0$$

as $N\to \infty$.

Thus, the maximum noise amplitude with a high probability does not exceed the universal threshold. In Refs. [20,21], additional considerations are also given, showing that it is reasonable not to consider $T>T_U$.

4. Statistical Estimate of the Mean Square Risk

Another approach to calculating the threshold value is based on minimizing the statistical risk estimate (2) constructed using Stein’s method [16]. If $h(y,T)$ is differentiable with respect to the variable y, then the following relation is valid:

$$\mathsf{E}{\left({\rho }_h(Y_{j,k},T)-{\mu }_{j,k}\right)}^2={\sigma }^2+\mathsf{E}{g}^2\left(Y_{j,k}\right)-2{\sigma }^2\mathsf{E}{g}^{\prime }\left(Y_{j,k}\right),$$

where $g\left(Y_{j,k}\right)=Y_{j,k}-{\rho }_h(Y_{j,k},T)$. That is, Stein’s method consists of the fact that, in (2), the expression containing unknown ${\mu }_{j,k}$ can be replaced by an expression containing only observed data.

Thus, the following expression can be used as an unbiased risk estimate:

$${\widehat{R}}_J\left(T\right)=\sum _{j=0}^{J-1}\sum _{k=0}^{2^j-1}F(Y_{j,k},T),$$

(3)

where

$$F(y,T)=\left\{\begin{array}{cc}{y}^2-{\sigma }^2\hfill & \mathrm{for} \left|y\right|\le T;\hfill \\ h^2(y,T)+{\sigma }^2-2{\sigma }^2{h}_y^{\prime }(y,T)\hfill & \mathrm{for} \left|y\right|>T.\hfill \end{array}\right.$$

(4)

Let $T_S$ denote the threshold that minimizes the estimate (3):

$${\widehat{R}}_J\left(T_S\right)=\underset{T\in [0,T_U]}{min}{\widehat{R}}_J\left(T\right).$$

This threshold mimics the theoretical “ideal” threshold $T_{min}$ that minimizes the theoretical risk (2):

$$R_J\left(T_{min}\right)=\underset{T\in [0,T_U]}{min}{R}_J\left(T\right).$$

Remark 1.

Since this paper studies the order of the expressions under consideration with respect to J, we will use the notations C, c and others for some positive constants that may depend on the model parameters but do not depend on J. In this case, in one chain of transformations, different constants may be designated by the same letters, and the transformations themselves are performed under the assumption that J is sufficiently large.

The following upper bound is valid for $R_J\left(T_{min}\right)$.

Lemma 2.

Let f be supported on some finite interval $[a,b]$ and be uniformly Lipschitz regular with some exponent $\gamma >0$. Then,

$$R_J\left(T_{min}\right)\le C{2}^{J/(2\gamma +1)}{J}^{(2\gamma +2)/(2\gamma +1)}.$$

(5)

Proof. 

Let the index $j^{*}$ be such that

$$|{\mu }_{j,k}|\le C{J}^{-1/2},j^{*}\le j\le J-1.$$

In this case, due to (1),

$$j^{*}\approx {\displaystyle \frac{J}{2\gamma +1}}+{\displaystyle \frac{1}{2\gamma +1}}{log}_{2}J+c.$$

Let us split the expression for risk into two sums:

$$R_J\left(T\right)=\sum _{j=0}^{J-1}\sum _{k=0}^{2^j-1}\mathsf{E}{({\rho }_h(Y_{j,k},T)-{\mu }_{j,k})}^2=$$

$$=\sum _{j=0}^{j^{*}-1}\sum _{k=0}^{2^j-1}\mathsf{E}{({\rho }_h(Y_{j,k},T)-{\mu }_{j,k})}^2+\sum _{j=j^{*}}^{J-1}\sum _{k=0}^{2^j-1}\mathsf{E}{({\rho }_h(Y_{j,k},T)-{\mu }_{j,k})}^2\equiv S_1+S_2.$$

Let us estimate the first sum from above. Given the restrictions on $h(y,T)$, we obtain

$$\mathsf{E}{({\rho }_h(Y_{j,k},T)-{\mu }_{j,k})}^2=\mathsf{E}((Y_{j,k}-h(Y_{j,k},T))\mathbb{I}\left(\right|Y_{j,k}|>T)-{\mu }_{j,k}{)}^2=$$

$$=\mathsf{E}(Y_{j,k}-h(Y_{j,k},T)-(Y_{j,k}-h(Y_{j,k},T))\mathbb{I}\left(\right|Y_{j,k}|\le T)-{\mu }_{j,k}{)}^2\le$$

$$\le 2\mathsf{E}{(Y_{j,k}-{\mu }_{j,k})}^2+2\mathsf{E}[h(Y_{j,k},T)+(Y_{j,k}-h(Y_{j,k},T))\mathbb{I}\left(\right|Y_{j,k}{|\le T)]}^2\le 2{\sigma }^2+C{T}^2.$$

Hence,

$$S_1=\sum _{j=0}^{j^{*}-1}\sum _{k=0}^{2^j-1}\mathsf{E}{({\rho }_h(Y_{j,k},T)-{\mu }_{j,k})}^2\le C{2}^{j^{*}}{T}^2\le C{T}^2{2}^{J/(2\gamma +1)}{J}^{1/(2\gamma +1)}.$$

Now, consider the second sum. Given the constraints on $h(y,T)$ and the fact that, for each term of the second sum $|{\mu }_{j,k}|\le C{J}^{-1/2}$ (and hence $|{\mu }_{j,k}|\le C/T$, since $T<T_U$), we obtain

$$\mathsf{E}{({\rho }_h(Y_{j,k},T)-{\mu }_{j,k})}^2=\mathsf{E}((Y_{j,k}-h(Y_{j,k},T))\mathbb{I}\left(\right|Y_{j,k}|>T)-{\mu }_{j,k}{)}^2\le$$

$$\le 2\mathsf{E}{(Y_{j,k}-h(Y_{j,k},T))}^2\mathbb{I}\left(\right|Y_{j,k}|>T)+2{\mu }_{j,k}^2\le$$

$$\le C\left(Texp\left\{-{\displaystyle \frac{{(T+{\mu }_{j,k})}^2}{2{\sigma }^2}}\right\}+Texp\left\{-{\displaystyle \frac{{(T-{\mu }_{j,k})}^2}{2{\sigma }^2}}\right\}\right)+2{\mu }_{j,k}^2\le CTexp\left\{-{\displaystyle \frac{T^2}{2{\sigma }^2}}\right\}+2{\mu }_{j,k}^2.$$

Taking into account (1), we have

$$S_2=\sum _{j=j^{*}}^{J-1}\sum _{k=0}^{2^j-1}\mathsf{E}{({\rho }_h(Y_{j,k},T)-{\mu }_{j,k})}^2\le C\left(Texp\left\{-{\displaystyle \frac{T^2}{2{\sigma }^2}}\right\}{2}^J+2^{J/(2\gamma +1)}{J}^{-{\displaystyle \frac{2\gamma }{2\gamma +1}}}\right).$$

Thus,

$$R_J\left(T\right)=S_1+S_2\le C\left(T^2{2}^{J/(2\gamma +1)}{J}^{1/(2\gamma +1)}+Texp\left\{-{\displaystyle \frac{T^2}{2{\sigma }^2}}\right\}{2}^J\right).$$

In this expression, the first term increases with the growing T, and the second one decreases (while, as noted, $T<T_U$). Repeating the reasoning of Ref. [22] and choosing the value

$$T^{*}=\sigma \sqrt{{\displaystyle \frac{4\gamma }{2\gamma +1}}ln{2}^J}-\sigma {\displaystyle \frac{2\gamma +3}{4\gamma +2}}\sqrt{{\displaystyle \frac{2\gamma +1}{4\gamma }}}{\displaystyle \frac{lnln{2}^J}{\sqrt{ln{2}^J}}}$$

to equalize the orders of these terms, we obtain the upper bound (5) (in this case, the value $T^{*}$ represents the lower bound for the “ideal” threshold $T_{min}$). The lemma is proven. □

For the threshold value $T_S$, the following statement holds, showing that it cannot be too small.

Lemma 3.

Let f be supported on some finite interval $[a,b]$ and be uniformly Lipschitz regular with some exponent $\gamma >0$. Let, for some positive θ,

$$|h_y^{\prime }(y,T)|\le C{T}^{\theta },|h_T^{\prime }(y,T)|\le C{T}^{\theta },\left|h_{yT}^{″}(y,T)\right|\le C{T}^{\theta },\forall y\in \mathbb{R}.$$

(6)

Then, for $0<\gamma \le 1/2$ and arbitrary $0<\alpha <{\displaystyle \frac{4\gamma }{2\gamma +1}}$,

$$\mathsf{P}(T_S\in [0,T_1])\le C{2}^{J}exp\left\{-c{J}^{-max\{3,\theta +1\}}{2}^{(1-\alpha )J}\right\}.$$

(7)

If $\gamma >1/2$, then, for an arbitrary ${\displaystyle \frac{2}{3}}<\alpha <{\displaystyle \frac{2}{3}}·{\displaystyle \frac{4\gamma }{2\gamma +1}}$,

$$\mathsf{P}(T_S\in [0,T_2])\le C{2}^{J}exp\left\{-c{J}^{-max\{5/2,\theta +1/2\}}{2}^{(1-\alpha )J}\right\},$$

(8)

where $T_1=\sigma \sqrt{\alpha ln{2}^J}$, $T_2=\sigma \sqrt{{\displaystyle \frac{3\alpha }{2}}ln{2}^J}$.

Proof. 

Let $0<\gamma \le 1/2$. Then, for $0<\alpha <{\displaystyle \frac{4\gamma }{2\gamma +1}}$, we obtain

$$\mathsf{P}\left(T_S\in [0,T_1]\right)\le \mathsf{P}\left(\underset{T\in [0,T_1]}{sup}({\widehat{R}}_J\left(T_{min}\right)-{\widehat{R}}_J\left(T\right))\ge 0\right)\le$$

$$\le \mathsf{P}\left(\underset{T\in [0,T_1]}{sup}\left|{\widehat{R}}_J\left(T_{min}\right)-{\widehat{R}}_J\left(T\right)-\mathsf{E}({\widehat{R}}_J\left(T_{min}\right)-{\widehat{R}}_J\left(T\right))\right|\ge \underset{T\in [0,T_1]}{inf}(R_J\left(T\right)-R_J\left(T_{min}\right))\right).$$

Next, since $0\le h(y,T)\le T$,

$$\mathsf{E}{({\rho }_h(Y_{J-1,k},T)-{\mu }_{J-1,k})}^2=\mathsf{E}((Y_{J-1,k}-h(Y_{J-1,k},T))\mathbb{I}\left(\right|Y_{J-1,k}|>T)-{\mu }_{J-1,k}{)}^2\ge$$

$$\ge \mathsf{E}\left(\right|Y_{J-1,k}{|-T)}^2\mathbb{I}\left(\right|Y_{J-1,k}|>T)-2|{\mu }_{J-1,k}\left|\mathsf{E}\right|Y_{J-1,k}\left|\mathbb{I}\right(|Y_{J-1,k}|>T)+{\mu }_{J-1,k}^2.$$

Using the inequality

$$1-\mathsf{\Phi }\left(y\right)>\phi \left(y\right)\left({\displaystyle \frac{1}{y}}-{\displaystyle \frac{1}{y^3}}\right),y>0,$$

where $\phi \left(y\right)$ and $\mathsf{\Phi }\left(y\right)$ are the density and distribution function of the standard normal law, and given that, by (1),

$$|{\mu }_{J-1,k}|\le {\displaystyle \frac{C}{2^{\gamma J}}},$$

we conclude that there is a constant $C_T>0$ such that, for $C_T<T\le T_1$,

$$\mathsf{E}\left(\right|Y_{J-1,k}{|-T)}^2\mathbb{I}\left(\right|Y_{J-1,k}|>T)-2|{\mu }_{J-1,k}\left|\mathsf{E}\right|Y_{J-1,k}\left|\mathbb{I}\right(|Y_{J-1,k}|>T)+{\mu }_{J-1,k}^2\ge$$

$$\ge {\displaystyle \frac{c}{T}}{e}^{-{\displaystyle \frac{T^2}{2{\sigma }^2}}}+{\mu }_{J-1,k}^2.$$

If $0<T\le C_T$, then, obviously, for some $\delta >0$, independent of J,

$$\mathsf{E}{({\rho }_h(Y_{J-1,k},T)-{\mu }_{J-1,k})}^2\ge \delta .$$

Thus,

$$\underset{T\in [0,T_1]}{inf}{R}_J\left(T\right)\ge \underset{T\in [0,T_1]}{inf}\sum _{k=0}^{2^{J-1}-1}\mathsf{E}{\left({\rho }_h(Y_{J-1,k},T)-{\mu }_{J-1,k}\right)}^2\ge {\displaystyle \frac{C{2}^J}{T_1}}{e}^{-{\displaystyle \frac{T_1^2}{2{\sigma }^2}}}={\displaystyle \frac{C{2}^{\frac{2-\alpha }{2}J}}{\sqrt{ln{2}^J}}}.$$

Therefore, taking into account (5),

$$\underset{T\in [0,T_1]}{inf}(R_J\left(T\right)-R_J\left(T_{min}\right))\ge {\displaystyle \frac{C{2}^{\frac{2-\alpha }{2}J}}{\sqrt{ln{2}^J}}}.$$

Denote

$$Z\left(T\right)={\widehat{R}}_J\left(T_{min}\right)-{\widehat{R}}_J\left(T\right)-\mathsf{E}({\widehat{R}}_J\left(T_{min}\right)-{\widehat{R}}_J\left(T\right)).$$

For $T<T^{\prime }$,

$$|Z\left(T\right)-Z\left(T^{\prime }\right)|=|{\widehat{R}}_J\left(T^{\prime }\right)-\mathsf{E}{\widehat{R}}_J\left(T^{\prime }\right)-{\widehat{R}}_J\left(T\right)+\mathsf{E}{\widehat{R}}_J\left(T\right)|.$$

Let $N(T,T^{\prime })=\#\{(j,k)|T<|Y_{j,k}|\le T^{\prime }\}$. Then, given the form of $F(y,T)$ and the constraints (6), we obtain the estimate

$$|Z\left(T\right)-Z\left(T^{\prime }\right)|\le C{T}^{\prime max\{2,\theta \}}N(T,T^{\prime })+C{T}^{\prime max\{2,\theta +1\}}(T^{\prime }-T)2^J.$$

Let us divide the interval $[0,T_1]$ into equal parts using points

$$T^l=l{\delta }_J,l=0,\dots ,2^J,{\delta }_J={\displaystyle \frac{T_1}{2^J}}.$$

Then,

$$A_J=\left\{\underset{T\in [0,T_1]}{sup}\left|{\widehat{R}}_J\left(T_{min}\right)-{\widehat{R}}_J\left(T\right)-\mathsf{E}({\widehat{R}}_J\left(T_{min}\right)-{\widehat{R}}_J\left(T\right))\right|\ge {\displaystyle \frac{C{2}^{\frac{2-\alpha }{2}J}}{\sqrt{ln{2}^J}}}\right\}\subset D_J\cup E_J,$$

where

$$D_J=\left\{\underset{l}{sup}\left|Z\left(T^l\right)\right|>{\displaystyle \frac{C{2}^{\frac{2-\alpha }{2}J}}{2\sqrt{ln{2}^J}}}\right\},E_J=\left\{\underset{l}{sup}\underset{T\in [T^l,T^{l+1})}{sup}\left|Z\left(T\right)-Z\left(T^l\right)\right|\ge {\displaystyle \frac{C{2}^{\frac{2-\alpha }{2}J}}{2\sqrt{ln{2}^J}}}\right\}.$$

Using Hoeffding’s inequality [23], we obtain

$$\mathsf{P}\left(D_J\right)\le \sum _{l=1}^{2^J}\mathsf{P}\left(|Z\left(T^l\right)|>{\displaystyle \frac{C{2}^{\frac{2-\alpha }{2}J}}{2\sqrt{ln{2}^J}}}\right)\le C{2}^{J}exp\left\{-c{\displaystyle \frac{2^{(1-\alpha )J}}{ln{2}^J{T}_{min}^{max\{4,2\theta \}}}}\right\}.$$

(9)

Next,

$$E_J\subset \left\{\underset{l}{sup}\underset{T\in [T^l,T^{l+1})}{sup}(C{T}^{max\{2,\theta \}}N(T^l,T)+C{T}^{max\{2,\theta +1\}}(T-T^l)2^J)\ge {\displaystyle \frac{C{2}^{\frac{2-\alpha }{2}J}}{2\sqrt{ln{2}^J}}}\right\}\subset$$

$$\subset \left\{\underset{l}{sup}C{\left(T^{l+1}\right)}^{max\{2,\theta \}}N(T^l,T^{l+1})\ge {\displaystyle \frac{C{2}^{\frac{2-\alpha }{2}J}}{2\sqrt{ln{2}^J}}}-C{\left(T_1\right)}^{max\{3,\theta +2\}}\right\}.$$

Since

$$\mathsf{E}N(T^l,T^{l+1})=\mathsf{E}\sum _{j=0}^{J-1}\sum _{k=0}^{2^j-1}\mathbb{I}(T^l<|Y_{j,k}|\le T^{l+1})\le C{T}_1,$$

we have

$$\left\{\underset{l}{sup}C{\left(T^{l+1}\right)}^{max\{2,\theta \}}N(T^l,T^{l+1})\ge {\displaystyle \frac{C{2}^{\frac{2-\alpha }{2}J}}{2\sqrt{ln{2}^J}}}-C{\left(T_1\right)}^{max\{3,\theta +2\}}\right\}\subset$$

$$\begin{gathered} \subset \left\{\underset{l}{sup}\left|C{\left(T^{l+1}\right)}^{max\{2,\theta \}}N(T^l,T^{l+1})-\mathsf{E}C{\left(T^{l+1}\right)}^{max\{2,\theta \}}N(T^l,T^{l+1})\right|\right.\ge \\ \left.\ge {\displaystyle \frac{C{2}^{\frac{2-\alpha }{2}J}}{2\sqrt{ln{2}^J}}}-C{\left(T_1\right)}^{max\{3,\theta +2\}}\right\}. \end{gathered}$$

Applying Hoeffding’s inequality one more time, we obtain

$$\mathsf{P}\left(E_n\right)\le C{2}^{J}exp\left\{-c{\displaystyle \frac{2^{(1-\alpha )J}}{ln{2}^J{T}_1^{max\{4,2\theta \}}}}\right\}.$$

(10)

Combining (9) and (10), we obtain

$$\mathsf{P}(T_S\in [0,T_1])\le \mathsf{P}\left(A_J\right)\le C{2}^{J}exp\left\{-c{\displaystyle \frac{2^{(1-\alpha )J}}{ln{2}^J{T}_1^{max\{4,2\theta \}}}}\right\}\le C{2}^{J}exp\left\{-c{J}^{-max\{3,\theta +1\}}{2}^{(1-\alpha )J}\right\}.$$

Now, let $\gamma >1/2$. Then, for $0<\alpha <1$, the estimate of the probability $\mathsf{P}(T_S\in [0,T_1])$ has exactly the same form. Next, let ${\displaystyle \frac{2}{3}}<\alpha <{\displaystyle \frac{2}{3}}·{\displaystyle \frac{4\gamma }{2\gamma +1}}$ and $\beta ={\displaystyle \frac{3\alpha }{2}}$. Then, repeating the previous reasoning, we obtain

$$\underset{T\in [T_1,T_2]}{inf}(R_J\left(T\right)-R_J\left(T_{min}\right))\ge {\displaystyle \frac{C{2}^{\frac{2-\beta }{2}J}}{\sqrt{ln{2}^J}}},$$

and

$$\mathsf{P}(T_S\in [T_1,T_2])\le \mathsf{P}\left(A_J^{\prime }\right),$$

where

$$A_J^{\prime }=\left\{\underset{T\in [T_1,T_2]}{sup}\left|{\widehat{R}}_J\left(T_{min}\right)-{\widehat{R}}_J\left(T\right)-\mathsf{E}({\widehat{R}}_J\left(T_{min}\right)-{\widehat{R}}_J\left(T\right))\right|\ge {\displaystyle \frac{C{2}^{\frac{2-\beta }{2}J}}{\sqrt{ln{2}^J}}}\right\}\subset D_J^{\prime }\cup E_J^{\prime },$$

$$D_J^{\prime }=\left\{\underset{l}{sup}\left|Z\left(T^l\right)\right|>{\displaystyle \frac{C{2}^{\frac{2-\beta }{2}J}}{2\sqrt{ln{2}^J}}}\right\},$$

$$E_J^{\prime }=\left\{\underset{l}{sup}\underset{T\in [T^l,T^{l+1})}{sup}\left|Z\left(T\right)-Z\left(T^l\right)\right|\ge {\displaystyle \frac{C{2}^{\frac{2-\beta }{2}J}}{2\sqrt{ln{2}^J}}}\right\},$$

$$T^l=T_1+l{\delta }_J,l=0,\dots ,2^J,{\delta }_J={\displaystyle \frac{T_2-T_1}{2^J}}.$$

It can also be shown that

$$\mathsf{D}({\widehat{R}}_J\left(T_{min}\right)-{\widehat{R}}_J\left(T\right))\le C{T}^{max\{3,2\theta -1\}}{2}^{{\displaystyle \frac{2-\alpha }{2}}J}.$$

Now, repeating the previous reasoning, to estimate the probabilities of the events $D_J^{\prime }$ and $E_J^{\prime }$, instead of Hoeffding’s inequality, we apply Bernstein’s inequality [23] and obtain

$$\mathsf{P}\left(D_J^{\prime }\right)\le C{2}^{J}exp\left\{-c{J}^{-max\{5/2,\theta +1/2\}}{2}^{(2-2\beta +\alpha )J/2}\right\},$$

$$\mathsf{P}\left(E_J^{\prime }\right)\le C{2}^{J}exp\left\{-c{J}^{-max\{5/2,\theta +1/2\}}{2}^{(2-2\beta +\alpha )J/2}\right\}.$$

Hence,

$$\mathsf{P}(T_S\in [T_1,T_2])\le C{2}^{J}exp\left\{-c{J}^{-max\{5/2,\theta +1/2\}}{2}^{(2-2\beta +\alpha )J/2}\right\}.$$

Combining this estimate with the estimate of the probability $\mathsf{P}(T_S\in [0,T_1])$ and taking into account that $\beta ={\displaystyle \frac{3\alpha }{2}}$, we obtain the statement of the lemma for $\gamma >1/2$.

The lemma is proven. □

5. Asymptotic Properties of the Mean Square Risk Estimate

In this section, we present statements describing the asymptotic properties of the risk estimate. For soft thresholding, similar results were obtained in [24,25].

Theorem 1.

Let f be supported on some finite interval $[a,b]$ and be uniformly Lipschitz regular with some exponent $\gamma >1/2$. Let $h(y,T)$ satisfy the conditions (6). Then,

$$\mathsf{P}\left({\displaystyle \frac{{\widehat{R}}_J\left(T_S\right)-R_J\left(T_{min}\right)}{{\sigma }^2\sqrt{2^{J+1}}}}<x\right)\to \mathsf{\Phi }\left(x\right)while J\to \infty .$$

(11)

Proof. 

Denote

$$U\left(T\right)={\widehat{R}}_J\left(T\right)-{\widehat{R}}_J\left(T_{min}\right)=\sum _{j=0}^{J-1}\sum _{k=0}^{2^j-1}(F(Y_{j,k},T)-F(Y_{j,k},T_{min})).$$

Let us write the numerator in the expression under the probability in (11) as

$${\widehat{R}}_J\left(T_S\right)-R_J\left(T_{min}\right)=U\left(T_S\right)+{\widehat{R}}_J\left(T_{min}\right)-R_J\left(T_{min}\right)$$

and show that

$${\displaystyle \frac{{\widehat{R}}_J\left(T_{min}\right)-R_J\left(T_{min}\right)}{{\sigma }^2\sqrt{2^{J+1}}}}⟹N(0,1),J\to \infty .$$

(12)

Let $V_J^2=\mathsf{D}({\widehat{R}}_J\left(T_{min}\right)-R_J\left(T_{min}\right))$. Given (1) and the fact that $\mathsf{D}{Y}_{j,k}=2{\sigma }^4+4{\sigma }^2{\mu }_{j,k}^2$, it can be shown that

$$\underset{J\to \infty }{lim}{\displaystyle \frac{V_J^2}{{\sigma }^4{2}^{J+1}}}=1.$$

Moreover, the Lindeberg condition is satisfied: for any $ϵ>0$, when $J\to \infty$,

$${\displaystyle \frac{1}{V_J^2}}\sum _{j=0}^{J-1}\sum _{k=0}^{2^j-1}\mathsf{E}\left[{(F(Y_{j,k},T_{min})-\mathsf{E}F(Y_{j,k},T_{min}))}^2\mathbb{I}\left(\right|F(Y_{j,k},T_{min})-\mathsf{E}F(Y_{j,k},T_{min})|>ϵV_J)\right]\to 0,$$

since, starting from some J, the indicators become equal to zero. Therefore, (12) holds. It remains to show that

$${\displaystyle \frac{U\left(T_S\right)}{\sqrt{2^J}}}\stackrel{\mathsf{P}}{⟶}0.$$

(13)

For ${\displaystyle \frac{2}{3}}<\alpha <{\displaystyle \frac{2}{3}}·{\displaystyle \frac{4\gamma }{2\gamma +1}}$, we denote $T_2=\sigma \sqrt{{\displaystyle \frac{3}{2}}\alpha ln{2}^J}$. Then, for an arbitrary $ϵ>0$,

$$\mathsf{P}\left({\displaystyle \frac{|U\left(T_S\right)|}{\sqrt{2^J}}}>ϵ\right)\le$$

$$\le \mathsf{P}\left(T_S\le T_2\right)+\mathsf{P}\left({\displaystyle \frac{\underset{T\in [T_2,T_U]}{sup}\left|U\left(T\right)-\mathsf{E}U\left(T\right)\right|+\underset{T\in [T_2,T_U]}{sup}\left|\mathsf{E}U\left(T\right)\right|}{\sqrt{2^J}}}>ϵ\right).$$

By Lemma 3, the first probability tends to zero. Moreover, repeating the reasoning from Lemma 3, we can verify that

$$\mathsf{P}\left({\displaystyle \frac{\underset{T\in [T_2,T_U]}{sup}\left|U\left(T\right)-\mathsf{E}U\left(T\right)\right|}{\sqrt{2^J}}}>ϵ\right)\le C{2}^{J}exp\left\{-c{J}^{-max\{1,\theta /2\}}{2}^{J/2}\right\}.$$

Finally, given that $\gamma >1/2$, it can be shown that

$${\displaystyle \frac{\underset{T\in [T_2,T_U]}{sup}\left|\mathsf{E}U\left(T\right)\right|}{\sqrt{2^J}}}\le \underset{T\in [T_2,T_U]}{sup}\left(C{T}^2{2}^{J/(2\gamma +1)}{J}^{1/(2\gamma +1)}+C{T}^{max\{1,\theta -1\}}{e}^{-{\displaystyle \frac{T^2}{2{\sigma }^2}}}{2}^J\right)\le$$

$$\le C{J}^{max\{1/2,(\theta -1)/2\}}{2}^{{\displaystyle \frac{2-3\alpha }{4}}J}\to 0.$$

Thus, (13) holds and the theorem is completely proven. □

Theorem 2.

Let f be supported on some finite interval $[a,b]$ and be uniformly Lipschitz regular with some exponent $\gamma >0$. Let $h(y,T)$ satisfy the conditions (6). Then,

$${\displaystyle \frac{{\widehat{R}}_J\left(T_S\right)-R_J\left(T_{min}\right)}{2^J}}\to 0a.s.when J\to \infty .$$

Proof. 

As in the previous theorem, we write

$${\widehat{R}}_J\left(T_S\right)-R_J\left(T_{min}\right)=U\left(T_S\right)+{\widehat{R}}_J\left(T_{min}\right)-R_J\left(T_{min}\right),$$

where

$$U\left(T\right)={\widehat{R}}_J\left(T\right)-{\widehat{R}}_J\left(T_{min}\right).$$

If $0<\gamma \le 1/2$, then, for $0<\alpha <{\displaystyle \frac{4\gamma }{2\gamma +1}}$, we denote $T_1=\sigma \sqrt{\alpha ln{2}^J}$. If $\gamma >1/2$, then, for ${\displaystyle \frac{2}{3}}<\alpha <{\displaystyle \frac{2}{3}}·{\displaystyle \frac{4\gamma }{2\gamma +1}}$, we denote $T_2=\sigma \sqrt{{\displaystyle \frac{3}{2}}\alpha ln{2}^J}$. In both situations, one of the estimates (7) or (8) is valid, the right-hand side of which we denote, in any case, by $p_J$.

Consider the case $0<\gamma \le 1/2$ (the case $\gamma >1/2$ is considered similarly). We have

$$\mathsf{P}\left({\displaystyle \frac{|U\left(T_S\right)|}{2^J}}>ϵ\right)\le$$

$$\le \mathsf{P}\left(T_S\le T_1\right)+\mathsf{P}\left({\displaystyle \frac{\underset{T\in [T_1,T_U]}{sup}\left|U\left(T\right)-\mathsf{E}U\left(T\right)\right|+\underset{T\in [T_1,T_U]}{sup}\left|\mathsf{E}U\left(T\right)\right|}{2^J}}>ϵ\right).$$

For the second probability, repeating the reasoning from Lemma 3, we obtain

$$q_J\equiv \mathsf{P}\left({\displaystyle \frac{\underset{T\in [T_1,T_U]}{sup}\left|U\left(T\right)-\mathsf{E}U\left(T\right)\right|}{2^J}}>ϵ\right)\le C{2}^{J}exp\left\{-c{J}^{-max\{1,\theta /2\}}{2}^J\right\},$$

$${\displaystyle \frac{\underset{T\in [T_1,T_U]}{sup}\left|\mathsf{E}U\left(T\right)\right|}{2^J}}\to 0.$$

Finally, using Hoeffding’s inequality, we obtain

$$r_J\equiv \mathsf{P}\left({\displaystyle \frac{{\widehat{R}}_J\left(T_{min}\right)-R_J\left(T_{min}\right)}{2^J}}>ϵ\right)\le Cexp\left\{-c{J}^{-max\{2,\theta \}}{2}^J\right\}.$$

Thus,

$$\sum _{J=1}^{\infty }(p_J+q_J+r_J)<\infty$$

and the statement of the theorem follows from the Borel–Cantelli lemma. □

6. Examples

Here, we provide some examples of thresholding functions that satisfy the conditions of Theorems 1 and 2 and have been proposed by various authors in recent years. Some of them depend on additional parameters.

• Thresholding function based on the hyperbolic tangent function [10]: for $\lambda >0$,

  $${\rho }_h(y,T)=\left\{\begin{array}{cc}\mathrm{sgn}\left(y\right)\left(\right|y|-T)+T\left({\displaystyle \frac{e^{\lambda \mathrm{sgn}\left(y\right)\left(\right|y|-T)/T}-e^{-\lambda \mathrm{sgn}\left(y\right)\left(\right|y|-T)/T}}{e^{\lambda \mathrm{sgn}\left(y\right)\left(\right|y|-T)/T}+e^{-\lambda \mathrm{sgn}\left(y\right)\left(\right|y|-T)/T}}}\right),\hfill & \left|y\right|>T;\hfill \\ 0,\hfill & \left|y\right|\le T.\hfill \end{array}\right.$$
• Hybrid thresholding function [11]: for $\lambda >1$,

  $${\rho }_h(y,T)=\left\{\begin{array}{cc}\mathrm{sgn}\left(y\right)\left(\right|y|-|y{|}^{1-\lambda }{T}^{\lambda }),\hfill & \left|y\right|>T;\hfill \\ 0,\hfill & \left|y\right|\le T.\hfill \end{array}\right.$$
• Sigmoid thresholding function [12]: for $\lambda \ge 0$,

  $${\rho }_h(y,T)=\left\{\begin{array}{cc}\mathrm{sgn}\left(y\right)\left(\right|y|-T)+{\displaystyle \frac{2T}{1+e^{-\lambda \mathrm{sgn}\left(y\right)\left(\right|y|-T)}}}-T,\hfill & \left|y\right|>T;\hfill \\ 0,\hfill & \left|y\right|\le T.\hfill \end{array}\right.$$
• Soft–hard thresholding function [13] (this function depends on two threshold values; threshold value $T^{\prime }$ is an additional parameter):

  $${\rho }_h(y,T)=\left\{\begin{array}{cc}0,\hfill & \left|y\right|\le T;\hfill \\ \mathrm{sgn}\left(y\right)T^{\prime }\left(\right|y|-T)/(T^{\prime }-T),\hfill & T<\left|y\right|\le T^{\prime };\hfill \\ y,\hfill & \left|y\right|>T^{\prime }.\hfill \end{array}\right.$$
• The smooth clipped absolute deviation (SCAD) penalty thresholding function [14]: for $\lambda >2$,

  $${\rho }_h(y,T)=\left\{\begin{array}{cc}0,\hfill & \left|y\right|\le T;\hfill \\ \mathrm{sgn}\left(y\right)\left(\right|y|-T),\hfill & T<\left|y\right|\le 2T;\hfill \\ \mathrm{sgn}\left(y\right)\left(\lambda \right|y|-|y|-\lambda T)/(\lambda -2),\hfill & 2T<\left|y\right|\le \lambda T;\hfill \\ y,\hfill & \left|y\right|>\lambda T.\hfill \end{array}\right.$$

Some other examples of thresholding functions can be found in the monograph [3] and the papers cited therein.

7. Conclusions

Wavelet analysis and threshold processing methods have found wide application in the signal analysis of noisy data. The advantages of these methods include high computational efficiency and the ability to adapt to the local features of the estimated function. Under some additional assumptions, it is possible to economically represent the signal function in such a way that the useful signal is concentrated in a relatively small number of coefficients with sufficiently large absolute values, after which threshold processing is performed to remove noise arising both for natural reasons (for example, due to background sources) and due to imperfections in the equipment used. Typically, the noise distribution is assumed to be additive and Gaussian. This model has been well studied in the literature. Various types of threshold processing and parameter selection strategies adapted to specific practical applications have been proposed. The analysis of errors arising from the application of threshold processing methods is an important practical task, since it allows one to evaluate the quality of both the methods themselves and the equipment used for processing. Most of the efforts in this area are devoted to studying the asymptotic order of the theoretical risk, which cannot be calculated in practice because it explicitly depends on the theoretical unobservable noise-free signal. However, a statistical risk estimate based only on the observed data can also be used to assess the quality of noise reduction methods.

In this paper, a model of a signal contaminated with additive Gaussian noise is considered, and the general formulation of the thresholding problem with threshold functions belonging to a certain special class is discussed. A number of examples of known threshold functions from the considered class are given, which are used in such problems as ECG and EEG signal processing, the processing of noisy images, electromyography, noise reduction in audio signals, seismic signal analysis, etc. The specific type of thresholding function for practical use depends on the application and the amount of available data. Tuning additional parameters of the specified thresholding methods improves their performance in some specific applications. For example, [10,11] present numerical studies demonstrating the effectiveness of hybrid thresholding and hyperbolic tangent-based thresholding in ECG signal processing.

This paper obtains lower bounds for thresholds that minimize the unbiased risk estimate. It also provides conditions under which this risk estimate is asymptotically normal and strongly consistent. The authors plan to conduct a number of further studies in this direction. In particular, we will seek to obtain results concerning the construction of confidence intervals and estimates of the convergence rate, which, in turn, will allow us to obtain specific error values in signal processing for a wide range of threshold processing methods. It is also planned to consider similar models for dependent sample components and non-Gaussian noise.