实验4DFT及频域滤波

1、实验名称：数字图像DFT及频域滤波一 实验目的1 了解数字图像各种正交变换的概念、原理和用途。2 深熟练掌握数字图像的 DFT/DCT 的原理、方法和实现流程，熟悉两种变换的性质，并能对数字图像 DFT 及 DCT 的结果进行必要解释。3 熟悉熟悉和掌握利用 MATLAB 工具进行数字图像 FFT 及 DCT 的基本步骤、MATLAB 函数使用及具体变换的处理流程，并能根据需要进行必要的频谱分析和可视化显示。4 熟悉利用空域滤波器构建对应的频域滤波器的方法和关键步骤。5 熟悉和掌握几种典型的频域低通滤波器及高通滤波器的原理、特性和作用。6 搞清空域图像处理与频域图像处理的异同，包括处理流程、各

2、自的优势等。掌握 频域滤波的基本原理和基本流程，并能编写出相应的程序代码。二 实验原理1、 FFT实验l 傅立叶变换：非周期函数表示为正弦和/或余弦乘以加权函数的积分。l 一维连续Fourier变换对函数f(x)进行傅立叶变换得到F(u) 逆变换，即将F(u)变换到f(x)为 l 一维离散Fourier变换正变换(DFT) 逆变换(IDFT) l 用幅值和相位表示傅立叶变换 l 用途：频域滤波，图像增强，频谱分析。三 实验步骤1. 模型图像的 FFT 实验图 Error! Main Document Only.1) 利用 MATLAB 程序自行生成如图（1）的（a）、（b）所示的二值图像，分别

3、对其分别进行离散傅立叶变换（DFT）计算。 （a） （b）2) 对变换结果做频谱中心化处理，并分别显示出其 2D 频谱图以及对应的 3D频谱图。3) 对以上两幅原始图像 FFT 后的频谱图进行分析，可以得出什么样的结论或验证了DFT的什么性质。图 22. 实际图像的 FFT 实验1) 读任意读取一幅灰度图像，对其进行 FFT 变换，分别画出变换前的原始图像及其 FFT 后频谱的 2D 平面图，要求变换结果做频谱中心化处理。2) 画出对应的3D中心化频谱。图 33. 数字图像的频域滤波处理1) 设定截止频率 D0 =100，试分别构建 256×256 的频域理想低通滤波器（ILPF）和

4、频域理想高通滤波器（IHPF）。编写 MATLAB 代码，分别画出它们的频域滤波器响应 3D 图及其对应的 2D 投影平面图。2) 任意读取一幅数字图像。编写 MATLAB 代码，分别利用理想低通滤波器（ILPF） 和理想高通滤波器（IHPF）对其进行频域滤波处理。具体可参考如下步骤：（1） 消除折叠现象的填充：P=2M，Q=2N，其中 M、N 为原图像的尺寸；（2） 原图 f (x, y) 的傅立叶变换，产生 F(u,v) ；（3） 频谱中心化处理（当然也可以不变换，则 H(u,v) 要改变）；（4） 用滤波器函数 H(u,v) 乘以 F(u,v) ；（5） 傅立叶反变换；（6） 取实数部分

5、，绝对值很小的虚数部分是浮点运算存在误差造成的；（7） 空间域中心还原变换（反中心化） 。3) 分别显示 ILPF 和 IHPF 两种滤波器滤波前的原始图像、滤波器频域响应（2D）及滤波结果图像。4) 对滤波结果进行必要分析，如振铃现象、图像模糊、图像变暗等说明原因，以及 尽可能克服这些现象的有效措施等。 图 4四 实验结果及分析图 Error! Main Document Only.1. 模型图像的 FFT 实验图 Error! Main Document Only.图 Error! Main Document Only.l 分析I1的白色矩形面积小与I，可以认为I1是I的缩小。l 而傅立叶

6、有如下的变换性质： 从图（5）和图（6）可以验证该性质。图 Error! Main Document Only.2. 实际图像的FFT实验图 Error! Main Document Only.图 Error! Main Document Only.3. 数字图像的频域滤波处理图 Error! Main Document Only.图 Error! Main Document Only.l 振铃现象：已知频域上的方波对应着空域的sinc波，频域：，傅立叶反变换到时域有：，频域相乘相当于时域作卷积。图像经过理想低通滤波器后，时域上相当于原始图像与sinc函数卷积，由于sinc函数振荡，则卷积后图

