数字图像处理实验报告.doc

数字图像处理试验报告
实验二：数字图像的空间滤波和频域滤波
姓名： XX学号： 2XXXXXXX实验日期：2017年4月26日
1. 实验目的
1. 掌握图像滤波的基本定义及目的。

2. 理解空间域滤波的基本原理及方法。

3. 掌握进行图像的空域滤波的方法。

4. 掌握傅立叶变换及逆变换的基本原理方法。

5. 理解频域滤波的基本原理及方法。

6. 掌握进行图像的频域滤波的方法。

2. 实验内容与要求
1. 平滑空间滤波：
1) 读出一幅图像，给这幅图像分别加入椒盐噪声和高斯噪声后并与前一张图显示在同一
图像窗口中。

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

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

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

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

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

2.锐化空间滤波
1) 读出一幅图像，采用3×3 的拉普拉斯算子 w = [ 1, 1, 1; 1 – 81;1,1, 1]
对其进行滤波。

2) 编写函数 w = genlaplacian(n) ，自动产生任一奇数尺寸n 的拉普拉斯算子，如 5
×5的拉普拉斯算子
w = [ 1 1 1 1 1
1 1 1 1 1
1 1 -24 1 1
1 1 1 1 1
1 1 1 1 1]
3) 分别采用5×5，9×9，15×15 和 25×25 大小的拉普拉斯算子对blurry_moon.tif
进
行锐化滤波，并利用式 g(x, y)
2 f (x, y) 完成图像的锐化增强，观察其有何
f (x, y)
不同，要求在同一窗口中显示。

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

5)自己设计锐化空间滤波器，并将其对噪声图像进行处理，显示处理后的图像；
3.傅立叶变换
1)读出一幅图像，对其进行快速傅立叶变换，分别显示其幅度图像和相位图像。

仅对
相位部分进行傅立叶反变换后查看结果图像。

2)仅对幅度部分进行傅立叶反变换后查看结果图像。

3)将图像的傅立叶变换 F 置为其共轭后进行反变换，比较新生成图像与原始图像的差
异。

4.平滑频域滤波
1)设计理想低通滤波器、巴特沃斯低通滤波器和高斯低通滤波器，截至频率自选，分
别给出各种滤波器的透视图。

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

( 提示 :1) 在频率域滤波同样要注意到填充问题；2）注意到 (-1);)
5.锐化频域滤波
1)设计理想高通滤波器、巴特沃斯高通滤波器和高斯高通滤波器，截至频率自选，分
别给出各种滤波器的透视图。

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

3.实验具体实现
1.平滑空间滤波：
（1） . 读出一幅图像，给这幅图像分别加入椒盐噪声和高斯噪声后并与前一张图显示
在同一图像窗口中。

