数字信号处理基于MATLAB的FFT算法设计课设说明书

1、课程设计说明书 目录 TOC o 1-3 h z u HYPERLINK l _Toc345629352 1引言 PAGEREF _Toc345629352 h 1 HYPERLINK l _Toc345629353 2课设要求 PAGEREF _Toc345629353 h 2 HYPERLINK l _Toc345629354 2.1课设题目 PAGEREF _Toc345629354 h 2 HYPERLINK l _Toc345629355 2.2设计内容及要求 PAGEREF _Toc345629355 h 2 HYPERLINK l _Toc345629356 2.3要求的设计成果

2、 PAGEREF _Toc345629356 h 2 HYPERLINK l _Toc345629357 3基于MATLAB的FFT算法实现 PAGEREF _Toc345629357 h 3 HYPERLINK l _Toc345629358 3.1系统总体流程图 PAGEREF _Toc345629358 h 3 HYPERLINK l _Toc345629359 3.2 DIT-FFT算法的基本原理 PAGEREF _Toc345629359 h 3 HYPERLINK l _Toc345629360 3.3 DIT-FFT算法的运算规律及编程思想 PAGEREF _Toc3456293

3、60 h 5 HYPERLINK l _Toc345629361 4 MATLAB实现程序 PAGEREF _Toc345629361 h 8 HYPERLINK l _Toc345629362 5用GUI界面实现运算 PAGEREF _Toc345629362 h 10 HYPERLINK l _Toc345629363 5.1 GUI简介 PAGEREF _Toc345629363 h 10 HYPERLINK l _Toc345629364 5.2界面设计 PAGEREF _Toc345629364 h 10 HYPERLINK l _Toc345629365 5.3 GUI实现程序 P

4、AGEREF _Toc345629365 h 12 HYPERLINK l _Toc345629366 5.4运行调试 PAGEREF _Toc345629366 h 13 HYPERLINK l _Toc345629367 6自编算法与内置算法结果比较 PAGEREF _Toc345629367 h 14 HYPERLINK l _Toc345629368 7总结 PAGEREF _Toc345629368 h 15 HYPERLINK l _Toc345629369 参考文献 PAGEREF _Toc345629369 h 16 HYPERLINK l _Toc345629370 附录 P

5、AGEREF _Toc345629370 h 17 HYPERLINK l _Toc345629371 附录 PAGEREF _Toc345629371 h 211引言MATLAB是矩阵实验室(Matrix Laboratory)的简称，是美国MathWorks公司出品的商数学软件，用于算法开发、数据可视化、数据分析以及数值计算的高级技术计算语言和交互式环境，主要包括MATLAB和Simulink两大部分。MATLAB 的应用范围非常广，包括信号和图像处理、通讯、控制系统设计、测试和测量、财务建模和分析以及计算生物学等众多应用领域。附加的工具箱（单独提供的专用 MATLAB 函数集）扩展了 M

6、ATLAB 环境，以解决这些应用领域内特定类型的问题。它以矩阵运算为基础，把计算、可视化、程序设计融合在一个简单易用的交互式工作环境中，是一款数据分析和处理功能都非常强大的工程适用软件。它可以将图片文件变换为离散的数据文件，然后利用其强大的矩阵运算能力处理数据，如数据滤波、傅立叶变换、时域和频域分析以及各种图的呈现等，它的信号处理与分析工具箱为图片分析提供了十分丰富的功能函数，利用这些功能函数可以快捷而又方便的完成图片信号的处理和分析以及信号的可视化。数字信号处理是MATLAB重要应用的领域之一。对于有限长序列x(n)，若要求其N点的傅里叶变换(DFT)需要经过次复数乘法运算和N*(N-1)次

7、复数加法运算。随着N的增加，运算量将急剧增加，而在实际问题中，N往往是较大的，如当N=1024时，完成复数乘法和复数加法的次数分别为百万以上，无论是用通用计算机还是用DSP芯片，都需要消耗大量的时间和机器内存，不能满足实时的要求。因此，DFT的这种运算只能进行理论上的计算，不适合对实时处理要求高的场合。因此，研究作为DSP的快速算法的FFT是相当必要的，快速傅里叶变换(FFT)是为提高DFT运算速度而采用的一种算法，快速算法的种类很多，而且目前仍在改进和提高，它是根据离散傅里叶变换的奇、偶、虚、实等特性，对离散傅立叶变换的算法进行改进获得的。基于本学期所学的DIT-FFT的运算规律和编程思想以

