基于matlab的图像分析

目录

1 引言 (1)

2 基于MATLAB的FFT算法实现 (2)

2.1系统总体流程图 (2)

2.2 FFT运算规律及编程思想 (2)

2.2.1图像信号的采集 (2)

2.2.2 DIT-FFT算法的基本原理 (3)

2.2.3 FFT算法的运算规律及编程思想 (5)

3 Matlab程序实现 (7)

3.1程序运行结果 (7)

3.2对比结果分析 (8)

4 系统人机对话界面 (9)

4.1 GUI简介 (9)

4.2 界面设计 (9)

4.3 运行调试 (10)

5 Matlab软件简介 (11)

6 心得体会 (12)

参考文献 (13)

附录Ⅰ (14)

附录Ⅱ (18)

1 引言

MATLAB是矩阵实验室（Matrix Laboratory）的简称，是美国MathWorks 公司出品的商数学软件，用于算法开发、数据可视化、数据分析以及数值计算的高级技术计算语言和交互式环境，主要包括MATLAB和Simulink两大部分。

MATLAB 的应用范围非常广，包括信号和图像处理、通讯、控制系统设计、测试和测量、财务建模和分析以及计算生物学等众多应用领域。附加的工具箱（单独提供的专用MATLAB 函数集）扩展了MATLAB 环境，以解决这些应用领域内特定类型的问题。它以矩阵运算为基础，把计算、可视化、程序设计融合在一个简单易用的交互式工作环境中，是一款数据分析和处理功能都非常强大的工程适用软件。它可以将声音文件变换为离散的数据文件，然后利用其强大的矩阵运算能力处理数据，如数据滤波、傅立叶变换、时域和频域分析、声音回放以及各种图的呈现等，它的信号处理与分析工具箱位语音信号分析提供了十分丰富的功能函数，利用这些功能函数可以快捷而又方便的完成语音信号的处理和分析以及信号的可视化。数字信号处理是MATLAB重要应用的领域之一。

对于有限长序列x(n),若要求其N点的傅里叶变换（DFT）需要经过2N次复数乘法运算和N*(N-1)次复数加法运算。随着N的增加，运算量将急剧增加，而在实际问题中，N往往是较大的，如当N=1024时，完成复数乘法和复数加法的次数分别为百万以上，无论是用通用计算机还是用DSP芯片，都需要消耗大量的时间和机器内存，不能满足实时的要求。因此，DFT的这种运算只能进行理论上的计算，不适合对实时处理要求高的场合。因此，研究作为DSP的快速算法的FFT是相当必要的，快速傅里叶变换（FFT）是为提高DFT运算速度而采用的一种算法，快速算法的种类很多，而且目前仍在改进和提高，它是根据离散傅里叶变换的奇、偶、虚、实等特性，对离散傅立叶变换的算法进行改进获得的。基于本学期所学的DIT-FFT的运算规律和编程思想以及Matlab的学习和使用，本课设要求在Matlab环境下编写基2 DIT-FFT算法实现对离散信号的快速傅里叶变换，再与Matlab软件自带的FFT函数实现对离散信号的傅里叶变换进行比较，如果得到的频谱相同，那么我们编写的程序就是正确的。本次课程设计是实现对选定图片进行FFT计算、还原（IFFT计算），并与系统FFT函数做对比，进行分析。如果有能力可以选做系统人机对话界面。用GUI界面完成人机交互方便使用的。本课程设计主要是对数字信号的分析。

2 基于MATLAB的FFT算法实现

2.1系统总体流程图

本设计要求找到一张明暗对比较大的图片；在Matlab环境下编写基FFT算法；利用自己编写的算法对已选择的图片信号进行计算，并显示出计算的结果，将计算的结果与Matlab数字信号处理工具箱中自带的fft函数进行对比研究，验证自编算法的正确性。系统的总体设计流程图如图2-1所示：

图2-1 系统的总体设计流程图

2.2 FFT运算规律及编程思想

2.2.1图像信号的采集

图像信号最好采用明暗对比比较大的灰度图像进行分析，这样实验结果对比比较明显。

在Matlab中用语句：

[filename, pathname]=uigetfile({'*.jpg;*.tif;*.bmp;*.gif' },'File Selector');

image=imread(strcat(pathname,filename));

用于读取图片的信号，Matlab 图像分析支持多种格式的图像信号，用上述语句时，在Matlab 中分析图像的时候可以系统自动检索所需分析的图片。

语句：image=rgb2gray(image);可以对图像进行灰度处理。 当我们要将图片显示出来的的时候只需要用语句： imshow(image)；