7、像也会振荡；或者说由于sinc函数有两个负边带，卷积后图像信号两侧出现“过冲现象”，而且能量不集中，即产生振铃效应。可以使用高斯低通滤波解决该现象。l 图像模糊：由于理想低通滤波器，滤去高通成分，削弱了图像的细节，从而使图像变得模糊l 图像变暗：由于高通滤波器会滤去低通成分，只突出图像细节，所以图像变暗，为改善这种情况，将高通滤波器的转移函数加上一常数量以便引入一些低频分量，也就是高频提升滤波器。 六.实验心得体会和建议l 心得体会：通过这次实验使我明白了可以利用傅立叶变换将空域滤波转换在频域滤波，更加符合人们的认知，而且便于在实际中实现。l 建议：可以让大家编写相关代码解决振铃现象，图像变暗

8、。七.程序源代码1. 模型图像的FFT实验fs=1000;I=zeros(256,256);%创建I矩阵I(98:158,98:158)=1;I1=zeros(256,256);%创建I1矩阵I1(118:138,118:138)=1;subplot(121),imshow(I);title('I');%分别在同一幅图中画出I和I1矩阵subplot(122),imshow(I1);title('I1');F=fft2(I);F_1=abs(fftshift(F);%分别对I和I1矩阵做二维傅立叶变换，并对结果中心化处理F1=fft2(I1);F1_1=abs(

9、fftshift(F1);figure,%分别显示2D频谱图subplot(121),imshow(log(1+F_1),);title('The 2D Spectrum of I');subplot(122),imshow(log(1+F1_1),);title('The 2D Spectrum of I1');figure,%分别显示3D频谱图subplot(121),mesh(F_1);title('The 3D Spectrum of I');xlabel('fx'),ylabel('fy'),zlabel

10、('amplitude');subplot(122),mesh(F1_1);title('The 3D Spectrum of I1');xlabel('fx'),ylabel('fy'),zlabel('amplitude');2. 实际图像的FFT实验clc,clear,close all;I=imread('barbara.jpg');I=double(I);I1=fft2(I);I2=abs(fftshift(I1);subplot(121),imshow(I,);title('原图

11、')subplot(122),imshow(log(1+I2),);title('频谱图')M,N=size(I2);m=71;n=71;Fx = (0:m-1)-floor(m/2);Fy = (0:n-1)-floor(n/2);figurex,y = meshgrid(Fx, Fy); surfl(x,y,I2(M/2-floor(m/2):M/2+floor(m/2),N/2-floor(n/2):N/2+floor(n/2)/max(I2(:);title('3D频域图'),xlabel('fx'),ylabel('fy

12、'),zlabel('振幅')3. 数字图像的频域滤波处理M=256;N=256;for u=1:2*M for v=1:2*N D(u,v)=sqrt(u-M)2+(v-N)2); endendd0=100;for u=1:2*M for v=1:2*N if D(u,v)>d0 H1(u,v)=0;%H1低通滤波 H2(u,v)=1;%H2高通滤波 else H1(u,v)=1; H2(u,v)=0; end endendfigure,%分别显示ILPF和IHPF的2D投影平面图subplot(121),imshow(H1,);title('2D IL

13、PF 频谱图')subplot(122),imshow(H2,);title('2D BLPF 频谱图')figure,%分别显示ILPF和IHPF的3D图subplot(121),mesh(H1);title('3D ILPF 频谱图')xlabel('u'),ylabel('v'),zlabel('振幅')subplot(122),mesh(H2);title('3D BLPF 频谱图')xlabel('u'),ylabel('v'),zlabel(

14、9;振幅')I=imread('cameraman.tif');%读入要处理的图像I=double(I);Ip=zeros(2*M,2*N);%预分配内存Ip(1:M,1:N)=I;%加0处理F=fft2(Ip);%2DfftF0=fftshift(F);%频域中心化处理F1=F0.*H1;F2=F0.*H2;Fa=real(ifft2(ifftshift(F1);%傅立叶反变换取实部Fb=real(ifft2(ifftshift(F2);figure,subplot(231),imshow(I,);title('原图');subplot(232),imshow(log(1+F1),);title('中心化频率谱');subplot(233),imshow(Fa(1:M,1:N),);title('ILPF滤波结果');subplot(235),imshow(log(1+F2),);title('中心化频率谱');subplot(236),imshow(Fb(1:M,1:N),);title(&