综合课设报告

东北大学秦皇岛分校计算机与通信工程学院综合课程设计设计题目 频率域滤波的Matlab 设计与实现专业名称电子信息工程班级学号4111321学生姓名张佳席指导教师韩光设计时间2014.6.20-2014.7.4东北大学秦皇岛分校计算机与通信工程学院综合课程设计任务书专业： 电子信息工程 班级： 41113 姓名： 张佳席 设计题目：频率域滤波的Matlab设计与实现一、设计实验条件计算机，数字信号处理实验室二、设计任务1、了解数字图像处理中关键问题；2、熟悉频率域滤波的产生背景及意义；3、熟悉频率域滤波的几种常用方法；4、掌握高斯低通滤波器算法的流程；5、用Matlab软件设计并实现数字图像的频率域滤波；三、设计说明书的内容1. 设计题目与设计任务（设计任务书）2. 前言（绪论）(设计的目的、意义等)3. 设计主体（各部分设计内容、分析、结论等）4. 结束语（设计的收获、体会等）5. 参考资料四、设计时间与设计时间安排1、设计时间： 2周2、设计时间安排： 熟悉实验设备、收集资料： 4 天设计图纸、实验、计算、程序编写调试： 2 天编写课程设计报告： 7 天答辩： 1 天目录1.设计题目与设计任务11.1设计题目11.2设计任务12.前言13.数字图像处理的关键问题14.频率域滤波的产生背景和意义24.1频率域滤波的产生背景24.2频率域滤波的意义35.频率域滤波的几种常用方法35.1低通滤波器35.2高通滤波器75.3带阻滤波器116.滤波器算法流程图147.程序仿真及结果147.1频率域滤波的matlab实现147.2低通滤波器仿真及结果157.3高通滤波器仿真及结果227.4带阻滤波器仿真及结果288.心得体会329.参考文献33341.设计题目与设计任务1.1设计题目频率域滤波的Matlab设计与实现1.2设计任务1、了解数字图像处理中关键问题；2、熟悉频率域滤波的产生背景及意义；3、熟悉频率域滤波的几种常用方法；4、掌握高斯低通滤波器算法的流程；5、用Matlab软件设计并实现数字图像的频率域滤波。2.前言数字图像处理是通过计算机对图像进行去除噪声、增强、复原、分割、提取特征等处理的方法和技术。它最早出现于20世纪50年代，当时的电子计算机已经发展到一定水平，人们开始利用计算机来处理图形和图像信息。早期的图像处理的目的是改善图像的质量，它以人为对象，以改善人的视觉效果为目的。现在图像处理还用于提取图像中所包含的某些特征或特殊信息图像。随着人类活动范围的不断扩大，图像处理的应用领域也将随之不断扩大。现广泛应用于航天和航空、生物医学工程、通信工程、工业和工程、军事公安、文化艺术、机器人视觉、视频和多媒体系统、科学可视化、电子商务等方面。在社会高度信息化的过程中，对人类获取信息中占有很大比例的视觉信息处理技术，毫无疑问将越来越受到人们的重视。在图像处理中，输入的是质量低的图像，输出的是改善质量后的图像，常用的图像处理方法有图像增强、复原、编码、压缩等。而频率域滤波则是实现图像处理的常用工具之一。3.数字图像处理的关键问题1）图像变换由于图像阵列很大，直接在空间域中进行处理，涉及计算量很大。因此，往往采用各种图像变换的方法，如傅立叶变换、沃尔什变换、离散余弦变换等间接处理技术，将空间域的处理转换为变换域处理，不仅可减少计算量，而且可获得更有效的处理（如傅立叶变换可在频域中进行数字滤波处理）。目前新兴研究的小波变换在时域和频域中都具有良好的局部化特性，它在图像处理中也有着广泛而有效的应用。2 ）图像编码压缩图像编码压缩技术可减少描述图像的数据量（即比特数），以便节省图像传输、处理时间和减少所占用的存储器容量。压缩可以在不失真的前提下获得，也可以在允许的失真条件下进行。编码是压缩技术中最重要的方法，它在图像处理技术中是发展最早且比较成熟的技术。3）图像增强和复原图像增强和复原的目的是为了提高图像的质量，如去除噪声，提高图像的清晰度的清晰度等。图像增强不考虑图像降质的原因，突出图像中所感兴趣的部分。如强化图像高频分量，可使图像中物体轮廓清晰，细节明显；如强调低频分量可减少图像中噪声影响。图像复原要求对图像降质的原因有一定的了解，一般讲根据降质过程建立“降质模型”，再采用某种滤波方法，恢复或重建原来的图像。4）图像分割图像分割是数字图像处理中的关键技术之一。图像分割时将图像中有意义的特征部分提取出来，包括图像中物体的边缘、区域等，这是进一步进行图像识别、分析和理解的基础。虽然目前已研究出不少边缘提取、区域分割的方法，但是还没有一种普遍使用于各种图像的有效方法。因此，对图像分割的研究还在不断深入之中，是目前图像处理中研究的热点之一。5）图像描述图像描述是图像识别和理解的必要前提。作为最简单的二值图像可采用其几何特性描述物体的特征。一般图像描述方法采用二维纹理特征描述。在计算机视觉中要进行物体描述的研究，因此，提出来体积描述、表面描述、广义圆柱体描述等方法。4.频率域滤波的产生背景和意义4.1频率域滤波的产生背景图像的频率是表征图像中灰度变化剧烈程度的指标，是灰度在平面空间上的梯度。如：大面积的沙漠在图像中是一片灰度变化缓慢的区域，对应的频率值很低；而对于地表属性变换剧烈的边缘区域在图像中是一片灰度变化剧烈的区域，对应的频率值较高。傅里叶变换在实际中有非常明显的物理意义，设f是一个能量有限的模拟信号，则其傅里叶变换就表示f的谱。从纯粹的数学意义上看，傅里叶变换是将一个函数转换为一系列周期函数来处理的。从物理效果看，傅里叶变换是将图像从空间域转换到频率域，其逆变换是将图像从频率域转换到空间域。换句话说，傅里叶变换的物理意义是将图像的灰度分布函数变换为图像的频率分布函数，傅里叶逆变换是将图像的频率分布函数变换为灰度分布函数。4.2频率域滤波的意义频率域滤波首先将输入的数字信号进行DFT变换，然后在频域进行各种有效的处理，然后进行DFT反变换，恢复为时域信号。这种用计算机对变换后的信号进行频域处理，比在时域中直接处理更加方便，计算量也大为减少，明显的提高了处理速度，并可采用二维数字滤波技术进行所需要的各种图像处理。5.频率域滤波的几种常用方法5.1低通滤波器在频谱中，低频主要对应图像在平滑区域的总体灰度级分布，而高频对应图像的细节部分，如边缘和噪声。因此，图像平滑可以通过衰减图像频谱中的高频部分来实现，这就建立了空间域图像与频域低通滤波之间的对应关系。1）理想低通滤波器我们最容易想到的衰减高频成分的方法就是在一个称为“截止频率”的位置“截断”所有的高频成分，将图像频谱中说有高于这一截止频率的频谱成分设为0，低于截止频率的成分设为保持不变，这样的滤波器我们称为理想低通滤波器。如果图像的宽度为M，高度为N，那么理想低通频域滤波器可形象化地描述为Hu,v=1, &(u-M2)2+(v-N2)2D00, &(u-M2)2+(v-N2)2D0其中D0表示理想低通滤波器的截止元素，滤波器的频率域原点在频谱图像的中心处，在一截止频率为半径的圆形区域之内的滤镜元素值都为1，而该圆之外的滤镜元素值全部为0。理想低通滤波器的频率特性在截止频率处十分陡峭，无法用硬件实现，这也是我们称为理想的原因，但其软件编程的模拟实现较为简单。理想低通滤波器在一定程度上可以去除噪声，但由此带来的也会带来较为明显的图像边缘和细节的模糊效应。图 1 理想低通滤波器曲面图Matlab实现：function out =lxdt(I,freq)%lxdt函数 构造理想低通滤波器%参数I 输入的灰度图像%参数freq 理想低通滤波器的截止频率M,N = size (I);out = ones (M,N);for i=1:M for j=1:N if(sqrt(i-M/2)2+(j-N/2)2) freq ) out(i,j)=0; end endendend2）高斯低通滤波器高斯低通滤波器的频率域二维形式由下式给出。Hu,v=e-u-M22+v-N22/22高斯函数具有相对简单的形式，而且它的傅里叶变换和傅里叶反函数都是实高斯函数。并且，当增大时，H(u,v)的图像倾于变宽，而其反函数的图像倾于变窄和变高。这体现了频率域和空间域的对应关系。频率域滤波器越窄，滤除的高频成分越多，图像就越平滑（模糊）；而在空间域，对应的滤波器就越宽，相应的卷积模板越平坦，平滑（模糊）效果就越明显。图 2 高斯低通滤波器曲面图MATLAB实现：function out = gsdt( I,sigma )%gsdt函数 构造频域高斯低通滤波器%参数I 输入的灰度图像%参数sigma 高斯函数的Sigma参数M,N = size(I);out = ones(M,N);for i=1:M for j=1:N out(i,j) = exp(-(i-M/2)2+(j-N/2)2)/(2*sigma2); endendend3）Butterworth低通滤波器n阶Butterworth低通滤波器的频率域二维形式由下式给出：Hu,v=11+1+Du,vD02n其中，D0为截止频率。当Du,vD0=1时，=0.5，它的特性是传递函数比较平滑，连续衰减，而不像理想滤波器那样陡峭变化，即明显的不连续。图 3 Butterworth低通滤波器曲面图Matlab实现：function out =butterworth_dt(I,freq)%lxdt函数 构造butterworth低通滤波器%参数I 输入的灰度图像%参数freq butterworth低通滤波器的截止频率M,N = size (I);n=2; %阶数设为2 m1=fix(M/2);n1=fix(N/2); for i=1:M for (j=1:N) f=sqrt(i-m1)2+(j-n1)2); out(i,j)=1/(1+(f/freq)(2*n); end end end5.2高通滤波器图像锐化可以通过衰减图像频谱中的低频成分来实现。这就建立了空间域图像锐化域频域高通滤波之间的对应关系1）高斯高通滤波器参考5.1中高斯低通滤波器中的H（u,v）函数，如果我们需要做相反的滤波操作，滤除低频成分而留下高频成分，则可以考虑简单地使用如下表达式来获得一个高斯高通滤波器。Hu,v=1-e-u-M22+v-N22/22图 4 高斯高通滤波器曲面图Matlab实现:function out = gsgt( I,sigma )%gsdt函数 构造频域高斯低通滤波器%参数I 输入的灰度图像%参数sigma 高斯函数的Sigma参数M,N = size(I);out = ones(M,N);for i=1:M for j=1:N out(i,j) = 1-exp(-(i-M/2)2+(j-N/2)2)/(2*sigma2); endendend2）频域拉普拉斯滤波器频域拉普拉斯算子的推导可以从一维开始，由傅里叶变换的性质可知FFTdnfxdxn=junF(u)因此拉普拉斯算子的傅里叶变换计算如下。FFT2fx,yx2+2fx,yy2=ju2Fu,v+jv2Fu,v =-u2+v2F(u,v)从而可得频域的拉普拉斯滤波器为Hu,v=-(u2+v2)根据频域图像频率的原点平移规律，将上式改写为Hu,v=-u-M22+v-N22M,N分别为图像的高和宽。图 5 拉普拉斯滤波器曲面图Matlab实现:function out = lpls( I)%gsdt函数 构造频域拉普拉斯滤波器%参数I 输入的灰度图像M,N = size(I);out = ones(M,N);for i=1:M for j=1:N out(i,j) = -(-(i-M/2)2+(j-N/2)2)/(2*sigma2); endendend3）Butterworth高通滤波器n阶Butterworth低通滤波器的频率域二维形式由下式给出：Hu,v=11+1+D0D(u,v)2n其中，D0为截止频率。图 6 Butterworth高通滤波器Matlab实现：function out =butterworth_gt(I,freq)%lxdt函数 构造butterworth高通滤波器%参数I 输入的灰度图像%参数freq butterworth高通滤波器的截止频率M,N = size (I);n=2; %阶数设为2m1=fix(M/2);n1=fix(N/2); for i=1:M for (j=1:N) f=sqrt(i-m1)2+(j-n1)2); out(i,j)=1/(1+(freq/f)(2*n); end end end5.3带阻滤波器“带阻”就是阻止频谱中某一频带范围的分量通过，其他频率成分不受影响。带阻滤波器常用来处理含有周期性噪声的图像。常用的带阻滤波器有理想带阻滤波器和高斯带阻滤波器。1）理想带阻滤波器理想带阻滤波器的表达式为H(u,v)=0, &|D-D0|W1, &|D-D0|W式中：D0阻塞频带中心频率到频率原点的距离；是阻塞频带宽度；是(u,v)点到频率原点的距离。Matlab实现：function out = lxdz( I,freq, width)%gsdt函数 构造频域理想带阻滤波器%参数I 输入的灰度图像%参数freq 阻带中心频率%参数width 阻带宽度M,N = size(I);out = ones(M,N);for i=1:M for j=1:N if(sqrt(i-M/2)2+(j-N/2)2)-freqwidth & freq-sqrt(i-M/2)2+(j-N/2)2) width) out (i,j) = 0; endendend图 7 理想带阻滤波器曲面图2）高斯带阻滤波器高斯带阻滤波器的表达式为Hu,v=1-e-12D2u,v-D02Du,vW2式中：D0是阻塞频带中心频率到频率原点的距离；W是阻塞频带宽度；D是(u,v)点到频率原点的距离。图 8 高斯带阻滤波器曲面图Matlab实现：function out = gsdz( I,freq, width)%gsdt函数 构造频域高斯带阻滤波器%参数I 输入的灰度图像%参数freq 阻带中心频率%参数width 阻带宽度M,N = size(I);out = ones(M,N);for i=1:M for j=1:N out(i,j) = 1-exp(-0.5*(i-M/2)2+(j-N/2)2)-freq2)/(sqrt(i.2+j.2)*width)2); endendend6.滤波器算法流程图图 9 频率域滤波流程图7.程序仿真及结果7.1频率域滤波的matlab实现function out = pylb( I,ff )%pylb函数 对灰度图像就行频域滤波%参数I 输入的空域图像%参数ff 应用的与原图像等大的频域滤镜%快速傅里叶变换f = fft2(double(I);%移动原点s = fftshift(f);%应用滤镜及反变换out = s.*ff;out = ifftshift(out);out = ifft2(out);%求模值out = abs(out);%归一化以便显示out = out/max(out(:);end7.2低通滤波器仿真及结果1）理想低通滤波器a.仿真程序：I = imread(C:ks.jpg); %打开图像I = I(:,:,1); %设定图像层数，使图像能够进行二维fft计算%显示原图像figure(1);subplot(2,3,1);imshow(I);title(原图像);%显示原图像频域temp = fft2(double(I);temp = log(1+abs(temp);temp = fftshift(temp);figure (1);subplot(2,3,4);imshow(temp, );title(原图像频域);%生成及应用理想低通滤波器滤镜，截止频率为60ff = lxdt(I, 60);out = pylb(I,ff);%显示理想低通滤波器滤波后的图像figure(1);subplot(2,3,2);imshow(out);title(理想低通滤波器滤波后图像，freq=60);%显示理想低通滤波器滤波后的图像频域temp = fft2(out);temp = log(1+abs(temp);temp = fftshift(temp);figure (1);subplot(2,3,5);imshow(temp, );title(理想低通滤波器滤波后频域，freq=60);%生成及应用理想低通滤波器滤镜，截止频率为80ff = lxdt(I, 80);out = pylb(I,ff);%显示理想低通滤波器滤波后的图像figure(1);subplot(2,3,3);imshow(out);title(理想低通滤波器滤波后图像，freq=80);%显示理想低通滤波器滤波后的图像频域temp = fft2(out);temp = log(1+abs(temp);temp = fftshift(temp);figure (1);subplot(2,3,6);imshow(temp, );title(理想低通滤波器滤波后频域，freq=80);b.仿真结果图 10 理想低通滤波器仿真图由图可见，当截止频率非常地时，只有非常靠近原点的低频成分能够通过，图像模糊严重；截止频率越高，通过的频率成分就越多，图像模糊的程度越小，所获得的图像也就越接近原图像。但可以看出，理想低通滤波器不能很好地兼顾噪声滤除和细节保留两个方面。2）高斯低通滤波器a. 仿真程序I = imread(C:ks.jpg); %打开图像I = I(:,:,1); %设定图像层数，使图像能够进行二维fft计算%显示原图像figure(1);subplot(2,3,1);imshow(I);title(原图像);%显示原图像频域temp = fft2(double(I);temp = log(1+abs(temp);temp = fftshift(temp);figure (1);subplot(2,3,4);imshow(temp, );title(原图像频域);%生成及应用高斯低通滤波器滤镜,sigma=60ff = gsdt(I, 60);out = pylb(I,ff);%显示高斯低通滤波器滤波后的图像figure(1);subplot(2,3,2);imshow(out);title(高斯低通滤波器滤波后图像,sigma=60);%显示高斯低通滤波器滤波后的图像频域temp = fft2(out);temp = log(1+abs(temp);temp = fftshift(temp);figure (1);subplot(2,3,5);imshow(temp, );title(高斯低通滤波器滤波后频域,sigma=60);%生成及应用高斯低通滤波器滤镜,sigma=100ff = gsdt(I, 100);out = pylb(I,ff);%显示高斯低通滤波器滤波后的图像figure(1);subplot(2,3,3);imshow(out);title(高斯低通滤波器滤波后图像,sigama=100);%显示高斯低通滤波器滤波后的图像频域temp = fft2(out);temp = log(1+abs(temp);temp = fftshift(temp);figure (1);subplot(2,3,6);imshow(temp, );title(高斯低通滤波器滤波后频域,sigma=100);图 11 高斯低通滤波器仿真图b. 仿真结果由图可以看出，相比于理想低通滤波器，高斯低通滤波器在有效抑制噪声的同时，模糊程度更低，对边缘带来的混叠程度更小。3）Butterworth低通滤波器仿真程序：I = imread(C:ks.jpg); %打开图像I = I(:,:,1); %设定图像层数，使图像能够进行二维fft计算%显示原图像figure(1);subplot(2,3,1);imshow(I);title(原图像);%显示原图像频域temp = fft2(double(I);temp = log(1+abs(temp);temp = fftshift(temp);figure (1);subplot(2,3,4);imshow(temp, );title(原图像频域);%生成及应用butterworth低通滤波器滤镜，截止频率为60ff = butterworth_dt(I, 60);out = pylb(I,ff);%显示butterworth低通滤波器滤波后的图像figure(1);subplot(2,3,2);imshow(out);title(butterworth低通滤波器滤波后图像，freq=60);%显示butterworth低通滤波器滤波后的图像频域temp = fft2(out);temp = log(1+abs(temp);temp = fftshift(temp);figure (1);subplot(2,3,5);imshow(temp, );title(butterworth低通滤波器滤波后频域，freq=60);%生成及应用butterworth低通滤波器滤镜，截止频率为80ff = butterworth_dt(I, 80);out = pylb(I,ff);%显示butterworth低通滤波器滤波后的图像figure(1);subplot(2,3,3);imshow(out);title(butterworth低通滤波器滤波后图像，freq=80);%显示butterworth低通滤波器滤波后的图像频域temp = fft2(out);temp = log(1+abs(temp);temp = fftshift(temp);figure (1);subplot(2,3,6);imshow(temp, );title(butterworth低通滤波器滤波后频域，freq=80);b. 仿真结果：图 12 Butterworth低通滤波器仿真图7.3高通滤波器仿真及结果1）高斯高通滤波器a. 仿真程序I = imread(C:ks.jpg); %打开图像I = I(:,:,1); %设定图像层数，使图像能够进行二维fft计算%显示原图像figure(1);subplot(2,3,1);imshow(I);title(原图像);%显示原图像频域temp = fft2(double(I);temp = log(1+abs(temp);temp = fftshift(temp);figure (1);subplot(2,3,4);imshow(temp, );title(原图像频域);%生成及应用高斯高通滤波器滤镜,sigma=20ff = gsgt(I, 20);out = pylb(I,ff);%显示高斯高通滤波器滤波后的图像figure(1);subplot(2,3,2);imshow(out);title(高斯高通滤波器滤波后图像,sigma=20);%显示高斯高通滤波器滤波后的图像频域temp = fft2(out);temp = log(1+abs(temp);temp = fftshift(temp);figure (1);subplot(2,3,5);imshow(temp, );title(高斯高通滤波器滤波后频域,sigma=20);%生成及应用高斯高通滤波器滤镜,sigma=40ff = gsgt(I, 40);out = pylb(I,ff);%显示高斯高通滤波器滤波后的图像figure(1);subplot(2,3,3);imshow(out);title(高斯高通滤波器滤波后图像,sigama=40);%显示高斯高通滤波器滤波后的图像频域temp = fft2(out);temp = log(1+abs(temp);temp = fftshift(temp);figure (1);subplot(2,3,6);imshow(temp, );title(高斯高通滤波器滤波后频域,sigma=40);b. 仿真结果图 13 高斯高通滤波器仿真图高斯高通滤波器可以较好地提取图像中的边缘信息，Sigma参数取值越小，边缘提取越不精确，会包含越多的非边缘信息；Sigma参数取值越大，边缘提取越精确，但可能包含不完整的边缘信息。2）拉普拉斯滤波器a. 仿真程序I = imread(C:ks.jpg); %打开图像I = I(:,:,1); %设定图像层数，使图像能够进行二维fft计算%显示原图像figure(1);subplot(2,2,1);imshow(I);title(原图像);%显示原图像频域temp = fft2(double(I);temp = log(1+abs(temp);temp = fftshift(temp);figure (1);subplot(2,2,3);imshow(temp, );title(原图像频域);%生成及应用拉普拉斯滤波器滤波器滤镜ff = gsgt(I, 20);out = pylb(I,ff);%显示拉普拉斯滤波器滤波后的图像figure(1);subplot(2,2,2);imshow(out);title(拉普拉斯滤波器滤波后图像);%显示拉普拉斯滤波器滤波后的图像频域temp = fft2(out);temp = log(1+abs(temp);temp = fftshift(temp);figure (1);subplot(2,2,4);imshow(temp, );title(拉普拉斯滤波器滤波后频域);b. 仿真结果图 14 拉普拉斯滤波器仿真图3）Butterworth低通滤波器a. 仿真程序I = imread(C:ks.jpg); %打开图像I = I(:,:,1); %设定图像层数，使图像能够进行二维fft计算%显示原图像figure(1);subplot(2,3,1);imshow(I);title(原图像);%显示原图像频域temp = fft2(double(I);temp = log(1+abs(temp);temp = fftshift(temp);figure (1);subplot(2,3,4);imshow(temp, );title(原图像频域);%生成及应用butterworth高通滤波器滤镜，截止频率为60ff = butterworth_gt(I, 60);out = pylb(I,ff);%显示butterworth高通滤波器滤波后的图像figure(1);subplot(2,3,2);imshow(out);title(butterworth高通滤波器滤波后图像，freq=60);%显示butterworth高通滤波器滤波后的图像频域temp = fft2(out);temp = log(1+abs(temp);temp = fftshift(temp);figure (1);subplot(2,3,5);imshow(temp, );title(butterworth高通滤波器滤波后频域，freq=60);%生成及应用butterworth高通滤波器滤镜，截止频率为80ff = butterworth_gt(I, 80);out = pylb(I,ff);%显示butterworth高通滤波器滤波后的图像figure(1);subplot(2,3,3);imshow(out);title(butterworth高通滤波器滤波后图像，freq=80);%显示butterworth高通滤波器滤波后的图像频域temp = fft2(out);temp = log(1+abs(temp);temp = fftshift(temp);figure (1);subplot(2,3,6);imshow(temp, );title(butterworth高通滤波器滤波后频域，freq=80);b. 仿真结果图 15 butterworth高通滤波器仿真图7.4带阻滤波器仿真及结果1）理想带阻滤波器a. 仿真程序I = imread(C:ks.jpg); %打开图像I = I(:,:,2); %设定图像层数，使图像能够进行二维fft计算%显示原图像figure(1);subplot(2,2,1);imshow(I);title(原图像);%显示原图像频域temp = fft2(double(I);temp = log(1+abs(temp);temp = fftshift(temp);figure (1);subplot(2,2,3);imshow(temp, );title(原图像频域);%生成及应用理想带阻滤波器滤镜,freq=300，width=10ff = lxdz(I, 300,40);out = pylb(I,ff);%显示理想带阻滤波器滤波后的图像figure(1);subplot(2,2,2);imshow(out);title(理想带阻滤波器滤波后图像,freq=300，width=10