本次课程设计就是分析灰度图像。通过用两种不同的方法对灰度图像的FFT 计算和IFFT 计算，来得到我们想要的结果。

采集到图像信号之后，就可以对图像信号进行分析和计算了。

2.2.2 DIT-FFT 算法的基本原理

快速傅里叶变换（FFT ）是为提高DFT 运算速度而采用的一种算法。 对一个有限长度序列x(n)的N 点的DFT 为：

所以，要求N 点的DFT,需要N2次的复数乘法运算,N*(N-1)次复数乘法运算算。随着N 的增加，运算量将急剧增加，而在实际问题中，N 往往是较大的，如当N=1024时，完成复数乘法和复数加法的次数分别为百万以上，无论是用通用计算机还是用DSP 芯片，都需要消耗大量的时间，不能满足实时的要求,，不适合于对实时处理要求高的场合。为了能实时处理DFT,要想减少DFT 的运算量可以有两个途径：第一是降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点DFT ，计算量可减小约一半；每一个N/2点序列按照同样的划分原则，可以划分为两个

)

()(N n k N n N k N kn N W W W ++==*)(*)(][][n k N n k N kn N W W W --==m

kn m

N kn N mkn mN kn N W W W W //,

==

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=logN ，运算流图有M 级蝶形，每一级都由N/2个蝶形运算构成，这样每一级蝶形运算都需要N/2次复数乘法和N 次复数加法。M 级运算共需要复数乘法次数为C=N/2*M,复数加法次数为C=N*M 。

()()()()()()(

)()()(

)1

N

021

N

N

0,2,4...1,3,5 (112)

2

212N

N

0,10,1112

2

2121

N

2

N

0,10,1

221N nk

n N N nk nk n n N N r k rk r r N N r k rk r r X k x n W x n W

x n W

x r W x r W x r W

x r W -=--==--+==--+====

+

=

+

+=

+

∑∑∑∑∑∑∑2j

2j 222

2

e

e

rk N rk

rk

rk

N

N N

W W π

π--===()()()1122

122

2

12,01N N rk k rk

N N N r r k N X k x r W W x r W X k W X k N --===+=+≤≤-∑∑（）（k）图2-2 DIT-FFT 蝶形运算流图符号

12

,

,1,0,)()12()()2(21-=?

??=+=N r r x r x r x r x

当N 值较大时，FFT 减少运算量的特点表现的越明显。

2.2.3 FFT 算法的运算规律及编程思想

为了编写DIT-FFT 算法的运算程序，首先要分析其运算规律，总结编程思想并绘出程序框图。

1. 原位计算

对M

N 2=点的FFT 共进行M 级运算，每级由N/2个蝶形运算组成。在同一级中，每个蝶的输入数据只对本蝶有用，且输出节点与输入节点在同一水平线上，这就意味着每算完一个蝶后，所得数据可立即存入原输入数据所占用的数组元素(存储单元)，这种原位(址)计算的方法可节省大量内存。

2. 蝶形运算

实现FFT 运算的核心是蝶形运算，找出蝶形运算的规律是编程的基础。蝶形运算是分级进行的；每级的蝶形运算可以按旋转因子的指数大小排序进行；如果指数大小一样则可从上往下依次蝶算。对M

N 2=点的FFT 共有M 级运算，用L 表示从左到右的运算级数(L=1，2，…，M )。第L 级共有12-=L B 个不同指数的旋转因子，用R 表示这些不同指数旋转因子从上到下的顺序(R=0，1，…，B-1)。第R 个旋转因子的指数R P L

M -=2

，旋转因子指数为P 的第一个蝶的第一节点标号

k 从R 开始，由于本级中旋转因子指数相同的蝶共有L

M -2个，且这些蝶的相邻间

距为L 2，故旋转因子指数为P 的最后一个蝶的第一节点标号k 为：

R N R L L L M +-=+?--22)12(，本级中各蝶的第二个节点与第一个节点都相距

B 点。

应用原位计算，蝶形运算可表示成如下形式：

L A (J)= 1-L A (J)+ 1-L A (J+B)* P

N W

L A (J+B)= 1-L A (J)-1-L A (J+B)* P

N W

总结上述运算规律，可采用如下运算方法进行DIT-FFT 运算。首先读入数据，根据数据长度确定运算级数M ，运算总点数M

N 2=，不足补0处理。然后对读入数据进行数据倒序操作。数据倒序后从第1级开始逐级进行，共进行M 级运算。在进行第L 级运算时，先算出该级不同旋转因子的个数1

2-=L B (也是该级中各个蝶形运算两输入数据的间距)，再从R=0开始按序计算，直到R=B-1结束。每个R 对应的旋转因子指数R P L

M -=2

，旋转因子指数相同的蝶从上往下依次逐个运算，

各个蝶的第一节点标号k 都是从R 开始，以L

2为步长，到R N L

+-2

(可简取极值

N-2)结束。考虑到蝶形运算有两个输出，且都要用到本级的两个输入数据，故第一个输出计算完毕后，输出数据不能立即存入输入地址，要等到第二个输出计算调用输入数据完毕后才能覆盖。这样数据倒序后的运算可用三重循环程序实现。整个蝶形运算流程图如图2-3所示：

图2-3 整个蝶形运算流程图

3 Matlab程序实现

3.1程序运行结果

用图片cj.jpg作为例子，运行调试程序。

图片cj.jpg如图3-1所示：

图3-1 原图像

程序运行开始的时候可以先建立一个M_file。生成M_file方法为：File->New->M-File，或者单击Matlab主界面的空白文档。如图3-2所示：

图3-2 生成M文件

将自行编写的程序写到文件中，然后保存，运行，程序执行之后会自动查询所要处理的图片。程序开始时要输入的窗口如图3-2所示：

图3-2 开始输入程序的界面

将已经编号的程序写到图3-2所示的窗口中保存之后就可以运行了。用自建的FFT和IFFT运行的结果如图3-3所示：

图3-3 程序的运行结果

3.2对比结果分析

如图3-3可以看出自编的FFT算法和系统自带的FFT算法运行结果有所不同，具体体现在三个方面。

一是：自建的FFT处理后的图片的大小发生了变化，宽度和内置的FFT处理后的结果相同；高度变小了。

二是：图像的亮度发生了变化，自建的FFT处理后的图片的亮度比系统自带的FFT处理后的图像的亮度高。

三是：程序的运行时间不同。自建的FFT算法运行的时间比系统自带的FFT 算法运行的时间要长。

但是不论何种算法实现FFT，它们的进行IFFT后得出的结果相同，和原图

片一样。

4 系统人机对话界面

4.1 GUI简介

图形用户界面（GUI），是一种提供人机交互的工具和方法。GUI是包含图形对象，如窗口、图标、菜单和文本等图文并茂的用户界面。

4.2 界面设计

用MATLAB图形用户界面开发环境设计GUI点的一般步骤是：

1.进行界面设计。

2.设计控件属性。

3.进行M语言编程。

以本设计要求为例介绍。

第一步，该选择本图形用户界面需要的控件：

七个推按钮(Push button)，其中六个用来运行所要显示的图片，第七个用来退出程序。

六个轴对象(axes)用来显示原图像，两种不同FFT处理后的图片以及进行IFFT后的图像。

完成人机设计界面如图4-1所示：

图4-1 人机界面设计图

第二步，设置控件属性：

双击组件可以设置文本框，推按钮的属性，如显示大小，名称和默认值等。

第三步，编写回调函数。

组件事件的发生是通过回调函数进行工作的。控件设置完成后保存，然后运行GUI（操作为ctrl+T），就会进入editor窗口，加入各个控件功能的函数代码。

完成后保存即可。

第四步，运行GUI。

运行editor窗口的程序后，会弹出已经激活的人机对话界面。

系统人机对话界面如图4-2所示。

图4-2 系统的人机对话界面

4.3 运行调试

运行GUI，已经弹出图4-2的系统人机对话界面，根据提示分别点击T1，T2，T3，T4，T5，T6可以分别运行出原图像，自编FFT处理后的图像，自编IFFT 后的图像，和原图像，内置FFT处理后的图像，内置IFFT处理后的图像。运行结果如图4-3所示。如果运行结束，可以通过点击“EXIT”推按钮退出该人机对话界面，返回Matlab。

图4-3 运行结果显示

5 Matlab软件简介

MATLAB产生的历史背景20世纪70年代中期，Cleve Moler博士和其同事在美国国家科学基金的资助下开发了调用EISPACK和LINPACK的FORTRAN 子程序库。EISPACK是特征值求解的FORTRAN程序库，LINPACK是解线性方程的程序库。在当时，这两个程序库代表矩阵运算的最高水平。在以后的数年里，MATLAB在多所大学里作为教学辅助软件使用，并作为面向大众的免费软件广为流传。1983年春天，Cleve Moler到Stanford大学讲学，MATLAB深深地吸引了工程师John Little。John Little敏锐地觉察到MATLAB在工程领域的广阔前景。同年，他和Cleve Moler、Sieve Bangert一起，用C语言开发了第二代专业版。这一代的MATLAB语言同时具备了数值计算和数据图示化的功能。1984年，Cleve Moler和John Lithe成立了MathWorks公司，正式把MATLAB推向市场，并继续进行MATLAB的研究和开发。在当今30多个数学类科技应用软件中，就软件数学处理的原始内核而言，可分为两大类。一类是数值计算型软件，如MATLAB、Xmath、Gauss等，这类软件长于数值计算，对处理大批数据效率高；另一类是数学分析型软件，如Mathematica、Maple等，这类软件以符号计算见长，能给出解析解和任意精度解，其缺点是处理大量数据时效率较低。

MathWorks公司顺应多功能需求之潮流，在其卓越数值计算和图示能力的基础上，又率先在专业水平上开拓了其符号计算、文字处理、可视化建模和实时控制能力，开发了适合多学科、多部门要求的新一代科技应用软件MATLAB。经过多年的国际竞争，MATLAB 已经占据了数值型软件市场的主导地位。在MATLAB进入市场前，国际上的许多应用软件包都是直接以FORTRAN和C语言等编程语言开发的。这种软件的缺点是使用面窄、接口简陋、程序结构不开放以及没有标准的基库，很难适应各学科的最新发展，因而很难推广。MATLAB 的出现，为各国科学家开发学科软件提供了新的基础。在MATLAB问世不久的20世纪80年代中期，原先控制领域里的一些软件包纷纷被淘汰或在MATLAB 上重建。时至今日，经过Math Works公司的不断完善，MATLAB已经发展成为适合多学科、多种工作平台的功能强劲的大型软件。在国外，MATLAB已经经受了多年考验。在欧美等高校，MATLAB已经成为线性代数、自动控制理论、数理统计、数字信号处理、时间序列分析、动态系统仿真等高级课程的基本教学工具；成为攻读学位的大学生、硕士生、博士生必须掌握的基本技能。在设计研究单位和工业部门，MATLAB被广泛用于科学研究和解决各种具体问题。

6 心得体会

本次实习的主要内容是通过用Matlab实现FFT的设计，可以实现对一张图片进行FFT运算和IFFT运算并且输出结果。把自己编写的FFT算法与Matlab 自带FFT算法进行比较。

在之前数字信号处理的学习以及完成实验的过程中，已经使用过Matlab，对其有了一些基础的了解和认识，通过这次的课程设计使我进一步了解了图像信号的分析方法，以及其中产生信号和绘制信号的基本命令和一些基础编程语言。让我感受到只有在了解课本知识的前提下，才能更好的应用这个工具，并且熟练的应用Matlab也可以很好的加深我对课程的理解，方便我的思维。这次课程设计使我了解了Matlab的使用方法，提高了自己的分析和动手实践能力。同时我相信，进一步加强对MATLAB的学习与研究对我今后的学习将会起到很大的帮助。

这次的课程设计是对本学期所学知识的一次重要巩固，使得在课堂上掌握的知识得到了真正的运用。在学习的过程中和同学讨论，更明白了理论知识与实践的联系。书到用时方恨少，有些知识学会是一回事，掌握是一回事，但应用起来，确实不是那么简单的，需要很多知识的融会贯通。

程序运行调试初期，曾经多次出现错误、不能产生图形等问题，但在我翻阅资料认真改正及老师同学的帮助下基本功能还是完成了，经过1个星期的上机实习，程序已得到一些完善，能完成基本的要求的功能。最后经过努力，又深入学习了图形用户界面（GUI），完成了选做要求的人机对话界面。

学习就是一个了解，疑惑，进而解惑的过程，这次实习就是提供了这样一个发现自己知识漏洞，与同学老师探讨进行解惑的的机会。

通过这次课程设计实习，我更深刻的了解了Matlab的运用，重新复习了FFT 中的重要的序列倒序和蝶形变换的程序，对课本上的知识有了更深的理解，使我对数字信号处理有了系统的认知。

在这里特别感谢魏老师和李老师，他们给了我们很大的发挥空间，让我们真正自己动手真正掌握了知识，感谢他们细心指导。也非常感谢我的同学，他们解开了我在实习中出现的诸多知识死角，谢谢大家！

参考文献

[1]范寿康DSP 技术与DSP芯片.北京：电子工业出版社

[2]程佩青.数字信号处理教程.北京：清华大学出版社出版，2001

[3]高西全丁玉美等.数字信号处理.北京：电子工业出版社，2009

[4] 李勇徐震.MA TLAB辅助现代工程数字信号处理.西安电子科技大学出版社

[5] 陈杰.Matlab宝典.电子工业出版社

[6] 苏金明张莲花刘波.MATLAB工具箱应用，电子工业出版社

附录Ⅰ

function image_process_FFT()

[filename, pathname]=uigetfile({'*.jpg;*.tif;*.bmp;*.gif' },'File Selector'); image=imread(strcat(pathname,filename)); %获取一个图像信号

if ndims(image)==3

image=rgb2gray(image);

end

scrsz=get(0,'ScreenSize');

figure('position',[0 0 scrsz(3)-1 scrsz(4)]);

set(gcf,'Name','快速傅里叶变换');

subplot(2,3,1);

imshow(image);

title('原始图像');

subplot(2,3,4);

imshow(image);

title('原始图像');

[r,c]=size(image);

%补0至最接近的2的整数次幂

array=image;

t=log2(r);

t1=floor(t);

t2=ceil(t);

if t1~=t2

array(2^t2,c)=0;

end

[r1,c1]=size(array);

t=log2(c1);

t3=floor(t);

t4=ceil(t);

if t3~=t4

array(r1,2^t4)=0;

end

[r1,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);

end

n=c1/2;

data_row=zeros(1,n,'double'); %按行方向计算时用到for m=1:n

data_row(m)=exp(-1i*2*pi*(m-1)/r1);

end

array=transform_fft2(array);

Ft=fftshift(array);

S1=log(1+abs(Ft));

subplot(2,3,2);

imshow(S1,[]);

title('自建FFT2函数结果');

array=transform_ifft2(array);

array=abs(array);

array=array(1:r,1:c);

subplot(2,3,3);

imshow(array,[]);

title('自建IFFT2结果');

F=fft2(image);

FC=fftshift(F);

S=log(1+abs(FC));

subplot(2,3,5)

imshow(S,[]);

title('内置FFT2结果');

array=ifft2(F);

array=round(abs(array));

subplot(2,3,6);

imshow(array,[]);

title('内置IFFT2结果');

return

function array=transform_fft2(array)

array=double(array);

[r1 c1]=size(array);

for j=1:r1

array(j,:)=transform_fft(array(j,:));

end

for j=1:c1

array(:,j)=transform_fft((array(:,j)));

end

N=length(array);

n=N/2;

w=zeros(1,n,'double');

for m=1:n

w(m)=exp(-1i*2*pi*(m-1)/N);

end

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);

