Wavelet threshold denoising based on MATLAB

Thank you so much:
http://blog.sina.com.cn/u/1861445474
Http://blog.chinaaet.com/detail/3083.html

The threshold denoising of signal is realized in MATLAB, which mainly includes two aspects: threshold denoising and threshold acquisition .

1. Threshold Acquisition

The functions implemented in MATLAB are ddencmp, Thselect, Wbmpen, and WWDCBM, and their usage is briefly explained below.

The DDENCMP format has the following three types of calls:
(1) [THR,SORH,KEEPAPP,CRIT]=DDENCMP (In1,in2,x)
(2) [thr,sorh,keepapp,crit]=ddencmp (IN1, ' WP ', X)
(3) [thr,sorh,keepapp,crit]=ddencmp (IN1, ' WV ', X)

The function ddencmp is used to obtain the default thresholds for the signal during noise cancellation or compression. The input parameter x is a one-dimensional or two-dimensional signal; The IN1 value is ' den ' or ' cmp ', ' den ' means de-noising, ' CMP ' means compression, ' WV ' or ' WP ', WV means select wavelet, WP means SELECT wavelet packet. The return value THR is the threshold value returned, SORH is the soft or hard threshold selection parameter, Keepapp is the low frequency signal, Crit is the entropy name (only used when selecting the Wavelet packet).
The calling format for the function Thselect is as follows:
Thr=thselect (X,TPTR);
Thr=thselect (X,TPTR) selects the adaptive threshold for signal X based on the threshold selection rules defined by the string tptr.

The selection rules for adaptive thresholds include the following four types:
tptr= ' rigrsure ', adaptive threshold selection uses Stein's unbiased risk estimation principle.
tptr= ' heursure ', using heuristic threshold selection.
tptr= ' Sqtwolog ', the threshold is equal to sqrt (2*log (Length (X))).
tptr= ' Minimaxi ', select the threshold value with the minimax principle.
The threshold selection rule is based on model y = f (t) + e,e is Gaussian white noise n (0,1).

The calling format for the function Wbmpen is as follows:
Thr=wbmpen (C,l,sigma,alpha);
Thr=wbmpen (C,l,sigma,alpha) returns the global threshold THR for de-noising. THR by a given wavelet coefficients selection rule calculation, the wavelet coefficients selection rule uses birge-massart penalty algorithm. {C,l] is a wavelet decomposition structure for denoising signals or images; Sigma is the standard deviation of the Gaussian white noise of the 0 mean; Alpha is the adjustment parameter used for punishment, it must be a real number greater than 1, generally take alpha=2.

Set t* so that Crit (t) =-sum (c (k) ^2,k<=t) + 2 * sigma^2 * t* (Alpha+log (n/t)) The minimum value, where C (k) is the absolute size of the large-to-small arrangement of the wavelet packet coefficients, n is the number of coefficients, then Thr=|c (t*) |.
Wbmpen (C,l,sigma,alpha,arg) calculates the threshold value and draws three curves.

2 * sigma^2 * t* (Alpha+log (n/t))
Sum (c (k) ^2, k<=t)
Crit (t)
The WDCBM format has the following two types of calls:
(1) [THR,NKEEP]=WDCBM (C,l,alpha);
(2) [THR,NKEEP]=WDCBM (C,L,ALPHA,M);
The function wdcbm is a threshold value that uses the Birge-massart algorithm to obtain a dimensional wavelet transformation. The return value THR is a scale-independent threshold, and Nkeep is the number of coefficients. [C,l] is a signal to be compressed or de-noising in the decomposition structure of the j=length (L)-2 layer, Lapha and M must be a real number greater than 1, THR is a vector about J, THR (i) is the threshold of layer I; Nkeep is also a vector of J, Nkeep (i) is the number of coefficients of layer I. In general compression, Alpha takes 1.5, and Alpha takes 3 when de-noising.

2. Threshold de-noising of signals
The functions of threshold denoising in MATLAB are Wden, wdencmp, Wthresh, Wthcoef, Wpthcoef and wpdencmp. The following is a brief introduction to their usage.

The function Wden has the following two types of invocation formats:
(1) [Xd,cxd,lxd]=wden (X,tptr,sorh,scal,n, ' wname ')
(2) [Xd,cxd,lxd]=wden (C,l,tptr,sorh,scal,n, ' wname ')

The function Wden is used for automatic noise cancellation of one-dimensional signals. X is the original signal, [c,l] is the wavelet decomposition of the signal, n is the number of layers of wavelet decomposition.
THR for threshold selection rules:
*tptr= ' rigrsure ', adaptive threshold selection uses Stein's unbiased risk estimation principle.
*tptr= ' heursure ', using heuristic threshold selection.
*tptr= ' Sqtwolog ', the threshold is equal to sqrt (2*log (Length (X))).
*tptr= ' Minimaxi ', select the threshold value with the minimax principle.
The SORH is the choice of soft or hard thresholds (corresponding to ' s ' and ' H ' respectively).
Scal refers to whether the threshold you are using needs to be re-adjusted to include the following three types of values:
*scal= ' One ' does not adjust;
*scal= ' SLN adjusts the threshold by estimating the noise layer based on the coefficients of the first layer.
*scal= ' mln ' adjusts thresholds based on different noise estimates.
XD is the signal after denoising, [Cxd,lxd] is the wavelet decomposition structure of the signal after denoising.

