# Supervised and Unsupervised Subband Adaptive Denoising Frameworks with Polynomial Threshold Function

Authors: Tierui Gong; Zhijia Yang; Gengshan Wang; Ping Jiao
Journal: Mathematical Problems in Engineering
 (2017)

---

## Abstract

Unlike inflexible structure of soft and hard threshold function, a unified linear matrix form with flexible structure for threshold function is proposed. Based on the unified linear flexible structure threshold function, both supervised and unsupervised subband adaptive denoising frameworks are established. To determine flexible coefficients, a direct mean-square error (MSE) minimization is conducted in supervised denoising while Stein’s unbiased risk estimate as a MSE estimate is minimized in unsupervised denoising. The SURE rule requires no hypotheses or a priori knowledge about clean signals. Furthermore, we discuss conditions to obtain optimal coefficients for both supervised and unsupervised subband adaptive denoising frameworks. Applying an Odd-Term Reserving Polynomial (OTRP) function as concrete threshold function, simulations for polynomial order, denoising performance, and noise effect are conducted. Proper polynomial order and noise effect are analyzed. Both proposed methods are compared with soft and hard based denoising technologies—VisuShrink, SureShrink, MiniMaxShrink, and BayesShrink—in denoising performance simulation. Results show that the proposed approaches perform better in both MSE and signal-to-noise ratio (SNR) sense.


---

## Body

## 1. Introduction

Signals are usually corrupted by noise in capturing and transmission stages due to environment disturbance and device error. Signal denoising has become an important research topic for a long time and a wide variety of denoising methods have been proposed. Due to their effectiveness and good performance, wavelet threshold methods have become a powerful tool for denoising problems since Donoho and several others’ fundamental works. The main purpose of these methods is to estimate a wide class of functions in some smoothness spaces from their corrupted versions [1]. The different energy distribution property between smoothness spaces’ functions and noise makes these methods effective: the energy of a function in smoothness is often concentrated on few coefficients while noise is spread on all coefficients.Of the various wavelet threshold schemes, soft and hard based threshold methods are the most popular technologies and have been theoretically verified by Donoho and Johnstone [2]. They gave a near optimal threshold value T=2σ2log⁡N in minimax sense, where the threshold is a function of noise variance σ and the number of samples N. This universal thresholding method is known as VisuShrink [3]. In MSE sense, better threshold methods have been proposed including SureShrink [4] and BayesShrink [5]. SureShrink gets the threshold via minimizing Stein’s unbiased risk estimate [6]. Assuming a generalized Gaussian distribution for wavelet coefficients, BayesShrink obtains the threshold through a Bayesian framework. Different a priori assumptions about the statistical distribution of wavelet coefficients, such as Gaussian scale mixture [7], Mixture of Laplacians [8], and alpha-stable [9], are also proposed, while the dependency on statistical distribution restricts their flexibility.Soft and hard based threshold schemes suffer their own flaws due to inherent defects of soft and hard threshold functions. For soft threshold schemes, systematical biased estimation could happen, while hard threshold schemes are less biased but less sensitive to small perturbations in the data. In addition, the more important drawback is that soft and hard threshold functions do not have continuous derivatives. Various improvements had been proposed by exploring new threshold functions [1, 10–17], but the nonnegative garrote-like functions [10–14] are still not differentiable. Zhang [1, 15], Nasri and Nezamabadi-pour [16], and Wu et al. [17], respectively, proposed a series of threshold functions with adjustable parameters. The differentiable property makes these threshold functions suitable for gradient based minimization problems.Most of mentioned methods above depend on a single parameter threshold, which makes those methods very sensitive to threshold value and lack of freedom. More flexible and convenient strategies had been proposed by Luisier et al. [18] and Smith et al. [19]. Both employ an idea of linear combination: the threshold function is formed by a linear combination of a set of parameters. The parameters are determined via minimizing a MSE problem. Compared to gradient based minimization, this minimization is linear and is easy to solve. In this paper, we adopt the linear combination idea and propose a unified matrix form with flexible structure for threshold function. Based on the unified flexible structure threshold function, both supervised and unsupervised subband adaptive denoising frameworks are established. In supervised denoising, a direct MSE minimization is conducted while in unsupervised denoising, we apply Stein’s unbiased risk estimate as MSE estimation. The SURE rule requires no hypotheses or a priori knowledge about clean signals. Furthermore, we discuss conditions for optimal solution for both supervised and unsupervised denoising frameworks.Our contributions can be summarized as follows:(1) a unified linear matrix form with flexible structure for threshold function and a concrete OTRP function are proposed; (2) both supervised and unsupervised denoising frameworks are established; (3) conditions for guaranteeing optimal solution for minimizing problems are discussed and provided.The paper is organized as follows. In Section2, both supervised and unsupervised denoising frameworks are introduced. In Section 3, frameworks employing proposed OTRP function are compared with soft and hard based denoising methods and comparison results are presented. Conclusions are made in the final.


## 2. Proposed Approaches

### 2.1. Problem Settings and Denoising Scheme

In time domain, it is assumed that the clean signalx is additively corrupted by noise e to produce the noisy signal y in the form of(1)y=x+e, where y=yii=1N, x=xii=1N, and e=eii=1N are discrete time series. Additive Gaussian white noise (AGWN) is only considered, with a zero mean and a σ2 variance; that is, ei~N0,σ2.We only consider orthonormal wavelet transform (OWT), which keeps energy conservation. In OWT domain, the AGWN is still Gaussian. Relationship of wavelet coefficients in vector form is(2)Y=X+E,where Y=Wy, X=Wx, and E=We. W is discrete wavelet transform (DWT) matrix. Each vector V(V=Y,X,E) is stacked by both detail and approximation wavelet coefficients as similarly shown in Figure 1. Each detail is a subband. Let Vj(1≤j≤J) denote the jth detail subband with a vector length of Kj. J represents the total number of detail subbands and also the number of decomposition level. There holds the relationship Kj=N/2j.Figure 1
A three-level wavelet decomposition (each color represents a level). There are three detail subbands each with a length ofKj(1≤j≤3) and one approximation subband with a length of K3.We measure the quality of denoising by mean-square error (MSE) defined as(3)MSEx,x^=εx-x^22,where ε· is the expectation operator. An estimate of MSE is given as(4)MSE^x,x^=1Nx-x^Tx-x^=1N∑i=1Nxi-x^i2.By considering OWT, energy conservation property guarantees that the total MSE is a weighted sum of MSE in each individual subband, which is in the form of(5)MSEx,x^=MSEX,X^=∑j=1Jωj×MSEXj,X^j,where ωj is a weight of corresponding subband. It only needs to minimize each individual subband MSE to get minimum of total MSE. In the rest of paper, we ignore subband index j and give a unified notation MSE(X,X^) to represent MSE in any subband.Wavelet domain denoising generally consists of a three-stage procedure. First, perform DWT on noisy signal. Then a threshold function is applied to wavelet domain coefficients. Finally, denoised signal is obtained through an inverse DWT (IDWT). We adopt this three-stage procedure to perform our proposed denoising methods. Each subband is adaptively denoised by supervised or unsupervised denoising method as shown in Figure2. That makes our approaches work in a subband adaptive manner. There are two input signals in Figure 2 where d represents the desired signal for supervised denoising and y denotes noisy signal. Corresponding uppercase letters of d and y denote their subbands. Both methods output an estimate for clean signal x.Figure 2
Proposed wavelet domain denoising schemes.