8、及MATLAB的学习和使用，本课设要求在MATLAB环境下编写基2 DIT-FFT算法实现对离散信号的快速傅里叶变换，再与MATLAB软件自带的FFT函数实现对离散信号的傅里叶变换进行比较，如果得到的频谱相同，那么我们编写的程序就是正确的。如果有能力可以选做系统人机对话界面，用GUI界面完成人机交互方便使用的。本课程设计主要是对数字信号的分析。2课设要求2.1课设题目基于MATLAB的FFT算法的设计2.2设计内容及要求1.设计内容：所设计的FFT算法应完成以下功能：(1) 在MATLAB环境下编写FFT算法（不调用系统现有函数）；(2) 实现对选定图片进行FFT计算、还原（IFFT计算），并

9、与系统FFT函数做对比，进行分析；(3) 设计GUI界面。2.设计要求：(1) 根据题目要求进行算法GUI总体设计； (2) 完成算法具体部分的设计，即算法原理图和算法原理说明；(3) 算法程序的设计，即对选定图片进行自编FFT计算与还原，并与自带函数进行对比与完整源程序；(4) 书写设计说明书。2.3要求的设计成果(1) 设计结果能正确仿真演示；(2) 设计说明书一份（包括总体设计、算法原理图及说明、系统GUI演示、源程序清单等）。3基于MATLAB的FFT算法实现3.1系统总体流程图本设计要求对一选定的图片进行FFT算法和IFFT算法分析。在MATLAB环境下编写基2 DIT-FFT算法，

10、利用自己编写的算法对图片进行频谱分析，并与MATLAB数字信号处理工具箱中的FFT函数进行对比研究，验证自编算法的正确性。所以得到系统总体流程图如图3-1所示。图片信号采集完成信号时域图完成信号频率响应编写FFT程序图实现输入信号的倒序实现一级中不同种蝶形算运实现一级中相同种蝶形运算与Matlab自带的FFT比较图3-1系统总体流程图3.2 DIT-FFT算法的基本原理快速傅里叶变换(FFT)是为提高DFT运算速度而采用的一种算法。对一个有限长度序列x(n)的N点的DFT为：所以，要求N点的DFT，需要N2次的复数乘法运算，N*(N-1)次复数乘法运算。随着N的增加，运算量将急剧增加，而在实际

11、问题中，N往往是较大的，如当N=1024时，完成复数乘法和复数加法的次数分别为百万以上，无论是用通用计算机还是用DSP芯片，都需要消耗大量的时间，不能满足实时的要求，不适合于对实时处理要求高的场合。为了能实时处理DFT，要想减少DFT的运算量可以有两个途径：第一是降N，N的值减小了，运算量就减少了；第二是利用旋转因子的周期性，对称性和可约性。利用这两个途径实现DFT的快速傅里叶变换(FFT)，FFT算法基本上可分为按时间抽取的FFT算法(DIT-FFT)和按频率抽取的FFT算法(DIF-FFT)。旋转因子的性质：(1) 周期性(2) 共轭对称性(3) 可约性本次课设要求用基2的按时间抽取的FF

12、T算法(DIT-FFT)实现FFT功能，设序列x(n)的长度为N，且N满足N=2M，M为正整数。若N不能满足上述关系，可以将序列x(n)补零实现。按时间抽取基2-FFT算法的基本思路是将N点序列按时间下标的奇偶分为两个N/2点序列，计算这两个N/2点序列的N/2点DFT，计算量可减小约一半；每一个N/2点序列按照同样的划分原则，可以划分为两个N/4点序列，最后，将原序列划分为多个2点序列，将计算量大大降低。按时间下标的奇偶将N点x(n)分别抽取组成两个N/2点序列，分别记为x1(n)和x2(n)，将x(n)的DFT转化为x1(n)和x2(n)的DFT的计算。利用旋转因子的可约性，即：用蝶形运算

