[人工智能] 图像傅里叶变换实验

开发: C++知识库 Java知识库 JavaScript Python PHP知识库人工智能区块链大数据移动开发嵌入式开发工具数据结构与算法开发测试游戏开发网络协议系统运维
教程: HTML教程 CSS教程 JavaScript教程 Go语言教程 JQuery教程 VUE教程 VUE3教程 Bootstrap教程 SQL数据库教程 C语言教程 C++教程 Java教程 Python教程 Python3教程 C#教程
数码: 电脑笔记本显卡显示器固态硬盘硬盘耳机手机 iphone vivo oppo 小米华为单反装机图拉丁

-> 人工智能 -> 图像傅里叶变换实验 -> 正文阅读

[人工智能]图像傅里叶变换实验

一、实验目的

1.掌握傅里叶变换及逆变换的基本原理方法。
2.理解频域滤波的基本原理及方法。
3.掌握进行图像的频域滤波的方法。

二、实验原理

1. 平滑滤波

(1) 线性平滑滤波器

线性低通平滑滤波器的所有系数都是正数，对3×3的模板来说，最简均值滤波器的是取所有系数为1，为了保持输出图像仍然在原图像的灰度值范围内，模板与像素邻域的乘积都要除以9。 MATLAB提供了fspecial函数生成滤波时所用的模板，并提供filter2函数用指定的滤波器模板对图像进行运算。MATLAB提供了一个函数imnoise来给图像增添噪声，请同学们通过MATLAB的help命令自行查阅fspecial、filter2及imnoise的用法。OpenCV也有类似的函数，请同学们自己查询掌握。

(2) 非线性平滑滤波器

中值滤波器是一种常用的非线性平滑滤波器，其滤波原理与均值滤波器方法类似，但不是计算加权求和，而是把邻域中的图像的像素按灰度级进行排序，然后选择该组的中间值作为输出像素值。 MATLAB提供了medfilt2函数来实现中值滤波。OpenCV也有类似的函数，请同学们自己查询。

2. 锐化滤波器

图像平滑往往使图像中的边界、轮廓变得模糊，为了减少这类不利效果的影响，需要利用图像锐化技术，使图像的边缘变得清晰。

(1) 线性锐化滤波器

线性高通滤波器是最常用的线性锐化滤波器。这种滤波器的中心系数都是正的，而周围的系数都是负的，所有的系数之和为0。

(2) 非线性锐化滤波

邻域平均可以模糊图像，因为平均对应积分，所以利用微分可以锐化图像。图像处理中最常用的微分方法是利用梯度。常用的空域非线性锐化滤波微分算子有sobel算子、prewitt算子、log算子等。

3. 频域增强

