1、小波域噪声抑制实验报告班级:姓名:学号:同做者: 2实验四 小波域噪声抑制一、实验目的在小波变换域实现图像噪声抑制。二、算法概要在分析突变信号和非平稳信号时,Fourier 分析显得无能为力,因此需要寻找新的分析工具,使之适合突变信号和非平稳信号。小波分析正是在这一背景下产生的。小波分析基本思想来源于可变窗口的伸缩和平移。Mallat 提出了实现小波变换的金字塔快速算法(Mallat 算法) ,使得小波变换真正进入了实际应用。小波变换是一种信号的时间-频率分析方法,它具有多分辨分析的特点,而且在时域和频域都具有表征信号局部特征的能力,是一种窗口大小固定不变但其形状可改变而且时间窗和频率窗都可以
2、改变的时频局部化分析方法。根据高频和低频使时间-频率变窄或变宽,即在低频部分具有较高的频率分辨率和较低的时间分辨率,在高频部分具有较高的时间分辨率和较低的频率分布率,很适合于探测正常信号中加带的瞬态反常现象并展示其成分,所以被誉为分析信号的显微镜。原则上讲,传统上使用傅立叶分析的地方,都可以用小波分析取代。小波分析优于傅立叶变换的地方是,它在时域和频域同时具有良好的局部化性质。1.连续小波变换设 ,其傅立叶变换为 ,当 满足允许条件(完全重构条件或恒2xLRw等分辨条件)(4-1)2RCd时,称 为一个基本小波。将 伸缩和平移后得xx(4-2 ),1,;0abbaR其中 为伸缩因子, 为平移因
3、子。a对于任意的函数 的连续小波变换和逆变换为2fxLR(4-3)12,2,1,xbf abaRfWf dC2.离散小波变换在实际应用中,尤其是在计算机上实现时,连续小波必须加以离散化。因此,有必3要讨论连续小波 和连续小波变换 的离散化。这一离散化都是针对连续尺度,abx,fWab参数 个连续平移参数 的,而不是针对变量 的。为了使离散后的函数组能覆盖整个a x所表示的平面,取 ,使得,b01,(4-4)00,jkabnajZ且将 改记为 ,即,abx,jkx(4-5 )220,000jj jjjkxkaxkba相应的离散小波变换为(4-6 ),f jkRCjkfxd取 ,则 ,就形成离散小
4、波变换中最常用的二进小波。02,1ab2,jjjkxx3.二维小波变换由于图像是二维信号,考虑二维尺度函数是可分离的情况,也就是:(4-7)(,)()xyy其中 是一维尺度函数。 是与之相对应的小波函数,那么就可以得到二维小()x波变换的三个基本小波:(4-8)123(,)(),xyy需要注意的是这里使用的上标只是索引而不是指数。 的函数 可以写为如下形式2()LRf(4-9 ), ,1,23(,)(,)()lkj jkkfxycxydxy由 Mallat 算法可知,二维图像的分解可以通过沿 方向和 方向分别进行一维滤波得到。于是图像的正交小波分解可以理解为一组独立的空间有向的频率通道上的信号
5、分解。每一尺度分解成四个子带 LL、HL、LH 和 HH,分别表征图像的低频信息及水平、垂直和斜方向上的细节。图 4.1 所示为图像的三层小波分解的结果。图 4.1 图像的三层小波分解示意图从小波变换的能量观点来看,高斯噪声的小波变换仍然是高斯分布,它是均匀分布在相空间的各个部分,而信号由于其带限性,它的小波变换系数仅仅集中在相空间的上的一4小部分。在小波域上,所有的小波系数都对噪声有贡献,也就是噪声的能量分布在所有的小波系数上,而只有一小部分小波系数对信号能量有贡献,所以可以把小波系数分为两类:第一类小波系数仅仅由噪声变换得到,这类小波系数幅值小,数目较多;第二类小波系数由信号变化得到,并包
6、含噪声的变换结果,这类小波系数幅值大,数目较少。从这点出发,可以通过系数幅值上的差异构造来达到去噪的效果。小波阈值去噪分三个步骤:利用小波变换对图像进行分解;设置一定的阈值对小波系数进行处理;小波逆变换得到空域结果图像。其中最关键的就是如何选择阈值,因为它控制图像噪声和信息的分离。对小波系数的阈值处理有硬门限和软门限两种方法,硬门限将小于阈值的小波系数置为零,大于阈值的保留;软门限是基于Donoho软门限思想的小波系数的非线性方法。 硬阈值方法:(4-10 ),0hWtTt软阈值方法(4-11 )sgn,0s ttTt其中 为小波系数, ,这是 Donoho 基于正交离散小波变换推倒出来的通用
7、W2logtN阈值方法,其中 为图像的长度, 为噪声标准方差。Donoho 提出的统一阈值计算简单,N已证明在高斯白噪声下对光滑信号是渐进最优的,故得到很广泛的应用。其中: 是像素点的坐标;,xy伸缩因子;a, 平移因子;bk小波分解低频子带系数;,lc小波分解高频子带系数;,jkd小波分解子带系数的统一表示;W阈值化处理后的各子带小波系数;各子带系数标准方差;降噪阈值。t三、算法步骤及流程图Step1:对原图像进行 3 级小波分解;Step2:通过式(4-10)对各子带系数进行阈值化处理;Step3:对阈值化处理的小波系数进行小波逆变换。5流程图如下:主程序流程图:硬阈值处理子程序流程图:
8、软阈值处理子程序流程图:输入原图像输出软阈值结果计算去噪阈值 t计算噪声标准方差 o利用公式 4-11进行阈值处理计算去噪阈值 t令噪声标准差 o=20输入图像,将图像像素点数记为图像大小N输出硬阈值结果利用公式 4-10进行阈值处理输入图像,将图像像素点数记为图像大小 N分别输出原图像、处理后图像、原图像直方图和去噪后图像直方图调用子程序对阈值化处理的小波系数进行小波逆变换调用子程序对各子代系数进行阈值处理调用子程序对原图像进行三级小波分解6四、实验程序主程序:tic;i=imread(SAR1.tif);i = i(:,:,1);figure;imshow(uint8(i);if isa(
9、i, double)i = im2double(i);endiSqmag=i.2;R=1000/max(max(iSqmag);iSqmag=iSqmag*R;A1,HL1,LH1,HH1=wt3det_spline(iSqmag,0);A2,HL2,LH2,HH2=wt3det_spline(A1,1);A3,HL3,LH3,HH3=wt3det_spline(A2,2);A4,HL4,LH4,HH4=wt3det_spline(A3,3);%Scale 3HH3_new=Threshold(HH3);HL3_new=Threshold(HL3);LH3_new=Threshold(LH3)
10、;%Scale 2HH2_new=Threshold(HH2);HL2_new=Threshold(HL2);LH2_new=Threshold(LH2);%Scale 1HH1_new=Threshold(HH1);HL1_new=Threshold(HL1);LH1_new=Threshold(LH1);% ReconstructionA2_new=iwt3det_spline(A3,HL3_new,LH3_new,HH3_new,2);A1_new=iwt3det_spline(A2_new,HL2_new,LH2_new,HH2_new,1);7A0_new=iwt3det_splin
11、e(A1_new,HL1_new,LH1_new,HH1_new,0);idespeck=abs(A0_new/R).0.5;idespeck_r=1*(idespeck-min(idespeck(:)/(max(idespeck(:)-min(idespeck(:);%Show the denoised image despeckledSAR=idespeck_r*255;figure(1);imshow(i*255,); title(-);figure(2);imshow(despeckledSAR,); title();figure(3);imhist(uint8(i*255);titl
12、e(-);figure(4);imhist(uint8(despeckledSAR);title();toc;硬阈值处理子程序:function current= priorshrinkG(former)a,b=size(former);N=a*b;o=std2(former);t=o*sqrt(2*log10(a*b);former(find(abs(former)t)=0; current=former;return;软阈值处理子程序:function current=Threshold(former)a,b=size(former);N=a*b;o=20;t=o*sqrt(2*log10(N);current(find(abs(former)t)=0;return; 8五、实验结果原图 硬阈值处理后图像软阈值处理后图像9六、算法综述在基于小波的图像去噪算法中,细节保持和平滑是总是相互冲突。图像中的强斑点噪声呈现脉冲特性,因此这些噪声像素点难以用线性平滑滤波器滤除。而中值滤波算法适合去除脉冲噪声并能很好地保持像边缘这样的细节结构。也可以使用局部几何结构来分类脉冲噪声或细节像素,中值滤波器滤除脉冲噪声而细节点保持不变。这样的非线性滤波器较好地保持了图像的细节同时抑制了斑点噪声。