基于MATLAB的FFT算法的设计

1、目 录1前言12FFT算法的基本原理22.1系统总体流程图22.2 FFT运算规律及编程思想2对图片的选择22.2.2 FFT算法的基本原理32.2.3 FFT算法的运算规律及编程思想43 软件简介53.1 Matlab简介53.1.1 Matlab软件概况53.1.2 Matlab的特点53.2 GUI简介6界面设计63.3对比结果与分析84心得体会10参考文献11附录 Matlab源程序12附录 GUI源程序161前言随着信息时代，数字时代的到来，数字信号处理已经成为一门极其重要的学科和技术领域。以DSP为核心芯片的处理系统日益变成了数字信号处理系统的主流。它广泛用于电子信息、通信、图像处

2、理、语音处理、生物医学、自动控制、地质探测等领域，受到工程设计和使用人员的青睐。MATLAB，它是美国Math Works公司推出的一种面向工程和科学计算的交互式计算软件。它以矩阵运算为基础，把计算、可视化、程序设计融合在一个简单易用的交互式工作环境中，是一款数据分析和处理功能都非常强大的工程适用软件。通过本次课设我们学会了分析和处理音频信号，首先要对图片信息进行采集，MATLAB的数据采集工具箱提供了一整套命令和函数，通过调用这些函数和命令，可直接控制图像进行数据采集。Window自带的程序也可驱动采集图片信息，并能保存该文件，供MATLAB相关函数直接读取写入。MATLAB语言是一种数据分

3、析和处理功能十分强大的计算机应用软件，它可以将图像文件变换位离散的数据文件，然后利用其强大的矩阵运算能力处理数据，如数据滤波、傅立叶变换、时域和频域分析、声音回放以及各种图的呈现等，它的信号处理与分析工具箱位语音信号分析提供了十分丰富的功能函数，利用这些功能函数可以快捷而又方便的完成图像信号的处理和分析以及信号的可视化，是人机交互更加便捷。信号处理是MATLAB重要应用的领域之一。对于有限长序列x(n),若要求其N点的傅里叶变换DFT需要经过次复数乘法运算和N*(N-1)次复数加法运算。随着N的增加，运算量将急剧增加，而在实际问题中，N往往是较大的，如当N=1024时，完成复数乘法和复数加法的

4、次数分别为百万以上，无论是用通用计算机还是用DSP芯片，都需要消耗大量的时间和机器内存，不能满足实时的要求。因此，DFT的这种运算只能进行理论上的计算，不适合对实时处理要求高的场合。因此，研究作为DSP的快速算法的FFT是相当必要的，快速傅里叶变换FFT是为提高DFT运算速度而采用的一种算法，快速算法的种类很多，而且目前仍在改进和提高，它是根据离散傅里叶变换的奇、偶、虚、实等特性，对离散傅立叶变换的算法进行改进获得的。基于本学期所学的DIT-FFT的运算规律和编程思想以及Matlab的学习和使用，本课设要求在Matlab环境下编写基2 DIT-FFT算法实现对离散信号的快速傅里叶变换，再与Ma

5、tlab软件自带的FFT函数实现对离散信号的傅里叶变换进行比较，如果得到的频谱相同，那么我们编写的程序就是正确的。用GUI界面完成人机交互方便使用的。本课程设计主要是对数字信号的分析。2FFT算法的基本原理2.1系统总体流程图本设计要求选择一个图片，并对该图像进行FFT变换。在Matlab环境下编写基2 DIT-FFT算法；利用自己编写的算法对选取的图片进行分析，并与Matlab数字信号处理工具箱中的fft函数进行对比研究，验证自编算法的正确性。所以得到系统总体流程图如下图所示。图片选择编写FFT程序，得到所需的图片对图片的路径选择对图片进行灰度变换实现选择图片的倒序一级中不同蝶形运算一级中相

