基于小波分析的信号去噪

Matlab与通信系统分析基于小波变换的信号去噪仿真院别

计算机与通信工程学院通信工程2012 年06月9日目录概述 1小波变换 1小波变换与信号去噪 1小波分析基本理论 2小波去噪原理 2小波去噪步骤 3小波去噪的阈值 4阈值函数 4阈值的选取 5三种阈值去噪形式 6小波去噪的MATLAB实现 7MATLAB中小波去噪的函数集合 7小波去噪的MATLAB仿真 95.结语 13参考文献 14人员分工 14概述小波变换J.Morlet1974A.CalderonHardy空间的原子分解和无条件基的深入研究J.O.Stromberg构造了历史上非常类似于现在的小1986S.Mallat立了构造小波基的统一方法——多尺度分析之后，小波分析才开始蓬勃发展起来，之I.Daubechies的《小波十讲》对小波的普及起了重要的推动作用。现在，小波变换已经在科技信息产业领域取得了令人瞩目的成就。小波变换具有Fourier变换不能解决的许多困难问题，因而小波变化被誉为“数学显微镜B超、CT、核磁共振成像的时间，提高分辨率等。小波变换与信号去噪效果了。1995年Donoho和Johnstone提出了小波收缩去噪的技术，他们研究的是在叠加性高1斯白噪声环境下检测出真实信号的情况，利用正交小波变换和高斯随机变量的性质对信号的小波分解系数做阈值量化，无失真的还原出真实信号。由小波变换的特性可知，高斯噪声的小波变换仍然是高斯分布的，它均匀分布在频率尺度空间的各部分，而信号由于其带限性，它的小波系数仅仅集中在频率尺度空间上的有限部分。相比于传统的傅里叶变换，小波分析方法是一种窗口大小（即窗口面积）固定但其形状可改变，时间窗和频率窗都可改变的时频局部化分析方法，它在低频部分具有较宽的时间窗和较高的频率分辨率，在高频部分具有较宽的频率窗和较高的时间分辨率。正是这种特性，使小波变换具有对信号的自适应性。本文对小波去噪方法进行了研究，验证了小波去噪的可靠性并对比了傅里叶去噪和小波去噪的效果。小波分析基本理论小波去噪原理叠加性高斯白噪声是最常见的噪声模型，受到叠加性高斯白噪声“污染”的观测信号可以表示为：yfi i

zi

i1,., (2.1)yz~N，i i i i为噪声水平，信号长度为n.为了从含噪信号yi

中还原出真实信号fi

，可以利用信号和噪声在小波变换下的不同的特性，通过对小波分解系数进行处理来达到信号和噪声分离的目的。在小波域，有效信号所产生的小波系数其模值往往较大；而噪声经过小波变换有白化趋势（高斯白噪声经过小波变换，在小波域仍然表现为很强的随机性，通常仍认为是高斯分布），其小波域对应系数模值很小。而且有用信号通常表现为低频信号或是一些比较平稳的信号，而噪声信号则通常表现为高频信号，所以我们可以先对含噪信号进行小波分解（如进行三层分解）：SCA1

CD1CA CD CD2 2 1

(2.2)CA3

CD1

CD2

CD32图2.1三层小波分解示意图其中CAi为分解的近似部分，CDii1,2,则噪声部分通常包含在CD1CD2CD3中，用门限阈值对小波系数进行处理（零），最后重构信号即可达到去噪的目的。小波去噪步骤总结去噪过程，可以分成以下三个步骤：对观测数据作小波分解变化：WyW f Wz (2.3)0 0 0yyy1 2

，„yn

fff1 2

，„fn

z是zz1 2

，„zn

，其中用到了小波分解变换是线性变换的性质。对小波系数W0

作门限阈值处理（根据具体情况可以使用软阈值处理或硬阈值22lognt n

