小波域噪声抑制.doc

实验二 小波域噪声抑制一、实验目的在小波变换域实现图像噪声抑制。二、算法概要在分析突变信号和非平稳信号时，Fourier分析显得无能为力，因此需要寻找新的分析工具，使之适合突变信号和非平稳信号。小波分析正是在这一背景下产生的。小波分析基本思想来源于可变窗口的伸缩和平移。Mallat提出了实现小波变换的金字塔快速算法（Mallat算法），使得小波变换真正进入了实际应用。小波变换是一种信号的时间-频率分析方法，它具有多分辨分析的特点，而且在时域和频域都具有表征信号局部特征的能力，是一种窗口大小固定不变但其形状可改变而且时间窗和频率窗都可以改变的时频局部化分析方法。根据高频和低频使时间-频率变窄或变宽，即在低频部分具有较高的频率分辨率和较低的时间分辨率，在高频部分具有较高的时间分辨率和较低的频率分布率，很适合于探测正常信号中加带的瞬态反常现象并展示其成分，所以被誉为分析信号的显微镜。原则上讲，传统上使用傅立叶分析的地方，都可以用小波分析取代。小波分析优于傅立叶变换的地方是，它在时域和频域同时具有良好的局部化性质。1.连续小波变换设，其傅立叶变换为，当满足允许条件（完全重构条件或恒等分辨条件） （4-1）时，称为一个基本小波。将伸缩和平移后得 （4-2）其中为伸缩因子，为平移因子。对于任意的函数的连续小波变换和逆变换为 （4-3）2.离散小波变换在实际应用中，尤其是在计算机上实现时，连续小波必须加以离散化。因此，有必要讨论连续小波和连续小波变换的离散化。这一离散化都是针对连续尺度参数个连续平移参数的，而不是针对变量的。为了使离散后的函数组能覆盖整个所表示的平面，取，使得 （4-4）且将改记为，即 （4-5）相应的离散小波变换为 （4-6）取，则，就形成离散小波变换中最常用的二进小波。3.二维小波变换由于图像是二维信号，考虑二维尺度函数是可分离的情况，也就是： （4-7）其中是一维尺度函数。是与之相对应的小波函数，那么就可以得到二维小波变换的三个基本小波： （4-8）需要注意的是这里使用的上标只是索引而不是指数。的函数可以写为如下形式 （4-9）由Mallat算法可知，二维图像的分解可以通过沿方向和方向分别进行一维滤波得到。于是图像的正交小波分解可以理解为一组独立的空间有向的频率通道上的信号分解。每一尺度分解成四个子带LL、HL、LH和HH，分别表征图像的低频信息及水平、垂直和斜方向上的细节。图4.1 所示为图像的三层小波分解的结果。图4.1 图像的三层小波分解示意图从小波变换的能量观点来看，高斯噪声的小波变换仍然是高斯分布，它是均匀分布在相空间的各个部分，而信号由于其带限性，它的小波变换系数仅仅集中在相空间的上的一小部分。在小波域上，所有的小波系数都对噪声有贡献，也就是噪声的能量分布在所有的小波系数上，而只有一小部分小波系数对信号能量有贡献，所以可以把小波系数分为两类：第一类小波系数仅仅由噪声变换得到，这类小波系数幅值小，数目较多；第二类小波系数由信号变化得到，并包含噪声的变换结果，这类小波系数幅值大，数目较少。从这点出发，可以通过系数幅值上的差异构造来达到去噪的效果。小波阈值去噪分三个步骤：利用小波变换对图像进行分解；设置一定的阈值对小波系数进行处理；小波逆变换得到空域结果图像。其中最关键的就是如何选择阈值，因为它控制图像噪声和信息的分离。对小波系数的阈值处理有硬门限和软门限两种方法，硬门限将小于阈值的小波系数置为零，大于阈值的保留；软门限是基于Donoho软门限思想的小波系数的非线性方法。 硬阈值方法： （4-10）软阈值方法 （4-11）其中为小波系数，这是Donoho基于正交离散小波变换推倒出来的通用阈值方法，其中为图像的长度，为噪声标准方差。Donoho提出的统一阈值计算简单，已证明在高斯白噪声下对光滑信号是渐进最优的，故得到很广泛的应用。三、变量说明 是像素点的坐标； 伸缩因子；, 平移因子； 小波分解低频子带系数；小波分解高频子带系数； 小波分解子带系数的统一表示； 阈值化处理后的各子带小波系数； 各子带系数标准方差； 降噪阈值。四、算法步骤Step1：对原图像进行3级小波分解；Step2：通过式（4-10）对各子带系数进行阈值化处理；Step3：对阈值化处理的小波系数进行小波逆变换。五、实验内容与要求（！）例程给出了一种冗余小波变换和逆变换的降噪算法，Threshold()函数实现小波系数的阈值化处理。（2）尝试对去噪前后的图像进行FCM聚类分割，对比结果。六，参考程序（1）冗余小波变换和逆变换的降噪算法 tic;i=imread(SAR1.tif);figure;imshow(uint8(i);%-从八位到十六位ifisa(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);%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,HL2_new,LH2_new,HH2_new,1);A0_new=iwt3det_spline(A1,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 imagedespeckledSAR=idespeck_r*255;figure(1);imshow(i*255,);title(原图像);figure(2);imshow(despeckledSAR,);title(去噪图像);figure(3);imhist(uint8(i*255);title(原图像直方图);figure(4);imhist(uint8(despeckledSAR);title(去噪图像直方图);fcm_noise(i*255);fcm_noise(despeckledSAR);toc;Threshold（）函数function y=Threshold(W)a,b=size(W);N=a*b;S=std2(W)t=S*sqrt(2*log10(N)/7;for i=1:a for j=1:bif abs(W(i,j)t y(i,j)=0;else %y(i,j)=W(i,j); %hard %y(i,j)=sign(y(i,j)*(abs(W(i,j)-t); %softend endend（2）对去噪前后的图像进行FCM聚类分割的程序function fcm(name)IM=name;%IM=despeckledSAR;IM=double(IM);figure(1);imshow(uint8(IM);maxX,maxY=size(IM);IMM=cat(3,IM,IM,IM); cc1=8;cc2=100;cc3=200;ttfcm=0;while(ttfcm15) ttfcm=ttfcm+1; c1=repmat(cc1,maxX,maxY); c2=repmat(cc2,maxX,maxY); c3=repmat(cc3,maxX,maxY); c=cat(3,c1,c2,c3); ree=repmat(0.000001,maxX,maxY); ree1=cat(3,ree,ree,ree); distance=IMM-c; distance=distance.*distance+ree1; daoShu=1./distance; daoShu2=daoShu(:,:,1)+daoShu(:,:,2)+daoShu(:,:,3); distance1=distance(:,:,1).*daoShu2; u1=1./distance1; distance2=distance(:,:,2).*daoShu2; u2=1./distance2; distance3=distance(:,:,3).*daoShu2; u3=1./distance3; ccc1=sum(sum(u1.*u1.*IM)/sum(sum(u1.*u1); ccc2=sum(sum(u2.*u2.*IM)/sum(sum(u2.*u2); ccc3=sum(sum(u3.*u3.*IM)/sum(sum(u3.*u3); tmpMatrix=abs(cc1-ccc1)/cc1,abs(cc2-ccc2)/cc2,abs(cc3-ccc3)/cc3; pp=cat(3,u1,u2,u3); for i=1:maxX for j=1:maxY if max(pp(i,j,:)=u1(i,j) ix2(i,j)=1; elseif max(pp(i,j,:)=u2(i,j) ix2(i,j)=2; else ix2(i,j)=3; end end end if max(tmpMatrix)0.0001 break; else cc1=ccc1; cc2=ccc2; cc3=ccc3; end for i=1:maxX for j