DIP实验报告DFT.docx
《DIP实验报告DFT.docx》由会员分享，可在线阅读，更多相关《DIP实验报告DFT.docx（15页珍藏版）》请在冰豆网上搜索。
DIP实验报告DFT
课程：
数字图像处理
课程作业实验报告
实验名称：
Two-DimensionalFourierTransformandFiltering
实验编号：
签名：
姓名：
学号：
截止提交日期：
年月日
摘要：
本次实验主要分为两个部分，第一部分：
编写一个2-DFFT程序包，实现滤波的中心变换、FFT变换、IFFT变换以及二维和三维显示幅度谱和相位谱；第二部分：
编写一个高斯低通滤波程序，实现图像的的滤波处理。
用原图减去低通滤波处理后的图像，与高通处理后的图像做比较，改变低通滤波器的参数，使其接近高通滤波处理的图像。
1、技术论述
1、二维离散傅里叶变换
设
表示一副大小为MxN的图像，其中x=0，1，2，…，M-1和y=0，1，2，…，N-1。
的二维离散傅里叶变换可表示为
，如下式所示：
（1.1）
其中u=0，1，2，…，M-1和v=0，1，2，…，N-1。
离散傅里叶逆变换由下式给出：
（1.2）
其中x=0，1，2，…，M-1和y=0，1，2，…，N-1。
令
和
分别表示
的实部和虚部，则傅里叶频谱定义为：
（1.3）
变换的相位角定义为：
（1.4）
使用我们熟悉的复数量的极坐标表示法，
表示如下：
（1.5）
在实际应用中，DFT及其逆变换可以通过使用快速傅里叶变换（FFT）算法来实现。
一个大小为MxN的图显数组f可以通过MATLAB工具箱中的函数fft2得到，函数fft2的简单语法为F=fft2（f）。
该函数返回一个大小仍为MxN的傅里叶变换，数据的原点在矩形的左上角，而四个四分之一周期交汇于频率矩形的中心。
在二维情况下，为了简化可视化程度，需要将F（0,0）处移动到（M/2，N/2）处。
利用傅里叶变换的平移性质，如下所示：
（1.6）
在式（1.6）中令
，可得如下表达式：
（1.7）
在本次实验中编写了一个中心变换的程序，在进行傅里叶变换前，需要将时域的信号乘以
，然后再进行FFT处理，此过程不可以倒置。
2、高斯低通滤波处理
二维高斯低通滤波器可由下式给出：
（2.1）
其中
是距频率矩形中心的距离，
是关于中心的扩展度的度量。
通过令
，可以得到本次实验所用到的滤波器表示方法：
（2.2）
其中
是截止频率，当
=
时，GLPF下降到其最大值的0.607处。
2、实验结果及讨论
1、实验结果
实验结果如图
（1）~图（9）所示：
图
（1）是中心化后的时域图及频谱图
图
（2）傅里叶变换后的二维幅度谱及相位谱
图（3）和图（4）是图
（2）中相应谱的三维图
图（5）傅里叶逆变换及中心化后和区实部后的图像
图（6）高斯低通滤波后的图像
图（7）和图（8）是高斯低通滤波器的三维显示及中心化后的图形
图（9）原图减去高斯低通滤波处理后图像的结果
图1中心变换后的输出图像
图2中心化后的二维频谱图
图3中心化后的三维幅度谱
图4中心化后的三维相位谱
图5傅里叶逆变换后中心化及取实部
图6高斯低通处理后图像
图7高斯低通滤波器中心化三维图像
图8高斯低通滤波器三维图像
图9原图减去低通处理后的图像
2、实验讨论
实验结果中，图
（1）是原图及其频谱与中心化后的图像与其频谱，图中原图的频谱图在四角有白色区（需要放大观察）、中心化后的频谱图比较易间（此图已经对数处理过），这与我们想要的结果一致。
可以直观的理解，将原图四角的白色区聚集到图形的中心。
图
（2）、（3）、（4）是其相应的二维与三位频谱图，三维更加直观。
图（5）是傅里叶逆变换的结果。
图（6）是高斯低通滤波器处理结果，输出图像变得模糊，体现了频率域低通滤波器的平滑效果。
图（7）和图（8）是高斯低通滤波器和中心化后的三维效果图，比较直观。
图（9）是原图减去高斯低通滤波处理后图像的结果，与高斯高通滤波处理后的图形很相似，通过改变低通滤波器的参数可以使其更加接近。
原理很简单，因为高斯高通滤波器就是一个全通滤波器减去高斯低通滤波器的结果。
3、程序附录
functionmain（）
%说明：
主函数包括三个部分：
1.中心化及DFT处理
%2.DFT逆变换
%3.高斯低通滤波
%第一部分：
中心化、DFT处理
ima=imread（'Fig4.41（a）.jpg'）;
imb=shifftmy（ima）;
imaf=fft2（ima）;
imbf=fft2（imb）;
imbf=log（1+abs（imbf））;
fa=fftshift（imaf）;%使用MATLAB自带的移频函数
fb=log（1+abs（fa））;
phase=angle（imaf）;
figure
subplot（2,2,1）;
imshow（ima）;
title（'Fig4.41（a）原图'）;
subplot（2,2,2）;
imshow（imb）;
title（'Fig4.41（a）时域*（-1）^（x+y）'）;
subplot（2,2,3）;
imshow（imaf,[]）;
title（'Fig4.41（a）频谱图'）;
subplot（2,2,4）;
imshow（imbf,[]）;
title（'中心化、对数变换后的频谱图'）;
figure
subplot（1,2,1）
imshow（fb,[]）;
title（'Fig4.41（a）中心化后的频谱图'）;
subplot（1,2,2）
imshow（phase,[]）;
title（'Fig4.41（a）相位图'）;
figure
mesh（fb）;
title（'Fig4.41（a）中心化后的三维频谱图'）;
figure
mesh（phase）;
title（'Fig4.41（a）三维相位图'）;
%第一部分结束
%第二部分：
DFT逆变换输出、中心化
imaif=ifft2（imaf）;
imaifsf=shifftmy（imaif）;
imar=real（imaifsf）;
figure
subplot（1,3,1）;
imshow（imaif）;
title（'FT逆变换结果'）;
subplot（1,3,2）;
imshow（imaifsf）;
title（'FT逆变换再*（-1）^（x+y）结果'）;
subplot（1,3,3）;
imshow（imar）;
title（'FT逆变换结果取其实部'）;
%第二部分结束
%第三部分：
高斯低通滤波处理
iml=ima;
PQ=paddedsize（size（iml））;
[U,V]=dftuv（PQ
（1）,PQ
（2））;
D0=0.05*PQ
（2）;
F=fft2（ima,PQ
（1）,PQ
（2））;
H=exp（-（U.^2+V.^2）/（2*（D0^2）））;
G=dftfilt（iml,H）;
figure
subplot（2,2,1）;
imshow（iml）;
title（'Fig4.41（a）原图'）;
subplot（2,2,2）;
HF=fftshift（H）;
imshow（HF,[]）;
title（'高斯低通滤波器'）;
subplot（2,2,3）;
F=log（1+abs（fftshift（F）））;
imshow（F,[]）;
title（'Fig4.41（a）频谱图'）;
subplot（2,2,4）;
imshow（G,[]）;
title（'低通滤波后'）;
figure
mesh（H）;
title（'高斯低通滤波器'）;
figure
mesh（HF）;
title（'高斯低通滤波器中心化'）;
imfz=iml-uint8（G）;
figure
imshow（imfz,[]）;
title（'原图减去滤波后的图像'）;
%第三部分结束
functionimb=shifftmy（ima）
%说明：
图像矩阵中心化处理
%输入参数：
%ima:
待处理图像
%输出参数：
%imb:
输出图像
%
[xs,ys]=size（ima）;
fori=1:
xs
forj=1:
ys
r=（-1）.^（i+j）;
imb（i,j）=ima（i,j）.*r;
end
end
end
functionPQ=paddedsize（AB,CD,PARAM）
%PADDEDSIZEComputespaddedsizesusefulforFFT-basedfiltering.
%PQ=PADDEDSIZE（AB）,whereABisatwo-elementsizevector,
%computesthetwo-elementsizevectorPQ=2*AB.
%
%PQ=PADDEDSIZE（AB,'PWR2'）computesthevectorPQsuchthat
%PQ
（1）=PQ
（2）=2^nextpow2（2*m）,wheremisMAX（AB）.
%
%PQ=PADDEDSIZE（AB,CD）,whereABandCDaretwo-elementsize
%vectors,computesthetwo-elementsizevectorPQ.Theelements
%ofPQarethesmallestevenintegersgreaterthanorequalto
%AB+CD-1.
%
%PQ=PADDEDSIZE（AB,CD,'PWR2'）computesthevectorPQsuchthat
%PQ
（1）=PQ
（2）=2^nextpow2（2*m）,wheremisMAX（[ABCD]）.
ifnargin==1
PQ=2*AB;
elseifnargin==2&~ischar（CD）
PQ=AB+CD-1;
PQ=2*ceil（PQ/2）;
elseifnargin==2
m=max（AB）;%Maximumdimension.
%Findpower-of-2atleasttwicem.
P=2^nextpow2（2*m）;
PQ=[P,P];
elseifnargin==3
m=max（[ABCD]）;%Maximumdimension.
P=2^nextpow2（2*m）;
PQ=[P,P];
else
error（'Wrongnumberofinputs.'）
end
function[U,V]=dftuv（M,N）
%DFTUV实现频域滤波器的网格函数
%Detailedexplanationgoeshere
u=0:
（M-1）;
v=0:
（N-1）;
idx=find（u>M/2）;%找大于M/2的数据
u（idx）=u（idx）-M;%将大于M/2的数据减去M
idy=find（v>N/2）;
v（idy）=v（idy）-N;
[V,U]=meshgrid（v,u）;
end
functiong=dftfilt（f,H）
%DFTFILTPerformsfrequencydomainfiltering.
%G=DFTFILT（F,H）filtersFinthefrequencydomainusingthe
%filtertransferfunctionH.Theoutput,G,isthefiltered
%image,whichhasthesamesizeasF.DFTFILTautomaticallypads
%FtobethesamesizeasH.FunctionPADDEDSIZEcanbeusedto
%determineanappropriatesizeforH.
%
%DFTFILTassumesthatFisrealandthatHisareal,uncentered
%circularly-symmetricfilterfunction.
%ObtaintheFFTofthepaddedinput.
F=fft2（f,size（H,1）,size（H,2））;
%Performfiltering.
g=real（ifft2（H.*F））;
%Croptooriginalsize.
g=g（1:
size（f,1）,1:
size（f,2））;
升级会员 