(2.4)门限阈值处理可以表示为 可以证明当n趋于无穷大时使用阈值公式(2.4)对小波tn系数作软阈值处理可以几乎完全去除观测数据中的噪声。3对处理过的小波系数作逆变换

1重构信号：0fW

0

Wtn

d (2.5)即可得到受污染采样信号去噪后的信号。小波去噪的阈值阈值函数在对小波系数作门限阈值处理操作时，可以使用软阈值处理方法或硬阈值处理方法，硬阈值处理只保留较大的小波系数并将较小的小波系数置零，表达式见式(3.1)。 (,t) ,

wwt

(3.1)软阈值处理将较小的小波系数置零但对较大的小波系数向零作了收缩，表达式见式(3.2)。wt,0, (w,t)0,Swt,

wtwt (3.2)wt直观形式见图3.1（图中取t=0.4）从图上我们可以看出软阈值处理是一种更为平滑小波系数连续性差,在某些点上会产生间断，可能会引起重构信号的振荡。软阈值处理实质上是在硬阈值的基础上将边界出现的不连续点收敛到零，对小波分解系数作了收缩，从而Donoho-Johnstone将这种去噪技术称之为小波收缩。这样可以有效避免间断，获得比较光滑的重建新号。针对硬阈值函数的不连续性和软阈值函数存在恒定偏差,以及阈值函数中参数选择的问题,提出了渐近半软阈值函数。半软阈值形式较为复杂，它将软阈值和硬阈值的优点结合起来，更有效地抑制了噪声，提高了信噪比，最大限度地提高了去噪的效果。由于篇幅所限，这里不多加介绍。4硬阈值信号10.50-0.5-10 10 20 30 40 50 60 70 80 90 100软阈值信号10.50-0.5-10 10 20 30 40 50 60 70 80 90 100阈值的选取

