多小波变换的源码，比单小波变换具有更好的性质，例如可同时具有对称性和正交性

flt='cl';fltap='clap'; %使用的多小波和预滤波器 deimg='c:\denoise\cameraman512.bmp'; %要处理的图像 img=imread(deimg); %从c:\noise目录下读入图像 %cc=0.0384;img1=imnoise(img,'gaussian',0,cc); %在某些论文中,都说加入了sigma=??的噪声,例如sigma=10等等, %见论文<基于小波变换的自适应多阀值图像去噪,中国图形图像学报,2005年第5期,表1和表2>, %我认为,应该用下面一行代码就可在MATLAB中实现,这是因为, %若数列Xn的标准差是1,则由D(cc*Xn)=cc^2*D(cc*Xn)=cc^2,得cc*Xn标准差是cc, %其中函数randn(m,n)生成均值是0,标准差是1的正态分布的mxn矩阵, %那么,它与使用imnoise(img,'gaussian',0,sigma_square)有什么关系? %imnoise中的sigma_square是方差,0=sigma_square=1,而cc是标准差,0=cc=255, %两者的关系是: sigma_square=(cc/255)^2 或者 sigma=cc/255, %例如:cc=10对应sigma=10/255, 也对应sigma_square=(10/255)^2 %另外,下面的这行程序已经变成一个此目录下的子程序,如下: % noise_img=addnoise(img,sigma) %上面的函数对于给定的图像矩阵img,加入均值为0,标准差为sigma的符合正态分布的噪声, %即高斯噪声,它与imnoise(img,'gaussian',0,sigma_square)区别在于,imnoise中的 %sigma_square是规范化后的方差,它在[0,1]之间,而addnoise的sigma是没有规范化的标准差, %它在[0,255]之间,其两者的转换关系见上面的说明,你可以这样理解这两个函数的不同, %对于取值于[0,255]的灰度图像矩阵img,函数imnoise先将img中的每个元素变化到[0,1]之间, %然后加入噪声,而函数addnoise则在img中直接加入噪声. cc=50; [m,n]=size(img); img1=double(img); img1=img1+cc*randn(m,n); %disp(rmse(img,img1)); pre=prep2D_appe(img1,fltap); %预滤波 trans_num=3;w=dec2D_pe(pre,flt,trans_num); %进行trans_num次多小波变换,第1次图像分成4x4块,每块(m/4)x(n/4) %第2次变换,左上角的4个图像块,被变换成4x4的块,每块(m/16)x(n/16) %下面对每次变换后的细节系数使用阀值算法修改块内的小波系数 au=0;bu=0;num=1;t=0.6;rrr=4; %disp(max(max(abs(w)))); while num=trans_num+1 for i=1:m/4:m-m/4+1 for j=n/2+1:n/4:n-n/4+1 mat=w(i:i+m/4-1,j:j+n/4-1);dd=dim(mat); xx=diffx(diffx(mat));yy=diffy(diffy(mat)); %fprintf('(%f,%f),(%f,%f)\n',i,j,i+m/4-1,j+n/4-1); [mm,nn]=size(mat); %ee=energy(mat);disp(ee); u0=(max(max(mat))-min(min(mat)))/rrr; qq=1/(mm*nn); for ii=1:mm for jj=1:nn %u=u0*exp(-sqrt(log10(1+t*abs(xx(ii,jj)+yy(ii,jj))))); u=u0*dd*exp(-qq*t*abs(xx(ii,jj)+yy(ii,jj))^2); %if num=10 fprintf('%f*',u);num=num+1;else fprintf('%f\n',u);;num=0;end if mat(ii,jj)>u mat(ii,jj)=mat(ii,jj)-u*t;au=au+1; elseif mat(ii,jj)<-u mat(ii,jj)=mat(ii,jj)+u*t;au=au+1; else mat(ii,jj)=0;bu=bu+1; end end end w(i:i+m/4-1,j:j+n/4-1)=mat; %pause; end end for i=m/2+1:m/4:m-m/4+1 for j=1:n/4:n/2-n/4+1 mat=w(i:i+m/4-1,j:j+n/4-1);dd=dim(mat); xx=diffx(diffx(mat));yy=diffy(diffy(mat)); %fprintf('(%f,%f),(%f,%f)\n',i,j,i+m/4-1,j+n/4-1); [mm,nn]=size(mat); %ee=energy(mat); u0=(max(max(mat))-min(min(mat)))/rrr; qq=1/(mm*nn); for ii=1:mm for jj=1:nn %u=u0*exp(-sqrt(log10(1+t*abs(xx(ii,jj)+yy(ii,jj))))); u=u0*dd*exp(-qq*t*abs(xx(ii,jj)+yy(ii,jj))^2); %if num=10 fprintf('%f*',u);num=num+1;else fprintf('%f\n',u);;num=0;end if mat(ii,jj)>u mat(ii,jj)=mat(ii,jj)-u*t;au=au+1; elseif mat(ii,jj)<-u mat(ii,jj)=mat(ii,jj)+u*t;au=au+1; else mat(ii,jj)=0;bu=bu+1; end end end w(i:i+m/4-1,j:j+n/4-1)=mat; %pause; end end m=m/2;n=n/2;num=num+1; end ww=rec2D_pe(w,flt,trans_num); www=postp2D_appe(ww,fltap); %www=ordfilt2(www,5,ones(3)); %3x3中值滤波 %www=www-0.01*t*(diffx(diffx(www))+diffy(diffy(www))); fprintf('%% 滤波前:snr(img,img1)=%f, snr(img,www)=%f\n',snr(img,img1),snr(img,www)); www=ordfilt2(www,5,ones(3)); %www=wiener2(www,[3,3]); fprintf('%% 下面是图像%s当t=%f,sigma=%f时的处理结果:\n',deimg,t,cc); fprintf('%% multiwavelets->%s, frefilter->%s, u=0, sigma=%f\n',flt,fltap,cc); fprintf('%% mse(img,img1)=%f, mse(img,www)=%f\n',mse(img,img1),mse(img,www)); fprintf('%% rmse(img,img1)=%f, rmse(img,www)=%f\n',rmse(img,img1),rmse(img,www)); fprintf('%% psnr(img,img1)=%f, psnr(img,www)=%f\n',psnr(img,img1),psnr(img,www)); fprintf('%% snr(img,img1)=%f, snr(img,www)=%f\n',snr(img,img1),snr(img,www)); fprintf('%% w(i,j>u:%f%%, w(i,j)=u:%f%%\n\n',au/(au+bu)*100,bu/(au+bu)*100); subplot(121);colormap(gray(256));image(img1);axis('off'); subplot(122);colormap(gray(256));image(www);axis('off');;