6、同蝶形运算与MATLAB自带的比较图2-1 系统总体流程图2.2 FFT运算规律及编程思想2.2.1对图片的选择保存一张图片，格式为bmp，jpg或者是gif中的任何一种。并将该图片文件保存在电脑的某个盘中。filename, pathname=uigetfile('*.jpg;*.tif;*.bmp;*.gif' ,'File Selector');image=imread(strcat(pathname,filename);这些代码实现了对图片的路径进行的选择，只要找到保存图片的位置即可打开图片。if ndims(image)=3 image=rgb2gra

7、y(image);end这些代码实现了对图片进行灰度变换，可以将原来是彩色的图片变换成黑白色，目的是为了取点时方便简单，而且速度更快，节省时间 。2.2.2 FFT算法的基本原理快速傅里叶变换FFT是为提高DFT运算速度而采用的一种算法。对一个有限长度序列x(n)的N点的DFT为：所以，要求N点的DFT,需要N2次的复数乘法运算,N*(N-1)次复数乘法运算算。随着N的增加，运算量将急剧增加，而在实际问题中，N往往是较大的，因此无论是用通用计算机还是用DSP芯片，都需要消耗大量的时间，不能满足实时的要求,，不适合于对实时处理要求高的场合。为了能实时处理DFT,要想减少DFT的运算量可以有两个途

8、径：第一是降N，N的值减小了，运算量就减少了；第二是利用旋转因子的周期性，对称性和可约性。利用这两个途径实现DFT的快速傅里叶变换FFT，FFT算法基本上可分为按时间抽取的FFT算法(DIT-FFT)和按频率抽取的FFT算法(DIF-FFT)。旋转因子的性质：(1)周期性(2)共轭对称性(3)可约性本次课设要求用用基2的按时间抽取的FFT算法（DIT-FFT）实现FFT功能，设序列x(n)的长度为N，且N满足N=2M,M为正整数。若N不能满足上述关系，可以将序列x(n)补零实现。按时间抽取基2-FFT算法的基本思路是将N点序列按时间下标的奇偶分为两个N/2点序列，计算这两个N/2点序列的N/2

9、点DFT；每一个N/2点序列按照同样的划分原则，可以划分为两个N/4点序列，最后，将原序列划分为多个2点序列，将计算量大大降低。按时间下标的奇偶将N点x(n)分别抽取组成两个N/2点序列，分别记为x1(n)和x2(n)，将x(n)的DFT转化为x1(n)和x2(n)的DFT的计算。利用旋转因子的可约性，即：用蝶形运算图可表示为：图2-2 蝶形运算图以此类推，还可以把x1(n)和x2(n)按n值得奇偶分为两个序列，这样就达到了降N得目的，从而减少了运算量。FFT对DFT的数学运算量改进：直接采用DFT进行计算，运算量为N2次复数乘法和N*(N-1)次复数乘法。当采用M次FFT时，由N=2M求得M

10、=logN，运算流图有M级蝶形，每一级都由N/2个蝶形运算构成，这样每一级蝶形运算都需要N/2次复数乘法和N次复数加法。M级运算共需要复数乘法次数为C=N/2*M,复数加法次数为C=N*M。当N值较大时，FFT减少运算量的特点表现的越明显。 FFT算法的运算规律及编程思想实现FFT运算的核心是蝶形运算，找出蝶形运算的规律是编程的基础。对点的FFT共有M级运算，用L表示从左到右的运算级数(L=1，2，M )。第L级共有个不同指数的旋转因子，用R表示这些不同指数旋转因子从上到下的顺序(R=0，1，B-1)。第R个旋转因子的指数，旋转因子指数为P的第一个蝶的第一节点标号k从R开始，由于本级中旋转因子

11、指数相同的蝶共有个，且这些蝶的相邻间距为，故旋转因子指数为P的最后一个蝶的第一节点标号k为：，本级中各蝶的第二个节点与第一个节点都相距B点。总结上述运算规律，可采用如下运算方法进行DIT-FFT运算。首先读入数据，根据数据长度确定运算级数M，运算总点数，不足补0处理。然后对读入数据进行数据倒序操作。运算时，先算出该级不同旋转因子的个数(也是该级中各个蝶形运算两输入数据的间距)，再从R=0开始按序计算，直到R=B-1结束。每个R对应的旋转因子指数。3 软件简介3.1 Matlab简介70年代后期,身为美国New Mexico大学计算机系系主任的Cleve Moler发现学生用

12、FORTRAN编写接口程序很费时间,于是他开始自己动手,利用业余时间为学生编写EISPACK和LINPACK的接口程序。Cleve Moler给这个接口程序取名为Matlab。1984年，为了推广Matlab在数值计算中的应用，Cleve Moler、Johon Little等正式成立了Math works公司，从而把Matlab推向市场，并开始了对Matlab工具相等的开发设计。 Matlab软件概况Matlab是MatrixLaboratory的缩写，意为矩阵实验室。它具有强大的矩阵处理功能和绘图功能，进还能进行文字处理，绘图，建模仿真等功能。随着版本的不断升级，它在数值计算及符

13、号计算功能上得到了进一步完善。Matlab已经发展成为多学科、多种工作平台的功能强大的大型软件。在欧美等高校，Matlab已经成为线性代数、自动控制理论、概率论及数理统计、数字信号处理、时间序列分析、动态系统仿真等高级课程的基本教学工具。Matlab软件打开后，其展现的窗口图如下：图3-1 运行Matlab窗口图 Matlab的特点Matlab有以下一些特点：Matlab的帮助功能很强大，自带有详细的帮助手册，基于HTML的完整的帮助功能，也可以用help命令来得到帮助信息。程序语法设计自由度大，方便我们编程。还有一个原因使Matlab受人们欢迎的，那就是Matlab源程序具有很大的开放性。除

14、了内部函数以外，所有Matlab的核心文件和工具箱文件都是可读可改的源文件，用户可通过对源文件的修改以及加入自己的文件构成新的工具箱。Matlab有强大的的图形绘制功能。在Matlab里，数据可视化的操作非常简单易用。Matlab还有较强的编辑图形界面的能力。可以用来声成图解和可视化的二维、三维图。Matlab还拥有功能强大的各种工具箱。其工具箱分为两类：功能性工具箱和学科性工具箱。功能性工具箱主要用来扩充其符号计算功能，图示建模仿真功能，文字处理功能以及与硬件实时交互功能。功能性工具箱用于多种学科。而学科性工具箱是专业性比较强的，如(control、signal proceessing 、c

15、ommumnication) toolbox等。这些工具箱都是由该领域内学术水平很高的专家编写的，所以用户无需编写自己学科范围内的基础程序，而直接进行高，精，尖的研究，能极大地促进我们的学习研究工作。虽然Matlab有很多优点，但它也有一些缺点，比如：由于Matlab的程序不用编译等预处理，也不生成可执行文件，程序为解释执行，所以速度较慢。3.2 GUI简介图形用户界面(GUI)，是一种提供人机交互的工具和方法。GUI是包含图形对象，如窗口、图标、菜单和文本等图文并茂的用户界面。一个设计非常优秀的GUI能够非常直观的使用户知道如何操作Matlab界面，了解设计者开发的意图。界面设计打开GUI后

16、未设置任何属性的界面如下图所示：图3-2 GUI界面用MATLAB图形用户界面开发环境设计GUI的一般步骤是：1.进行界面设计。2.设计控件属性。3.进行M语言编程。以本设计要求为例介绍。第一步，该选择本图形用户界面需要的控件：八个推按钮(Push button)，用来运行和退出。六个轴对象(axes)用来显示信号时域波形，频谱图，和两种不同FFT实现的频谱图，即运用Matlab自带的。第二步，设置各个对象的属性。如果要设置各个对象在表面上的寓意，只要改变它的string即可。Tag是它的本质属性，是唯一的。控制属性完成后，图形用户界面如图所示。图3-3 控件属性设置完成后的图形用户界面第三步

17、，编写回调函数。组件事件的发生是通过回调函数进行工作的。回调函数是用户指定组件需要完成一定动作的函数，如本次设计中单击按钮后的功能就是由编写回调函数实现的。一个组件中可以包括多个回调函数，但是几乎每一个组件都有的属性是Callback属性。GUIDE可创建一个M回调函数文件的构架，文件将自动处理并将所有文件的句柄传递到handles结构数组中。Handles是GUI中一个特殊的结构数组，它是GUI中所有组件共用的一个结构。在handles结构数组中还可以添加用户自定义的，在不同的回调函数之间共享的数据。回调函数的声明格式为：Function object_Callback(hObject，ev

18、entData，handles)声明格式各部分为：object为发生事件控制的Tag属性字符串。hObject为发生事件的控制句柄。eventData为保留字段。本设计中各个按钮回调函数完成3个任务：1得到用户输入图片。例如当点击原图按钮时，就会得到图片路径选择路径，得到图片。图3-4 得到图片2计算数据。3建立图形。3.3对比结果与分析得到原图的两个图形对比如下：图3-5 原图对比由以上两图对比可以看出，自建的原图和自带的原图得到的结果并不是完全一样的，在色彩上有差别，明暗程度不一样，自建的比自带的颜色更鲜明，更深。图3-6 FFT对比图由以上两图对比可以看出，自建的fft和系统自带的iff

19、t也是有细微的差别的，自建的得到的图片比系统自带的图像颜色更甚而且由中间汇聚的光束可以看出，自建的函数得到的结果光束比较大，发散的范围更广，而系统自带的函数得到的结果相对来说比较窄。图3-7 IFFT对比图由以上两图对比可以看出，自建的IFFT和系统自带的IFFT总体上是一样的，但是有细微的差别，自建的相对来说颜色较暗，更深而自带的得到的图像颜色更亮更鲜明。4心得体会通过本课程设计使我了解了数字图像的基本概念，掌握数字图像处理的基本内容，如图像点运算、几何变换、增强处理、图像复原等的基本原理和Matlab实现方法。通过本次课程设计，让我们掌握如何学习一门语言，如何进行资料查阅搜集，如何自己解决

20、问题等方法，养成良好的学习习惯。扩展理论知识，培养学生的综合设计能力。本次课设是通过用Matlab实现FFT的设计，对图像信号进行分析，并画出采样信号的时域与频域图。把自己编写的FFT算法与Matlab自带FFT算法进行比较。程序运行调试时，自己选择输入要采样的图像信号，采样像素以及要变换的范围，可以实现对不同信号的信号采样和进行不同点的FFT运算。在之前数字信号处理的学习以及完成实验的过程中，已经使用过Matlab，对其有了一些基础的了解和认识，通过这次的课程设计使我进一步了解了信号的产生，采样及频谱分析的方法。让我感受到只有在了解课本知识的前提下，才能更好的应用这个工具，并且熟练的应用Ma

21、tlab也可以很好的加深我对课程的理解，方便我的思维。这次课程设计使我了解了Matlab的使用方法，提高了自己的分析和动手实践能力。同时我相信，进一步加强对MATLAB的学习与研究对我今后的学习将会起到很大的帮助。这次的课程设计是对本学期所学知识的一次重要巩固，使得在课堂上掌握的知识得到了真正的运用。在学习的过程中和同学讨论，更明白了理论知识与实践的联系。书到用时方恨少，有些知识学会是一回事，掌握是一回事，但应用起来，确实不是那么简单的，需要很多知识的融会贯通。程序运行调试初期，曾经多次出现错误、不能产生图形等问题，但在我翻阅资料认真改正及老师同学的帮助下基本功能还是完成了，经过1个星期的上机

22、实习，程序已得到一些完善，能完成基本的要求的功能。最后经过努力，又深入学习了图形用户界面(GUI)，完成了选做要求的人机对话界面。学习就是一个了解，疑惑，进而解惑的过程，这次课设就是提供了这样一个发现自己知识漏洞，与同学老师探讨进行解惑的的机会。通过这次课程设计，我更深刻的了解了Matlab的运用，重新复习了FFT 中的重要的变换的程序，对课本上的知识有了更深的理解，使我对数字信号处理有了系统的认知。在这里特别感谢魏老师和我身边的同学们，他们给了我们很大的发挥空间，让我们真正自己动手真正掌握了知识，感谢他们细心指导，他们解开了我在实习中出现的诸多知识死角，谢谢大家！参考文献1楼顺天,李博菡.基

23、于MATLAB的系统分析与设计信号处理.西安电子科技大学出版社，19982奥本海姆.离散时间信号处理.科学出版社，20003宗孔德,胡广书.数字信号处理.清华大学出版社，19974程佩青.数字信号处理教程.北京：清华大学出版社出版，20015高西全 丁玉美等.数字信号处理.北京：电子工业出版社，20096陈杰.Matlab宝典.电子工业出版社7苏金明 张莲花 刘波.MATLAB工具箱应用，电子工业出版社附录 Matlab源程序function image_process_FFT()filename, pathname=uigetfile('*.jpg;*.tif;*.bmp;*.gif

24、' ,'File Selector');image=imread(strcat(pathname,filename); if ndims(image)=3 image=rgb2gray(image);endscrsz=get(0,'ScreenSize');figure('position',0 0 scrsz(3)-1 scrsz(4);set(gcf,'Name','快速傅里叶变换');subplot(2,3,1);imshow(image);title('原始图像');subplot(

25、2,3,4);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 array(r1,2t4)=0;end r1,c1=size(array);n=r1/2;data_col=zeros(1,n,'double'); for m=1:n data_col(m)=exp

26、(-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,3,2);imshow(S1,);title('自建FFT结果');array=transform_ifft2(array);array=abs(array);array=array(1:r,1:c);s

27、ubplot(2,3,3);imshow(array,);title('自建IFFT结果');F=fft2(image);FC=fftshift(F);S=log(1+abs(FC);subplot(2,3,5)imshow(S,);title('内置FFT结果');array=ifft2(F);array=round(abs(array);subplot(2,3,6);imshow(array,);title('内置IFFT结果');returnfunction array=transform_fft2(array)array=double(a

28、rray);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');fo

29、r 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); array(k*t1*2+j+1)=data1+data2; array(k*

30、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(array(:,i);endarray=array/(r1*c1);附录 GUI

31、源程序function varargout = untitled3(varargin)gui_Singleton = 1;gui_State = struct('gui_Name', mfilename, . 'gui_Singleton', gui_Singleton, . 'gui_OpeningFcn', untitled3_OpeningFcn, . 'gui_OutputFcn', untitled3_OutputFcn, . 'gui_LayoutFcn', , . 'gui_Callback&

32、#39;, );if nargin && ischar(varargin1) gui_State.gui_Callback = str2func(varargin1);end if nargout varargout1:nargout = gui_mainfcn(gui_State, varargin:);else gui_mainfcn(gui_State, varargin:);endfunction untitled3_OpeningFcn(hObject, eventdata, handles, varargin)handles.output = hObject;gui

33、data(hObject, handles); function varargout = untitled3_OutputFcn(hObject, eventdata, handles) varargout1 = handles.output;function pushbutton1_Callback(hObject, eventdata, handles) filename, pathname=uigetfile('*.jpg;*.tif;*.bmp;*.gif' ,'File Selector');image=imread(strcat(pathname,f

34、ilename); if ndims(image)=3 image=rgb2gray(image);endscrsz=get(0,'ScreenSize');set(gcf,'Name','快速傅里叶变换');axes(handles.axes1);imshow(image);title('原图); function pushbutton2_Callback(hObject, eventdata, handles) filename, pathname=uigetfile('*.jpg;*.tif;*.bmp;*.gif'

35、 ,'File Selector');image=imread(strcat(pathname,filename); if ndims(image)=3 image=rgb2gray(image);endscrsz=get(0,'ScreenSize');set(gcf,'Name','快速傅里叶变换');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(c

36、1);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=transform_fft2(array);Ft=fftshift(ar

37、ray);S1=log(1+abs(Ft);axes(handles.axes2);imshow(S1,);title('自建FFT');function pushbutton3_Callback(hObject, eventdata, handles) filename, pathname=uigetfile('*.jpg;*.tif;*.bmp;*.gif' ,'File Selector');image=imread(strcat(pathname,filename); if ndims(image)=3 image=rgb2gray(im

38、age);end scrsz=get(0,'ScreenSize');set(gcf,'Name','快速傅里叶变换');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 array(r1,2t4)=0;end r1,c1=size(array);n=r1/2;data_col=zeros(1,n,&#

39、39;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); array=transform_ifft2(array);array=abs(array);array=array(1:r,1:c);axes(handl

40、es.axes3);imshow(array,);title('自建IFFT);function pushbutton4_Callback(hObject, eventdata, handles)filename, pathname=uigetfile('*.jpg;*.tif;*.bmp;*.gif' ,'File Selector');image=imread(strcat(pathname,filename);if ndims(image)=3 image=rgb2gray(image);endscrsz=get(0,'ScreenSize

41、');set(gcf,'Name','快速傅里叶变换');axes(handles.axes4);imshow(image);title('原图');function pushbutton5_Callback(hObject, eventdata, handles) filename, pathname=uigetfile('*.jpg;*.tif;*.bmp;*.gif' ,'File Selector');image=imread(strcat(pathname,filename); if ndims(

42、image)=3 image=rgb2gray(image);end scrsz=get(0,'ScreenSize');set(gcf,'Name','快速傅里叶变换');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 array(r1,2t4)=0;end r1,c1=size(array);

43、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); F=fft2(image);FC=fftshift(F);S=log(1+abs(FC);axes

44、(handles.axes5);imshow(S,);title('自带FFT');function pushbutton6_Callback(hObject, eventdata, handles) filename, pathname=uigetfile('*.jpg;*.tif;*.bmp;*.gif' ,'File Selector');image=imread(strcat(pathname,filename);if ndims(image)=3 image=rgb2gray(image);endscrszset(gcf,'Name','快速傅里叶变换');r,c=size(image); array=image;t=log2(r);t1=floor(t);t2=ceil(t);if t1=t2 arra