高频阈值处理matlab,[MATLAB图像处理] 帮忙看看小波阈值图像去噪程序的问题

帮忙看看小波阈值图像去噪程序的问题

clear all;

close all;

I=imread('lena.bmp');

X=rgb2gray(I);

[height width]=size(X);

X=im2double(X);

A=imnoise(X,'gaussian',0,0.01);

wname='sym4';

n=3;

[c,s]=wavedec2(A,n,wname);

for i=1:3

step(i)=s((i+1),1)*s((i+1),2);

end

num(1,1)=s(1,1)*s(1,2)+1; %获取各层各高频分量在c向量中的坐标

num(1,2)=num(1,1)+s(2,1)*s(2,2);

num(1,3)=num(1,2)+s(2,1)*s(2,2);

num(2,1)=num(1,3)+s(2,1)*s(2,2);

num(2,2)=num(2,1)+s(3,1)*s(3,2);

num(2,3)=num(2,2)+s(3,1)*s(3,2);

num(3,1)=num(2,3)+s(3,1)*s(3,2);

num(3,2)=num(3,1)+s(4,1)*s(4,2);

num(3,3)=num(3,2)+s(4,1)*s(4,2);

m=0.02;

C=c;

for i=1:3

[H,V,D]=detcoef2('a',c,s,i);%提取第i层各高频系数

B=[H V D];

[L,T]=size(B);

for k=1:L

for w=1:T

sigma=median(abs(B(k,w)))/0.6745;

end

end

th=sigma*sqrt(2*log10(L*T));%计算阈值

ch=c(1,num(4-i,1):num(4-i,3)+step(4-i)-1);%对各高频系数进行阈值处理

chl=length(ch);

for j=1:chl

if abs(ch(j))>=th

ch(j)=sign(ch(j))*(abs(ch(j))-2*th/(1+exp(m*(ch(j)^2))));%阈值处理函数

else

ch(j)=0;

end

end

C(1,num(4-i,1):num(4-i,3)+step(4-i)-1)=ch(1,1:chl);

end

X0=waverec2(C,s,wname);

disp('原图像信噪比为：'),sny=10*log10(sum(sum(double(X).^2))/sum(sum((double(X)-double(A)).^2)))

disp('去噪后图像信噪比为：'),sndz=10*log10(sum(sum(double(X).^2))/sum(sum((double(X)-double(X0)).^2)))

figure(1)

imshow(X)

title('原图像')

figure(2)

imshow(A)

title('添加噪声后的图像')

figure(3)

imshow(X0,[])

title('去噪后图像')