频域增强是利用图像变换方法将原来的图像空间中的图像以某种形式转换到其他空间中，然后利用该空间的特有性质方便地进行图像处理，最后再转换回原来的图像空间中，从而得到处理后的图像。
频域增强的主要步骤是：
(1）选择变换方法，将输入图像变换到频域空间。
(2）在频域空间中，根据处理目的设计一个传递函数，并进行处理。
(3）将所得结果用逆变换得到增强的图像。
(4）常用的频域增强方法有低通滤波和高通滤波。

4. 低通滤波

图像的能量大部分集中在幅度谱的低频和中频部分，而图像的边缘和噪声对应于高频部分。因此能降低高频成分幅度的滤波器就能减弱噪声的影响。由卷积定理，在频域实现低通滤波的数学表达式：

(1) 理想低通滤波器（ILPF）

(2) 巴特沃斯低通滤波器（BLPF）

(3) 指数型低通滤波器（ELPF）

5. 高通滤波

由于图像中的细节部分与其高频分量相对应，所以高通滤波可以对图像进行锐化处理。高通滤波与低通滤波相反，它是高频分量顺利通过，使低频分量受到削弱。高通滤波器和低通滤波器相似，其转移函数分别为:

(1) 理想高通滤波器（IHPF）

(2) 巴特沃斯高通滤波器（BLPF）

(3) 指数型高通滤波器（ELPF）

图像经过高通滤波处理后，会丢失许多低频信息，所以图像的平滑区基本上会消失。所以，可以采用高频加强滤波来弥补。高频加强滤波就是在设计滤波传递函数时，加上一个大于0小于1的常数，即：

三、实验内容步骤（记录实验主要步骤，且在调试成功后，将结果截屏或拍照保存）

要求使用MATLAB或OpenCV来完成以下实验。

1. 平滑空域滤波

（1) 读出一幅灰度图像，给这幅图像分别加入椒盐噪声和高斯噪声后并与前一张图显示在同一图像窗口中。

lena = imread('lena.jpg'); %读入原始图像
lena_gray = rgb2gray(lena); %原始图像灰度化

figure, subplot(1, 3, 1);
imshow(lena_gray);
title('原始图像');
lena_pepper = imnoise(lena_gray, 'salt & pepper', 0.02); %加入椒盐噪声

subplot(1, 3, 2);
imshow(lena_pepper);
title('椒盐噪声图像');

lena_gauss = imnoise(lena_gray, "gaussian", 0.02); %加入高斯噪声
subplot(1, 3, 3);
imshow(lena_gauss);
title('高斯噪声图像');

（2) 对加入噪声图像选用不同的平滑（低通）模板做运算，对比不同模板所形成的效果，要求在同一窗口中显示。

lena = imread('lena.jpg'); %读入原始图像

lena_gray = rgb2gray(lena); %原始图像灰度化
figure, subplot(1, 3, 1);
imshow(lena_gray);
title('原始图像');

lena_pepper = imnoise(lena_gray, 'salt & pepper', 0.02); %加入椒盐噪声
subplot(1, 3, 2);
imshow(lena_pepper);
title('椒盐噪声图像');

lena_gauss = imnoise(lena_gray, "gaussian", 0.02); %加入高斯噪声
subplot(1, 3, 3);
imshow(lena_gauss);
title('高斯噪声图像');

%对椒盐噪声图像进行滤波处理
h = fspecial('average', 3);
I1 = filter2(h, lena_pepper) / 255; %二维数字滤波器
I2 = medfilt2(lena_pepper, [3 3]); %二维中位数滤波
figure, subplot(2, 2, 1);
imshow(lena);
title('原始图像');

subplot(2, 2, 2);
imshow(lena_pepper);
title('椒盐噪声图像');

subplot(2, 2, 3);
imshow(I1);
title('3*3 均值滤波图像');

subplot(2, 2, 4);
imshow(I2);
title('3*3 中值滤波图像');

%对高斯噪声图像进行滤波处理
G1 = filter2(h, lena_gauss) / 255; %二维数字滤波器
G2 = medfilt2(lena_gauss, [3 3]); %二维中位数滤波

figure, subplot(2, 2, 1);
imshow(lena);
title('原始图像');

subplot(2, 2, 2);
imshow(lena_gauss);
title('高斯噪声图像');

subplot(2, 2, 3);
imshow(G1);
title('3*3 均值滤波图像');

subplot(2, 2, 4);
imshow(G2);
title('3*3 中值滤波图像');

（3) 如果使用MATLAB，使用函数imfilter时，分别采用不同的填充方法（或边界选项，如零填充、’replicate’、’symmetric’、’circular’）进行低通滤波，显示处理后的图像。

close all
lena=imread('lena.jpg');    %将原始图像读入
lena_gray=rgb2gray(lena);   %将图像灰度化
h=fspecial('motion',50,45);  %创建预定义的二维滤波器
filteredimg=imfilter(lena_gray,h);    % 对图像进行滤波,默认通过使用相关来完成
boundaryReplicate=imfilter(lena_gray,h,"replicate") ;    %图像大小通过复制外边界的值来扩展
boundary0=imfilter(lena_gray,h,0);   % 输入图像的边界通过用值X（无引号）来填充扩展，其默认值为0
boundarysymmetric=imfilter(lena_gray,h,"symmetric");  %图像大小通过镜像反射其边界来扩展
boundarycircular=imfilter(lena_gray,h,'circular');   %图像大小通过将图像看成是一个二维周期函数的一个周期来扩展

figure,subplot(3,2,1);%为各个分块创建分区
imshow(lena);   
title('原始图像');
subplot(3,2,2);
imshow(filteredimg);
title('Motion Blurred');
subplot(3,2,3);
imshow(boundaryReplicate);
title('Replicate');
subplot(3,2,4);
imshow(boundary0);
title('0-Padding');
subplot(3,2,5);
imshow(boundarysymmetric);
title('symmetric');
subplot(3,2,6);
imshow(boundarycircular);
title('cicular');

（4) 运用for循环，将加有椒盐噪声的图像进行10次、20次均值滤波，查看其特点，显示均值处理后的图像（提示:利用fspecial函数的’average’类型生成均值滤波器）。