13、可表示为如图3-2所示：图3-2 DIT-FFT蝶形运算流图符号以此类推，还可以把x1(n)和x2(n)按n值得奇偶分为两个序列，这样就达到了降N得目的，从而减少了运算量。FFT对DFT的数学运算量改进：直接采用DFT进行计算，运算量为N2次复数乘法和N*(N-1)次复数乘法。当采用M次FFT时，由N=2M求得M=logN，运算流图有M级蝶形，每一级都由N/2个蝶形运算构成，这样每一级蝶形运算都需要N/2次复数乘法和N次复数加法。M级运算共需要复数乘法次数为C=N/2*M，复数加法次数为C=N*M。当N值较大时，FFT减少运算量的特点表现的越明显。3.3 DIT-FFT算法的运算规律及编程思想

14、为了编写DIT-FFT算法的运算程序，首先要分析其运算规律，总结编程思想并绘出程序框图。1.原位计算对点的FFT共进行M级运算，每级由N/2个蝶形运算组成。在同一级中，每个蝶的输入数据只对本蝶有用，且输出节点与输入节点在同一水平线上，这就意味着每算完一个蝶后，所得数据可立即存入原输入数据所占用的数组元素（存储单元）中，这种原位（址）计算的方法可节省大量内存。2.蝶形运算实现FFT运算的核心是蝶形运算，找出蝶形运算的规律是编程的基础。蝶形运算是分级进行的，每级的蝶形运算可以按旋转因子的指数大小排序进行，如果指数大小一样则可从上往下依次蝶算。对点的FFT共有M级运算，用L表示从左到右的运算级数(L

15、=1，2，M )。第L级共有个不同指数的旋转因子，用R表示这些不同指数旋转因子从上到下的顺序(R=0，1，B-1)。第R个旋转因子的指数，旋转因子指数为P的第一个蝶的第一节点标号k从R开始，由于本级中旋转因子指数相同的蝶共有个，且这些蝶的相邻间距为，故旋转因子指数为P的最后一个蝶的第一节点标号k为：，本级中各蝶的第二个节点与第一个节点都相距B点。应用原位计算，蝶形运算可表示成如下形式： (J)= (J)+ (J+B)* (J+B)= (J)-(J+B)*总结上述运算规律，可采用如下运算方法进行DIT-FFT运算。首先读入数据，根据数据长度确定运算级数M，运算总点数，不足补0处理。然后对读入数据

16、进行数据倒序操作。数据倒序后从第1级开始逐级进行，共进行M级运算。在进行第L级运算时，先算出该级不同旋转因子的个数（也是该级中各个蝶形运算两输入数据的间距），再从R=0开始按序计算，直到R=B-1结束。每个R对应的旋转因子指数，旋转因子指数相同的蝶从上往下依次逐个运算，各个蝶的第一节点标号k都是从R开始，以为步长，到（可简取极值N-2）结束。考虑到蝶形运算有两个输出，且都要用到本级的两个输入数据，故第一个输出计算完毕后，输出数据不能立即存入输入地址，要等到第二个输出计算调用输入数据完毕后才能覆盖。这样数据倒序后的运算可用三重循环程序实现。整个蝶形运算流程图如图3-3所示。送入x(n),MN=2

17、M倒序L=1:MB=2(L-1)J=0:B-1P=J*2(M-L) K=J:2L:NT=A(K)+A(K+B)*WNPA(K+B)=A(K)-A(K+B)*WNPA(K)=T输出开始 结束束束图3-3整个蝶形运算流程图3.序列倒序为了保证运算输出的X(k)按顺序排列，要求序列x(n)倒序输入，即在运算前要先对输入的序列进行位序颠倒。如果总点数为的x(n)的顺序数是用M位二进制数表示，则倒序数只需将顺序数的二进制位倒置即可，按照这一规律用硬件电路和汇编语言很容易产生倒序数。但用MATLAB等高级语言实现倒序时，直接倒置二进制数位的方法不可取，还须找出产生倒序的十进制规律。将十进制顺序数用I表示，

18、与之对应的二进制数用IB表示。十进制倒序数用J表示，与之对应的二进制数用JB表示。JB是IB的位倒置结果，十进制顺序数I增加1，相当于IB最低位加1且逢2向高位进1，即相当于JB最高位加1且逢2向低位进1。JB的变化规律反映到J的变化分二种情况：如果JB的最高位是0，则直接由加1得到下一个倒序值；如果JB的最高位是1，则要先将最高位变0，再在次高位加1。但次高位加1时，同样要判断0、1值，如果是0 ，则直接加1，否则要先将次高位变0，再判断下一位。依此类推，直到完成最高位加1，逢2向右进位的运算。利用这一算法可按顺序数I的递增顺序，依次求得与之对应的倒序数J。为了节省内存，数据倒序可原址进行，

