基于MATLAB小波去噪方法的研究

基金项目:中科院优秀青年学者奖基金项目收稿日期:2002-06-17　　第20卷　第7期计　算　机　仿　真2003年7月　　文章编号:1006-9348(2003)07-0119-04基于MATLAB小波去噪方法的研究杜浩藩,丛爽(中国科学技术大学自动化系,安徽合肥230027)摘要:在实际的计算机控制系统中,采样信号不可避免的受到各种噪声和干扰的污染,使得由辨识采样信号得到的系统模型存在偏差而妨碍了系统控制精度的提高。Donoho和Johnstone提出的小波收缩去噪算法对去除叠加性高斯白噪声非常有效。我们对此在MATLAB环境下做了详尽的仿真研究,并利用滤波器、傅立叶变换/反变换和小波收缩三种方法对实际电机随动系统的辨识采样数据进行了去噪对比实验,总结了应用小波去噪的一些实际经验。关键词:小波变化;滤波;去噪中图分类号:TP301.6　　文献标识码:A1　引言在计算机控制中,需要从控制系统中获得采样信号(观测数据)用于系统的模型辨识或实时控制,而这些采样信号将不可避免的受到噪声的污染。如何从这些受噪声干扰的信号中估计得到“纯净”的信号是建立系统高精度模型和实现高性能控制的关键。滤波器去噪是实际应用最广泛的一种方法,但时常在滤除噪声的同时导致了有用信号的失真,它是从纯频域的角度来分析应该消除哪些频率范围内的噪声。1995年Donoho和Johnstone[1]、[2]提出了小波收缩去噪的技术,他们研究的是在叠加性高斯白噪声环境下检测出真实信号的情况,利用正交小波变换和高斯随机变量的性质对信号的小波分解系数做阈值量化(thresholding),无失真的还原出真实信号。文章[3]中对Donoho-Johnstone的去噪方法做了总结推广,研究了在高斯白噪声情况下选择不同的小波变换形式(如非正交小波变换)的去噪效果,并公式化了实际数据中的几种更复杂的噪声模型。本文对Donoho-Johnstone的收缩去噪方法在MATLAB环境下做了较为详尽的研究,总结了应用小波收缩去噪的一些实际经验。文章组织如下:第一节是小波去噪原理;第二节是小波去噪过程中两种门限阈值处理方式和几种阈值选取方法;第三节是MATLAB中的小波去噪对比实验研究;最后是小结。2　小波去噪原理叠加性高斯白噪声是最常见的噪声模型,受到叠加性高斯白噪声“污染”的观测信号可以表示为:di=fi+ε·zi,i=1,…,N(▲)其中di为含噪信号,fi为“纯净”采样信号,zi为独立同分布的高斯白噪声N(0,1)(用符号可表示为:zi～iidN(0,1)),ε为噪声水平,信号长度为N。为了从含噪信号di中还原出真实信号fi,可以利用信号和噪声在小波变换下的不同的特性,通过对小波分解系数进行处理来达到信号和噪声分离的目的。在实际工程应用中,有用信号通常表现为低频信号或是一些比较平稳的信号,而噪声信号则通常表现为高频信号,所以我们可以先对含噪信号进行小波分解(如进行三层分解):S=cA1+cD1=cA2+cD2+cD1=cA3+cD3+cD2+cD1(1)其中cAi为分解的近似部分,cDi为分解的细节部分,i=1,2,3,则噪声部分通常包含在cD1,cD2,cD3中,用门限阈值对小波系数进行处理,重构信号即可达到去噪的目的[4]。总结去噪过程,可以分成以下三个步骤:1)对观测数据作小波分解变化[3]:W0d=W0f+ε·W0z(2)其中d表示观测数据向量f1,f2,…fN,f是真实信号向量f1,f2,…,fN,z是高斯随机向量z1,z2,…,zN,其中用到了小波分解变换是线性变换的性质(小波分解变换本质上是一种积分变换)。2)对小波系数W0d作门限阈值处理(根据具体情况可以使用软阈值处理或硬阈值处理,而且可以选择不同的阈值形式,这将在后面作详细讨论),比如选取最著名的阈值形式:tN=ε2logN(3)　　门限阈值处理可以表示为ηtNW0d,可以证明当N趋于无穷大时使用阈值形式(3)对小波系数作软阈值处理可以几乎完全去除观测数据中的噪声[5]。3)对处理过的小波系数作逆变换W-10重构信号:—119—f＊=W-10ηtNW0d(4)即可得到受污染采样信号去噪后的信号。3　阈值的选取和阈值量化Donoho-Johnstone的小波收缩去噪方法的关键步骤是如何选择阈值和如何进行门限阈值处理,在这我们将作较为详细的讨论。3.1　软阈值和硬阈值在对小波系数作门限阈值处理操作时,可以使用软阈值处理方法或硬阈值处理方法,硬阈值处理只保留较大的小波系数并将较小的小波系数置零:ηH(w,t)=w,|w|≥t0,|w|t(5)　　软阈值处理将较小的小波系数置零但对较大的小波系数向零作了收缩:ηS(w,t)=w-t,w≥t0,|w|tw+t,w≤t(6)　　直观形式见图1(图中取t=0.4)。从图上我们可以看出软阈值处理是一种更为平滑的形式,在去噪后能产生更为光滑的结果,而硬阈值处理能够更多的保留真实信号中的尖峰等特征。软阈值处理实质上是对小波分解系数作了收缩,从而Donoho-Johnstone将这种去噪技术称之为小波收缩(WaveletShrinkage)。图1　硬阈值和软阈值3.2　阈值的几种形式阈值的选取有多种形式,选取规则都是基于含噪信号模型(▲)中信号水平为1的情况,对于噪声水平未知或非白噪声的情况可以在去噪时重新调整得到的阈值。在Matlab中有4种阈值函数形式可以选用:1)sqtwolog:采用固定的阈值形式,如式(3),Donoho-Johnstone称之为VisuShrink,因为这种阈值形式在软门限阈值处理中能够得到直观意义上很好的去噪效果。2)minimaxi:采用极大极小原理选择的阈值,和sqtwolog一样也是一种固定的阈值,它产生一个最小均方误差的极值,计算公式为:t=0,N≤320.3936+0.1829＊log2N(7)　　3)rigrsure:采用史坦(Stein)的无偏似然估计原理进行阈值选择,首先得到一个给定阈值t的风险估计,选择风险最小的阈值t＊作为最终选择。4)heursure:选择启发式阈值,它是sqtwolog和rigrsure的综合,当信噪比很小时rigrsure估计有很大的噪声,这时heur-sure采用固定阈值sqtwolog。4　MATLAB中的小波去噪MATLAB中的小波工具包提供了全面的小波变化及其应用的各种功能,其中小波去噪方面实现了Donoho-Johnstone等的去噪算法,而且可以选择使用图形界面操作工具或者去噪函数集合两种形式,图形界面操作工具直观易用,而利用函数集合可以实现更灵活强大的功能。我们利用小波去噪函数集合在MATLAB中作了一系列实验,充分体会到了小波去噪的强大功能。4.1　小波去噪函数集合下面是几个最为常用的小波去噪函数:1)x=wnoise(fun,n):产生Donoho-Johnstone设计的6种用于测试小波去噪效果的典型测试数据,函数根据输入参数fun的值输出名为`blocks',`bumps',`heavysine”,`doppler',`quadchirp'或`mishmash'的6种函数数据,数据长度为2N。这6种测试数据在验证算法和仿真实验时非常有用。2)[xd,cxd,lxd]=wden(x,tptr,sorh,scal,level,`wname':最主要的一维小波去噪函数。其中输入参数x为输入需要去噪的信号,tptr为2.2节中的4种阈值形式,sorh设定为`s'或`h'表示用软门限阈值或硬门限阈值处理。2.2节中说过4种阈值形式是基于信号水平为1的高斯白噪声模型推导得到的,当噪声不是白噪声时,必须在小波分解的不同层次估计噪声水平,scal=`one'不进行重新估计,scal=`sln'只根据第一层小波分解系数估计噪声水平,scal=`mln'在每个不同的小波分解层次估计噪声水平,根据scal参数的设定,wden()函数决定最终应用于每一个小波分解层次的阈值数值。最后两个参数level和wname表示利用名为wname的小波对信号作level层分解。输出去噪后的数据xd及xd的附加小波分解结构[cxd,lxd]。还有功能更强大的用于一维或二维小波去噪或压缩的函数wdencmp()。3)thr=thselect(x,tptr):去噪阈值选择函数。4)y=wthresh(x,sorh,t):对信号x作阈值为t的门限阈值处理。4.2　小波去噪对比实验我们对小波去噪时各种参数设置作了详尽的对比研究,实验信号是由wnoise()函数产生的含标准的高斯白噪声信噪比为3的`heavysine'信号,用wden()函数进行去噪处理,对—120—在各种可选参数下的去噪效果作了对比研究,得出了以下结论。1)四种阈值选取方式的对比(tptr的设置)在MATLAB下运行如下程序段:tptr=[`rigrsure';`heursure';`sqtwolog';`minimaxi'];forI=1:4,xd=wden(x,TPTR(I,:),`s',`sln',5,`sym8');subplot(2,2,I);plot(xd);title(TPTR(I,:));end;得到四个去噪效果图(见图2),可以看出固定阈值形式(sqtwolog)和启发式阈值形式(heursure)的去噪更彻底,而由于rigrsure和minimaxi阈值选取规则较为保守(阈值较小,导致只有部分系数置零)噪声去除不彻底,但在信号的高频成分和噪声有重叠时利用这两种阈值可以把弱小信号从噪声中分离出来。图2　四种阈值选取方式对比2)软门限阈值和硬门限阈值处理比较(sorh的设置)实验表明软门限阈值处理方式一般能够取得更为平滑和理想的去噪效果。3)scal的比较由于我们的实验信号包含的是标准高斯白噪声,此时scal的设定对去噪效果影响不大。4)分解层次的比较(level)在Matlab下运行下面程序段:forlevel=1:10,xd=wden(x,`heursure',`s',`one',level,`sym8');subplot(2,5,level);plot(xd);title([`level'num2str(level)]);end;可以得到不同分解层次的去噪效果图,从中可以看出随着小波分解层次的增加去噪效果变好,但是分解层次增加到5层以上后去噪效果改善已经不明显,反而增加了计算代价,实际应用中分解层次取3～5层即可。5)小波函数的选取(wname)小波函数可以选取一个正交小波,如Daubechies(dbN),symlets(symN),coiflets(coifN),具体选择可以根据实际需要决定,在我们的实验中选择的是sym8。4.3　实际电机辨识数据去噪对比实验我们对一个实际的计算机控制的直流电机随动系统的模型辨识数据用滤波器,并采用傅立叶变换/反变换和小波收缩方法作了去噪对比实验。首先我们给实际随动系统加上辨识输入信号:ControlU(k)=500+PRBS75(k)即在系统工作点500个数字量输入信号的基础上加上幅值为75个数字量的伪随机二进制信号,读取系统的输出,由于系统传感器量测噪声和各种电磁干扰的影响,使得直接由量测数据辨识得到的系统模型阶次过高(3阶以上),而且这些较高的阶次所辨识的是噪声的模型,难以得到一个清晰简洁的被控对象的模型,给进一步对各种控制算法的设计和控制精度的提高造成了障碍。我们分别用滤波器,傅立叶变换/反变换和小波收缩对对辨识信号中的一段进行去噪处理,见图3,其中a)是原始辨识数据。最简单去噪方法是设计滤波器对量测数据进行滤波,我们使用了一个2阶截止频率为15Hz的FIR滤波器,滤除了大多数噪声,见图3b);直接对含噪信号作傅立叶变换(FFT)(在作FFT前先去除了信号中的直流成分)得到信号的频谱,硬性的将频谱中的200～800Hz频段的增益置为0,再用傅立叶反变换(IFFT)还原得到去噪后的数据,从图上可以看出去噪效果和滤波器近似,见图3c),但是傅立叶变换/反变换方法在允许的情况下可以对信号的