close all
lena=imread('lena.jpg');    %将原始图像读入
lena_gray=rgb2gray(lena);   %将图像灰度化
lena_pepper=imnoise(lena_gray,'salt & pepper',0.02);    %加入椒盐噪声
h=fspecial('average');      %生成均值滤波器
for i=1:10  %对椒盐噪声图像进行10次均值滤波
    J1=imfilter(lena_pepper,h); end
for j=1:20  %对椒盐噪声图像进行20次均值滤波
    J2=imfilter(lena_pepper,h); end;
figure,subplot(1,3,1);
imshow(lena_pepper);
title('椒盐噪声图像');
subplot(1,3,2);
imshow(J1);
title('10次均值滤波');
subplot(1,3,3);
imshow(J2);
title('20次均值滤波');

（5) 对加入椒盐噪声的图像分别采用均值滤波法，和中值滤波法对有噪声的图像做处理，要求在同一窗口中显示结果。

close all
lena=imread('lena.jpg');    %读入原始图像
lena_gray=rgb2gray(lena);   %将图像灰度化
lena_pepper=imnoise(lena_gray,'salt & pepper',0.02);    %加入椒盐噪声
h=fspecial('average');  %生成均值滤波器
J1=imfilter(lena_pepper,h);
J2=medfilt2(lena_pepper);
figure,subplot(1,3,1);
imshow(lena_pepper);
title('椒盐噪声图像');
subplot(1,3,2);
imshow(J1);
title('均值滤波法');
subplot(1,3,3);
imshow(J2);
title('中值滤波法');

（6) 自己设计平滑空间滤波器，并将其对噪声图像进行处理，显示处理后的图像。

lena=imread('lena.jpg');    %读入原始图像
lena_gray=rgb2gray(lena);   %将图像灰度化
lena_pepper=imnoise(lena_gray,'salt & pepper',0.02);    %加入椒盐噪声
[m,n]=size(lena_pepper);
figure,subplot(1,2,1);
imshow(lena_pepper);
title('椒盐噪声图');
s=zeros(1,9);
for i=2:1:m-1
    for j=2:1:n-1
        h=1;
        for p=i-1:1:i+1
            for q=j-1:1:j+1
                s(h)=lena_pepper(p,q);
                h=h+1;
            end
        end
        s=sort(s);
       I(i,j)=s(3,3);
    end
end

subplot(1,2,2);
imshow(I);
title('噪声处理后的图像');

2. 锐化空域滤波

(1) 读出一幅灰度图像，采用的拉普拉斯算子对其进行滤波。

close all
lena=imread('lena.jpg');    %读入原始图像
lena_gray=rgb2gray(lena);   %将图像灰度化
lena_double=im2double(lena_gray) %将图像转换为双精度值
w=[1,1,1;
    1,-8,1;
    1,1,1]