### 2.2. Proposed Threshold Function

The general problem can be formulated as a construction of threshold function that minimizes the MSE. By employing linear combination idea, an estimate of a clean signal can be represented by a linear combination of atoms; atoms mean the columns of a matrix. Involving wavelet transform, an estimate of clean wavelet coefficients can be represented by a linear combination of atoms which are decided by noisy wavelet coefficients. A general threshold function can be demonstrated as(6)θ=Ψa,where θ is an estimate for clean wavelet coefficients. Ψ=φ1φ2⋯φM is a matrix composed of M atoms and each atom φm=φm(Y1)φm(Y2)⋯φm(YK)T is decided by a function of noisy wavelet coefficients. Dimension of an atom is decided by subband wavelet coefficients number Kj=N/2j: that is, K=Kj when denoising subband j. Besides, atom dimension K does not have direct relationship with atom number M. A more deep connection is discussed in Section 2.5. Furthermore, atom number M has nothing to do with detail subband number or decomposition level J. a=a1a2⋯aMT is a weight vector with each element having a weight coefficient for corresponding atom, so θ is a weighted sum of atoms and can be denoted as θ=θ(Y1)θ(Y2)⋯θ(YK)T.In this paper, the adopted concrete threshold function is OTRP function with orderM, denoted as(7)θY=∑m=1MamY2m-1.This polynomial function could fit almost any curves along with varying polynomial order. The most important advantages are its differentiability and its linear flexible structure over soft and hard threshold functions. Clearly, (7) conforms with (6) where atom is φm=Y12m-1Y22m-1⋯YK2m-1T,(1≤m≤M). A comparison among soft and hard threshold function and a realization of a three-order θ(Y) are shown in Figure 3.Figure 3
Visual comparison among soft and hard threshold function and a realization of a three-orderθ(Y) with coefficients a3a2a1=0.000640.0160.2.


### 2.3. Supervised Denoising

In supervised denoising, a desired signald is available; thus, a direct MSE minimization problem can be carried out. As (5) indicates, total MSE is a weighted sum of MSE in each individual subband. It only needs to minimize a subband MSE to determine its flexible parameter a in (6). With prior knowledge about desired signal d, a direct MSE can be written as(8)MSED,θ=εD-θ22=εD-θTD-θ=εDTD-DTθ-θTD+θTθ=εDTD-2εθTD+εθTθ=εDTD-2aTεΨTD+aTεΨTΨa.The optimal parameter aopt can be obtained through ∂MSE(D,θ)/∂a=0 on condition that ΨTΨ is positive semidefinite and invertible. These conditions are discussed in Section 2.5. We first derive the optimal solution form as(9)∂MSED,θ∂a=2εΨTΨa-2εΨTD=0⟹(10)aopt=εΨTΨ-1εΨTD.As seen from (10), the optimal parameter aopt is affected by matrix Ψ and desired signal D. The estimates of εΨTΨ and εΨTD are obtained via an estimate of expectation operator in (4): (11)ε^ΨTΨ=ε^φ1Tφ1ε^φ1Tφ2⋯ε^φ1TφMε^φ2Tφ1ε^φ2Tφ2⋯ε^φ2TφM⋮⋮⋱⋮ε^φMTφ1ε^φMTφ2⋯ε^φMTφM=1Kφ1Tφ1φ1Tφ2⋯φ1TφMφ2Tφ1φ2Tφ2⋯φ2TφM⋮⋮⋱⋮φMTφ1φMTφ2⋯φMTφM=1KΨTΨ,ε^ΨTD=ε^φ1TDε^φ2TD⋮ε^φMTD=1Kφ1TDφ2TD⋮φMTD=1KΨTD.ε^{ΨTΨ}-1=K(ΨTΨ)-1 is obtained from ε^{ΨTΨ}=ΨTΨ/K; thus, coefficient 1/K is eliminated by K. So, the estimate of aopt is obtained as(12)aoptest=ΨTΨ-1ΨTD.


### 2.4. Unsupervised Denoising

In practice, when a desired signal is not available, a direct MSE minimization is impossible. Constructing an estimate of MSE seems a reasonable choice. A practical approach is Stein’s unbiased risk estimate [6]. This part mainly introduces a SURE based unsupervised denoising.


#### 2.4.1. Stein’s Unbiased Risk Estimate (SURE)

This section mainly introduces Stein’s theorem stated in [6, 20] and its tailored version.Theorem 1.
Y
=
Y
1
Y
2
⋯
Y
K is a normal random vector with mean X=X1X2⋯XK and identity as covariance matrix. Let Y+h(Y) be an estimate of X, where function(13)hY=h1Yh2Y⋯hKY:RK⟶RK,hkY:RK⟶R1,∀j=1,2,…,kis (weakly) differentiable. If(14)∇·hY=∑k=1K∂hk∂Yk,ε∇·hY<∞,then(15)R=εY+hY-X2=K+εhY2+2∇·hY.Stein’s unbiased risk estimate of MSEY+hY,X is(16)R^=K+hY2+2∇·hY.The proposed structure of threshold function does not strictly satisfy conditions ofh(Y) in Theorem 1. Because each element of θ in (6) is R→R mapping while each element in h(Y) is RK→R1 mapping. Therefore, a tailored version of SURE for θ is stated in the following theorem.Theorem 2.
Y
=
Y
1
Y
2
⋯
Y
K is a normal random vector with mean X=X1X2⋯XK and σ2I is covariance matrix. Let θ=θY1θY2⋯θYK be an estimate of X, where θy:R→R is (weakly) differentiable. If(17)∇·θ=∑k=1Kθ′Yk,ε∇·θ<∞,then (18)MSEX,θ=εXTX-2εθTY+2σ2ε∇·θ+εθTθ. The unbiased estimate of MSE(X,θ) is(19)R^=1KXTX+1K-2θTY+2σ21Tθ′+θTθ.Proof.
Combining (2) and (8), the MSE is represented as(20)MSEX,θ=εXTX-2εθTX+εθTθ=εXTX-2εθTY-E+εθTθ=εXTX-2εθTY+2εθTE+εθTθ.The only unknown part in the above equation is εθTE while equality εEkθYk=σ2εθ′Yk, 1≤k≤K proved in [18] ensures that εθTE can be transformed to a known part. Concrete deduction formula is in the following equation:(21)εθTE=ε∑k=1KEkθYk=∑k=1KεEkθYk=σ2∑k=1Kεθ′Yk=σ2ε∑k=1Kθ′Yk=σ2ε∇·θ.So, MSE can be rewritten as(22)MSEX,θ=εXTX-2εθTY+2σ2ε∇·θ+εθTθ.Estimate the expectation of εXTX, εθTY, ε∇·θ, and εθTθ via(23)ε^XTX=1KXTX,ε^θTY=1KθTY,ε^∇·θ=1K1Tθ′,ε^θTθ=1KθTθ.It is obvious that (19) is an unbiased estimate of MSE(X,θ).
Proof is completed.


#### 2.4.2. Coefficients Determination