The format (1) returns the wavelet decomposition structure of the de-noising signal xd and the signal xd after the signal x is processed by the wavelet coefficients of n-layer decomposition [cxd,lxd]. Format (2) The return parameter is the same as the format (1), but its structure is obtained by the threshold processing of the signal's wavelet decomposition structure [c,l] directly.

The function wdencmp has the following three types of invocation formats:
(1) [XC,CXC,LXC,PERF0,PERFL2]=WDENCMP (' GBL ', X, ' Wname ', N,thtr,sorh,keepapp);
(2) [XC,CXC,LXC,PERF0,PERFL2]=WDENCMP (' Lvd ', X, ' Wname ', N,THTR,SORH);
(3) [XC,CXC,LXC,PERF0,PERFL2]=WDENCMP (' Lvd ', c,l, ' Wname ', N,THTR,SORH);

The function wdencmp is used for denoising or compressing one-dimensional or two-dimensional signals. Wname is the wavelet function used, GBL (global abbreviation) means that each layer is treated with the same threshold, LVD means that each layer is processed with different thresholds, n represents the number of layers of the wavelets decomposition, the thr is a threshold vector, for the format (2) and (3) each layer requires a threshold value, Therefore, the length of the threshold vector thr is n,sorh to select the soft or hard threshold value (respectively, the value is ' s ' and ' H '), the parameter Keepapp value is 1 o'clock, the low frequency coefficient does not carry out threshold quantization, conversely, the low frequency coefficient to be quantified threshold. XC is a signal to be de-noising or compressed, [CXC,LXC] is the wavelet decomposition structure of XC, PERF0 and PERFL2 are the norm percentages for recovering or compressing l^2. If [c,l] is the wavelet decomposition structure of x, then perfl2=100* (norm of the CXC vector/norm of the C vector) ^2; if X is a one-dimensional signal, the wavelet wname is an orthogonal wavelet, then perfl2=100| | xc| | ^2/| | x| | ^2.

The calling format for the function Wthresh is as follows:
Y=wthresh (x,sorh,t)
Y=wthresh (x,sorh,t) returns the input vector or matrix X after a soft threshold (if sorh= ' s ') or a hard threshold value (if sorh= ' h ') is processed by the signal. T is the threshold value.
Y=wthresh (x, ' s ', T) returns the Y=sig (x) * (| X|-T) +, which compares the absolute value of the signal to the threshold, the point less than or equal to the threshold becomes 0, and the point greater than the threshold value is the difference between the point value and the threshold value.
Y=wthresh (X, ' h ', T) returns the Y=x*1 (| X|>T), which compares the absolute value of the signal to the threshold, the point less than or equal to the threshold becomes 0, and the point greater than the threshold remains unchanged. In general, the signal treated with a hard threshold is more coarse than the signal processed with a soft threshold value.

There are four types of call formats for function Wthcoef:
(1) Nc=wthcoef (' d ', c,l,n,p)
(2) Nc=wthcoef (' d ', c,l,n)
(3) Nc=wthcoef (' A ', c,l)
(4) Nc=wthcoef (' t ', C,L,N,T,SORH)
The function Wthcoef is used for threshold processing of one-dimensional signal wavelet coefficients.

Format (1) returns the wavelet decomposition structure [c,l] The new wavelet decomposition vector, which is processed by the compression rate defined by the vector N and P, nc,[nc,l] constitutes a new wavelet decomposition structure. n contains the compressed detail vector, and p is a vector that puts the percentage information of the smaller factor at 0. The length of n and P must be the same, vector n must meet 1<=n (i) <=length (L)-2.
Format (2) returns the wavelet decomposition structure [c,l] after the detail coefficient specified in vector n is set 0 after the wavelet decomposition vector nc.
Format (3) returns the wavelet decomposition structure [c,l] after the approximate coefficient of 0 after the wavelet decomposition vector nc.
Format (4) returns the wavelet decomposition structure [c,l] after the vector n is processed by the threshold value of the wavelet decomposition vector nc. If sorh= ' s ', then the soft threshold value, or if sorh= ' H ' is a hard threshold. n the scale vector containing the details, and T is the threshold vector corresponding to the n phase. The lengths of N and T must be equal.

The function wpdencmp has the following two types of invocation formats:
(1) [XD,TREED,PERF0,PERFL2]=WPDENCMP (X,sorh,n, ' wname ', Crit,par,keepapp)
(2) [XD,TREED,PERF0,PERFL2]=WPDENCMP (Tree,sorh,crit,par,keepapp)
The function wpdencmp is used to compress or de-noising the signal using Wavelet packet transform.

Format (1) returns the input signal x (one-dimensional or two-dimensional) of the de-noising or compressed signal xd. The output parameter treed is the best wavelet packet Decomposition Tree for XD, and PERFL2 and PERF0 are the percentage of energy recovered and compressed L2. perfl2=100* (wavelet packet coefficients of x's wavelet packet coefficient norm/x) ^2; if X is a one-dimensional signal, the wavelet wname is an orthogonal wavelet, then perfl2=100*| | xd| | ^2/| | x| | ^2. The value of the Sorh is ' s ' or ' h ', which indicates a soft or hard threshold value.

The input parameter n is the decomposition layer of the wavelet packet, and Wname is a string containing the wavelet name. The function uses the entropy and threshold parameter par defined by the string crit to achieve optimal decomposition. If keepapp=1, the wavelet coefficients of the approximate signal are not quantified by the threshold value, otherwise, the threshold is quantified.
The format (2) is the same as the output parameter of format (1), and the input option is the same, except that it is de-noising or compressed from the wavelet packet Decomposition Tree of the signal.