k=conv2(lena_double,w,"same");
figure,subplot(1,2,1);      
imshow(lena);
title('原始图像');
subplot(1,2,2);      
imshow(k);
title('滤波处理后的图像');

(2) 写函数w = genlaplacian(n)，自动产生任一奇数尺寸的拉普拉斯算子，如的拉普拉斯算子

在这里插入图片描述

num=input('请输入数字n：');
n=num;
W=ones(n,n);
for i =1:n
    for j=1:n
        if(i==fix(n/2)+1&&j==fix(n/2)+1)
            W(i,j)=n*n-1;
        end
    end
end
display(W)

(3) 分别采用5×5 、9×9、15×15 和 25×25大小的拉普拉斯算子对图像进行锐化滤波，并利用式 .完成图像的锐化增强，观察其有何不同，要求在同一窗口中显示。

function [w]=lapulasi(num)
n=num;
w=ones(n);
x=fix(n/2)+1;
w(x,x)=-(n*n-1);
end

function [w]=lapulasi(num)
n=num;
w=ones(n);
x=fix(n/2)+1;
w(x,x)=-(n*n-1);
end
lena=imread("lena.jpg");    %读入原始图像
lena_double=im2double(lena);    %将图像转换为双精度值
figure,subplot(2,3,1);
imshow(lena_double);
title('原始图像');
w0=lapulasi(3);
w1=lapulasi(5);
w2=lapulasi(9);
w3=lapulasi(15);
w4=lapulasi(25);
lena_0=lena_double-imfilter(lena_double,w0,"replicate");
subplot(2,3,2);
imshow(lena_0);
title('3*3 拉普拉斯');
lena_1=lena_double-imfilter(lena_double,w1,"replicate");
subplot(2,3,3);
imshow(lena_1);
title('5*5 拉普拉斯');
lena_2=lena_double-imfilter(lena_double,w2,"replicate");
subplot(2,3,2);
imshow(lena_2);
title('9*9 拉普拉斯');
lena_3=lena_double-imfilter(lena_double,w3,"replicate");
subplot(2,3,5);
imshow(lena_3);
title('15*15 拉普拉斯');
lena_4=lena_double-imfilter(lena_double,w4,"replicate");
subplot(2,3,6);
imshow(lena_4);
title('25*25 拉普拉斯');

（4) 采用不同的梯度算子对图像进行锐化滤波，并比较其效果。

[I,map]=imread('lena_gray.jpg'); 
I=double(I); 
subplot(2,3,1) 
imshow(I,map); 
title('原始图像'); 
[Gx,Gy]=gradient(I); % 梯度计算 
G=sqrt(Gx.*Gx+Gy.*Gy); % 矩阵 
J1=G; % 梯度1 
subplot(2,3,2)
imshow(J1,map); 
title(' 操作1图像'); 
J2=I; % 梯度2 
K=find(G>=7); 
J2(K)=G(K); 
subplot(2,3,3) 
imshow(J2,map); 
title('操作2图像'); 
J3=I; % 梯度3 
K=find(G>=7); 
J3(K)=255; 
subplot(2,3,4) 
imshow(J3,map); 
title('操作3图像'); 
J4=I; % 梯度4 
K=find(G<=7); 
J4(K)=255; 
subplot(2,3,5) 
imshow(J4,map); 
title('操作4图像'); 
J5=I; %梯度5 
K=find(G<=7); 
J5(K)=0; 
Q=find(G>=7); 
J5(Q)=255; 
subplot(2,3,6) 
imshow(J5,map); 
title('操作5图像');

（5) 自己设计锐化空间滤波器，并将其对噪声图像进行处理，显示处理后的图像。