From (6), we first obtain element derivative of θ in the form of(24)θ′Yk=aTφ1′Ykφ2′Yk⋯φM′YkT.So, the formula of ∇·θ can be deduced as(25)∇·θ=∑k=1Kθ′Yk=aT∑k=1Kφ1′Yk∑k=1Kφ2′Yk⋯∑k=1KφM′YkT=aT1Tφ1′1Tφ2′⋯1TφM′T=aT1TΨ′T.Thus,(26)ε∇·θ=aTε1TΨ′T.Combining (6) and (18) and (26), MSE(X,θ) can be rewritten as(27)MSEX,θ=εXTX-2aTεΨTY+2σ2aTε1TΨ′T+aTεΨTΨa.aopt is obtained through ∂MSE(X,θ)/∂a=0 on condition that ΨTΨ is positive semidefinite and invertible. These conditions the same as in supervised denoising are discussed in Section 2.5. Here, we first derive the optimal solution form as(28)∂MSEX,θ∂a=2εΨTΨa-2εΨTY+2σ2ε1TΨ′T=0⟹(29)aopt=εΨTΨ-1εΨTY-σ2ε1TΨ′T.An estimate of aopt is obtained via Theorem 2. From (6), we first obtain pointwise derivative of θ that can be denoted as(30)θ′=Ψ′a.Then (19) can be rewritten as(31)R^=1KXTX+1K-2aTΨTY+2σ21TΨ′a+aTΨTΨa=1KXTX+1K-2aTΨTY+2σ2aT1TΨ′T+aTΨTΨa.The result can be found by setting the derivative over a to zero. Still positive semidefinite and invertible must be satisfied for ΨTΨ. Estimate for optimal solution is achieved by(32)∂R^∂a=1K-2ΨTY+2σ21TΨ′T+2ΨTΨa=0⟹(33)aoptest=ΨTΨ-1ΨTY-σ21TΨ′T.


#### 2.4.3. Noise Variance Estimate

As it can be seen from (29) and (33), noise variance σ is needed in unsupervised denoising. In practice, there are various methods that can provide estimate for AGWN variance, such as the most popular median absolute deviation (MAD) [3], SVD-based [21], and block-based [22, 23] methods. MAD implements simply and achieves accurate estimate. For simplicity, we adopt MAD to estimate noise variance. The estimate is denoted as(34)σ^=MAD0.6745,MAD=MedianY1-Y-,…,YK-Y-,where Y- is the mean value in corresponding subband.


### 2.5. Optimal Solution Guarantee Discussion

Optimal solution assurance conditions for both supervised and unsupervised denoising are discussed in this section. To guarantee optimal solution ofa by setting ∂MSE(X,θ)/∂a to zero, the matrix ΨTΨ must be positive semidefinite. According to the definition of positive semidefinite matrix, A is positive semidefinite when xTAx≥0 for any vector x. So, it is easy to verify positive semidefinite of ΨTΨ; that is, for any vector x, there always exists(35)xTΨTΨx=Ψx,Ψx=Ψx22≥0.Thus, it can be concluded that ΨTΨ is always positive semidefinite. Positive semidefinite property guarantees that MSE(X,θ) is convex while invertible of ΨTΨ makes the optimal solution of a in the form of (10), (12) and (29), (33). It must be full rank to ensure invertible ΨTΨ, so we can deduce that Ψ must be full row or column rank because of rank(Ψ)=rank(ΨTΨ). In practice, we do not use too many atoms to form matrix Ψ and the number of noisy coefficients K in a subband is much greater than atom number M(i.e.,K≫M). In other words, with row number far greater than column number, Ψ is a very “thin” matrix. Thus, matrix Ψ tends to be full column rank and matrix ΨTΨ tends to be full rank.


## 2.1. Problem Settings and Denoising Scheme

In time domain, it is assumed that the clean signalx is additively corrupted by noise e to produce the noisy signal y in the form of(1)y=x+e, where y=yii=1N, x=xii=1N, and e=eii=1N are discrete time series. Additive Gaussian white noise (AGWN) is only considered, with a zero mean and a σ2 variance; that is, ei~N0,σ2.We only consider orthonormal wavelet transform (OWT), which keeps energy conservation. In OWT domain, the AGWN is still Gaussian. Relationship of wavelet coefficients in vector form is(2)Y=X+E,where Y=Wy, X=Wx, and E=We. W is discrete wavelet transform (DWT) matrix. Each vector V(V=Y,X,E) is stacked by both detail and approximation wavelet coefficients as similarly shown in Figure 1. Each detail is a subband. Let Vj(1≤j≤J) denote the jth detail subband with a vector length of Kj. J represents the total number of detail subbands and also the number of decomposition level. There holds the relationship Kj=N/2j.Figure 1
A three-level wavelet decomposition (each color represents a level). There are three detail subbands each with a length ofKj(1≤j≤3) and one approximation subband with a length of K3.We measure the quality of denoising by mean-square error (MSE) defined as(3)MSEx,x^=εx-x^22,where ε· is the expectation operator. An estimate of MSE is given as(4)MSE^x,x^=1Nx-x^Tx-x^=1N∑i=1Nxi-x^i2.By considering OWT, energy conservation property guarantees that the total MSE is a weighted sum of MSE in each individual subband, which is in the form of(5)MSEx,x^=MSEX,X^=∑j=1Jωj×MSEXj,X^j,where ωj is a weight of corresponding subband. It only needs to minimize each individual subband MSE to get minimum of total MSE. In the rest of paper, we ignore subband index j and give a unified notation MSE(X,X^) to represent MSE in any subband.Wavelet domain denoising generally consists of a three-stage procedure. First, perform DWT on noisy signal. Then a threshold function is applied to wavelet domain coefficients. Finally, denoised signal is obtained through an inverse DWT (IDWT). We adopt this three-stage procedure to perform our proposed denoising methods. Each subband is adaptively denoised by supervised or unsupervised denoising method as shown in Figure2. That makes our approaches work in a subband adaptive manner. There are two input signals in Figure 2 where d represents the desired signal for supervised denoising and y denotes noisy signal. Corresponding uppercase letters of d and y denote their subbands. Both methods output an estimate for clean signal x.Figure 2
Proposed wavelet domain denoising schemes.


## 2.2. Proposed Threshold Function

The general problem can be formulated as a construction of threshold function that minimizes the MSE. By employing linear combination idea, an estimate of a clean signal can be represented by a linear combination of atoms; atoms mean the columns of a matrix. Involving wavelet transform, an estimate of clean wavelet coefficients can be represented by a linear combination of atoms which are decided by noisy wavelet coefficients. A general threshold function can be demonstrated as(6)θ=Ψa,where θ is an estimate for clean wavelet coefficients. Ψ=φ1φ2⋯φM is a matrix composed of M atoms and each atom φm=φm(Y1)φm(Y2)⋯φm(YK)T is decided by a function of noisy wavelet coefficients. Dimension of an atom is decided by subband wavelet coefficients number Kj=N/2j: that is, K=Kj when denoising subband j. Besides, atom dimension K does not have direct relationship with atom number M. A more deep connection is discussed in Section 2.5. Furthermore, atom number M has nothing to do with detail subband number or decomposition level J. a=a1a2⋯aMT is a weight vector with each element having a weight coefficient for corresponding atom, so θ is a weighted sum of atoms and can be denoted as θ=θ(Y1)θ(Y2)⋯θ(YK)T.In this paper, the adopted concrete threshold function is OTRP function with orderM, denoted as(7)θY=∑m=1MamY2m-1.This polynomial function could fit almost any curves along with varying polynomial order. The most important advantages are its differentiability and its linear flexible structure over soft and hard threshold functions. Clearly, (7) conforms with (6) where atom is φm=Y12m-1Y22m-1⋯YK2m-1T,(1≤m≤M). A comparison among soft and hard threshold function and a realization of a three-order θ(Y) are shown in Figure 3.Figure 3
Visual comparison among soft and hard threshold function and a realization of a three-orderθ(Y) with coefficients a3a2a1=0.000640.0160.2.