img=imread('lena.png')
figure,subplot(1,3,1); imshow(img);title('
原始图像 '); img2=imnoise(img,'salt &
pepper',0.02); subplot(1,3,2);
imshow(img2); title('椒盐噪声图像
');
img3=imnoise(img,'gaussian',0.02);
subplot(1,3,3),imshow(img3);
title('高斯噪声图像 ');
实验结果如下：
（2） . 对加入噪声图像选用不同的平滑（低通）模板做运算，对比不同模板所形成的效果，要求在同一窗口中显示。

平滑滤波是低频增强的空间域滤波技术。

它的目的有两个，一是模
糊，二是消除噪声。

将空间域低通滤波按线性和非线性特点有：线性、非线性平滑滤波器，
线性平滑滤波器包括均值滤波器，非线性的平滑滤波器有最大值滤波器，中值滤波器，最小值滤波器。

代码如下：
img=imread('lena.png') img=rgb2gray(img);
figure,subplot(1,3,1);
imshow(img);title('原始图像');
img2=imnoise(img,'salt & pepper',0.02);
subplot(1,3,2);imshow(img2);title(' 椒盐噪声图像 ');
img3=imnoise(img,'gaussian',0.02);
subplot(1,3,3),imshow(img3); title(' 高斯噪声图像
');
%对椒盐噪声图像进行滤波处理
h=fspecial('average',3);
I1=filter2(h,img2)/255;
I2=medfilt2(img2,[3 3]);
figure,subplot(2,2,1),imshow(img),title(' 原图像 ');
subplot(2,2,2),imshow(img2),title(' 椒盐噪声图 ');
subplot(2,2,3),imshow(I1),title('3*3 均值滤波图 ');
subplot(2,2,4),imshow(I2),title('3*3 中值滤波图 ');
%对高斯噪声图像进行滤波处理
G1=filter2(h,img3)/255;
G2=medfilt2(img3,[3 3]);
figure,subplot(2,2,1),imshow(img),title(' 原图像 ');
subplot(2,2,2),imshow(img3),title(' 高斯噪声图 ');
subplot(2,2,3),imshow(G1),title('3*3 均值滤波图 ');
subplot(2,2,4),imshow(G2),title('3*3 中值滤波图 ');
实验结果如下：
(3).使用函数充、’ replicate imfilter ’、’时，分别采用不同
symmetric ’、’c ircular
的填充方法（或边界选项，如零填
’）进行低通滤波，显示处理后的图像。

g = imfilter(f, w, filtering_mode, boundary_options, size_options) ，其中， f 为输入图像， w 为滤波掩模， g 为滤波后图像。

h=fspecial('motion',50,45); %创建一个运动模糊滤波器
filteredimg=imfilter(img,h);
boundaryReplicate=imfilter(img,h,'replicate');
boundary0=imfilter(img,h,0);
boundarysymmetric=imfilter(img,h,'symmetric');
boundarycircular=imfilter(img,h,'circular');
figure,subplot(3,2,1),imshow(img),title('Original Image');
subplot(3,2,2),imshow(filteredimg),title('Motion Blurred Image');
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('circular');
实验结果如下：
（4） . 运用 for循环，将加有椒盐噪声的图像进行10 次， 20 次均值滤波，查看其特点
示均值处理后的图像（提示: 利用 fspecial函数的’ ave rage’类型生成均值滤波器）。

, 显
代码如下：
h=fspecial('average');
J1=imfilter(img2,h);
end
for j=1:20
J2=imfilter(img2,h);
end
figure,subplot(1,3,1),imshow(img2),title('salt & pepper Noise');
subplot(1,3,2),imshow(J1),title('10 Average Filtering');
subplot(1,3,3),imshow(J2),title('20 Average Filtering');
实验结果 :
（5） . 对加入椒盐噪声的图像分别采用均值滤波法，和中值滤波法对有噪声的图像做处理，要求在同一窗口中显示结果。

代码如下：
h1=fspecial('average');
J=imfilter(img2,h1);
J2=medfilt2(img2);
figure,subplot(1,3,1),imshow(img2),title('salt & pepper Noise');
subplot(1,3,2),imshow(J),title('Averaging Filtering');
subplot(1,3,3),imshow(J2),title('Median Filtering');
实验结果为：
（6） . 自己设计平滑空间滤波器，并将其对噪声图像进行处理，显示处理后的图像。

代码如下：
[m n]=size(img2);
figure,subplot(1,2,1),imshow(img2);
s=zeros(1,9);
for j=2:1:n-1
h=1;
for p=i-1:1:i+1
for q=j-1:1:j+1
s(h)=img2(p,q);
h=h+1;
end
end
s=sort(s);
I(i,j)=s(5);
end
end
subplot(1,2,2),imshow(I);
实验结果：
2. 锐化空间滤波
（1）读出一幅图像，采用3×3的拉普拉斯算子w = [ 1, 1, 1; 1 –81;1,1,1] 对其进行滤波。

代码如下：
img=imread('lena.png');
img=rgb2gray(img);
img=im2double(img);
w=[1,1,1;
1,-8,1;
1,1,1]
k=conv2(img,w,'same');
imshow(k);
实验结果为：
（2）编写函数 w = genlaplacian(n) ，自动产生任一奇数尺寸n 的拉普拉斯算子，如 5×5 的拉普拉斯算子
w = [ 1 1 1 1 1
1 1 1 1 1
1 1 -24 1 1
1 1 1 1 1
1 1 1 1 1]
代码如下：
num=input('please enter any num:');
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 大小的拉普拉斯算子对blurry_moon.tif
进行锐化滤波，并利用式完成图像的锐化增强，观察其有何不同，要求在同一窗口中显示。

代码如下：基于上一题要求形
成一函数：
function [W]=lapulasi(num)
n=num,W=ones(n),x=fix(n/2)+1;
W(x,x)=-(n*n-1);
其他代码：
f=imread('moon.tif');
f=im2double(f);
figure,subplot(2,3,1),imshow(f),title('Original Image'); w0=lapulasi(3),w1=lapulasi(5),w2=lapulasi(9);
w3=lapulasi(15),w4=lapulasi(25);
f0=f-imfilter(f,w0,'replicate');
subplot(2,3,2),imshow(f0),title('3*3 lapulasi');
f1=f-imfilter(f,w1,'replicate');
subplot(2,3,3),imshow(f1),title('5*5 lapulasi');
f2=f-imfilter(f,w2,'replicate');
subplot(2,3,4),imshow(f2),title('9*9 lapulasi');
f3=f-imfilter(f,w3,'replicate');
subplot(2,3,5),imshow(f3),title('15*15 lapulasi');
f4=f-imfilter(f,w4,'replicate');
subplot(2,3,6),imshow(f4),title('25*25 lapulasi');
实验结果为：
（4）采用不同的梯度算子对该幅图像进行锐化滤波，并比较其效果。

代码如下：
[I,map]=imread('moon.tif');
I=double(I);
figure,subplot(2,3,1),imshow(I,map),title('Original Image');
[Gx,Gy]=gradient(I);
G=sqrt(Gx.*Gx+Gy.*Gy),J1=G;
subplot(2,3,2),imshow(J1,map),title('Operator1 Image');
J2=I;K=find(G>=7);
J2(K)=G(K);
subplot(2,3,3),imshow(J2,map);
title('Operator2 Image');
J3=I;
K=find(G>=7);
J3(K)=255;
subplot(2,3,4),imshow(J3,map),title('Operator3 Image');
J4=I;
K=find(G<=7);
J4(K)=255;
subplot(2,3,5),imshow(J4,map),title('Operator4 Image');
J5=I;
K=find(G<=7);
J5(K)=0;
Q=find(G>=7);
J5(Q)=255;
subplot(2,3,6),imshow(J5,map),title('Operator5 Image');
实验效果如下：
（5）自己设计锐化空间滤波器，并将其对噪声图像进行处理，显示处理后的图像；
代码如下：
I=imread('lena.png');
I=rgb2gray(I);
h=fspecial('sobel');
h1=h'*0.5;
h2=h';
h3=h'*1.5;
z1=imfilter(I,h1);
z2=imfilter(I,h2);
z3=imfilter(I,h3);
figure,subplot(2,2,1),imshow(I),title('Original Image');
subplot(2,2,2),imshow(z1);title('Vertical filtering1');
subplot(2,2,3),imshow(z2);title('Vertical filtering2');
subplot(2,2,4),imshow(z3);title('Vertical filtering3')
运行结果如下：
3.傅立叶变换
（1） . 读出一幅图像，对其进行快速傅立叶变换，分别显示其幅度图像和相位图像。

仅对相位部分进行傅立叶反变换后查看结果图像。

代码为：
img=imread('lena.png');
img=rgb2gray(img);
f1=fft2(img); %快速傅里叶变换
f2=log(1+abs(f1)); %振幅谱
f3=fftshift(f1);
f4=angle(f1); %相位谱
figure,subplot(1,3,1),imshow(img),title('Original Image');
subplot(1,3,2),imshow(log(1+abs(f3)),[]),title('amplitude spectrum');
subplot(1,3,3),imshow(f4),title('phase spectrum');
实验结果为：
代码为：
f=ifft2(abs(f1));
figure,subplot(1,3,1),imshow(img),title('Original Image');
subplot(1,3,2),imshow(log(1+abs(f3)),[]),title('amplitude spectrum');
subplot(1,3,3),imshow(log(1+abs(f)),[]),title('absamplitude spectrum');
（2） . 仅对幅度部分进行傅立叶反变换后查看结果图
像。

实验结果为：
（3） . 将图像的傅立叶变换 F 置为其共轭后进行反变换，比较新生成图像与原始图像的差异。

代码
f1=fft2(img);
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(img),title('Original Image');
subplot(1,2,2),imshow(real(f7),[]),title('inverse fourier transform');
实验效果：
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));%D(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('理想低通滤波透视图') ;
%%%%%%%%%%%%2阶巴特沃斯低通滤波透视图%%%%%%%%%%%%%%%%55 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('n=2 Butterworth lowess filter')
%%%%%%%%%%%%%%%%%%%%%%高斯低通滤波%%%%%%%%%%%%%%%%%%%%% 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('Gaussian lowess filter');
实验结果：
由上图可知，当频带中心宽度相同时，理想低通滤波器为圆柱形图像，二阶巴特沃斯
低通滤波器的面线比较紧凑，高斯滤波图像最为平滑。

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

( 提示 :1) 在频率域滤波同样要注意到填充问题； 2）注意到 (-1)x+y;)
%理想低通滤波
img=imread('lena.png');
img=rgb2gray(img);
f=double(img);
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(img),title('Original Image');
subplot(2,2,2),imshow(J2),title('d0=15 lowpss filter');
%%%%%%%%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 lowpss filter');
%%%%%%%%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 lowpss filter');
实验结果：
当截止频率d0=15 时，滤波后的图像比较模糊，振铃现象也很明显；当d0=30 时，图像模糊程度减弱，振铃现象仍存在。

当 d0=100 时，滤波后的图像比较清晰，但高频分量损失后，图
像边沿仍然存在一点振铃现象。

2阶巴特沃斯低通滤波
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(img),title('Original Image');
subplot(2,2,2),imshow(J2),title('d0=15 Butterworth lowpss filter');
%%%%%%%%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 Butterworth lowpss filter');
%%%%%%%%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 Butterworth lowpss filter');
实验效果为：
采用截止频率分别为15 ， 30,100 进行实验，滤波后的图像逐渐清晰，整体的振铃现象
没有理想低通滤波器的严重。

%高斯低通滤波
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(img),title('Original Image');
subplot(2,2,2),imshow(J2),title('d0=15 Gaussian filter');
%%%%%%%%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
Gaussian filter');
%%%%%%%%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 Gaussian filter');
实验结果为：
从实验对比
可知，高斯滤波器的平滑效果不如二阶巴特沃斯滤波器，但用高斯滤波后的图像无振铃现象产生。

5.锐化频域滤波
（1）设计理想高通滤波器、巴特沃斯高通滤波器和高斯高通滤波器，截至频率自选，分别
给出各种滤波器的透视图。

%%%%%%%%%%%%%%%%理想高通滤波器%%%%%%%%%%%%%%%%%%%%%%%%%%% a=100;
b=100;
U=0:a;
V=0:b;
M=length(U);N=length(V);
D0=15; %D0 是频带的中心半径;
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)的值if(a>=D0)
H(u,v)=1;
else
H(u,v)=0;
end
end
end
figure;
surf(U,V,H),title('理想高通滤波透视图') ;
高通滤波器
%%%%%%%%%%%%%%%巴特沃斯高通滤波器%%%%%%%%%%%%%%%%%%%%%%%% a=100;
b=100;
U=0:a;
V=0:b;
M=length(U);N=length(V);
D0=15; %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));%D(u,v)的值
b=1+(a/D0)^2*n;
H(u,v)=-1/b;
end
end
figure,surf(U,V,H),title('n=2 Butterworth high filter')
巴特沃斯高通滤波器
%%%%%%%%%%%%%%高斯高通滤波
a=100;
b=100;
U=0:a;
V=0:b;
M=length(U);N=length(V);
D0=15; %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));
figure,surf(U,V,H),title('Gaussian high filter');
高斯高通滤波器
（2）读出一幅图像，分别采用理想高通滤波器、巴特沃斯高通滤波器和高斯高通滤波器对其进行滤波（截至频率自选），再做反变换，观察不同的截止频率下采用不同高通滤波器得到的图像与原图像的区别。

