18组实验四王栩灿肖凌峰张皓.docx
《18组实验四王栩灿肖凌峰张皓.docx》由会员分享，可在线阅读，更多相关《18组实验四王栩灿肖凌峰张皓.docx（29页珍藏版）》请在冰点文库上搜索。
18组实验四王栩灿肖凌峰张皓
实验四傅里叶变换及频率滤波
时间：
2017.11.3组员：
王栩灿/肖凌峰/张皓第18组
一、实验目的
1.了解傅里叶变换在图像处理中的应用；
2.掌握频率滤波的原理和特点；
3.掌握傅里叶变换的基本性质；
4.通过实验了解二维频谱的分布特点。
二、实验内容
1.构造同时包含50Hz和120Hz的正弦信号进行傅里叶变换并分析。
2.以高斯光束为例（自己选的例子亦可），利用matlab实现二维图像的傅里叶变换，选取不同的抽样间隔，观察抽样间隔对结果的影响并分析。
3.对某一图像（如lena.jpg）
（1）编程实现图像的理想低通、高斯低通、布特沃斯低通滤波，分析不同的截止频率对滤波后图像的影响。
（2）编程实现图像的理想高通、高斯高通、布特沃斯高通滤波，分析不同的截止频率对滤波后图像的影响。
（3）在空域中实现高斯低通滤波，并与频域的结果进行比较。
4.利用FFT实现二维高斯光束的菲涅尔衍射和夫琅禾费衍射，并分析光束在不同传输距离上的衍射图样（传输距离自定，至少三个数值说明即可）。
利用一幅图像进行菲涅尔衍射和夫琅禾费衍射后会是什么样子？
三、实验原理
①快速傅里叶变换
有限长序列可以通过离散傅里叶变换（DFT）将其频域也离散化成有限长序列。
但其计算量太大，很难实时地处理问题，因此引出了快速傅里叶变换（FFT）.1965年，Cooley和Tukey提出了计算离散傅里叶变换（DFT）的快速算法，将DFT的运算量减少了几个数量级。
从此，对快速傅里叶变换（FFT）算法的研究便不断深入，数字信号处理这门新兴学科也随FFT的出现和发展而迅速发展。
根据对序列分解与选取方法的不同而产生了FFT的多种算法，基本算法是基2DIT和基2DIF。
FFT在离散傅里叶反变换、线性卷积和线性相关等方面也有重要应用。
快速傅立叶变换（FFT），是离散傅氏变换的快速算法，它是根据离散傅氏变换的奇、偶、虚、实等特性，对离散傅立叶变换的算法进行改进获得的。
FFT的基本思想是把原始的N点序列，依次分解成一系列的短序列。
充分利用DFT计算式中指数因子所具有的对称性质和周期性质，进而求出这些短序列相应的DFT并进行适当组合，达到删除重复计算，减少乘法运算和简化结构的目的。
此后，在这思想基础上又开发了高基和分裂基等快速算法，随着数字技术的高速发展，1976年出现建立在数论和多项式理论基础上的维诺格勒傅里叶变换算法（WFTA）和素因子傅里叶变换算法。
它们的共同特点是，当N是素数时，可以将DFT算转化为求循环卷积，从而更进一步减少乘法次数，提高速度。
②理想低通/高通滤波器
是一种假想的低通/高通滤波器，其对于高于/低于截止频率的信号完全截止，而对于低于/高于截止频率的信号完全无失真传输。
③布特沃斯低通/高通滤波器
截止频率位于距原点D0的n阶布特沃斯滤波器（BLPF）的传递函数定义为
与IFPF不同，BLPF并没有字通过频率和滤除频率之间给出明显截止的急剧不连续性。
且空间域的一阶布特沃斯滤波器没有振铃现象。
④高斯低通/高通滤波器
高斯滤波器的二维形式由下式给出：
⑤Fresnel衍射与Fraunhofer衍射
在光学里，菲涅耳衍射（Fresneldiffraction）指的是光波在近场区域的衍射，即光源或衍射的图样的屏与衍射孔（障碍物）的距离是有限的。
菲涅耳衍射积分式可以用来计算光波在近场区域的传播，因法国物理学者奥古斯丁·菲涅耳而命名，是基尔霍夫衍射公式的近似。
夫琅禾费衍射是指把单色点光源放在透镜的焦点上，经过透镜后的单色平行光垂直照射衍射屏时，在屏后面不同距离上会观察到一些衍射现象，其中当屏远离到足够大的距离后，光斑中心出现一个较大的亮斑，外围是一些较弱的明暗相间的同心圆环，此后再往外移动，衍射花样出现稳定分布，中心处总是亮的，只是半径不断扩大而已，这种衍射称为夫琅禾费衍射，又称远场衍射（而当距离较近时为菲涅尔衍射，较近时没有夫琅禾费衍射）。
4、实验程序与结果
构造同时包含50Hz和120Hz的正弦信号进行傅里叶变换并分析。
主程序
f1=fun_GenSinWave（1000,1000,50）;
%plot（f1）;
f2=fun_GenSinWave（1000,1000,120）;
%figure;plot（f2）;
f=f1+f2;
y=fft（f）;
figure;plot（y）;
正弦波生成函数
function[output]=fun_GenSinWave（fs,N,f）
%partone,generate50Hzand100HzsinsignalanddoFourierchange
fsm=fs;%采样频率
Ns=N;
n=0:
N-1;
t=n/fs;
f1=f;
x=sin（2*pi*f1*t）;
可以看到不同的频率在做傅里叶变换后被区分开来。
②以高斯光束为例（自己选的例子亦可），利用matlab实现二维图像的傅里叶变换，选取不同的抽样间隔，观察抽样间隔对结果的影响并分析。
采用lena图像进行二维傅里叶变换
主程序：
%二维傅里叶变换
f=imread（'lena512color.tiff'）;%以lena图为例
f=rgb2gray（f）;%转为灰度图
f=double（f）;%转为double值
F1=fft2（f）;%使用MATLAB库傅里叶变换函数
F1=abs（F1）;
F1=fftshift（F1）;%浮点移动到中央位置
F1=log10（F1）;
F1=uint8（F1）;%取对数，恢复至无符号八位整形
imshow（F1,[]）%显示图像。
imwrite（F1,'lenaDFT1.jpg'）;
F2=downsample（f,2）;%使用downsample函数进行降采样
F2=fft2（f）;%使用MATLAB库傅里叶变换函数
F2=abs（F2）;
F2=fftshift（F2）;%浮点移动到中央位置
F2=log10（F2）;
F2=uint8（F2）;%取对数，恢复至无符号八位整形
figure,imshow（F2,[]）%显示图像。
imwrite（F2,'lenaDFT2.jpg'）;
原图像
二维傅里叶变换后图像二维傅里叶变换后图像（降采样）
采用了matlab库函数中的快速傅里叶变换函数（FFT）进行变换，虽然从结果图上不易看到降采样产生的影响，但是在查看图像的矩阵中数值时可以观察到像元的值是不一样的。
③对某一图像（如lena.jpg）
（1）编程实现图像的理想低通、高斯低通、布特沃斯低通滤波，分析不同的截止频率对滤波后图像的影响。
（2）编程实现图像的理想高通、高斯高通、布特沃斯高通滤波，分析不同的截止频率对滤波后图像的影响。
主程序：
clearall
clc
f=imread（'图3-5lena.jpg'）;
R1=30;
R2=90;
R3=160;
f=double（f）;
F=fft2（f）;
F=fftshift（F）;%原图像的傅里叶逆变换
%------------------------理想低通滤波--------------------------------------
Lowpass_Filters1=fun_Lowpass_Filters（F,R1）;
imwrite（Lowpass_Filters1,'R=30的理想低通滤波.bmp'）;
Lowpass_Filters2=fun_Lowpass_Filters（F,R2）;
imwrite（Lowpass_Filters2,'R=90的理想低通滤波.bmp'）;
Lowpass_Filters3=fun_Lowpass_Filters（F,R3）;
imwrite（Lowpass_Filters3,'R=160的理想低通滤波.bmp'）;
%------------------------理想高通滤波--------------------------------------
Highpass_Filters1=fun_Highpass_Filters（F,R1）;
imwrite（Highpass_Filters1,'R=30的理想高通滤波.bmp'）;
Highpass_Filters2=fun_Highpass_Filters（F,R2）;
imwrite（Highpass_Filters2,'R=90的理想高通滤波.bmp'）;
Highpass_Filters3=fun_Highpass_Filters（F,R3）;
imwrite（Highpass_Filters3,'R=160的理想高通滤波.bmp'）;
%------------------------高斯低通滤波--------------------------------------
Gaussian_Lowpass_Filters1=fun_Gaussian_Lowpass_Filters（F,R1）;
imwrite（Gaussian_Lowpass_Filters1,'R=30的高斯低通滤波.bmp'）;
Gaussian_Lowpass_Filters2=fun_Gaussian_Lowpass_Filters（F,R2）;
imwrite（Gaussian_Lowpass_Filters2,'R=90的高斯低通滤波.bmp'）;
Gaussian_Lowpass_Filters3=fun_Gaussian_Lowpass_Filters（F,R3）;
imwrite（Gaussian_Lowpass_Filters3,'R=160的高斯低通滤波.bmp'）;
%------------------------高斯高通滤波--------------------------------------
Gaussian_Highpass_Filters1=fun_Gaussian_Highpass_Filters（F,R1）;
imwrite（Gaussian_Highpass_Filters1,'R=30的高斯高通滤波.bmp'）;
Gaussian_Highpass_Filters2=fun_Gaussian_Highpass_Filters（F,R2）;
imwrite（Gaussian_Highpass_Filters2,'R=90的高斯高通滤波.bmp'）;
Gaussian_Highpass_Filters3=fun_Gaussian_Highpass_Filters（F,R3）;
imwrite（Gaussian_Highpass_Filters3,'R=160的高斯高通滤波.bmp'）;
%------------------------布特沃斯低通滤波-----------------------------------
Butterworth_Lowpass_Filters1=fun_Butterworth_Lowpass_Filters（F,R1）;
imwrite（Butterworth_Lowpass_Filters1,'R=30的布特沃斯低通滤波.bmp'）;
Butterworth_Lowpass_Filters2=fun_Butterworth_Lowpass_Filters（F,R2）;
imwrite（Butterworth_Lowpass_Filters1,'R=90的布特沃斯低通滤波.bmp'）;
Butterworth_Lowpass_Filters3=fun_Butterworth_Lowpass_Filters（F,R3）;
imwrite（Butterworth_Lowpass_Filters1,'R=160的布特沃斯低通滤波.bmp'）;
%------------------------布特沃斯高通滤波-----------------------------------
Butterworth_Highpass_Filters1=fun_Butterworth_Highpass_Filters（F,R1）;
imwrite（Butterworth_Highpass_Filters1,'R=30的布特沃斯高通滤波.bmp'）;
Butterworth_Highpass_Filters2=fun_Butterworth_Highpass_Filters（F,R2）;
imwrite（Butterworth_Highpass_Filters2,'R=90的布特沃斯高通滤波.bmp'）;
Butterworth_Highpass_Filters3=fun_Butterworth_Highpass_Filters（F,R3）;
imwrite（Butterworth_Highpass_Filters3,'R=160的布特沃斯高通滤波.bmp'）;
使用的自定义函数
理想低通/高通
function[output]=fun_Highpass_Filters（H_Image,R）
%------------------------变量------------------------------
%Image读取的图像
%R理想低通滤波的半径
Low=H_Image;
r=R;
[M,N]=size（Low）;
h=zeros（M,N）;
m=round（M/2）;
n=round（N/2）;
fori=1:
M
forj=1:
N
S=sqrt（（i-m）^2+（j-n）^2）;
ifS<=r
h（i,j）=1;
else
h（i,j）=0;
end
end
end
Lowpass_Filters=h.*Low;%理想低通滤波
Lowpass_Filters=uint8（real（ifft2（ifftshift（Lowpass_Filters））））;
figure（'Name','理想低通滤波'）;
imshow（Lowpass_Filters）;
output=Lowpass_Filters;
function[output]=fun_Highpass_Filters（H_Image,R）
%------------------------变量------------------------------
%Image读取的图像
%R理想高通滤波的半径
High=H_Image;
r=R;
[M,N]=size（High）;
h=zeros（M,N）;
m=round（M/2）;
n=round（N/2）;
fori=1:
M
forj=1:
N
S=sqrt（（i-m）^2+（j-n）^2）;
ifS<=r
h（i,j）=0;
else
h（i,j）=1;
end
end
end
Highpass_Filters=h.*High;%理想低通滤波
Highpass_Filters=uint8（real（ifft2（ifftshift（Highpass_Filters））））;
figure（'Name','理想高通滤波'）;
imshow（Highpass_Filters）;
output=Highpass_Filters;
高斯低通/高通
function[output]=fun_Gaussian_Lowpass_Filters（G_Image,R）
%------------------------------变量------------------------------
%H_Image经过傅里叶变换处理的图像
%R高斯低通滤波的半径
Gaussian_Lowpass=G_Image;
[M,N]=size（Gaussian_Lowpass）;
m1=fix（M/2）;
n1=fix（N/2）;
r=R;
fori=1:
M
forj=1:
N
d=sqrt（（i-m1）^2+（j-n1）^2）;
G_H（i,j）=exp（-d^2/2/r^2）;
end
end
Gaussian_Lowpass_Filters=Gaussian_Lowpass.*G_H;
Gaussian_Lowpass_Filters=ifftshift（Gaussian_Lowpass_Filters）;
Gaussian_Lowpass_Filters=ifft2（Gaussian_Lowpass_Filters）;
Gaussian_Lowpass_Filters=mat2gray（real（Gaussian_Lowpass_Filters））;
figure（'Name','高斯低通滤波'）;
imshow（Gaussian_Lowpass_Filters）;
output=Gaussian_Lowpass_Filters;
function[output]=fun_Gaussian_Highpass_Filters（G_Image,R）
%------------------------变量------------------------------
%H_Image经过傅里叶变换处理的图像
%R高斯高通滤波的半径
Gaussian_Highpass=G_Image;
[M,N]=size（Gaussian_Highpass）;
m1=fix（M/2）;
n1=fix（N/2）;
r=R;
fori=1:
M
forj=1:
N
d=sqrt（（i-m1）^2+（j-n1）^2）;
G_Hh（i,j）=1-exp（-d^2/2/r^2）;
end
end
Gaussian_Highpass_Filters=Gaussian_Highpass.*G_Hh;
Gaussian_Highpass_Filters=ifftshift（Gaussian_Highpass_Filters）;
Gaussian_Highpass_Filters=ifft2（Gaussian_Highpass_Filters）;
Gaussian_Highpass_Filters=mat2gray（real（Gaussian_Highpass_Filters））;
figure（'Name','高斯高通滤波'）;
imshow（Gaussian_Highpass_Filters）;
output=Gaussian_Highpass_Filters;
布特沃斯低通/高通
function[output]=fun_Butterworth_Lowpass_Filters（B_Image,R）
%------------------------变量------------------------------
%B_Image经过傅里叶变换处理的图像
%R布特沃斯低通滤波的半径
Butterworth_Lowpass=B_Image;
[M,N]=size（Butterworth_Lowpass）;
r=R;
m2=fix（M/2）;
n2=fix（N/2）;
fori=1:
M
forj=1:
N
d=sqrt（（i-m2）^2+（j-n2）^2）;
B_Hl=1/（1+（d/r）^（2*2））;
resultBl（i,j）=B_Hl*Butterworth_Lowpass（i,j）;
end
end
Butterworth_Lowpass_Filters=ifftshift（resultBl）;
Butterworth_Lowpass_Filters=ifft2（Butterworth_Lowpass_Filters）;
Butterworth_Lowpass_Filters=uint8（real（Butterworth_Lowpass_Filters））;
figure（'Name','布特沃斯低通滤波'）;
imshow（Butterworth_Lowpass_Filters）;
output=Butterworth_Lowpass_Filters;
function[output]=fun_Butterworth_Highpass_Filters（B_Image,R）
%------------------------变量------------------------------
%B_Image经过傅里叶变换处理的图像
%R布特沃斯高通滤波的半径
Butterworth_Highpass=B_Image;
[M,N]=size（Butterworth_Highpass）;
r=R;
m2=fix（M/2）;
n2=fix（N/2）;
fori=1:
M
forj=1:
N
d=sqrt（（i-m2）^2+（j-n2）^2）;
B_HH=1/（1+（r/d）^（2*2））;
resultBH（i,j）=B_HH*Butterworth_Highpass（i,j）;
end
end
Butterworth_Highpass_Filters=ifftshift（resultBH）;
Butterworth_Highpass_Filters=ifft2（Butterworth_Highpass_Filters）;
Butterworth_Highpass_Filters=uint8（real（Butterworth_Highpass_Filters））;
figure（'Name','布特沃斯高通滤波'）;
imshow（Butterworth_Highpass_Filters）;
output=Butterworth_Highpass_Filters;
结果图：
理想低通滤波
半径R
30
90
160
图片
理想高通滤波
半径R
30
90
160
图片
布特沃斯低通滤波
半径R
30
90
160
图片
布特沃斯高通滤波
半径R
30
90
160
图片
高斯低通滤波
半径R
30
90
160
图片
高斯高通滤波
半径R
30
90
160
图片
（3）在空域中实现高斯低通滤波，并与频域的结果进行比较。
空域高斯滤波（R=90）
比较可以得出，相对于频率域高斯滤波，空域高斯滤波得到的图像灰度值较高一些。
4.利用FFT实现二维高斯光束的菲涅尔衍射和夫琅禾费衍射，并分析光束在不同传输距离上的衍射图样（传输距离自定，至少三个数值说明即可）。
实验主程序
由于版本不同及兼容问题汉字部分会变成乱码，源程序可看。
%===========·ÆÄù¶ûÑÜÉäÓë·òÀź̷ÑÑÜÉä========
%
%
%
%
%========================================
clearall;
clc;
升级会员 