MATLAB中的小波阈值去噪

MATLAB中实现了信号的阈值去噪，主要包括阈值去噪和阈值获取两方面。

# b; ~9 {/ e5 v* i) m' F) K2 h1.阈值获取
' G4 S+ g$ e( q7 cMATLAB中实现阈值获取的函数有ddencmp、thselect、wbmpen和wwdcbm，下面对它们的用法进行简单的说明。

ddencmp的调用格式有以下三种：
: j( M( ~4 D: U（1）[THR,SORH,KEEPAPP,CRIT]=ddencmp(IN1,IN2,X)
  y3 m& D  Z5 u) g4 R& ?0 t（2）[THR,SORH,KEEPAPP,CRIT]=ddencmp(IN1,'wp',X)
& l7 M6 z  F% W* ?& j0 X! P2 B（3）[THR,SORH,KEEPAPP,CRIT]=ddencmp(IN1,'wv',X)
1 ~, y7 h9 q# R: k9 {函数ddencmp用于获取信号在消噪或压缩过程中的默认阈值。输入参数X为一维或二维信号；IN1取值为'den'或'cmp'，'den'表示进行去噪，'cmp'表示进行压缩；IN2取值为'wv'或'wp'，wv表示选择小波，wp表示选择小波包。返回值THR是返回的阈值；SORH是软阈值或硬阈值选择参数；KEEPAPP表示保存低频信号；CRIT是熵名（只在选择小波包时使用）。

函数thselect的调用格式如下：/ D6 M& m+ Z, K) W. M
THR=thselect(X,TPTR);
; E1 \* M. J. xTHR=thselect(X,TPTR)根据字符串TPTR定义的阈值选择规则来选择信号X的自适应阈值。
4 t$ y) }$ `: i4 O% \$ Z2 K自适应阈值的选择规则包括以下四种：
( P# ^( s  G, |*TPTR='rigrsure'，自适应阈值选择使用Stein的无偏风险估计原理。6 s  V4 m2 |  U( K/ o. S2 `
*TPTR='heursure'，使用启发式阈值选择。  N4 ]3 p/ j' F2 ?- o$ J
*TPTR='sqtwolog'，阈值等于sqrt(2*log(length(X))).
3 V) \" W5 V7 D0 {$ ]( Q7 \*TPTR='minimaxi'，用极大极小原理选择阈值。& z/ j, V8 [5 ^) c7 ]
阈值选择规则基于模型 y = f(t) + e，e是高斯白噪声N(0,1)。

函数wbmpen的调用格式如下：
0 S1 _9 a. l4 D' ?: R, V) j' ETHR=wbmpen(C,L,SIGMA,ALPHA);( T3 x* H. I5 b: |$ J7 k
THR=wbmpen(C,L,SIGMA,ALPHA)返回去噪的全局阈值THR。THR通过给定的一种小波系数选择规则计算得到，小波系数选择规则使用Birge-Massart的处罚算法。{C,L]是进行去噪的信号或图像的小波分解结构；SIGMA是零均值的高斯白噪声的标准偏差；ALPHA是用于处罚的调整参数，它必须是一个大于1的实数，一般取ALPHA=2。
: ^* N7 o' J- J: ]设t*使crit(t)=-sum(c(k)^2,k<=t) + 2 * SIGMA^2 * t*(ALPHA+log(n/t))的最小值，其中c(k)是按绝对值从大到小排列的小波包系数，n是系数的个数，则THR=|c(t*)|。
: o- r/ @7 n% X5 _! ~wbmpen(C,L,SIGMA,ALPHA,ARG)计算阈值并画出三条曲线。
/ g) d, t0 i, L0 n0 Y- C* {& ~2 * SIGMA^2 * t*(ALPHA+log(n/t))
. i" L4 P- E9 e( tsum(c(k)^2, k<=t)
* F0 u' l. T9 d- k9 ~0 R; r: G" kcrit(t)

wdcbm的调用格式有以下两种：
; I# ?8 l9 g6 e, l, _& l（1）[THR,NKEEP]=wdcbm(C,L,ALPHA);
4 a  \5 S1 ?% }% I（2）[THR,NKEEP]=wdcbm(C,L,ALPHA,M);
3 |2 Y4 W5 j, L2 S0 G函数wdcbm是使用Birge-Massart算法获取一维小波变换的阈值。返回值THR是与尺度无关的阈值，NKEEP是系数的个数。[C,L]是要进行压缩或消噪的信号在j=length(L)-2层的分解结构；LAPHA和M必须是大于1的实数；THR是关于j的向量，THR（i)是第i层的阈值；NKEEP也是关于j的向量，NKEEP(i)是第i层的系数个数。一般压缩时ALPHA取1.5，去噪时ALPHA取3.