19、当I = J时不需要交换，当I J时需要交换数据。另外，为了避免再次调换前面已经调换过的一对数据，只对IGUIDE”；(2) 在“Start”菜单中选择“MATLAB”下的“GUIDE”命令；(3) 在工具栏中单击“GUIDE”的图标，进入GUI默认窗口界面，如图5-1所示。这里选择空白界面类型(Black GUI)，单击“OK”按钮，MATLAB将启动GUIDE，如图5-2所示。图5-1 GUIDE的启动界面图5-2具有空白界面的GUIDE在GUI空白界面中，位于中央的深灰色部分为绘制控件的画布，拖动右下角小黑点可以调整画布尺寸的大小。在GUI界面的左侧为MATLAB控制面板。控制面板的外观

20、可以通过设置GUIDE的属性进行简要的修改，选择GUIDE中的“File”菜单下的“Preference”命令，在弹出的对话框中选择“Show names in Component palette”复选框，如图5-3所示。单击“OK”按钮后，控制面板中在不同的空间旁边会显示相应空间的名称，如图5-4所示。 图5-3 GUIDE的属性对话框 图5-4显示控件名称控制面板上有各种控件，其中有推按钮(Push Button)、单按钮(Radio Button)、复选框(Checkbox)、文本框(Edit Text)、静态文本框(Static Text)、下拉框(Popup Menu)、列表框(Li

21、stbox)、单选按钮(Toggle Button)、轴对象(axes)等。从控制面板中选择上述控件，按照一定的界面布局用鼠标将控件拖动到画布上。5.3 GUI实现程序以本设计要求为例介绍。第一步，该选择本图形用户界面需要的控件：九个推按钮(Push button)，分别为原图，灰度图，内置FFT图，内置IFFT图，原图，灰度图，自编FFT图，自编IFFT图，退出。八个轴对象(axes)用来显示原图，灰度图，内置FFT图，内置IFFT图，原图，灰度图，自编FFT图，自编IFFT图。完成的GUI界面如图5-5所示：图5-5完成的GUI界面第二步，设置控件属性：双击组件可以设置推按钮的属性，如显示

22、大小，名称和默认值等。上图就设置了推按钮的名称，这样可以更清楚的明白每个推按钮的功能。第三步，编写回调函数。组件事件的发生是通过回调函数进行工作的。控件设置完成后保存，然后运行GUI（操作为Ctrl+T），就会进入Editor窗口，加入各个控件功能的函数代码。完成后保存即可。第四步，运行GUI。运行Editor窗口的程序后，会弹出已经激活的GUI界面。点击代表各个函数操作的按钮，就会出现进行了相应函数操作的图片。如图5-6所示：图5-6运行GUI后的结果5.4运行调试当点击推按钮，有时相应的图片将不会出现，提示出现错误，这时极有可能是图片路径设置错误，则将回调函数中图片路径修改一下，就可出现正

23、确的结果。6自编算法与内置算法结果比较我们知道MATLAB软件自带FFT算法和IFFT算法，我们可以通过比较自编算法运行结果与内置算法运行结果来检验自编算法的正确性。通过观察图4-3和图5-6的运行结果可知，经过自编FFT运算得到的图形和内置FFT得到的图形大体一致，经过自编IFFT运算得出的图形和内置IFFT得出的图形一致。只是内置的FFT得到的图片比自编得到的FFT图片更清晰，其他相差不多，大致可认为两种运算结果是一致的。IFFT算法是在对图像进行FFT算法处理的基础上对图片进行的快速傅里叶逆变换，因此可以观察到图片进行IFFT算法处理后，又恢复为灰度图，而且较之前的灰度图清晰度较差。本次

