exp_6.m

先由原始图像（任选）产生待恢复的图像；(产生方法如下：冲激 函数为

%-----------------------------------------------------------%
%------------------------图像恢复---------------------------%
%-----------------------------------------------------------%
f=double(imread('Cameraman.bmp'));%输入实验图像 灰度图像
%-----------------------------------------------------------%
%------------------------图像模糊---------------------------%
%-----------------------------------------------------------%
[A B]=size(f);%f是A×B尺寸
%冲激函数模板 使用9*9模板 即h是9×9尺寸，设为C×D
for i=1:9
    for j=1:9
          h(i,j)=exp(sqrt((i-5)^2+(j-5)^2)/240);
    end
end
h=h/sum(sum(h));
[C D]=size(h);
%对f和h进行延拓，周期延拓为M=A+C-1，N=B+D-1
M=A+C-1;
N=B+D-1;
fe=f;
fe(M,:)=0;
fe(:,N)=0;
he=h;
he(M,:)=0;
he(:,N)=0;
%将f和h的延拓结果进行卷积得到模糊图像
g=conv2(fe,he);
ge=g(1:M,1:N);
%添加高斯随机噪声
ne=randn([M,N]);
ge_8=ge+sqrt(8)*ne;%均值为0，方差为8的高斯随机噪声
ge_16=ge+sqrt(16)*ne;%均值为0，方差为16的高斯随机噪声
ge_32=ge+sqrt(32)*ne;%均值为0，方差为32的高斯随机噪声

figure(1)
subplot(2,2,1);
imshow(f,[]);
title('Cameraman.bmp');
subplot(2,2,2);
imshow(ge_8((M-A)/2+1:(M+A)/2,(N-B)/2+1:(N+B)/2),[]);
title('模糊并带有方差为8高斯随机噪声的图像');
subplot(2,2,3);
imshow(ge_16((M-A)/2+1:(M+A)/2,(N-B)/2+1:(N+B)/2),[]);
title('模糊并带有方差为16高斯随机噪声的图像');
subplot(2,2,4);
imshow(ge_32((M-A)/2+1:(M+A)/2,(N-B)/2+1:(N+B)/2),[]);
title('模糊并带有方差为32高斯随机噪声的图像');
%-----------------------------------------------------------%
%-------------------------逆滤波----------------------------%
%-----------------------------------------------------------%
%取频域
G_8=fft2(ge_8);
G_16=fft2(ge_16);
G_32=fft2(ge_32);

H=fft2(he);

%为了防止H(u,v)在UV平面上取0或很小且消除振铃效应
%如果  H(u,v)<=d,   W(u,v)=k
%      H(u,v)>d,    W(u,v)=1/H(u,v)
% 其中k和d为小于1的常数，且选的较小为好
d=0.1;
k=0.1;
for u=1:M
    for v=1:N
        if abs(H(u,v))<=d
            W(u,v)=k;
        else
            W(u,v)=1/H(u,v);
        end
    end
end

f_8_1=abs(ifft2(G_8.*W));
f_16_1=abs(ifft2(G_16.*W));
f_32_1=abs(ifft2(G_32.*W));

figure(2)
subplot(2,2,1);
imshow(f,[]);
title('Cameraman.bmp');
subplot(2,2,2);
imshow(f_8_1(1:A,1:B),[]);
title('逆滤波恢复结果(噪声方差8)');
subplot(2,2,3);
imshow(f_16_1(1:A,1:B),[]);
title('逆滤波恢复结果(噪声方差16)');
subplot(2,2,4);
imshow(f_32_1(1:A,1:B),[]);
title('逆滤波恢复结果(噪声方差32)');
%-----------------------------------------------------------%
%------------------------维纳滤波---------------------------%
%-----------------------------------------------------------%
F=fft2(fe);

%高斯随机噪声的FFT
N_8=fft2(sqrt(8)*ne); 
N_16=fft2(sqrt(16)*ne);
N_32=fft2(sqrt(32)*ne);

F_8=conj(H)./(abs(H).^2+abs(N_8).^2./abs(F).^2).*G_8;
F_16=conj(H)./(abs(H).^2+abs(N_16).^2./abs(F).^2).*G_16;
F_32=conj(H)./(abs(H).^2+abs(N_32).^2./abs(F).^2).*G_32;

f_8_2=abs(ifft2(F_8));
f_16_2=abs(ifft2(F_16));
f_32_2=abs(ifft2(F_32));

figure(3)
subplot(2,2,1);
imshow(f,[]);
title('Cameraman.bmp');
subplot(2,2,2);
imshow(f_8_2(1:A,1:B),[]);
title('维纳滤波恢复结果(噪声方差8)');
subplot(2,2,3);
imshow(f_16_2(1:A,1:B),[]);
title('维纳滤波恢复结果(噪声方差16)');
subplot(2,2,4);
imshow(f_32_2(1:A,1:B),[]);
title('维纳滤波恢复结果(噪声方差32)');