## 2.3. Supervised Denoising

In supervised denoising, a desired signald is available; thus, a direct MSE minimization problem can be carried out. As (5) indicates, total MSE is a weighted sum of MSE in each individual subband. It only needs to minimize a subband MSE to determine its flexible parameter a in (6). With prior knowledge about desired signal d, a direct MSE can be written as(8)MSED,θ=εD-θ22=εD-θTD-θ=εDTD-DTθ-θTD+θTθ=εDTD-2εθTD+εθTθ=εDTD-2aTεΨTD+aTεΨTΨa.The optimal parameter aopt can be obtained through ∂MSE(D,θ)/∂a=0 on condition that ΨTΨ is positive semidefinite and invertible. These conditions are discussed in Section 2.5. We first derive the optimal solution form as(9)∂MSED,θ∂a=2εΨTΨa-2εΨTD=0⟹(10)aopt=εΨTΨ-1εΨTD.As seen from (10), the optimal parameter aopt is affected by matrix Ψ and desired signal D. The estimates of εΨTΨ and εΨTD are obtained via an estimate of expectation operator in (4): (11)ε^ΨTΨ=ε^φ1Tφ1ε^φ1Tφ2⋯ε^φ1TφMε^φ2Tφ1ε^φ2Tφ2⋯ε^φ2TφM⋮⋮⋱⋮ε^φMTφ1ε^φMTφ2⋯ε^φMTφM=1Kφ1Tφ1φ1Tφ2⋯φ1TφMφ2Tφ1φ2Tφ2⋯φ2TφM⋮⋮⋱⋮φMTφ1φMTφ2⋯φMTφM=1KΨTΨ,ε^ΨTD=ε^φ1TDε^φ2TD⋮ε^φMTD=1Kφ1TDφ2TD⋮φMTD=1KΨTD.ε^{ΨTΨ}-1=K(ΨTΨ)-1 is obtained from ε^{ΨTΨ}=ΨTΨ/K; thus, coefficient 1/K is eliminated by K. So, the estimate of aopt is obtained as(12)aoptest=ΨTΨ-1ΨTD.


## 2.4. Unsupervised Denoising

In practice, when a desired signal is not available, a direct MSE minimization is impossible. Constructing an estimate of MSE seems a reasonable choice. A practical approach is Stein’s unbiased risk estimate [6]. This part mainly introduces a SURE based unsupervised denoising.


### 2.4.1. Stein’s Unbiased Risk Estimate (SURE)

This section mainly introduces Stein’s theorem stated in [6, 20] and its tailored version.Theorem 1.
Y
=
Y
1
Y
2
⋯
Y
K is a normal random vector with mean X=X1X2⋯XK and identity as covariance matrix. Let Y+h(Y) be an estimate of X, where function(13)hY=h1Yh2Y⋯hKY:RK⟶RK,hkY:RK⟶R1,∀j=1,2,…,kis (weakly) differentiable. If(14)∇·hY=∑k=1K∂hk∂Yk,ε∇·hY<∞,then(15)R=εY+hY-X2=K+εhY2+2∇·hY.Stein’s unbiased risk estimate of MSEY+hY,X is(16)R^=K+hY2+2∇·hY.The proposed structure of threshold function does not strictly satisfy conditions ofh(Y) in Theorem 1. Because each element of θ in (6) is R→R mapping while each element in h(Y) is RK→R1 mapping. Therefore, a tailored version of SURE for θ is stated in the following theorem.Theorem 2.
Y
=
Y
1
Y
2
⋯
Y
K is a normal random vector with mean X=X1X2⋯XK and σ2I is covariance matrix. Let θ=θY1θY2⋯θYK be an estimate of X, where θy:R→R is (weakly) differentiable. If(17)∇·θ=∑k=1Kθ′Yk,ε∇·θ<∞,then (18)MSEX,θ=εXTX-2εθTY+2σ2ε∇·θ+εθTθ. The unbiased estimate of MSE(X,θ) is(19)R^=1KXTX+1K-2θTY+2σ21Tθ′+θTθ.Proof.
Combining (2) and (8), the MSE is represented as(20)MSEX,θ=εXTX-2εθTX+εθTθ=εXTX-2εθTY-E+εθTθ=εXTX-2εθTY+2εθTE+εθTθ.The only unknown part in the above equation is εθTE while equality εEkθYk=σ2εθ′Yk, 1≤k≤K proved in [18] ensures that εθTE can be transformed to a known part. Concrete deduction formula is in the following equation:(21)εθTE=ε∑k=1KEkθYk=∑k=1KεEkθYk=σ2∑k=1Kεθ′Yk=σ2ε∑k=1Kθ′Yk=σ2ε∇·θ.So, MSE can be rewritten as(22)MSEX,θ=εXTX-2εθTY+2σ2ε∇·θ+εθTθ.Estimate the expectation of εXTX, εθTY, ε∇·θ, and εθTθ via(23)ε^XTX=1KXTX,ε^θTY=1KθTY,ε^∇·θ=1K1Tθ′,ε^θTθ=1KθTθ.It is obvious that (19) is an unbiased estimate of MSE(X,θ).
Proof is completed.


### 2.4.2. Coefficients Determination

From (6), we first obtain element derivative of θ in the form of(24)θ′Yk=aTφ1′Ykφ2′Yk⋯φM′YkT.So, the formula of ∇·θ can be deduced as(25)∇·θ=∑k=1Kθ′Yk=aT∑k=1Kφ1′Yk∑k=1Kφ2′Yk⋯∑k=1KφM′YkT=aT1Tφ1′1Tφ2′⋯1TφM′T=aT1TΨ′T.Thus,(26)ε∇·θ=aTε1TΨ′T.Combining (6) and (18) and (26), MSE(X,θ) can be rewritten as(27)MSEX,θ=εXTX-2aTεΨTY+2σ2aTε1TΨ′T+aTεΨTΨa.aopt is obtained through ∂MSE(X,θ)/∂a=0 on condition that ΨTΨ is positive semidefinite and invertible. These conditions the same as in supervised denoising are discussed in Section 2.5. Here, we first derive the optimal solution form as(28)∂MSEX,θ∂a=2εΨTΨa-2εΨTY+2σ2ε1TΨ′T=0⟹(29)aopt=εΨTΨ-1εΨTY-σ2ε1TΨ′T.An estimate of aopt is obtained via Theorem 2. From (6), we first obtain pointwise derivative of θ that can be denoted as(30)θ′=Ψ′a.Then (19) can be rewritten as(31)R^=1KXTX+1K-2aTΨTY+2σ21TΨ′a+aTΨTΨa=1KXTX+1K-2aTΨTY+2σ2aT1TΨ′T+aTΨTΨa.The result can be found by setting the derivative over a to zero. Still positive semidefinite and invertible must be satisfied for ΨTΨ. Estimate for optimal solution is achieved by(32)∂R^∂a=1K-2ΨTY+2σ21TΨ′T+2ΨTΨa=0⟹(33)aoptest=ΨTΨ-1ΨTY-σ21TΨ′T.