%理想高通滤波
img=imread('lena.png');
img=rgb2gray(img);
f=double(img);
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(img),title('Original Image');
subplot(2,2,2),imshow(J2),title('d0=15 high filter');
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 high filter');
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 high filter');
理想高通滤波器得到的图片效果为：
当 d0=15 时，滤波后的图像灰度变化部分基本保留；当 d0=30 时，图像的轮廓还比较清楚；
当 d0=80 时，只能看到细微的图像轮廓；
%巴特沃斯高通滤波
img=imread('lena.png');
img=rgb2gray(img);
f=double(img);
g=fft2(f); %傅立叶变换
g=fftshift(g);
[M,N]=size(g);
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(img),title('Original Image');
subplot(2,2,2),imshow(J2),title('d0=15 Butterworth high filter');
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 Butterworth high filter');
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 Butterworth high filter');
n=2 时的巴特沃斯高通滤波效果为：
滤波结果比理想高通滤波器更加平滑，边缘失真情况较小。

%高斯高通滤波器
img=imread('lena.png');
img=rgb2gray(img);
f=double(img);
g=fft2(f); %傅立叶变换
g=fftshift(g);
[M,N]=size(g);
n=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-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(img),title('Original Image');
subplot(2,2,2),imshow(J2),title('d0=15 Gaussian filter');
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-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 Gaussian filter');
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);
h=1-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=80 Gaussian filter');
高斯高通滤波效果为：
高斯高通的滤波效果在这 3 种滤波器中最为平滑，过滤效果比较清晰。