图3.1硬阀值和软阀值目前常见的阈值选取规则主要有固定阈值估计，极值阈值估计，无偏似然估计和启发式阈值估计几种。sqtwolog:2log(Length(2log(Length(N))其中，N为信号长度。minimaxisqtwolog量是实现在最坏条件下最大均方误差最小的任选量。0, n32t[0.39360.1829log2n],n32

(3.5)rigrsure:采用史坦的无偏似然估计原理进行阈值选择，首先得到一个给定阈值5t的思然估计，再将非似然的t最小化，就得到了所选的阈值。heursure：选择启发式阈值，它是sqtwolog和rigrsure的综合，是最优预测变量heursure采用固定阈值sqtwolog一般来讲，极值阈值估计和无偏似然估计方法比较保守，当噪声在信号的高频段分布较少时，这两种阈值估计方法效果较好可以将微弱的信号提取出来。而固定阈值估计和启发式阈值估计法去噪比较彻底，在去噪时显得更为有效，但是也容易把有用的信号误认为噪声而去掉。三种阈值去噪形式在实际应用中，处理信号和图像常借助计算机工具。MATLAB相应的为我们提供了小波分析工具箱，对于利用小波变换方法对数据的处理提供了方便。利用小波分析工具箱中的函数对信号和图像进行阈值去噪大致可分为三种：比较平滑，但容易丢失信号的有用成分。默认阈值去噪处理。该方法利用''ddencmp''用''wdencmp''函数进行去噪处理。给定软(或硬)阈值去噪处理。在实际的去噪处理中，阈值往往可以通过经验公''wthresh''函数进行。在三种方法中，强制去噪处理方法较简单，去噪后的信号较平滑，但容易丢失有用信号。对于第三种方法，阈值需要通过经验公式获得。因此此方法本文也不考虑。对于初学者来说，自己设定阈值是一件较困难的事情，为此，对于阈值需要根据经验设定的选项不予以考虑。所以本文中采取默认阈值方法消噪。另外在利用阈值进行小波系数的处理时，可以选择全局阈值处理和部分阈值处理。全局阈值即在每层分解的小波系数进行处理时采用相同的阈值。局部阈值即在每层分解的小波系数进行处理时根据不同层的噪声估计值不断调整阈值，采取不同的阈值。全局阈值对于不同特征的信号区分度差，不论信号的低频域分量小，还是高频分量小，都是由默认阈值来处理，有时效果会较差。比如对于含高频分量较大，而低频分量较小的信号，采用全局阈值算法消噪往往会丢失大量的高频部分信息，从而引起偏差。6但是全局阈值较为简单，且可由''ddencmp''直接获得，若采用局部阈值处理，仍然需要按照经验获取阈值向量，本文不予考虑。综上，本文采用全局默认阈值的处理方法。小波去噪的Matlab实现MATLAB噪方面实现Donoho-Johnstone等的去噪算法，而且可以选择使用图形界面操作工具或者去噪函数集合两种形式，这里利用小波去噪函数集合编写MATLAB程序，完成仿真。4.1Matlab中小波去噪的函数集合下面是几个最为常用的与小波去噪有关的函数：wnoise函数：调用方式：x=wnoise(fun,n,snr)；作用：产生Donoho-Johnstone设计的6种用于测试小波去噪效果的典型测试数据，函数根据输入参数fun的值输出名为''blocks''，''bumps''，''heavy''，''doppler''，''quadchirp''或''mishmash''的6种函数数据，数据长度为2n。信噪比为''snr''。这6种测试数据在验证和仿真实验时非常有用。wden函数：调用方式：[XD,CXD,LXD]=wden(X,TPTR,SORH,SCAL,N,''wname'')；[XD,CXD,LXD]=wden(C,L,TPTR,SORH,SCAL,N,''wname'')；作用：''wden''是最主要的一维小波去噪函数。其中输入参数''x''为输入需要的信号''TPTR''为3.2节中4种阀值形式，''SORH''设定为''s''表示用软门限阈值处理表示用硬门限阈值处理3.2节中提到的4种阈值形式是基于信号水平为1的高斯白噪声模型推导得到的，当噪声不是白噪声时，必须在小波分解的不同层次估计噪声水平，SCAL=''one''不进行重新估计，SCAL=''sln''只根据第一层小波分解系数估计噪声水平，SCAL=''mln''在每个不同的小波分解层次估计噪声水平根据''SCAL''参数的设定''wden''函数决定最终应用于每一个小波分解层次的阈值函数。最后两个参数 ''N''和''wname''表示利用名为''wname''的小波对信号进行N层分解除''wden''外还有功能更强大的用于一维或二维小波去噪或压缩的函数''wdencmp''。ddencmp函数：7调用方式：[THR,SORH,KEEPAPP,CRIT]=ddencmp(IN1,IN2,X)；[THR,SORH,KEEPAPP,CRIT]=ddencmp(IN1,''wp'',X)；[THR,SORH,KEEPAPP,CRIT]=ddencmp(IN1,''wv'',X)；作用：函数''ddencmp''用于获取信号在消噪或压缩过程中的默认阈值。输入参数''X''为一维或二维信号；''IN1''取值为''den''或''cmp''，''den''表示进行去噪，''cmp''表示进行压缩；''IN2''取值为''wv''或''wp''，''wv''表示选择小波，''wp''表示选择小波包。返回值''THR''是返回的阈值；''SORH''是软阈值或硬阈值选择参数；''KEEPAPP''表示保存低频信号；''CRIT''是熵名（只在选择小波包时使用）。wdencmp函数：调用方式：[XC,CXC,LXC]=wdencmp(''gbl'',X,''wname'',N,THTR,SORH,KEEPAPP)；[XC,CXC,LXC]=wdencmp(''lvd'',X,''wname'',N,THTR,SORH)；[XC,CXC,LXC]=wdencmp(''lvd'',C,L,''wname'',N,THTR,SORH)；''wdencmp''''wname''是所用的小波函数，''gbl''(''global''的缩写)表示每层采用不同的阈值进行处理，''N''表示小波分解的层数，''THR''为阈值向量，对于后面两种格式，''THR''的长度为N，''SORH''表示选择软阈值或硬（分别取值为''s''和参数''KEEPAPP''取值为1''XC''是''XC''的小波分解结构。MATLAB中与小波去噪有关的函数总结见表4.1和表4.2。表4.1小波分解和重构函数分类分类函数名appcoefdetcoefdwtwavedecupcoefupwlevidwtwaverecwrcoef说明多尺度一维小波分解一维小波系数的单尺度重构单尺度一维离散小波逆变换多尺度一维小波重构8表4.2压缩和消噪函数分类名称说明阈值获取函数ddencmp获取在消噪和压缩中的默认阈值thselect去噪的阈值选择wbmpen获取一维小波去噪阈值wdcbm用Birge-Massart算法获取小波变换阈值去噪函数wden用小波变换对一维信号自动消噪wdencmp用小波进行消噪或压缩阈值处理函数wthcoef一维信号小波系数的阈值处理wthresh软阈值或硬阈值处理4.2小波去噪的MATLAB仿真以下为仿真源代码：clc;clearsnr=4; %设置信噪比%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%MATLAB中用''wnoise''产生测试信号%原始信号为xref，含高斯白噪声的信号为x%信号类型为blocks（由函数中参数1决定）%长度均为2^11（由函数中的参数11决定）%信噪比snr=4（由函数中的参数snr决定）%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%[xref,x]=wnoise(1,11,snr);xref=xref(1:2000); %2000点x=x(1:2000); %2000点%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%用全局默认阈值进行去噪处理%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%[thr,sorh,keepapp]=ddencmp('den','wv',x); %获取全局默认阈值xd=wdencmp('gbl',x,'sym8',3,thr,sorh,keepapp); %利用全局默认阈值对信号去噪%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%9%下面是作图函数，作出原始信号和含噪声信号的图%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%figuresubplot(211);plot(xref); %title('原始信号');subplot(212);plot(x);title('含噪声信号'); %画出含噪声信号的图%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%下面用傅里叶变换进行原信号和噪声信号的频谱分析%通过在频域上的分析选择合适的滤波器对信号进行去噪%最后将这一经典的去噪结果与小波去噪结果进行比较%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%dt=1/(2^11); %时域分辨率Fs=1/dt; %计算频域分辨率df=Fs/2000;xxref=fft(xref); %对原始信号做快速傅里叶变换xxref=fftshift(xxref); %将频谱图平移xxref=abs(xxref); %取傅里叶变换的幅值xx=fft(x); %对含噪声信号做快速傅里叶变换xx=fftshift(xx); %将频谱搬移absxx=abs(xx); %取傅里叶变换的幅值ff=-1000*df:df:1000*df-df; %figure(2)subplot(211);plot(ff,xxref); %title('原始信号的频谱图');subplot(212);plot(ff,absxx);title('含信号噪声的频谱图'); %画出含噪声信号的频谱图%进行低通滤波，滤波频率为0~200的相对频率indd2=1:800; %0频左边高频率系数置零10xx(indd2)=zeros(size(indd2));indd2=1201:2000;xx(indd2)=zeros(size(indd2)); %0频右边高频系数置零xden=ifft(xx); %滤波后的信号作傅里叶逆变换xden=abs(xden); %figure(3)subplot(211);plot(xd); %title('小波去除噪后的信号');subplot(212);plot(xden); %title('傅里叶分析去噪的信号');以下是仿真结果：原始信号和含高斯白噪声的信号波形见图4.1。原始信号和含高斯白噪声的信号经4.24.3。原始信号151050-50 200 400 600 800 1000 1200 1400 1600 1800 2000含噪声信号20100-100 200 400 600 800 1000 1200 1400 1600 1800 2000图4.1原始信号和含高斯白噪声的信号11100008000600040002000

原始信号的频谱图0-1500100005000