### 2.4.3. Noise Variance Estimate

As it can be seen from (29) and (33), noise variance σ is needed in unsupervised denoising. In practice, there are various methods that can provide estimate for AGWN variance, such as the most popular median absolute deviation (MAD) [3], SVD-based [21], and block-based [22, 23] methods. MAD implements simply and achieves accurate estimate. For simplicity, we adopt MAD to estimate noise variance. The estimate is denoted as(34)σ^=MAD0.6745,MAD=MedianY1-Y-,…,YK-Y-,where Y- is the mean value in corresponding subband.


## 2.4.1. Stein’s Unbiased Risk Estimate (SURE)

This section mainly introduces Stein’s theorem stated in [6, 20] and its tailored version.Theorem 1.
Y
=
Y
1
Y
2
⋯
Y
K is a normal random vector with mean X=X1X2⋯XK and identity as covariance matrix. Let Y+h(Y) be an estimate of X, where function(13)hY=h1Yh2Y⋯hKY:RK⟶RK,hkY:RK⟶R1,∀j=1,2,…,kis (weakly) differentiable. If(14)∇·hY=∑k=1K∂hk∂Yk,ε∇·hY<∞,then(15)R=εY+hY-X2=K+εhY2+2∇·hY.Stein’s unbiased risk estimate of MSEY+hY,X is(16)R^=K+hY2+2∇·hY.The proposed structure of threshold function does not strictly satisfy conditions ofh(Y) in Theorem 1. Because each element of θ in (6) is R→R mapping while each element in h(Y) is RK→R1 mapping. Therefore, a tailored version of SURE for θ is stated in the following theorem.Theorem 2.
Y
=
Y
1
Y
2
⋯
Y
K is a normal random vector with mean X=X1X2⋯XK and σ2I is covariance matrix. Let θ=θY1θY2⋯θYK be an estimate of X, where θy:R→R is (weakly) differentiable. If(17)∇·θ=∑k=1Kθ′Yk,ε∇·θ<∞,then (18)MSEX,θ=εXTX-2εθTY+2σ2ε∇·θ+εθTθ. The unbiased estimate of MSE(X,θ) is(19)R^=1KXTX+1K-2θTY+2σ21Tθ′+θTθ.Proof.
Combining (2) and (8), the MSE is represented as(20)MSEX,θ=εXTX-2εθTX+εθTθ=εXTX-2εθTY-E+εθTθ=εXTX-2εθTY+2εθTE+εθTθ.The only unknown part in the above equation is εθTE while equality εEkθYk=σ2εθ′Yk, 1≤k≤K proved in [18] ensures that εθTE can be transformed to a known part. Concrete deduction formula is in the following equation:(21)εθTE=ε∑k=1KEkθYk=∑k=1KεEkθYk=σ2∑k=1Kεθ′Yk=σ2ε∑k=1Kθ′Yk=σ2ε∇·θ.So, MSE can be rewritten as(22)MSEX,θ=εXTX-2εθTY+2σ2ε∇·θ+εθTθ.Estimate the expectation of εXTX, εθTY, ε∇·θ, and εθTθ via(23)ε^XTX=1KXTX,ε^θTY=1KθTY,ε^∇·θ=1K1Tθ′,ε^θTθ=1KθTθ.It is obvious that (19) is an unbiased estimate of MSE(X,θ).
Proof is completed.


## 2.4.2. Coefficients Determination

From (6), we first obtain element derivative of θ in the form of(24)θ′Yk=aTφ1′Ykφ2′Yk⋯φM′YkT.So, the formula of ∇·θ can be deduced as(25)∇·θ=∑k=1Kθ′Yk=aT∑k=1Kφ1′Yk∑k=1Kφ2′Yk⋯∑k=1KφM′YkT=aT1Tφ1′1Tφ2′⋯1TφM′T=aT1TΨ′T.Thus,(26)ε∇·θ=aTε1TΨ′T.Combining (6) and (18) and (26), MSE(X,θ) can be rewritten as(27)MSEX,θ=εXTX-2aTεΨTY+2σ2aTε1TΨ′T+aTεΨTΨa.aopt is obtained through ∂MSE(X,θ)/∂a=0 on condition that ΨTΨ is positive semidefinite and invertible. These conditions the same as in supervised denoising are discussed in Section 2.5. Here, we first derive the optimal solution form as(28)∂MSEX,θ∂a=2εΨTΨa-2εΨTY+2σ2ε1TΨ′T=0⟹(29)aopt=εΨTΨ-1εΨTY-σ2ε1TΨ′T.An estimate of aopt is obtained via Theorem 2. From (6), we first obtain pointwise derivative of θ that can be denoted as(30)θ′=Ψ′a.Then (19) can be rewritten as(31)R^=1KXTX+1K-2aTΨTY+2σ21TΨ′a+aTΨTΨa=1KXTX+1K-2aTΨTY+2σ2aT1TΨ′T+aTΨTΨa.The result can be found by setting the derivative over a to zero. Still positive semidefinite and invertible must be satisfied for ΨTΨ. Estimate for optimal solution is achieved by(32)∂R^∂a=1K-2ΨTY+2σ21TΨ′T+2ΨTΨa=0⟹(33)aoptest=ΨTΨ-1ΨTY-σ21TΨ′T.


## 2.4.3. Noise Variance Estimate

As it can be seen from (29) and (33), noise variance σ is needed in unsupervised denoising. In practice, there are various methods that can provide estimate for AGWN variance, such as the most popular median absolute deviation (MAD) [3], SVD-based [21], and block-based [22, 23] methods. MAD implements simply and achieves accurate estimate. For simplicity, we adopt MAD to estimate noise variance. The estimate is denoted as(34)σ^=MAD0.6745,MAD=MedianY1-Y-,…,YK-Y-,where Y- is the mean value in corresponding subband.


## 2.5. Optimal Solution Guarantee Discussion

Optimal solution assurance conditions for both supervised and unsupervised denoising are discussed in this section. To guarantee optimal solution ofa by setting ∂MSE(X,θ)/∂a to zero, the matrix ΨTΨ must be positive semidefinite. According to the definition of positive semidefinite matrix, A is positive semidefinite when xTAx≥0 for any vector x. So, it is easy to verify positive semidefinite of ΨTΨ; that is, for any vector x, there always exists(35)xTΨTΨx=Ψx,Ψx=Ψx22≥0.Thus, it can be concluded that ΨTΨ is always positive semidefinite. Positive semidefinite property guarantees that MSE(X,θ) is convex while invertible of ΨTΨ makes the optimal solution of a in the form of (10), (12) and (29), (33). It must be full rank to ensure invertible ΨTΨ, so we can deduce that Ψ must be full row or column rank because of rank(Ψ)=rank(ΨTΨ). In practice, we do not use too many atoms to form matrix Ψ and the number of noisy coefficients K in a subband is much greater than atom number M(i.e.,K≫M). In other words, with row number far greater than column number, Ψ is a very “thin” matrix. Thus, matrix Ψ tends to be full column rank and matrix ΨTΨ tends to be full rank.


## 3. Experiments and Results

### 3.1. Experiment Settings