24、课程设计编写的程序严格按程序框图编写，思路清晰、容易理解，程序的运行过程在命令窗中一目了然。通过与自编函数运算的结果比对，虽存在一点小误差，但在允许误差之内，所以程序编写正确。7总结本次实习的主要内容是通过用MATLAB实现FFT的设计，可以实现对选定图片信号进行分析，并做出经FFT算法及IFFT算法后的图片结果。把自己编写的FFT算法与MATLAB自带FFT算法进行比较。程序运行调试时，自己选择输入要采样的图片，可以实现对不同图片的FFT运算。在之前数字信号处理的学习以及完成实验的过程中，已经使用过MATLAB，对其有了一些基础的了解和认识，通过这次的课程设计使我进一步了解了信号的产生，采样

25、及频谱分析的方法，以及其中产生信号和绘制信号的基本命令和一些基础编程语言。让我感受到只有在了解课本知识的前提下，才能更好的应用这个工具，并且熟练的应用MATLAB也可以很好的加深我对课程的理解，方便我的思维。这次课程设计使我了解了MATLAB的使用方法，提高了自己的分析和动手实践能力。同时我相信，进一步加强对MATLAB的学习与研究对我今后的学习将会起到很大的帮助。这次的课程设计是对本学期所学知识的一次重要巩固，使得在课堂上掌握的知识得到了真正的运用。在学习的过程中和同学讨论，更明白了理论知识与实践的联系。书到用时方恨少，有些知识学会是一回事，掌握是一回事，但应用起来，确实不是那么简单的，需要

26、很多知识的融会贯通。程序运行调试初期，曾经多次出现错误、不能产生图形等问题，但在我翻阅资料认真改正及老师同学的帮助下基本功能还是完成了，经过1个星期的上机实习，程序已得到一些完善，能完成基本的要求的功能。最后经过努力，又深入学习了图形用户界面(GUI)，完成了选做要求的GUI界面。学习就是一个了解，疑惑，进而解惑的过程，这次实习就是提供了这样一个发现自己知识漏洞，与同学老师探讨进行解惑的的机会。通过这次课程设计实习，我更深刻的了解了MATLAB的运用，重新复习了FFT 中的重要的序列倒序和蝶形变换的程序，对课本上的知识有了更深的理解，使我对数字信号处理有了系统的认知。在这里特别感谢老师，他们给

27、了我们很大的发挥空间，让我们真正自己动手真正掌握了知识，感谢他们细心指导。也非常感谢我的同学，他们解开了我在实习中出现的诸多知识死角，谢谢大家！参考文献1 范寿康.DSP 技术与DSP芯片.北京电子工业出版社，20082 程佩青.数字信号处理教程.北京清华大学出版社出版，20013 高西全，丁玉美等.数字信号处理.北京电子工业出版社，20094 余成波，陶红艳.数字信号处理及MATLAB 实现.北京清华大学出版社，2008 5 曹弋，赵阳.MATLAB 实用教程.北京电子工业出版社，2007 6 奥本海姆.离散时间信号处理.科学出版社，20007 宗孔德，胡广书.数字信号处理.清华大学出版社，