9 @8 ~& o# W& c! d# e. ?# [
2.信号的阈值去噪
1 \0 n+ E) y, q" ]MATLAB中实现信号的阈值去噪的函数有wden、wdencmp、wthresh、wthcoef、wpthcoef以及wpdencmp。下面对它们的用法作简单的介绍。

函数wden的调用格式有以下两种：
& Y: F. T, F6 g" T/ M& o; o（1）[XD,CXD,LXD]=wden(X,TPTR,SORH,SCAL,N,'wname')
0 {8 I, \% x7 G) N9 T（2）[XD,CXD,LXD]=wden(C,L,TPTR,SORH,SCAL,N,'wname')
% Z( c4 ^! v6 o3 R9 m0 V2 E函数wden用于一维信号的自动消噪。X为原始信号，[C,L]为信号的小波分解，N为小波分解的层数。
1 E8 @: L- L, {" T! S0 P" g+ aTHR为阈值选择规则：2 P/ O, D- i, R! D' H! d
*TPTR='rigrsure'，自适应阈值选择使用Stein的无偏风险估计原理。
0 ~4 k2 n. o# l7 K' x  z" S*TPTR='heursure'，使用启发式阈值选择。, [+ }/ s9 R* g$ y. _  Z, Y7 N3 V
*TPTR='sqtwolog'，阈值等于sqrt(2*log(length(X))).
+ {/ A$ {6 B' ~& T/ I*TPTR='minimaxi'，用极大极小原理选择阈值。
2 N; a$ v# ?/ A4 oSORH是软阈值或硬阈值的选择（分别对应's'和'h'）。  ^5 a4 m+ v) u3 O6 V& |* d
SCAL指所使用的阈值是否需要重新调整，包含下面三种：( J3 `2 D# `5 d+ L
*SCAL='one' 不调整；
1 Z! U: k) d5 h* v*SCAL='sln' 根据第一层的系数进行噪声层的估计来调整阈值。
% I1 j1 S) _! X4 t. K% Z*SCAL='mln' 根据不同的噪声估计来调整阈值。) N! j1 R  e  Q" E0 ^/ u" e
XD为消噪后的信号，[CXD,LXD]为消噪后信号的小波分解结构。格式（1）返回对信号X经过N层分解后的小波系数进行阈值处理后的消噪信号XD和信号XD的小波分解结构[CXD,LXD]。格式（2）返回参数与格式（1）相同，但其结构是由直接对信号的小波分解结构[C,L]进行阈值处理得到的。

函数wdencmp的调用格式有以下三种：
$ d( e5 |: [$ l2 |, `6 N(1)[XC,CXC,LXC,PERF0,PERFL2]=wdencmp('gbl',X,'wname',N,THTR,SORH,KEEPAPP);4 w: }9 l, l7 W# g% c! E
(2)[XC,CXC,LXC,PERF0,PERFL2]=wdencmp('lvd',X,'wname',N,THTR,SORH);) z0 q4 u4 {9 P5 z
(3)[XC,CXC,LXC,PERF0,PERFL2]=wdencmp('lvd',C,L,'wname',N,THTR,SORH);3 ?$ M3 |* ?4 ?3 |- o, W
函数wdencmp用于一维或二维信号的消噪或压缩。wname是所用的小波函数，gbl(global的缩写)表示每一层都采用同一个阈值进行处理，lvd表示每层采用不同的阈值进行处理，N表示小波分解的层数，THR为阈值向量，对于格式（2）和（3）每层都要求有一个阈值，因此阈值向量THR的长度为N，SORH表示选择软阈值或硬阈值（分别取值为's'和'h'），参数KEEPAPP取值为1时，则低频系数不进行阈值量化，反之，低频系数要进行阈值量化。XC是要进行消噪或压缩的信号，[CXC,LXC]是XC的小波分解结构，PERF0和PERFL2是恢复或压缩L^2的范数百分比。如果[C,L]是X的小波分解结构，则PERFL2=100*(CXC向量的范数/C向量的范数)^2；如果X是一维信号，小波wname是一个正交小波，则PERFL2=100||XC||^2/||X||^2。

函数wthresh的调用格式如下：
; I" f+ ?$ F) A: D, oY=wthresh(X,SORH,T)! ^- T  k1 W  ^; t2 I" {5 ~
Y=wthresh(X,SORH,T) 返回输入向量或矩阵X经过软阈值（如果SORH='s'）或硬阈值（如果SORH='h'）处理后的信号。T是阈值。8 j6 n: v/ L/ F6 O) {# q7 ^9 C
Y=wthresh(X,'s',T)返回的是Y=SIG(X)*(|X|-T)+，即把信号的绝对值与阈值进行比较，小于或等于阈值的点变为零，大于阈值的点为该点值与阈值的差值。( l% `3 G3 ~1 ^5 y! t
Y=wthresh(X,'h',T)返回的是Y=X*1(|X|>T)，即把信号的绝对值和阈值进行比较，小于或等于阈值的点变为零，大于阈值的点保持不变。一般来说，用硬阈值处理后的信号比用软阈值处理后的信号更粗糙。

函数wthcoef的调用格式下面四种：1 q7 X6 Z' m1 {' V- j! S
（1）NC=wthcoef('d',C,L,N,P); T7 G- s* e4 a; M
（2）NC=wthcoef('d',C,L,N)
3 l- a- |; n: `3 ?4 D1 g  E2 Y（3）NC=wthcoef('a',C,L)" S& a, q2 w- Y- w2 V. a/ C( B" G  o
（4）NC=wthcoef('t',C,L,N,T,SORH)7 m$ p8 R  Q1 g5 k" O6 {/ t
函数wthcoef用于一维信号小波系数的阈值处理。: X0 r! R* T. B, K& m
格式（1）返回小波分解结构[C,L]经向量N和P定义的压缩率处理后的新的小波分解向量NC，[NC,L]构成一个新的小波分解结构。N包含被压缩的细节向量，P是把较小系数置0的百分比信息的向量。N和P的长度必须相同，向量N必须满足1<=N(i)<=length(L)-2。; F$ M$ J6 M- y5 K- o
格式（2）返回小波分解结构[C,L]经过向量N中指定的细节系数置0后的小波分解向量NC。& C0 c6 N: d$ F; P! e: i
格式（3）返回小波分解结构[C,L]经过近似系数置0后的小波分解向量NC。: D5 ^0 j' S2 @$ T9 R# B! M/ j6 s; y
格式（4）返回小波分解结构[C,L]经过将向量N作阈值处理后的小波分解向量NC。如果SORH=’s‘，则为软阈值；如果SORH='h'则为硬阈值。N包含细节的尺度向量，T是N相对应的阈值向量。N和T的长度必须相等。

函数wpdencmp的调用格式有以下两种：
0 _& Z; }2 d. @! V8 v8 X（1）[XD,TREED,PERF0,PERFL2]=wpdencmp(X,SORH,N,'wname',CRIT,PAR,KEEPAPP)1 d6 `* ^' I3 o
（2）[XD,TREED,PERF0,PERFL2]=wpdencmp(TREE,SORH,CRIT,PAR,KEEPAPP)8 b1 N9 O: l, `  `7 f9 O- P" U" m% z* Q
函数wpdencmp用于使用小波包变换进行信号的压缩或去噪。) S7 e" ?2 }: w7 j! S( C" ^4 S, I
格式（1）返回输入信号X（一维或二维）的去噪或压缩后的信号XD。输出参数TREED是XD的最佳小波包分解树；PERFL2和PERF0是恢复和压缩L2的能量百分比。PERFL2=100*(X的小波包系数范数/X的小波包系数)^2；如果X是一维信号，小波wname是一个正交小波，则PERFL2=100*||XD||^2/||X||^2。SORH的取值为's'或'h'，表示的是软阈值或硬阈值。
+ U- x$ J' m7 u& ~# X' w+ J输入参数N是小波包的分解层数，wname是包含小波名的字符串。函数使用由字符串CRIT定义的熵和阈值参数PAR实现最佳分解。如果KEEPAPP=1，则近似信号的小波系数不进行阈值量化；否则，进行阈值量化。
/ a. D' U" j" W6 _格式（2）与格式（1）的输出参数相同，输入选项也相同，只是它从信号的小波包分解树TREE进行去噪或压缩。
7 j5 j4 t) d% [  `0 S& ]" s