Standard experiment settings for latter simulations are given. The adopted signals in simulations are “blocks,” “quadchirp,” “multitone,” and “multiband” with AGWN. In order to explain the superiority, the proposed methods are also tested by real audio signal. All those clean signals are shown in Figure4 with each signal having length of 8192 for the first four signals and a length of 65536 for audio signal. The adopted orthogonal wavelet is “sym8” and decomposition level is J=5. So the jth detail subband has a vector length of N/2j(1≤j≤5), where N can be 8192 or 65536.Figure 4
Clean “blocks,” “quadchirp,” “multitone,” and “multiband,” each signal with length of 8192 and audio signal with a length of 65536.In supervised denoising, coefficienta is determined by matrix Ψ and desired subband D in each individual subband. Desired subbands come from desired signal by OWT. Desired signal can be clean or noisy. In practice, it may only have access to noisy signal with known noise level; thus, we adopt it as desired signal to determine coefficient a and then use a to derive the denoised version of noisy signal with unknown noise level. A remarkable notice is that the desired signal and noisy signal are, respectively, the same clean signal corrupted by known and unknown noise level. In unsupervised denoising, coefficient a is determined by matrix Ψ, noisy signal subband Y, and noise variance. There are no requirements for clean signal or desired signal. For simplicity, an estimate for noise variance is provided by MAD.Applying the proposed OTRP function as concrete threshold function, supervised and unsupervised denoising methods are compared with several other available techniques, that is, VisuShrink, SureShrink, MiniMaxShrink, and BayesShrink. Wavelab 850 is applied to perform the other four techniques; their noise variance estimates are also provided by MAD. Denoising quality is measured by both MSE and SNR, which are defined as(36)MSE=1N∑i=1Nx^i-xi2,SNR=10log⁡∑i=1Nxi2∑i=1Nx^i-xi2.


### 3.2. Polynomial Order Simulation

Before carrying out the final denoising, proper polynomial order needs be explored. Due to the fact that polynomial order is equal to atom number of matrixΨ, it needs to select proper polynomial order to guarantee full column rank of Ψ and thus full rank of ΨTΨ. Although it has already been discussed in Section 2.5 that subband wavelet coefficients number is far greater than polynomial order and ΨTΨ tends to be full rank, simulations are still needed. Whether ΨTΨ is full rank or not is determined by polynomial order and subband wavelet coefficients. Thus, a probability statistical table is made through simulations. The achieved probabilities make up Table 1 where PO represents polynomial order and SB denotes subband. The probabilities reveal how often ΨTΨ is full rank in different PO for different signals corrupted with different noise level.Table 1
Ψ
T
Ψ full rank probabilities statistical table.
Signals
σ
PO1
PO2
PO3
PO4
PO5
PO6
SB1
SB2
SB3
SB4
SB5
SB1
SB2
SB3
SB4
SB5
SB1
SB2
SB3
SB4
SB5
SB1
SB2
SB3
SB4
SB5
SB1
SB2
SB3
SB4
SB5
SB1
SB2
SB3
SB4
SB5
Blocks
0.6
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
0
100
100
100
100
0
100
100
100
54
0
0.8
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
0
100
100
100
100
0
100
100
97
51
0
1.0
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
0
100
100
100
100
0
100
100
95
50
0
1.2
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
1
100
100
100
100
0
100
100
89
34
0
1.4
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
1
100
100
100
100
0
96
97
80
22
0
1.6
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
2
100
100
100
97
0
62
69
44
18
0
1.8
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
1
100
100
100
95
0
1
8
7
4
0
2.0
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
1
100
100
100
93
0
0
0
0
0
0
Quadchirp
0.6
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
98
13
0.8
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
87
10
1.0
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
99
64
5
1.2
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
99
100
98
79
38
3
1.4
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
99
91
68
47
11
1
1.6
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
99
90
42
23
15
3
1
1.8
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
98
84
0
0
0
0
1
2.0
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
99
91
74
0
0
0
0
0
Multitone
0.6
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
94
30
100
0.8
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
76
0
100
1.0
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
99
39
0
98
1.2
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
97
2
0
94
1.4
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
99
100
99
66
0
0
85
1.6
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
96
95
100
69
9
0
0
73
1.8
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
94
73
98
6
0
0
0
42
2.0
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
85
41
93
0
0
0
0
17
Multiband
0.6
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
91
95
40
0.8
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
74
59
23
1.0
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
99
100
100
100
99
56
10
14
1.2
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
99
100
99
100
70
15
3
4
1.4
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
99
100
99
91
9
1
0
4
1.6
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
95
99
93
43
0
0
0
0
1.8
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
99
99
93
90
0
0
0
0
0
2.0
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
98
85
84
85
0
0
0
0
0Six orders of PO are tested in Table1. Due to 5-level wavelet decomposition, there are five detail SBs to be processed. Each SB corresponds to a ΨTΨ and different PO corresponds to different dimension of ΨTΨ. So, there are 30 combinations shown in Table 1 header where each PO contains five SBs. In x-axis direction, Table 1 reflects how PO influences matrix full rank probabilities. In y-axis direction of Table 1, four signals are tested with each corrupted by different noise variance σ. This reflects how noise level influences matrix full rank probabilities. For example, the italic number shown in Table 1 means matrix full rank probability is 100%. This matrix is determined by PO2 and “blocks” signal’s SB3 in σ=1.6 noise level. All probabilities are obtained by counting matrix full rank times in a 100-trail.From Table1, PO1–PO3 always ensure 100% matrix full rank for four signals in all tested noise level. Saying 100% matrix full rank for a signal we mean 100% matrix full rank for every subband of this signal. PO4 only cannot ensure “blocks” due to almost 0% probability for SB5 (bold underlined numbers, bold numbers indicate the probabilities under 100%). For PO5 and PO6, there are more and more bold numbers than PO4. In horizontal direction, the trend is higher PO lower probability for satisfying matrix full rank. In vertical direction, higher noise level tends to lead lower probability for full rank, but PO is a dominant factor. From statistical analysis, we restrict our option from PO1 to PO3.Visual comparisons of PO1 to PO3 in supervised denoising are shown in Figure5. As it can be seen, in both MSE and SNR sense PO1 is inferior to PO2 and PO3 in all four signals. PO2 is inferior to PO3 in blocks and almost equals PO3 in the three remaining signals. In unsupervised denoising, similar conclusions can be made and are not repetitive state. In the end, we choose PO2 for final denoising simulation.Figure 5
PO1 to PO3 comparisons in supervised denoising condition.


### 3.3. Denoising Performance Comparison