lena=imread("lena.jpg");
lena_gray=rgb2gray(lena);   %转换为灰度图像
h=fspecial("sobel");
h1=h'*0.5;
h2=h';
h3=h'*1.5;
z1=imfilter(lena_gray,h1);
z2=imfilter(lena_gray,h2);
z3=imfilter(lena_gray,h3);
figure,subplot(2,2,1);
imshow(lena_gray);
title('原始图像');
subplot(2,2,2);
imshow(z1);
title('梯度算子1');
subplot(2,2,3);
imshow(z2);
title('梯度算子2');
subplot(2,2,4);
imshow(z3);
title('梯度算子3');

3. 傅里叶变换

（1) 读出一幅灰度图像，对其进行快速傅里叶变换，分别显示其幅度图像和相位图像。

lena=imread("lena.jpg");
lena_gray=rgb2gray(lena);
f1=fft2(lena_gray); %快速傅里叶变换
f2=log(1+abs(f1));  %振幅谱
f3=fftshift(f1);
f4=angle(f1);   %相位谱
figure,subplot(1,3,1);
imshow(lena_gray);
title('原始图像');
subplot(1,3,2);
imshow(log(1+abs(f3)),[]);
title('幅度图像');
subplot(1,3,3);
imshow(f4);
title('相位图像');

（2) 仅对相位部分进行傅里叶逆变换后查看结果图像。

lena=imread("lena.jpg");
lena_gray=rgb2gray(lena);
f1=fft2(lena_gray); %快速傅里叶变换
f=ifft(abs(f1));
figure,subplot(1,3,1);
imshow(lena_gray);
title('原始图像');
subplot(1,3,2);
imshow(log(1+abs(f3)),[]);
title('振幅图像');
subplot(1,3,3);
imshow(log(1+abs(f)),[]);
title('相位图像');

（3) 仅对幅度部分进行傅里叶逆变换后查看结果图像。

lena=imread("lena.jpg");
lena_gray=rgb2gray(lena);
f1=fft2(lena_gray); %快速傅里叶变换
f2=log(1+abs(f1));
f3=fftshift(f1);
f4=angle(f1);
f5=-f4;
f6=double(f3*exp(f4));   %傅里叶变换的复共轭
f7=ifft2(f6);   %反傅里叶变换
figure,subplot(1,2,1);
imshow(lena_gray);
title('原始图像');
subplot(1,2,2);
imshow(real(f7),[]);
title('逆傅里叶变换');

（4) 将图像的傅里叶变换置为其共轭后进行逆变换，比较新生成图像与原始图像的差异。

lena=imread("lena.jpg");
lena_gray=rgb2gray(lena);
f1=fft2(lena_gray); %快速傅里叶变换
f2=log(1+abs(f1));
f3=fftshift(f1);
f4=angle(f1);
f5=f4;
f6=double(f3*exp(f4));  %傅里叶变换的复共轭
figure,subplot(1,2,1);
imshow(lena_gray);
title('原始图像');
subplot(1,2,2);
imshow(real(f7),[]);
title('逆傅里叶变换');

4. 平滑频域滤波

（1) 设计理想低通滤波器、巴特沃斯低通滤波器和高斯低通滤波器，截止频率自选。

%理想低通滤波器的透视图
a=100;
b=100;
U=0:a;
V=0:b;
M=length(U);
N=length(V);
D0=10;  %D0是频带的中心半径
x1=50;
y1=50;
x0=-50;
y0=-50;
m=fix(M/2);
n=fix(N/2);
H=zeros(M,N);
n=2;
for u=1:M
    for v=1:N
        a=sqrt((U(u)-50).*(U(u)-50)+(V(v)-50).*(V(v)-50));  %(u,v)的值
        if(a<=D0)   %理想滤波器
            H(u,v)=1;
        else
            H(u,v)=0;
        end
    end
end
figure,subplot(1,3,1);
surf(U,V,H);
title('理想低通滤波透视图');
%巴特沃斯低通滤波器透视图
a=100;    
b=100;    
U=0:a;    
V=0:b;    
M=length(U);N=length(V);    
D0=10;%W=200;%D0 是频带的中心半径;W 是频带的宽度
x1=50;y1=50;    
x0=-50;y0=-50;    
m=fix(M/2); 
n=fix(N/2);    
H=zeros(M,N);    
n=2;    
for u=1:M    
    for v=1:N    
        a=sqrt((U(u) - 50) .* (U(u)-50) + (V(v) - 50) .* (V(v) - 50));%D(u,v)的值
        b=1+(a/D0)^2*n;    
        H(u,v)=1/b;  
    end    