28、1997附录以下是用Matlab实现的程序编码：function image_process_FFT()filename, pathname=uigetfile(*.jpg;*.tif;*.bmp;*.gif ,File Selector);image=imread(strcat(pathname,filename);scrsz=get(0,ScreenSize);figure(position,0 0 scrsz(3)-1 scrsz(4);set(gcf,Name,快速傅里叶变换);subplot(2,4,1);imshow(image);title(原始图像);subplot(2,4,5

29、);imshow(image);title(原始图像);if ndims(image)=3 image=rgb2gray(image);endsubplot(2,4,2);imshow(image);title(灰度图像);subplot(2,4,6);imshow(image);title(灰度图像)r,c=size(image);array=image;t=log2(r);t1=floor(t);t2=ceil(t);if t1=t2 array(2t2,c)=0;endr1,c1=size(array);t=log2(c1);t3=floor(t);t4=ceil(t);if t3=t4

30、 array(r1,2t4)=0;endr1,c1=size(array);n=r1/2;data_col=zeros(1,n,double);for m=1:n data_col(m)=exp(-1i*2*pi*(m-1)/r1);endn=c1/2;data_row=zeros(1,n,double);for m=1:n data_row(m)=exp(-1i*2*pi*(m-1)/r1);endarray=transform_fft2(array);Ft=fftshift(array);S1=log(1+abs(Ft);subplot(2,4,7);imshow(S1,);title(自

31、建FFT2函数结果);array=transform_ifft2(array);array=abs(array);array=array(1:r,1:c);subplot(2,4,8);imshow(array,);title(自建IFFT2函数结果);F=fft2(image);FC=fftshift(F);S=log(1+abs(FC);subplot(2,4,3)imshow(S,);title(内置FFT2函数结果);array=ifft2(F);array=round(abs(array);subplot(2,4,4);imshow(array,);title(内置IFFT2函数结果

32、);returnfunction array=transform_fft2(array)array=double(array);r1 c1=size(array);for j=1:r1 array(j,:)=transform_fft(array(j,:);endfor j=1:c1 array(:,j)=transform_fft(array(:,j);endfunction array1=transform_fft(array)N=length(array);n=N/2;w=zeros(1,n,double); for m=1:n w(m)=exp(-1i*2*pi*(m-1)/N);en

33、d p=log2(N);array1=zeros(1,N,double);for q=1:p t1=2(q-1); t2=2(p-1); for k=0:(t2/t1-1) for j=0:(t1-1) if mod(q,2)=1 data1=array(k*t1+j+1); data2=array(k*t1+j+t2+1); array1(k*t1*2+j+1)=data1+data2; array1(k*t1*2+j+t1+1)=(data1-data2)*w(k*t1+1); else data1=array1(k*t1+j+1); data2=array1(k*t1+j+t2+1);

34、array(k*t1*2+j+1)=data1+data2; array(k*t1*2+j+t1+1)=(data1-data2)*w(k*t1+1); end end endendif mod(p,2)=1 return else array1=array; returnendfunction array=transform_ifft2(array) array=conj(array);r1,c1=size(array);for i=1:r1 array(i,:)=transform_fft(array(i,:);endfor i=1:c1 array(:,i)=transform_fft(

35、array(:,i);endarray=array/(r1*c1);附录以下是用GUI界面实现的程序编码：function pushbutton1_Callback(hObject, eventdata, handles)axes(handles.axes1);x=imread(sxy.jpg); imshow(x);title(原图);function pushbutton2_Callback(hObject, eventdata, handles)axes(handles.axes2);x=imread(sxy.jpg);f=rgb2gray(x);imshow(f);title(灰度图)

36、;function pushbutton3_Callback(hObject, eventdata, handles)axes(handles.axes3);img=imread(sxy.jpg);f=rgb2gray(img);F=fft2(f);FS=fftshift(F);S=log(1+abs(FS);imshow(S,);title(FFT图);function pushbutton4_Callback(hObject, eventdata, handles)axes(handles.axes4);img=imread(sxy.jpg);f=rgb2gray(img);F=fft2(

37、f);FS=fftshift(F); S=log(1+abs(FS);fr=real(ifft2(ifftshift(FS); ret=im2uint8(mat2gray(fr);imshow(ret,);title(IFFT图);function pushbutton5_Callback(hObject, eventdata, handles)axes(handles.axes5);x=imread(sxy.jpg); imshow(x);title(原图);function pushbutton6_Callback(hObject, eventdata, handles)axes(hand

38、les.axes6);image=imread(sxy.jpg);if ndims(image)=3 image=rgb2gray(image);endimshow(image);title(灰度图);function pushbutton7_Callback(hObject, eventdata, handles)axes(handles.axes7);image=imread(sxy.jpg);if ndims(image)=3 image=rgb2gray(image);endr,c=size(image);array=image;t=log2(r);t1=floor(t);t2=cei

39、l(t);if t1=t2 array(2t2,c)=0;endr1,c1=size(array);t=log2(c1);t3=floor(t);t4=ceil(t);if t3=t4 array(r1,2t4)=0;endr1,c1=size(array);n=r1/2;data_col=zeros(1,n,double);for m=1:n data_col(m)=exp(-1i*2*pi*(m-1)/r1);endn=c1/2;data_row=zeros(1,n,double);for m=1:n data_row(m)=exp(-1i*2*pi*(m-1)/r1);endarray=

40、transform_fft2(array);Ft=fftshift(array);S1=log(1+abs(Ft);imshow(S1,);title(FFT图);function array=transform_fft2(array)array=double(array);r1 c1=size(array); for j=1:r1 array(j,:)=transform_fft(array(j,:);endfor j=1:c1 array(:,j)=transform_fft(array(:,j);endfunction array1=transform_fft(array)N=lengt

41、h(array);n=N/2;w=zeros(1,n,double); for m=1:n w(m)=exp(-1i*2*pi*(m-1)/N);endp=log2(N);array1=zeros(1,N,double);for q=1:p t1=2(q-1); t2=2(p-1); for k=0:(t2/t1-1) for j=0:(t1-1) if mod(q,2)=1 data1=array(k*t1+j+1); data2=array(k*t1+j+t2+1); array1(k*t1*2+j+1)=data1+data2; array1(k*t1*2+j+t1+1)=(data1-

42、data2)*w(k*t1+1); else data1=array1(k*t1+j+1); data2=array1(k*t1+j+t2+1); array(k*t1*2+j+1)=data1+data2; array(k*t1*2+j+t1+1)=(data1-data2)*w(k*t1+1); end end endendif mod(p,2)=1 return else array1=array; returnendfunction pushbutton8_Callback(hObject, eventdata, handles)axes(handles.axes8);image=im

43、read(sxy.jpg);if ndims(image)=3 image=rgb2gray(image);endr,c=size(image);array=image;t=log2(r);t1=floor(t);t2=ceil(t);if t1=t2 array(2t2,c)=0;endr1,c1=size(array);t=log2(c1);t3=floor(t);t4=ceil(t);if t3=t4 array(r1,2t4)=0;endr1,c1=size(array);n=r1/2;data_col=zeros(1,n,double);for m=1:n data_col(m)=e

44、xp(-1i*2*pi*(m-1)/r1);endn=c1/2;data_row=zeros(1,n,double); for m=1:n data_row(m)=exp(-1i*2*pi*(m-1)/r1);endarray=transform_fft2(array);Ft=fftshift(array);S1=log(1+abs(Ft);array=transform_ifft2(array);array=abs(array);array=array(1:r,1:c);imshow(array,);title(IFFT图);function array=transform_fft2(arr

45、ay)array=double(array);r1 c1=size(array)for j=1:r1 array(j,:)=transform_fft(array(j,:);endfor j=1:c1 array(:,j)=transform_fft(array(:,j);endfunction array1=transform_fft(array)N=length(array);n=N/2;w=zeros(1,n,double); for m=1:n w(m)=exp(-1i*2*pi*(m-1)/N);endp=log2(N);array1=zeros(1,N,double);for q=

46、1:p t1=2(q-1); t2=2(p-1); for k=0:(t2/t1-1) for j=0:(t1-1) if mod(q,2)=1 data1=array(k*t1+j+1); data2=array(k*t1+j+t2+1); array1(k*t1*2+j+1)=data1+data2; array1(k*t1*2+j+t1+1)=(data1-data2)*w(k*t1+1); else data1=array1(k*t1+j+1); data2=array1(k*t1+j+t2+1); array(k*t1*2+j+1)=data1+data2; array(k*t1*2

47、+j+t1+1)=(data1-data2)*w(k*t1+1); end end endendif mod(p,2)=1 return else array1=array; returnendfunction array=transform_ifft2(array)array=conj(array);r1,c1=size(array);for i=1:r1 array(i,:)=transform_fft(array(i,:);endfor i=1:c1 array(:,i)=transform_fft(array(:,i);endarray=array/(r1*c1);function p

48、ushbutton9_Callback(hObject, eventdata, handles)clc;close all;close(gcf);附录资料：不需要的可以自行删除各类滤波器的MATLAB程序理想低通滤波器IA=imread(lena.bmp);f1,f2=freqspace(size(IA),meshgrid);Hd=ones(size(IA);r=sqrt(f1.2+f2.2);Hd(r0.2)=0;Y=fft2(double(IA);Y=fftshift(Y);Ya=Y.*Hd;Ya=ifftshift(Ya);Ia=ifft2(Ya);figuresubplot(2,2,1

49、),imshow(uint8(IA);subplot(2,2,2),imshow(uint8(Ia);figuresurf(Hd,Facecolor,interp,Edgecolor,none,Facelighting,phong); 二、理想高通滤波器IA=imread(lena.bmp);f1,f2=freqspace(size(IA),meshgrid);Hd=ones(size(IA);r=sqrt(f1.2+f2.2);Hd(r0.2)=0;Y=fft2(double(IA);Y=fftshift(Y);Ya=Y.*Hd;Ya=ifftshift(Ya);Ia=real(ifft2(

50、Ya);figuresubplot(2,2,1),imshow(uint8(IA);subplot(2,2,2),imshow(uint8(Ia);figuresurf(Hd,Facecolor,interp,Edgecolor,none,Facelighting,phong); Butterworth低通滤波器IA=imread(lena.bmp);f1,f2=freqspace(size(IA),meshgrid);D=0.3;r=f1.2+f2.2;n=4;for i=1:size(IA,1) for j=1:size(IA,2) t=r(i,j)/(D*D); Hd(i,j)=1/(t

51、n+1); endendY=fft2(double(IA);Y=fftshift(Y);Ya=Y.*Hd;Ya=ifftshift(Ya);Ia=real(ifft2(Ya);figuresubplot(2,2,1),imshow(uint8(IA);subplot(2,2,2),imshow(uint8(Ia);figuresurf(Hd,Facecolor,interp,Edgecolor,none,Facelighting,phong); Butterworth高通滤波器IA=imread(lena.bmp);f1,f2=freqspace(size(IA),meshgrid);D=0.

52、3;r=f1.2+f2.2;n=4;for i=1:size(IA,1) for j=1:size(IA,2) t=(D*D)/r(i,j); Hd(i,j)=1/(tn+1); endendY=fft2(double(IA);Y=fftshift(Y);Ya=Y.*Hd;Ya=ifftshift(Ya);Ia=real(ifft2(Ya);figuresubplot(2,2,1),imshow(uint8(IA);subplot(2,2,2),imshow(uint8(Ia);figuresurf(Hd,Facecolor,interp,Edgecolor,none,Facelighting

53、,phong); 高斯低通滤波器IA=imread(lena.bmp);IB=imread(babarra.bmp);f1,f2=freqspace(size(IA),meshgrid);D=100/size(IA,1);r=f1.2+f2.2;Hd=ones(size(IA);for i=1:size(IA,1) for j=1:size(IA,2) t=r(i,j)/(D*D); Hd(i,j)=exp(-t); endendY=fft2(double(IA);Y=fftshift(Y);Ya=Y.*Hd;Ya=ifftshift(Ya);Ia=real(ifft2(Ya);figures

54、ubplot(2,2,1),imshow(uint8(IA);subplot(2,2,2),imshow(uint8(Ia);figuresurf(Hd,Facecolor,interp,Edgecolor,none,Facelighting,phong); 高斯高通滤波器IA=imread(lena.bmp);IB=imread(babarra.bmp);f1,f2=freqspace(size(IA),meshgrid);%D=100/size(IA,1);D=0.3;r=f1.2+f2.2;for i=1:size(IA,1) for j=1:size(IA,2) t=r(i,j)/(D

55、*D); Hd(i,j)=1-exp(-t); endendY=fft2(double(IA);Y=fftshift(Y);Ya=Y.*Hd;Ya=ifftshift(Ya);Ia=real(ifft2(Ya);figuresubplot(2,2,1),imshow(uint8(IA);subplot(2,2,2),imshow(uint8(Ia);figuresurf(Hd,Facecolor,interp,Edgecolor,none,Facelighting,phong); 梯形低通滤波器IA=imread(lena.bmp);IB=imread(babarra.bmp);f1,f2=f

56、reqspace(size(IA),meshgrid);%D=100/size(IA,1);D0=0.1;D1=0.4;r=sqrt(f1.2+f2.2);Hd=zeros(size(IA);Hd(r=D0 & r(i,j)=D1 Hd(i,j)=(D1-r(i,j)/(D1-D0); end endendY=fft2(double(IA);Y=fftshift(Y);Ya=Y.*Hd;Ya=ifftshift(Ya);Ia=real(ifft2(Ya);figuresubplot(2,2,1),imshow(uint8(IA);subplot(2,2,2),imshow(uint8(Ia);figuresurf(Hd,Facecolor,interp,Edgecolor,none,Facelighting,phong); 梯形高通滤波器IA=imread(lena.bmp);IB=imread(babarra.bm