As mentioned above, comparisons are made among supervised and unsupervised denoising methods and VisuShrink, SureShrink, MiniMaxShrink, and BayesShrink. In both MSE and SNR sense, denoising performance comparisons in global and local views with different noise varianceσ are presented in Figure 6. σ changes in range of 0.5–2.0 with step 0.1 and each data point plotted in Figure 6 is a 50-average result. A remarkable notice is that we adopt a clean signal corrupted by noise in σ=0.5 level as desired signal for supervised denoising. Moreover, MSEs in Figure 6 are plotted in log scale and SNRs in Figure 6 are plotted in linear scale.Figure 6
Comparison of both MSE and SNR changing with noise varianceσ for “blocks,” “quadchirp,” “multitone,” and “multiband” denoising.From global and local views of Figure6, supervised denoising method achieves better MSE and SNR quality than unsupervised denoising method and BayesShrink for all four noisy signals. All these three methods have better MSE and SNR performance than VisuShrink, SureShrink, and MiniMaxShrink. In noisy “blocks” processing, BayesShrink is slightly better than unsupervised denoising method when σ is greater than 0.9 while poorer when σ is less than 0.9 in both MSE and SNR sense. However, unsupervised denoising method has better MSE and SNR performance than BayesShrink in noisy quadchirp, multitone, and multiband processing. Thus, a conclusion can be made that the proposed supervised and unsupervised denoising methods show their advantages.Our approaches are also validated for real-world signals. We test supervised and unsupervised denoising approaches by an AGWN corrupted real audio signal. In both MSE and SNR sense, comparisons with changingσ in range of 0.1–0.5 are shown in Figure 7 and each data point plotted is still a 50-average result. We adopt clean audio signal corrupted by noise in σ=0.5 level as desired signal for supervised denoising. In addition, MSEs are plotted in log scale and SNRs are plotted in linear scale. From global and local views of Figure 7, the proposed supervised and unsupervised denoising methods have smaller MSEs and higher SNRs than VisuShrink, SureShrink, MiniMaxShrink, and BayesShrink. A conclusion can still be made that the proposed supervised and unsupervised denoising methods show their advantages in real-world condition.Figure 7
Denoising performance comparisons in both MSE and SNR sense with changing noise varianceσ for real audio signal denoising.Remark. All denoising methods except for supervised denoising involved in Figures 6 and 7 require noise variance σ, and estimates of it are achieved by MAD for all those methods. So, proposed unsupervised method and the other four techniques work under the same condition, which reflects superiority of the proposed method. In practice, MAD with simple implementation and estimate accuracy do satisfy our requirements.


### 3.4. Noise Effect Analysis

In this section, different noise levels of desired signal for supervised denoising are analyzed. It only aims at real audio signal. In simulation, desired signal noise levels are set in rangeσ=0.1 to σ=2.5 with equal step 0.4. Both MSE and SNR results are shown in Figure 8 and each data point plotted in Figure 8 is averaged by 50 times. It must be noticed that σ on x-axis represents noise level of noisy signal to be processed. All desired signal noise levels are compared with BayesShrink denoising mothed. From Figure 8 it can be seen that all σ=0.1 to σ=2.5 desired signal noise levels have better denoising performance than BayesShrink in MSE and SNR sense when σ≥0.18 on x-axis. While when σ<0.18 on x-axis BayesShrink does better than σ=2.1 and σ=2.5 and worse than σ=0.1 to σ=1.7 in MSE and SNR sense. In addition, denoising performance becomes worse when desired signal noise level gets higher. In conclusion, if some noisy desired signal (even higher noise level than the signal to be processed) is available at hand, it is still possible to achieve good denoising performance.Figure 8
Different desired signal noise level and BayesShrink denoising performance comparisons in both MSE and SNR sense with changing noise varianceσ for real audio signal denoising.


## 3.1. Experiment Settings

Standard experiment settings for latter simulations are given. The adopted signals in simulations are “blocks,” “quadchirp,” “multitone,” and “multiband” with AGWN. In order to explain the superiority, the proposed methods are also tested by real audio signal. All those clean signals are shown in Figure4 with each signal having length of 8192 for the first four signals and a length of 65536 for audio signal. The adopted orthogonal wavelet is “sym8” and decomposition level is J=5. So the jth detail subband has a vector length of N/2j(1≤j≤5), where N can be 8192 or 65536.Figure 4
Clean “blocks,” “quadchirp,” “multitone,” and “multiband,” each signal with length of 8192 and audio signal with a length of 65536.In supervised denoising, coefficienta is determined by matrix Ψ and desired subband D in each individual subband. Desired subbands come from desired signal by OWT. Desired signal can be clean or noisy. In practice, it may only have access to noisy signal with known noise level; thus, we adopt it as desired signal to determine coefficient a and then use a to derive the denoised version of noisy signal with unknown noise level. A remarkable notice is that the desired signal and noisy signal are, respectively, the same clean signal corrupted by known and unknown noise level. In unsupervised denoising, coefficient a is determined by matrix Ψ, noisy signal subband Y, and noise variance. There are no requirements for clean signal or desired signal. For simplicity, an estimate for noise variance is provided by MAD.Applying the proposed OTRP function as concrete threshold function, supervised and unsupervised denoising methods are compared with several other available techniques, that is, VisuShrink, SureShrink, MiniMaxShrink, and BayesShrink. Wavelab 850 is applied to perform the other four techniques; their noise variance estimates are also provided by MAD. Denoising quality is measured by both MSE and SNR, which are defined as(36)MSE=1N∑i=1Nx^i-xi2,SNR=10log⁡∑i=1Nxi2∑i=1Nx^i-xi2.


## 3.2. Polynomial Order Simulation

Before carrying out the final denoising, proper polynomial order needs be explored. Due to the fact that polynomial order is equal to atom number of matrixΨ, it needs to select proper polynomial order to guarantee full column rank of Ψ and thus full rank of ΨTΨ. Although it has already been discussed in Section 2.5 that subband wavelet coefficients number is far greater than polynomial order and ΨTΨ tends to be full rank, simulations are still needed. Whether ΨTΨ is full rank or not is determined by polynomial order and subband wavelet coefficients. Thus, a probability statistical table is made through simulations. The achieved probabilities make up Table 1 where PO represents polynomial order and SB denotes subband. The probabilities reveal how often ΨTΨ is full rank in different PO for different signals corrupted with different noise level.Table 1
Ψ
T
Ψ full rank probabilities statistical table.
Signals
σ
PO1
PO2
PO3
PO4
PO5
PO6
SB1
SB2
SB3
SB4
SB5
SB1
SB2
SB3
SB4
SB5
SB1
SB2
SB3
SB4
SB5
SB1
SB2
SB3
SB4
SB5
SB1
SB2
SB3
SB4
SB5
SB1
SB2
SB3
SB4
SB5
Blocks
0.6
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
0
100
100
100
100
0
100
100
100
54
0
0.8
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
0
100
100
100
100
0
100
100
97
51
0
1.0
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
0
100
100
100
100
0
100
100
95
50
0
1.2
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
1
100
100
100
100
0
100
100
89
34
0
1.4
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
1
100
100
100
100
0
96
97
80
22
0
1.6
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
2
100
100
100
97
0
62
69
44
18
0
1.8
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
1
100
100
100
95
0
1
8
7
4
0
2.0
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
1
100
100
100
93
0
0
0
0
0
0
Quadchirp
0.6
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
98
13
0.8
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
87
10
1.0
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
99
64
5
1.2
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
99
100
98
79
38
3
1.4
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
99
91
68
47
11
1
1.6
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
99
90
42
23
15
3
1
1.8
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
98
84
0
0
0
0
1
2.0
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
99
91
74
0
0
0
0
0
Multitone
0.6
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
94
30
100
0.8
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
76
0
100
1.0
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
99
39
0
98
1.2
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
97
2
0
94
1.4
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
99
100
99
66
0
0
85
1.6
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
96
95
100
69
9
0
0
73
1.8
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
94
73
98
6
0
0
0
42
2.0
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
85
41
93
0
0
0
0
17
Multiband
0.6
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
91
95
40
0.8
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
74
59
23
1.0
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
99
100
100
100
99
56
10
14
1.2
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
99
100
99
100
70
15
3
4
1.4
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
99
100
99
91
9
1
0
4
1.6
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
95
99
93
43
0
0
0
0
1.8
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
99
99
93
90
0
0
0
0
0
2.0
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
100
98
85
84
85
0
0
0
0
0Six orders of PO are tested in Table1. Due to 5-level wavelet decomposition, there are five detail SBs to be processed. Each SB corresponds to a ΨTΨ and different PO corresponds to different dimension of ΨTΨ. So, there are 30 combinations shown in Table 1 header where each PO contains five SBs. In x-axis direction, Table 1 reflects how PO influences matrix full rank probabilities. In y-axis direction of Table 1, four signals are tested with each corrupted by different noise variance σ. This reflects how noise level influences matrix full rank probabilities. For example, the italic number shown in Table 1 means matrix full rank probability is 100%. This matrix is determined by PO2 and “blocks” signal’s SB3 in σ=1.6 noise level. All probabilities are obtained by counting matrix full rank times in a 100-trail.From Table1, PO1–PO3 always ensure 100% matrix full rank for four signals in all tested noise level. Saying 100% matrix full rank for a signal we mean 100% matrix full rank for every subband of this signal. PO4 only cannot ensure “blocks” due to almost 0% probability for SB5 (bold underlined numbers, bold numbers indicate the probabilities under 100%). For PO5 and PO6, there are more and more bold numbers than PO4. In horizontal direction, the trend is higher PO lower probability for satisfying matrix full rank. In vertical direction, higher noise level tends to lead lower probability for full rank, but PO is a dominant factor. From statistical analysis, we restrict our option from PO1 to PO3.Visual comparisons of PO1 to PO3 in supervised denoising are shown in Figure5. As it can be seen, in both MSE and SNR sense PO1 is inferior to PO2 and PO3 in all four signals. PO2 is inferior to PO3 in blocks and almost equals PO3 in the three remaining signals. In unsupervised denoising, similar conclusions can be made and are not repetitive state. In the end, we choose PO2 for final denoising simulation.Figure 5
PO1 to PO3 comparisons in supervised denoising condition.