end    
subplot(1,3,2);
surf(U,V,H);  
title('巴特沃斯低通滤波器透视图');
%高斯低通滤波
a=100;    
b=100;    
U=0:a;    
V=0:b;    
M=length(U);
N=length(V);    
D0=10; %D0 是频带的中心半径
x1=50;
y1=50;    
x0=-50;
y0=-50;    
m=fix(M/2);
n=fix(N/2);    
H=zeros(M,N);          
for u=1:M    
    for v=1:N    
        D1=((u-m-x0)^2+(v-n-y0).^2)^0.5;    
        D2=((u-m+x0)^2+(v-n+y0).^2)^0.5;    
        D11=((u-m-x1)^2+(v-n-y1).^2)^0.5;    
        D21=((u-m+x1)^2+(v-n+y1).^2)^0.5;    
        H(u,v) = (U(u) - 50) .* (U(u)-50) + (V(v) - 50) .* (V(v)  - 50);   
    end    
end   
S=50;    
H = -H/(2*S);    
H = exp(H) / (sqrt(2*pi) * sqrt(S));    
subplot(1,3,3),surf(U,V,H),title('高斯低通滤波');

（2) 读出一幅灰度图像，分别采用理想低通滤波器、巴特沃斯低通滤波器和高斯低通滤波器对其进行滤波（截止频率自选），再做逆变换，观察不同的截止频率下采用不同低通滤波器得到的图像与原图像的区别，特别注意振铃效应。

%d0=15 理想低通滤波
lena=imread("lena.jpg");
lena_gray=rgb2gray(lena);
f=double(lena_gray);
g=fft2(f);  %傅里叶变换
g=fftshift(g);
[M,N]=size(g);
d0=15;
m=fix(M/2);
n=fix(N/2);
for i=1:M
    for j=1:N
        d=sqrt((i-m)^2+(j-n)^2);
        if(d<d0)
            h=1;
        else
            h=0;
        end
        result(i,j)=h*g(i,j);
    end
end
result=ifftshift(result);
J1=ifft2(result); 
J2=uint8(real(J1));  
figure,subplot(2,2,1);
imshow(lena);
title('原始图像'); 
subplot(2,2,2);
imshow(J2);
title('d0=15 理想低通滤波器');
%d0=30 的理想低通滤波
d0=30;  
m=fix(M/2); 
n=fix(N/2); 
for i=1:M  
    for j=1:N  
        d=sqrt((i-m)^2+(j-n)^2);  
        if(d<=d0)  
            h=1;  
        else   
            h=0;  
        end  
        result(i,j)=h*g(i,j);  
    end 
end  
result=ifftshift(result); 
J1=ifft2(result); 
J2=uint8(real(J1));  
subplot(2,2,3);
imshow(J2);
title('d0=30 理想低通滤波器');
%d0=100 的理想低通滤波
d0=100;  
m=fix(M/2); 
n=fix(N/2); 
for i=1:M  
    for j=1:N  
        d=sqrt((i-m)^2+(j-n)^2);  
        if(d<=d0)  
            h=1;  
        else   
            h=0;  
        end  
        result(i,j)=h*g(i,j);  
    end 
end  
result=ifftshift(result); 
J1=ifft2(result); 
J2=uint8(real(J1));  
 subplot(2,2,4);
 imshow(J2);
 title('d0=100 理想低通滤波器');

当截止频率 d0=15 时，滤波后的图像比较模糊，振铃现象也很明显；当 d0=30 时，图像模糊程度减弱，振铃现象仍存在。当d0=100 时，滤波后的图像比较清晰，但高频分量损失后，图像边沿仍然存在一点振铃现象。