array(k*t1*2+j+1)=data1+data2;

array(k*t1*2+j+t1+1)=(data1-data2)*w(k*t1+1);

end

end

end

end

if mod(p,2)==1

return

else

array1=array;

return

end

function array=transform_ifft2(array) array=conj(array);

[r1,c1]=size(array);

for i=1:r1

array(i,:)=transform_fft(array(i,:)); end

for i=1:c1

array(:,i)=transform_fft(array(:,i)); end

array=array/(r1*c1);

附录Ⅱ

function varargout = image_fft(varargin)

gui_Singleton = 1;

gui_State = struct('gui_Name', mfilename, ...

'gui_Singleton', gui_Singleton, ...

'gui_OpeningFcn', @image_fft_OpeningFcn, ...

'gui_OutputFcn', @image_fft_OutputFcn, ...

'gui_LayoutFcn', [] , ...

'gui_Callback', []);

if nargin && ischar(varargin{1})

gui_State.gui_Callback = str2func(varargin{1});

end

if nargout

[varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:}); else

gui_mainfcn(gui_State, varargin{:});

end

function image_fft_OpeningFcn(hObject, eventdata, handles, varargin) handles.output = hObject;

guidata(hObject, handles);

function varargout = image_fft_OutputFcn(hObject, eventdata, handles) varargout{1} = handles.output;

function pushbutton1_Callback(hObject, eventdata, handles)

image=imread('cj.jpg');

if ndims(image)==3

image=rgb2gray(image);

end

axes(handles.one_axes)

imshow(image);

title('原始图像');

function pushbutton2_Callback(hObject, eventdata, handles)

image=imread('cj.jpg');

if ndims(image)==3

image=rgb2gray(image);

end

axes(handles.two_axes)

[r,c]=size(image);

array=image;

t=log2(r);

t1=floor(t);

t2=ceil(t);

if t1~=t2

array(2^t2,c)=0;

end

[r1,c1]=size(array);

t=log2(c1);

t3=floor(t);

t4=ceil(t);

if t3~=t4

array(r1,2^t4)=0;

end

[r1,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);

end

n=c1/2;

data_row=zeros(1,n,'double'); %按行方向计算时用到for m=1:n

data_row(m)=exp(-1i*2*pi*(m-1)/r1);

end

%正变换

array=transform_fft2(array);

Ft=fftshift(array);

S1=log(1+abs(Ft));

imshow(S1,[]);

title('自建FFT2函数结果');