## 3.3. Denoising Performance Comparison

As mentioned above, comparisons are made among supervised and unsupervised denoising methods and VisuShrink, SureShrink, MiniMaxShrink, and BayesShrink. In both MSE and SNR sense, denoising performance comparisons in global and local views with different noise varianceσ are presented in Figure 6. σ changes in range of 0.5–2.0 with step 0.1 and each data point plotted in Figure 6 is a 50-average result. A remarkable notice is that we adopt a clean signal corrupted by noise in σ=0.5 level as desired signal for supervised denoising. Moreover, MSEs in Figure 6 are plotted in log scale and SNRs in Figure 6 are plotted in linear scale.Figure 6
Comparison of both MSE and SNR changing with noise varianceσ for “blocks,” “quadchirp,” “multitone,” and “multiband” denoising.From global and local views of Figure6, supervised denoising method achieves better MSE and SNR quality than unsupervised denoising method and BayesShrink for all four noisy signals. All these three methods have better MSE and SNR performance than VisuShrink, SureShrink, and MiniMaxShrink. In noisy “blocks” processing, BayesShrink is slightly better than unsupervised denoising method when σ is greater than 0.9 while poorer when σ is less than 0.9 in both MSE and SNR sense. However, unsupervised denoising method has better MSE and SNR performance than BayesShrink in noisy quadchirp, multitone, and multiband processing. Thus, a conclusion can be made that the proposed supervised and unsupervised denoising methods show their advantages.Our approaches are also validated for real-world signals. We test supervised and unsupervised denoising approaches by an AGWN corrupted real audio signal. In both MSE and SNR sense, comparisons with changingσ in range of 0.1–0.5 are shown in Figure 7 and each data point plotted is still a 50-average result. We adopt clean audio signal corrupted by noise in σ=0.5 level as desired signal for supervised denoising. In addition, MSEs are plotted in log scale and SNRs are plotted in linear scale. From global and local views of Figure 7, the proposed supervised and unsupervised denoising methods have smaller MSEs and higher SNRs than VisuShrink, SureShrink, MiniMaxShrink, and BayesShrink. A conclusion can still be made that the proposed supervised and unsupervised denoising methods show their advantages in real-world condition.Figure 7
Denoising performance comparisons in both MSE and SNR sense with changing noise varianceσ for real audio signal denoising.Remark. All denoising methods except for supervised denoising involved in Figures 6 and 7 require noise variance σ, and estimates of it are achieved by MAD for all those methods. So, proposed unsupervised method and the other four techniques work under the same condition, which reflects superiority of the proposed method. In practice, MAD with simple implementation and estimate accuracy do satisfy our requirements.


## 3.4. Noise Effect Analysis

In this section, different noise levels of desired signal for supervised denoising are analyzed. It only aims at real audio signal. In simulation, desired signal noise levels are set in rangeσ=0.1 to σ=2.5 with equal step 0.4. Both MSE and SNR results are shown in Figure 8 and each data point plotted in Figure 8 is averaged by 50 times. It must be noticed that σ on x-axis represents noise level of noisy signal to be processed. All desired signal noise levels are compared with BayesShrink denoising mothed. From Figure 8 it can be seen that all σ=0.1 to σ=2.5 desired signal noise levels have better denoising performance than BayesShrink in MSE and SNR sense when σ≥0.18 on x-axis. While when σ<0.18 on x-axis BayesShrink does better than σ=2.1 and σ=2.5 and worse than σ=0.1 to σ=1.7 in MSE and SNR sense. In addition, denoising performance becomes worse when desired signal noise level gets higher. In conclusion, if some noisy desired signal (even higher noise level than the signal to be processed) is available at hand, it is still possible to achieve good denoising performance.Figure 8
Different desired signal noise level and BayesShrink denoising performance comparisons in both MSE and SNR sense with changing noise varianceσ for real audio signal denoising.


## 4. Conclusions

In order to utilize the effectiveness of wavelet domain denoising, two subband adaptive denoising schemes were proposed in this paper. For any subband, a unified linear matrix form with flexible structure of threshold function was proposed and OTRP function was proposed for a concrete realization. Based on the unified linear flexible structure threshold function, both supervised and unsupervised subband adaptive denoising frameworks were established. A direct MSE minimization was conducted in supervised denoising while Stein’s unbiased risk estimate as MSE estimate was minimized in unsupervised denoising for flexible coefficients determination. Conditions to obtain optimal coefficients were further discussed.Applying the concrete OTRP function, simulations for polynomial order, denoising performance, and noise effect were conducted. A full rank probability statistical table was generated in polynomial order simulation. The table reflects that PO1–PO3 always ensure full rank ofΨTΨ. By choosing PO2, denoising performance was compared among proposed supervised and unsupervised denoising methods and VisuShrink, SureShrink, MiniMaxShrink, and BayesShrink. Comparison results demonstrate that the proposed methods have advantages over those soft and hard threshold function based denoising schemes in both simulation and real-world condition. Final simulation of noise effect analysis for supervised denoising demonstrates that supervised denoising method is robust for different desired signal noise level.



---

*Source: 5203214-2017-03-09.xml*