5. 锐化频域滤波

（1) 设计理想高通滤波器、巴特沃斯高通滤波器和高斯高通滤波器，截止频率自选。

lena=imread("lena.jpg");
lena_gray=rgb2gray(lena);
f=double(lena_gray);
g=fft2(f);  %傅里叶变换
g=fftshift(g);
[M,N]=size(g);

%d0=15的巴特沃斯高通滤波器
nn=2;
d0=15;
m=fix(M/2);
n=fix(N/2); 
for i=1:M  
    for j=1:N  
        d=sqrt((i-m)^2+(j-n)^2);  
        h=1/(1+0.414*(d/d0)^(2*nn));  
        result(i,j)=h*g(i,j);  
    end 
end  
result=ifftshift(result); 
J1=ifft2(result); 
J2=uint8(real(J1));

figure,subplot(2,2,1);
imshow(lena);
title('原始图像'); 
subplot(2,2,2);
imshow(J2);
title('d0=15 巴特沃斯高通滤波器')
%d0=30的巴特沃斯高通滤波器
d0=30;  
m=fix(M/2);
n=fix(N/2); 
for i=1:M  
    for j=1:N  
        d=sqrt((i-m)^2+(j-n)^2);  
        h=1/(1+0.414*(d/d0)^(2*nn));  
        result(i,j)=h*g(i,j);  
    end 
end  
result=ifftshift(result); 
J1=ifft2(result);
J2=uint8(real(J1));   
subplot(2,2,3);
imshow(J2);
title('d0=30 巴特沃斯高通滤波器');
%d0=100的巴特沃斯高通滤波器
d0=100;  
m=fix(M/2);
n=fix(N/2); 
for i=1:M  
    for j=1:N  
        d=sqrt((i-m)^2+(j-n)^2);  
        h=1/(1+0.414*(d/d0)^(2*nn));  
        result(i,j)=h*g(i,j);  
    end 
end  
result=ifftshift(result); 
J1=ifft2(result);
J2=uint8(real(J1));   
subplot(2,2,4);
imshow(J2);
title('d0=100 巴特沃斯高通滤波器');
%d0=15的高斯低通滤波
d0=15;  
m=fix(M/2);
n=fix(N/2); 
for i=1:M  
    for j=1:N  
        d=sqrt((i-m)^2+(j-n)^2);  
        h=exp(-(d.^2)./(2*(d0^2)));  
        result(i,j)=h*g(i,j);  
    end 
end  
result=ifftshift(result); 
J1=ifft2(result);
J2=uint8(real(J1)); 
figure,subplot(2,2,1);
imshow(lena);
title('原始图像'); 
subplot(2,2,2);
imshow(J2);
title('d0=15 高斯高通滤');  
%d0=30的高斯低通滤波
d0=30;  
m=fix(M/2);
n=fix(N/2); 
for i=1:M  
    for j=1:N  
        d=sqrt((i-m)^2+(j-n)^2);  
        h=exp(-(d.^2)./(2*(d0^2)));  
        result(i,j)=h*g(i,j);  
    end 
end  
result=ifftshift(result);
J1=ifft2(result);
J2=uint8(real(J1)); 
subplot(2,2,3);
imshow(J2);
title('d0=30 高斯高通滤波');  
%d0=100的高斯低通滤波 
d0=100;  
m=fix(M/2);
n=fix(N/2); 
for i=1:M  
    for j=1:N  
        d=sqrt((i-m)^2+(j-n)^2);  
        h=exp(-(d.^2)./(2*(d0^2)));  
        result(i,j)=h*g(i,j);  
    end 
end  
result=ifftshift(result); 
J1=ifft2(result);
J2=uint8(real(J1));   
subplot(2,2,4);
imshow(J2);
title('d0=100 高斯高通滤波');

(2) 读出一幅灰度图像，分别采用理想高通滤波器、巴特沃斯高通滤波器和高斯高通滤波器对其进行滤波（截止频率自选），再做逆变换，观察不同的截止频率下采用不同高通滤波器得到的图像与原图像的区别。

close all
lena=imread('lena.jpg');    %将原始图像读入
lena_gray=rgb2gray(lena);   %将图像灰度化
lena_double=double(lena_gray);
g=fft2(f);  %傅里叶变换
g=fftshift(g);
[M,N]=size(g);
%d0=15的理想高通滤波器
d0=15;
m=fix(M/2);
n=fix(N/2); 
for i=1:M  
    for j=1:N  
        d=sqrt((i-m)^2+(j-n)^2);
        if(d>=d0)
            h=1;
        else
            h=0;
        end
        result(i,j)=h*g(i,j);
    end
end
result=ifftshift(result);
J1=ifft2(result);
J2=uint8(real(J1));
figure,subplot(2,2,1);
imshow(lena);
title('原始图像');
subplot(2,2,2);
imshow(J2);
title('d0=15 理想高通滤波器');
%d0=30的理想高通滤波器
d0=30;
m=fix(M/2);
n=fix(N/2); 
for i=1:M  
    for j=1:N  
        d=sqrt((i-m)^2+(j-n)^2);
        if(d>=d0)
            h=1;
        else
            h=0;
        end
        result(i,j)=h*g(i,j);
    end
end
result=ifftshift(result);
J1=ifft2(result);
J2=uint8(real(J1));
subplot(2,2,3);
imshow(J2);
title('d0=30 理想高通滤波器');
%d0=80的理想高通滤波器
d0=80;
m=fix(M/2);
n=fix(N/2); 
for i=1:M  
    for j=1:N  
        d=sqrt((i-m)^2+(j-n)^2);
        if(d>=d0)
            h=1;
        else
            h=0;
        end
        result(i,j)=h*g(i,j);
    end
end
result=ifftshift(result);
J1=ifft2(result);
J2=uint8(real(J1));
subplot(2,2,4);
imshow(J2);
title('d0=80 理想高通滤波器');
%d0=15的巴特沃斯高通滤波器
nn=2;
d0=15;
m=fix(M/2);
n=fix(N/2);
for i=1:M
    for j=1:N
        d=sqrt((i-m)^2+(j-n)^2);
        if(d==0)
            h=0;
        else
            h=1/(1+0.414*(d0/d)^(2*nn));
        end
        result(i,j)=h*g(i,j);
    end
end
result=ifftshift(result);
J1=ifft2(result);
J2=uint8(real(J1));
figure,subplot(2,2,1);
imshow(lena);
title('原始图像');
subplot(2,2,2);
imshow(J2);
title('d0=15 巴特沃斯高通滤波')
%d0=30的巴特沃斯高通滤波器
d0=30;
m=fix(M/2);
n=fix(N/2);
for i=1:M
    for j=1:N
        d=sqrt((i-m)^2+(j-n)^2);
        if(d==0)
            h=0;
        else
            h=1/(1+0.414*(d0/d)^(2*nn));
        end
        result(i,j)=h*g(i,j);
    end
end
result=ifftshift(result);
J1=ifft2(result);
J2=uint8(real(J1));

subplot(2,2,3);
imshow(J2);
title('d0=30 巴特沃斯高通滤波')
%d0=80的巴特沃斯高通滤波器
d0=80;
m=fix(M/2);
n=fix(N/2);
for i=1:M
    for j=1:N
        d=sqrt((i-m)^2+(j-n)^2);
        if(d==0)
            h=0;
        else
            h=1/(1+0.414*(d0/d)^(2*nn));
        end
        result(i,j)=h*g(i,j);
    end
end
result=ifftshift(result);
J1=ifft2(result);
J2=uint8(real(J1));
subplot(2,2,4);
imshow(J2);
title('d0=80 巴特沃斯高通滤波')