基于MATLAB的FFT算法的设计

基于MATLAB的FFT算法的设计
基于MATLAB的FFT算法的设计

目录

1前言 (1)

2FFT算法的基本原理 (2)

2.1系统总体流程图 (2)

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

2.2.1对图片的选择 (2)

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

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

3 软件简介 (5)

3.1 Matlab简介 (5)

3.1.1 Matlab软件概况 (5)

3.1.2 Matlab的特点 (5)

3.2 GUI简介 (6)

3.2.1界面设计 (6)

3.3对比结果与分析 (8)

4心得体会 (10)

参考文献 (11)

附录ⅠMatlab源程序 (12)

附录ⅡGUI源程序 (16)

1前言

随着信息时代,数字时代的到来,数字信号处理已经成为一门极其重要的学科和技术领域。以DSP为核心芯片的处理系统日益变成了数字信号处理系统的主流。它广泛用于电子信息、通信、图像处理、语音处理、生物医学、自动控制、地质探测等领域,受到工程设计和使用人员的青睐。

MATLAB,它是美国Math Works公司推出的一种面向工程和科学计算的交互式计算软件。它以矩阵运算为基础,把计算、可视化、程序设计融合在一个简单易用的交互式工作环境中,是一款数据分析和处理功能都非常强大的工程适用软件。通过本次课设我们学会了分析和处理音频信号,首先要对图片信息进行采集,MATLAB的数据采集工具箱提供了一整套命令和函数,通过调用这些函数和命令,可直接控制图像进行数据采集。Window自带的程序也可驱动采集图片信息,并能保存该文件,供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函数实现对离散信号的傅里叶变换进行比较,如果得到的频谱相同,那么我们编写的程序就是正确的。用GUI界面完成人机交互方便使用的。本课程设计主要是对数字信号的分析。

2FFT 算法的基本原理

2.1系统总体流程图

本设计要求选择一个图片,并对该图像进行FFT 变换。在Matlab 环境下编写基2 DIT-FFT 算法;利用自己编写的算法对选取的图片进行分析,并与Matlab 数字信号处理工具箱中的fft 函数进行对比研究,验证自编算法的正确性。所以得到系统总体流程图如下图所示。

图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=rgb2gray(image); end

这些代码实现了对图片进行灰度变换,可以将原来是彩色的图片变换成黑白色,目的是为了取点时方便简单,而且速度更快,节省时间 。

图片选择

编写FFT 程序,得到所需的图片

对图片的路径选择

对图片进行灰度变换

实现选择图片的倒序

一级中不同蝶形运算

一级中相同蝶形运算

与M A T L A B 自带的 比较

2.2.2 FFT 算法的基本原理

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

所以,要求N 点的DFT,需要N 2次的复数乘法运算,N*(N-1)次复数乘法运算算。随着N 的增加,运算量将急剧增加,而在实际问题中,N 往往是较大的,因此无论是用通用计算机还是用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/4点序列,最后,将原序列划分为多个2点序列,将计算量大大降低。

按时间下标的奇偶将N 点x(n)分别抽取组成两个N/2点序列,分别记为x1(n)和x2(n),将x(n)的DFT 转化为x1(n)和x2(n)的DFT 的计算。

()

()()()1

1

,0,1,2,...,1

1,0,1,2,...,1

N n k N

n N n k

N k X

k x n W

k N x n X

k W n N N

-=--==

=-=

=-∑∑

)

()(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 //,

==1

2

,

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

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

利用旋转因子的可约性,即:

用蝶形运算图可表示为:

以此类推,还可以把x1(n)和x2(n)按n 值得奇偶分为两个序列,这样就达到了降N 得目的,从而减少了运算量。FFT 对DFT 的数学运算量改进:

直接采用DFT 进行计算,运算量为N 2次复数乘法和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 减少运算量的特点表现的越明显。

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

实现FFT 运算的核心是蝶形运算,找出蝶形运算的规律是编程的基础。对

M

N 2

=点的FFT 共有M 级运算,用L 表示从左到右的运算级数(L=1,2,…,M )。

第L 级共有1

2-=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 点。

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

N

2

=,不足补0处理。然后对读入

数据进行数据倒序操作。运算时,先算出该级不同旋转因子的个数1

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

M -=2

2j

2j

222

2

e

e

rk

N rk

rk

rk N

N N

W W ππ--===()()()1

1

2

2

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 蝶形运算图

3 软件简介

3.1 Matlab简介

70年代后期,身为美国New Mexico大学计算机系系主任的Cleve Moler发现学生用FORTRAN编写接口程序很费时间,于是他开始自己动手,利用业余时间为学生编写EISPACK和LINPACK的接口程序。Cleve Moler给这个接口程序取名为Matlab。1984年,为了推广Matlab在数值计算中的应用,Cleve Moler、Johon Little等正式成立了Math works公司,从而把Matlab推向市场,并开始了对Matlab 工具相等的开发设计。

3.1.1 Matlab软件概况

Matlab是Matrix Laboratory的缩写,意为矩阵实验室。它具有强大的矩阵处理功能和绘图功能,进还能进行文字处理,绘图,建模仿真等功能。随着版本的不断升级,它在数值计算及符号计算功能上得到了进一步完善。Matlab已经发展成为多学科、多种工作平台的功能强大的大型软件。在欧美等高校,Matlab 已经成为线性代数、自动控制理论、概率论及数理统计、数字信号处理、时间序列分析、动态系统仿真等高级课程的基本教学工具。

Matlab软件打开后,其展现的窗口图如下:

图3-1运行Matlab窗口图

3.1.2 Matlab的特点

Matlab有以下一些特点:

Matlab的帮助功能很强大,自带有详细的帮助手册,基于HTML的完整的帮助功能,也可以用help命令来得到帮助信息。

程序语法设计自由度大,方便我们编程。

还有一个原因使Matlab受人们欢迎的,那就是Matlab源程序具有很大的开

放性。除了内部函数以外,所有Matlab的核心文件和工具箱文件都是可读可改的源文件,用户可通过对源文件的修改以及加入自己的文件构成新的工具箱。

Matlab有强大的的图形绘制功能。在Matlab里,数据可视化的操作非常简单易用。Matlab还有较强的编辑图形界面的能力。可以用来声成图解和可视化的二维、三维图。

Matlab还拥有功能强大的各种工具箱。其工具箱分为两类:功能性工具箱和学科性工具箱。功能性工具箱主要用来扩充其符号计算功能,图示建模仿真功能,文字处理功能以及与硬件实时交互功能。功能性工具箱用于多种学科。而学科性工具箱是专业性比较强的,如(control、signal proceessing 、commumnication) toolbox等。这些工具箱都是由该领域内学术水平很高的专家编写的,所以用户无需编写自己学科范围内的基础程序,而直接进行高,精,尖的研究,能极大地促进我们的学习研究工作。

虽然Matlab有很多优点,但它也有一些缺点,比如:由于Matlab的程序不用编译等预处理,也不生成可执行文件,程序为解释执行,所以速度较慢。3.2 GUI简介

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

一个设计非常优秀的GUI能够非常直观的使用户知道如何操作Matlab界面,了解设计者开发的意图。

3.2.1界面设计

打开GUI后未设置任何属性的界面如下图所示:

图3-2GUI界面

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

1.进行界面设计。

2.设计控件属性。

3.进行M语言编程。

以本设计要求为例介绍。

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

八个推按钮(Push button),用来运行和退出。

六个轴对象(axes)用来显示信号时域波形,频谱图,和两种不同FFT实现的频谱图,即运用Matlab自带的。

第二步,设置各个对象的属性。

如果要设置各个对象在表面上的寓意,只要改变它的string即可。Tag是它的本质属性,是唯一的。控制属性完成后,图形用户界面如图所示。

图3-3控件属性设置完成后的图形用户界面

第三步,编写回调函数。

组件事件的发生是通过回调函数进行工作的。回调函数是用户指定组件需要完成一定动作的函数,如本次设计中单击按钮后的功能就是由编写回调函数实现的。一个组件中可以包括多个回调函数,但是几乎每一个组件都有的属性是Callback属性。

GUIDE可创建一个M回调函数文件的构架,文件将自动处理并将所有文件的句柄传递到handles结构数组中。Handles是GUI中一个特殊的结构数组,它是GUI中所有组件共用的一个结构。在handles结构数组中还可以添加用户自定义的,在不同的回调函数之间共享的数据。

回调函数的声明格式为:

Function object_Callback(hObject,eventData,handles)

声明格式各部分为:

object为发生事件控制的Tag属性字符串。

hObject为发生事件的控制句柄。

eventData为保留字段。

本设计中各个按钮回调函数完成3个任务:

1得到用户输入图片。

例如当点击原图按钮时,就会得到图片路径选择路径,得到图片。

图3-4得到图片

2计算数据。

3建立图形。

3.3对比结果与分析

得到原图的两个图形对比如下:

图3-5原图对比

由以上两图对比可以看出,自建的原图和自带的原图得到的结果并不是完全一样的,在色彩上有差别,明暗程度不一样,自建的比自带的颜色更鲜明,更深。

图3-6 FFT对比图

由以上两图对比可以看出,自建的fft和系统自带的ifft也是有细微的差别的,自建的得到的图片比系统自带的图像颜色更甚而且由中间汇聚的光束可以看出,自建的函数得到的结果光束比较大,发散的范围更广,而系统自带的函数得到的结果相对来说比较窄。

图3-7 IFFT对比图

由以上两图对比可以看出,自建的IFFT和系统自带的IFFT总体上是一样的,但是有细微的差别,自建的相对来说颜色较暗,更深而自带的得到的图像颜色更亮更鲜明。

4心得体会

通过本课程设计使我了解了数字图像的基本概念,掌握数字图像处理的基本内容,如图像点运算、几何变换、增强处理、图像复原等的基本原理和Matlab 实现方法。

通过本次课程设计,让我们掌握如何学习一门语言,如何进行资料查阅搜集,如何自己解决问题等方法,养成良好的学习习惯。扩展理论知识,培养学生的综合设计能力。

本次课设是通过用Matlab实现FFT的设计,对图像信号进行分析,并画出采样信号的时域与频域图。把自己编写的FFT算法与Matlab自带FFT算法进行比较。程序运行调试时,自己选择输入要采样的图像信号,采样像素以及要变换的范围,可以实现对不同信号的信号采样和进行不同点的FFT运算。

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

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

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

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

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

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

参考文献

[1]楼顺天,李博菡.基于MA TLAB的系统分析与设计—信号处理.西安电子科技大学出版社,1998

[2]奥本海姆.离散时间信号处理.科学出版社,2000

[3]宗孔德,胡广书.数字信号处理.清华大学出版社,1997

[4]程佩青.数字信号处理教程.北京:清华大学出版社出版,2001

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

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

[7]苏金明张莲花刘波.MA TLAB工具箱应用,电子工业出版社

附录Ⅰ 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);

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('自建FFT结果');

array=transform_ifft2(array);

array=abs(array);

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

subplot(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结果');

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

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

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

附录Ⅱ GUI源程序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', []);

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 untitled3_OpeningFcn(hObject, eventdata, handles, varargin) handles.output = hObject;

guidata(hObject, handles);

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

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

end

scrsz=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' },'File Selector'); image=imread(strcat(pathname,filename));

if ndims(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(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));

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

array=transform_ifft2(array);

array=abs(array);

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

axes(handles.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);

end

scrsz=get(0,'ScreenSize');

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(image)==3

image=rgb2gray(image);

end

scrsz=get(0,'ScreenSize');

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

[r,c]=size(image);

Matlab中的FFT使用说明

FFT是Fast Fourier Transform(快速傅里叶变换)的简称,FFT算法在MATLAB 中实现的函数是Y=fft(x,n)。刚接触频谱分析用到FFT时,几乎都会对MATLAB 的fft函数产生一些疑惑,下面以看一个例子(根据MATLA帮助修改)。 Fs = 2000; % 设置采样频率 T = 1/Fs; % 得到采用时间 L = 1000; % 设置信号点数,长度1 秒 t = (0:L-1)*T; % 计算离散时间, % 两个正弦波叠加 f1 = 80; A1 = 0.5; % 第一个正弦波100Hz,幅度0.5 f2 = 150; A2 = 1.0 ; % 第2个正弦波150Hz,幅度 1.0 A3 = 0.5; % 白噪声幅度; x = A1*sin(2*pi*f1*t) + A2*sin(2*pi*f2*t); % 产生离散时间信号; y = x + A3*randn(size(t)); % 叠加噪声; % 时域波形图 subplot(2,1,1) plot(Fs*t(1:50),x(1:50)) title('Sinusoids Signal') xlabel('time (milliseconds)') subplot(2,1,2) plot(Fs*t(1:50),y(1:50)) title('Signal Corrupted with Zero-Mean Random Noise') xlabel('time (milliseconds)') NFFT = 2A nextpow2(L); % 设置FFT点数,一般为2 的N次方,如1024,512 等Y = fft(y,NFFT)/L; % 计算频域信号, f = Fs/2*linspace(0,1,NFFT/2+1); %频率离散化,fft后对应的频率是-Fs/2到Fs/2,由NFFT个离散频点表示 % 这里只画出正频率; % Plot single-sided amplitude spectrum. figure; plot(f,2*abs(Y(1:NFFT/2+1))); % fft 后含幅度和相位,一般观察幅度谱,并把负频率加上去, title('Single-Sided Amplitude Spectrum of y(t)') xlabel('Frequency (Hz)')

Simulink下的频谱分析方法及matlab的FFT编程

Simulink下的频谱分析方法 实现功能: 信号发生器一个信号输入,实时显示其频谱分析 调用模块: 信号源(Signal Processing Blockset -> Signal Processing Sources -> Sine Wave) Tip 1:不能用连续的信号源 频谱观察窗(Signal Processing Blockset -> Signal Processing Sources -> Spectrum Scope)Tip 2: 不能用普通的观察窗 Tip 3:必须构上设置中的Buffer input. Buffer size 越大越精细。 Tip 4: 剩下的tips读帮助。 连接关系: 如下图所示 原理框图实验结果:

输出示意图------------------------------ ------------------------------ 实现功能: 从Workspace读取一组数,进行频谱分析 调用模块: From Workspace Tip 1: 采样时间不能用0,即必须使用离散模式 Tip 2: 从其他模型中Scope保存出来的“Structure with time”的数据可以直接用频谱观察窗(同上一功能) ------------------------------ ------------------------------ 实现功能: 从dSPACE读取一组数,进行频谱分析 实现方法:

1. 从dSPACE读数保存成文件,数据导入Workspace(过程略) 2. 采用从其他模型的Scope保存数据为“Structure with time”的方式构建一个结构变量ScopeData1 3. 使用以下代码将dSPACE数据dscapture拷贝到结构变量ScopeData1中 %% =[0::]; %纯粹为占位,19157为dSPACE保存数据长度 for i=1:19157 end %% 4. 采用下图中的模型进行频谱分析 实验结果: 通过以上方法对单轴压电加速度传感器进行灵敏度分析,下图分别为采用dSPACE和直接利用示波器分析的结果对比。

按时间抽取的基2FFT算法分析与MATLAB实现

按时间抽取的基2FFT 算法分析及MATLAB 实现 一、DIT-FFT 算法的基本原理 基2FFT 算法的基本思想是把原始的N 点序列依次分解成一系列短序列,充分利用旋转因子的周期性和对称性,分别求出这些短序列对应的DFT ,再进行适当的组合,得到原N 点序列的DFT ,最终达到减少运算次数,提高运算速度的目的。 按时间抽取的基2FFT 算法,先是将N 点输入序列x(n)在时域按奇偶次序分解成2个N/2点序列x1(n)和x2(n),再分别进行DFT 运算,求出与之对应的X1(k)和X2(k),然后利用图1所示的运算流程进行蝶形运算,得到原N 点序列的DFT 。只要N 是2的整数次幂,这种分解就可一直进行下去,直到其DFT 就是本身的1点时域序列。 图1 DIT-FFT 蝶形运算流图 二、DIT-FFT 算法的运算规律及编程思想 1.原位计算 对N=M 2点的FFT 共进行M 级运算,每级由N/2个蝶形运算组成。在同一级中,每个蝶的输入数据只对本蝶有用,且输出节点与输入节点在同一水平线上,这就意味着每算完一个蝶后,所得数据可立即存入原输入数据所占用的数组元素(存储单元),经过M 级运算后,原来存放输入序列数据的N 个存储单元中可依次存放X(k)的N 个值,这种原位(址)计算的方法可节省大量内存。 2.旋转因子的变化规律 N 点DIT ―FFT 运算流图中,每个蝶形都要乘以旋转因子p W N ,p 称为旋转因子的指数。例如N =8 =3 2 时各级的旋转因子: 第一级:L=1, 有1个旋转因子:p W N =J /4W N =J 2L W J=0 第二级:L=2,有2个旋转因子:p W N =J /2W N =J 2L W J=0,1 第三级:L=3,有4个旋转因子:p W N =J W N =J 2L W J=0,1,2,3 对于N =M 2的一般情况,第L 级共有1 -L 2个不同的旋转因子: p W N =J 2L W J=0,1,2,… ,1 -L 2-1 L 2=M 2×M -L 2 = N ·M -L 2 故: 按照上面两式可以确定第L 级运算的旋转因子

Matlab编程实现FFT变换.

Matlab编程实现FFT变换及频谱分析的程序代码 内容 1.用Matlab产生正弦波,矩形波,以及白噪声信号,并显示各自时域波形图 2.进行FFT变换,显示各自频谱图,其中采样率,频率、数据长度自选 3.做出上述三种信号的均方根图谱,功率图谱,以及对数均方根图谱 4.用IFFT傅立叶反变换恢复信号,并显示恢复的正弦信号时域波形图 源程序 %*************************************************************** **********% % FFT实践及频谱分析% %*************************************************************** **********% %*************************************************************** **********% %***************1.正弦波****************% fs=100;%设定采样频率 N=128; n=0:N-1; t=n/fs; f0=10;%设定正弦信号频率 %生成正弦信号 x=sin(2*pi*f0*t); figure(1); subplot(231); plot(t,x);%作正弦信号的时域波形 xlabel('t'); ylabel('y'); title('正弦信号y=2*pi*10t时域波形'); grid; %进行FFT变换并做频谱图 y=fft(x,N);%进行fft变换 mag=abs(y);%求幅值 f=(0:length(y)-1)'*fs/length(y);%进行对应的频率转换 figure(1); subplot(232); plot(f,mag);%做频谱图 axis([0,100,0,80]); xlabel('频率(Hz)'); ylabel('幅值'); title('正弦信号y=2*pi*10t幅频谱图N=128'); grid; %求均方根谱

MATLAB中FFT结果的物理意义

FFT结果的物理意义 最近正在做一个音频处理方面的项目,以前没有学过fft,只是知道有这么个东西,最近这一用才发现原来欠缺这么多,最基本的,连fft的输入和输出各自代表什么都不知道了,终于在网上查到这样的一点资料,得好好保存了,也欢迎大家分享。 FFT是离散傅立叶变换的快速算法,可以将一个信号变换到频域。有些信号在时域上是很难看出什么特征的,但是如果变换到频域之后,就很容易看出特征了。这就是很多信号分析采用FFT变换的原因。另外,FFT可以将一个信号的频谱提取出来,这在频谱分析方面也是经常用的。虽然很多人都知道FFT是什么,可以用来做什么,怎么去做,但是却不知道FFT之后的结果是什意思、如何决定要使用多少点来做FFT。现在圈圈就根据实际经验来说说FFT结果的具体物理意义。一个模拟信号,经过ADC采样之后,就变成了数字信号。采样定理告诉我们,采样频率要大于信号频率的两倍,这些我就不在此罗嗦了。 采样得到的数字信号,就可以做FFT变换了。N个采样点,经过FFT之后,就可以得到N个点的FFT结果。为了方便进行FFT运算,通常N取2的整数次方。假设采样频率为Fs,信号频率F,采样点数为N。那么FFT之后结果就是一个为N点的复数。每一个点就对应着一个频率点。这个点的模值,就是该频率值下的幅度特性。具体跟原始信号的幅度有什么关系呢?假设原始信号的峰值为A,那么FFT的结果的每个点(除了第一个点直流分量之外)的模值就是A的N/2倍。而第一个点就是直流分量,它的模值就是直流分量的N倍。而每个点的相位呢,就是在该频率下的信号的相位。第一个点表示直流分量(即0Hz),而最后一个点N的再下一个点(实际上这个点是不存在的,这里是假设的第N+1个点,也可以看做是将第一个点分做两半分,另一半移到最后)则表示采样频率Fs,这中间被N-1个点平均分成N等份,每个点的频率依次增加。例如某点n所表示的频率为:Fn=(n-1)*Fs/N(ps:横坐标第n个点对应的频率值Fn的计算公式。整个横坐标代表了采样频率Fs,被分为N点。故其频率分辨率为Fs/N)。由上面的公式可以看出,Fn所能分辨到频率为为Fs/N,如果采样频率Fs为1024Hz,采样点数为1024点,则可以分辨到1Hz。1024Hz的采样率采样1024点,刚好是1秒,也就是说,采样1秒时间的信号并做FFT,则结果可以分析到1Hz,如果采样2秒时间的信号并做FFT,则结果可以分析到0.5Hz。如果要提高频率分辨力,则必须增加采样点数,也即采样时间。频率分辨率和采样时间是倒数关系。 假设FFT之后某点n用复数a+bi表示,那么这个复数的模就是An=根号a*a+b*b,相位就是Pn=atan2(b,a)。根据以上的结果,就可以计算出n点(n≠1,且n<=N/2)对应的信号的表达式为:An/(N/2)*cos(2*pi*Fn*t+Pn),即2*An/N*cos(2*pi*Fn*t+Pn)。对于n=1点的信号,是直流分量,幅度即为A1/N。由于FFT结果的对称性,通常我们只使用前半部分的结果,即小于采样频率一半的结果。 好了,说了半天,看着公式也晕,下面圈圈以一个实际的信号来做说明。假设我们有一个信号,它含有2V的直流分量,频率为50Hz、相位为-30度、幅度为3V的交流信号,以及一个频率为75Hz、相位为90度、幅度为1.5V的交流信号。用数学表达式就是如下:S=2+3*cos(2*pi*50*t-pi*30/180)+1.5*cos(2*pi*75*t+pi*90/180)式中cos参数为弧度,所以-30度和90度要分别换算成弧度。我们以256Hz的采样率对这个信号进行采样,总共采样256点。按照我们上面的分析,Fn=(n-1)*Fs/N,我们可以知道,每两个点之间的间距就是1Hz,第n个点的频率就是n-1。我们的信号有3个频率:0Hz、50Hz、75Hz,应该分别在第1个点、第51个点、第76个点上出现峰值,其它各点应该接近0。实际情况如何呢?我们来看看FFT的结果的模值如图所示。

MATLAB中FFT的使用方法

MATLAB中FFT的使用方法 一.调用方法 X=FFT(x); X=FFT(x,N); x=IFFT(X); x=IFFT(X,N) 用MATLAB进行谱分析时注意: (1)函数FFT返回值的数据结构具有对称性。 例: N=8; n=0:N-1; xn=[4 3 2 6 7 8 9 0]; Xk=fft(xn) →Xk = 39.0000 -10.7782 + 6.2929i 0 - 5.0000i 4.7782 - 7.7071i 5.0000 4.7782 + 7.7071i 0 + 5.0000i -10.7782 - 6.2929i Xk与xn的维数相同,共有8个元素。Xk的第一个数对应于直流分量,即频率值为0。 (2)做FFT分析时,幅值大小与FFT选择的点数有关,但不影响分析结果。在IFFT时已经做了处理。要得到真实的振幅值的大小,只要将得到的变换后结果乘以2除以N即可。 二.FFT应用举例 例1:x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t)。采样频率fs=100Hz,分别绘制N=128、1024点幅频图。 clf; fs=100;N=128; %采样频率和数据点数

n=0:N-1;t=n/fs; %时间序列 x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t); %信号 y=fft(x,N); %对信号进行快速Fourier变换 mag=abs(y); %求得Fourier变换后的振幅 f=n*fs/N; %频率序列 subplot(2,2,1),plot(f,mag); %绘出随频率变化的振幅 xlabel('频率/Hz'); ylabel('振幅');title('N=128');grid on; subplot(2,2,2),plot(f(1:N/2),mag(1:N/2)); %绘出Nyquist频率之前随频率变化的振幅xlabel('频率/Hz'); ylabel('振幅');title('N=128');grid on; %对信号采样数据为1024点的处理 fs=100;N=1024;n=0:N-1;t=n/fs; x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t); %信号 y=fft(x,N); %对信号进行快速Fourier变换 mag=abs(y); %求取Fourier变换的振幅 f=n*fs/N; subplot(2,2,3),plot(f,mag); %绘出随频率变化的振幅 xlabel('频率/Hz'); ylabel('振幅');title('N=1024');grid on; subplot(2,2,4) plot(f(1:N/2),mag(1:N/2)); %绘出Nyquist频率之前随频率变化的振幅 xlabel('频率/Hz'); ylabel('振幅');title('N=1024');grid on; 运行结果:

MATLAB关于FFT频谱分析的程序

MATLAB关于FFT频谱分析的程序 %***************1.正弦波****************% fs=100;%设定采样频率 N=128; n=0:N-1; t=n/fs; f0=10;%设定正弦信号频率 %生成正弦信号 x=sin(2*pi*f0*t); figure(1); subplot(231); plot(t,x);%作正弦信号的时域波形 xlabel('t'); ylabel('y'); title('正弦信号y=2*pi*10t时域波形'); grid; %进行FFT变换并做频谱图 y=fft(x,N);%进行fft变换 mag=abs(y);%求幅值 f=(0:length(y)-1)'*fs/length(y);%进行对应的频率转换 figure(1); subplot(232); plot(f,mag);%做频谱图 axis([0,100,0,80]); xlabel('频率(Hz)'); ylabel('幅值');

title('正弦信号y=2*pi*10t幅频谱图N=128'); grid; %求均方根谱 sq=abs(y); figure(1); subplot(233); plot(f,sq); xlabel('频率(Hz)'); ylabel('均方根谱'); title('正弦信号y=2*pi*10t均方根谱'); grid; %求功率谱 power=sq.^2; figure(1); subplot(234); plot(f,power); xlabel('频率(Hz)'); ylabel('功率谱'); title('正弦信号y=2*pi*10t功率谱'); grid; %求对数谱 ln=log(sq); figure(1); subplot(235); plot(f,ln);

利用MATLAB实现信号DFT的计算

07级电信(2)班 刘坤洋 24 实验一 利用MATLAB 实现信号DFT 的计算 一、实验目的: 1、熟悉利用MATLAB 计算信号DFT 的方法 2、掌握利用MATLAB 实现由DFT 计算线性卷积的方法 二、实验设备:电脑、matlab 软件 三、实验内容: 1、练习用matlab 中提供的内部函数用于计算DFT (1) fft (x ),fft (x ,N ),ifft (x ),ifft (x ,N )的含义及用法 (2) 在进行DFT 时选取合适的时域样本点数N 请举例,并编程实现 题目: 源程序: >> N=30; %数据的长度 >>L=512; %DFT 的点数 >>f1=100; f2=120; >>fs=600; %抽样频率 >>T=1/fs; %抽样间隔 >>ws=2*pi*fs; >>t=(0:N-1)*T; >>f=cos(4*pi*f1*t)+cos(4*pi*f2*t); >>F=fftshift(fft(f,L)); >>w=(-ws/2+(0:L-1)*ws/L)/(2*pi); >>hd=plot(w,abs(F)); >>ylabel('幅度谱') >> xlabel('频率/Hz') 的频谱 分析利用)π4cos()π4cos()(DFT 21t f t f t x +=Hz 600,Hz 120,Hz 10021===s f f f

>> title('my picture') 结果图: (3) 在对信号进行DFT 时选择hamming 窗增加频率分辨率 请举例,并编程实现 题目: 源程序:>> N=50; %数据的长度 >>L=512; %DFT 的点数 >>f1=100;f2=150; >>fs=600; %抽样频率 >>T=1/fs; %抽样间隔 >>ws=2*pi*fs; >>t=(0:N-1)*T; >>f=cos(4*pi*f1*t)+0.15*cos(4*pi*f2*t); 的频谱 分析利用)π4cos(15.0)π4cos()(DFT 21t f t f t x +=Hz 600,Hz 150,Hz 10021===s f f f

关于使用Matlab里Powergui的FFTTool分析的问题及解决办法

首先设置 POWERLIB—》powergui,将该模块拖入模型中即可 在需要进行频谱分析的地方连接一示波器 示波器参数设定: Parameters—》Data history—》Save data to workspace; Format—》Structure with time. 运行一次后,双击powergui—》FFT Analysis. 1. 问题1及解决办法 仿真完成后,采用Powergui分析FFT,有时会发生错误:"simulation time of the signals is not enough long for the given fundamental frequency". 很多论坛说是仿真时间短了,可能这也是原因,不过更有可能是这样: FFT的数据来自于示波器SCOPE,在SCOPE PARAMETERS/GENERAL选项卡/SAMPLING 中,有DECIMATION和SAMPLE TIME两项,DECIMATION的意思是 The Decimation parameter allows you to write data at every nth sample, where n is the decimation factor. The default decimation, 1, writes data at every time step. 所以,如果选择DECIMATION,记录数据的时刻为第N个采样点,采样点间的时间间隔为采样步长,而在MATLAB Simulink中,如果采用变步长仿真,采样周期就是变化的,这样就很难对采样的数据进行FFT分析,或许软件只认可采样周期一定的数据,所以会出现文首的错误。 如果选择sample time,那么采样周期固定(与仿真步长无关),这样就可以进行FFT 分析了。所以如果遇到文首的错误,可以尝试将示波器的SAMPLing改为sample time,并设定采样周期,Sampling time

实验二 FFT算法的MATLAB实现

班级:学号:姓名 实验二FFT算法的MATLAB实现 (一)实验目的: (1)掌握用matlab进行FFT在数字信号处理中的高效率应用。 (2)学习用FFT对连续信号和时域离散信号进行谱分析。 (二)实验内容及运行结果: 题1:若x(n)=cos(nπ/6)是一个N=12的有限序列,利用MATLAB计算它的DFT 并进行IDFT变换同时将原图与IDFT变换后的图形进行对比。当求解IFFT变换中,采样点数少于12时,会产生什么问题。 程序代码: N=12; n=0:11; Xn=cos(n*pi/6); k=0:11; nk=n'*k; WN=exp(-j*2*pi/N) WNnk=WN.^nk XK=Xn*WNnk; figure(1) stem(Xn) figure(2) stem(abs(XK)) 运行结果:

IFFT变换中,当采样点数少于12时图像如下图显示:

分析:由图像可以看出,当采样点数小于12时,x(n)的频谱不变,周期为6,而XK 的频谱图发生改变。 题2:对以下序列进行谱分析 132()()103()8470x n R n n n x n n n =+≤≤?? =-≤≤??? 其他n 选择FFT 的变换区间N 为8和16点两种情况进行频谱分析,分别打印其幅频特 性曲线并进行对比、分析和讨论。 ㈠ 程序代码: x=ones(1,3);nx=0:2; x1k8=fft(x,8); F=(0:length(x1k8)-1)'*2/length(x1k8); %进行对应的频率转换 stem(f,abs(x1k8));%8点FFT title('8点FFTx_1(n)'); xlabel('w/pi'); ylabel('幅度'); N=8时:

MATLAB中FFT使用详解

MATLAB中FFT使用详解 一.调用方法 X=FFT(x); X=FFT(x,N); x=IFFT(X); x=IFFT(X,N) 用MA TLAB进行谱分析时注意: (1)函数FFT返回值的数据结构具有对称性。 例: N=8; n=0:N-1; xn=[4 3 2 6 7 8 9 0]; Xk=fft(xn) → Xk = 39.0000 -10.7782 + 6.2929i 0 - 5.0000i 4.7782 - 7.7071i 5.0000 4.7782 + 7.7071i 0 + 5.0000i -10.7782 - 6.2929i Xk与xn的维数相同,共有8个元素。Xk的第一个数对应于直流分量,即频率值为0。 (2)做FFT分析时,幅值大小与FFT选择的点数有关,但不影响分析结果。在IFFT时已经做了处理。要得到真实的振幅值的大小,只要将得到的变换后结果乘以2除以N即可。 二.FFT应用举例 例1:x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t)。采样频率fs=100Hz,分别绘制N=128、1024点幅频图。 clf; fs=100;N=128; %采样频率和数据点数 n=0:N-1;t=n/fs; %时间序列 x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t); %信号 y=fft(x,N); %对信号进行快速Fourier变换 mag=abs(y); %求得Fourier变换后的振幅 f=n*fs/N; %频率序列 subplot(2,2,1),plot(f,mag); %绘出随频率变化的振幅

xlabel('频率/Hz'); ylabel('振幅');title('N=128');grid on; subplot(2,2,2),plot(f(1:N/2),mag(1:N/2)); %绘出Nyquist频率之前随频率变化的振幅 xlabel('频率/Hz'); ylabel('振幅');title('N=128');grid on; %对信号采样数据为1024点的处理 fs=100;N=1024;n=0:N-1;t=n/fs; x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t); %信号 y=fft(x,N); %对信号进行快速Fourier变换 mag=abs(y); %求取Fourier变换的振幅 f=n*fs/N; subplot(2,2,3),plot(f,mag); %绘出随频率变化的振幅 xlabel('频率/Hz'); ylabel('振幅');title('N=1024');grid on; subplot(2,2,4) plot(f(1:N/2),mag(1:N/2)); %绘出Nyquist频率之前随频率变化的振幅 xlabel('频率/Hz'); ylabel('振幅');title('N=1024');grid on; 运行结果: fs=100Hz,Nyquist频率为fs/2=50Hz。整个频谱图是以Nyquist频率为对称轴的。并且可以明显识别出信号中含有两种频率成分:15Hz和40Hz。由此可以知道FFT变换数据的对称性。因此用FFT对信号做谱分析,只需考察0~Nyquist频率范围内的福频特性。若没有给出采样频率和采样间隔,则分析通常对归一化频率0~1进行。另外,振幅的大小与所用采样点数有关,采用128点和1024点的相同频率的振幅是有不同的表现值,但在同一幅图中,40Hz

用MATLAB进行FFT频谱分析

用MATLAB 进行FFT 频谱分析 假设一信号: ()()292.7/2cos 1.0996.2/2sin 1.06.0+++=t t R ππ 画出其频谱图。 分析: 首先,连续周期信号截断对频谱的影响。 DFT 变换频谱泄漏的根本原因是信号的截断。即时域加窗,对应为频域卷积,因此,窗函数的主瓣宽度等就会影响到频谱。 实验表明,连续周期信号截断时持续时间与信号周期呈整数倍关系时,利用DFT 变换可以得到精确的模拟信号频谱。举一个简单的例子: ()ππ2.0100cos +=t Y 其周期为0.02。截断时不同的持续时间影响如图一.1:(对应程序shiyan1ex1.m ) 图 错误!文档中没有指定样式的文字。.1 140.0160.0180.02 截断时,时间间期为周期整数倍,频谱图 0.0250.03 20 40 60 80 100 截断时,时间间期不为周期整数倍,频谱图

其次,采样频率的确定。 根据Shannon 采样定理,采样带限信号采样频率为截止频率的两倍以上,给定信号的采样频率应>1/7.92,取16。 再次,DFT 算法包括时域采样和频域采样两步,频域采样长度M 和时域采样长度N 的关系要符合M ≧N 时,从频谱X(k)才可完全重建原信号。 实验中信号R 经采样后的离散信号不是周期信号,但是它又是一个无限长的信号,因此处理时时域窗函数尽量取得宽一些已接近实际信号。 实验结果如图一.2:其中,0点位置的冲激项为直流分量0.6造成(对应程序为shiyan1.m ) 图 错误!文档中没有指定样式的文字。.2 ?ARMA (Auto Recursive Moving Average )模型: 将平稳随机信号x(n)看作是零均值,方差为σu 2的白噪声u(n)经过线性非移变系统H(z)后的输出,模型的传递函数为 020406080100120140160180200 0.4 0.50.60.7 0.800.050.10.150.20.250.30.350.40.450.5 50100 150

MATLAB中FFT的使用方法

MATLAB中FFT的使用方法 调用方法 X=FFT(x); X=FFT(x,N);%N为FFT后的数据点数,如果实际信号的数据点数小于N的话,则需要在FFT变换时增加采样点数,或者通过采用频率细分法在原数据后面补充一定数量的0,从而满足N个数据点 X=IFFT(X); X=IFFT(X,N) 一、用MATLAB进行谱分析时注意: (1)函数FFT返回值的数据结构具有对称性。 例: N=8; n=0:N-1; xn=[4 3 2 6 7 8 9 0]; Xk=fft(xn) Xk = 39.0000 -10.7782 + 6.2929i 0 - 5.0000i 4.7782 - 7.7071i 5.0000 4.7782 + 7.7071i 0 + 5.0000i -10.7782 - 6.2929i Xk与xn的维数相同,共有8个元素。Xk的第一个数对应于直流分量,即频率值为0。 (2)做FFT分析时,幅值大小与FFT选择的点数有关,但不影响分析结果。

在IFFT时已经做了处理。要得到真实的振幅值的大小,只要将得到的变换后结果乘以2除以N即可。 二、FFT应用举例 例1:x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t)。采样频率fs=100Hz,分别绘制N=128、1024点幅频图。 clf; fs=100;N=128; %采样频率和数据点数 n=0:N-1;t=n/fs; %时间序列 x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t); %信号 y=fft(x,N); %对信号进行快速Fourier变换 mag=abs(y); %求得Fourier变换后的振幅 f=n*fs/N; %频率序列 subplot(2,2,1),plot(f,mag); %绘出随频率变化的振幅 xlabel('频率/Hz'); ylabel('振幅');title('N=128');grid on; subplot(2,2,2),plot(f(1:N/2),mag(1:N/2)); %绘出Nyquist频率之前随频率变化的振幅 xlabel('频率/Hz'); ylabel('振幅');title('N=128');grid on; %对信号采样数据为1024点的处理 fs=100; N=1024;

基于matlab的FFT算法程序设计

数字通信课程设计报告书 课题名称 基于matlab 的FFT 算法程序设计 姓 名 学 号 院 系 物理与电信工程系 专 业 电子信息工程 指导教师 2010年 01 月15日 ※※※※※※※※※ ※ ※ ※※ ※※ ※※ ※※※※※ ※※ 2007级数字通信 课程设计

基于matlab的FFT算法程序设计 0712401-36 李晔 (湖南城市学院物理与电信工程系通信工程专业,益阳,413000) 一、设计目的 1.通过该设计,进一步了解MATLAB软件。 2.通过该设计,进一步熟悉MATLAB的语法规则和编辑方式。 3.通过该设计,掌握傅里叶变换的含义和方法。 二、设计的主要要求 掌握Fourier变换,解了关于MATLAB软件在数字信号处理方面的应用,熟悉MATLAB的语法规则和编程。用MATLAB实现快速Fourier变换。 三、整体设计方案 对信号x=sin(2*pi*f0*t)进行频谱分析,用MATLAB仿真。选取抽样频率为fs=100Hz,依照下列条件用MATLAB软件对信号xt进行傅里叶变换y=fft(xt,N)并绘制频谱图,观察所产生的六幅频谱图进行对比,并进行分析。 四、程序设计 fs=100;%设定采样频率 N=128; n=0:N-1; t=n/fs; f0=10;%设定正弦信号频率

%生成正弦信号 x=sin(2*pi*f0*t); figure(1); subplot(321); plot(t,x);%作正弦信号的时域波形 xlabel('t'); ylabel('y'); title('正弦信号y=2*pi*10t时域波形'); grid; %进行FFT变换并做频谱图 y=fft(x,N);%进行fft变换 mag=abs(y);%求幅值 m=length(y); f=(0:m/2-1)'*fs/m;%进行对应的频率转换 figure(1); subplot(322); plot(f,mag(1:m/2));%做频谱图 axis([0,100,0,80]); xlabel('频率(Hz)'); ylabel('幅值'); title('正弦信号y=2*pi*10t幅频谱图N=128'); grid; %求均方根谱 sq=abs(y); figure(1); subplot(323); plot(f,sq(1:m/2)); xlabel('频率(Hz)'); ylabel('均方根谱');

利用MATLAB编写FFT快速傅里叶变换

一、实验目的 1.利用MATLAB 编写FFT 快速傅里叶变换。 2.比较编写的myfft 程序运算结果与MATLAB 中的FFT 的有无误差。 二、实验条件 PC 机,MATLAB7.0 三、实验原理 1. FFT (快速傅里叶变换)原理: 将一个N 点的计算分解为两个N/2点的计算,每个N/2点的计算再进一步分解为N/4点的计算,以此类推。根据DFT 的定义式,将信号x[n]根据采样号n 分解为偶采样点和奇采样点。设偶采样序列为y[n]=x[2n],奇采样序列为z[n]=x[2n+1]。 上式中的k N W -为旋转因子N k j e /2π-。下式则为y[n]与z[n]的表达式: 2. 蝶形变换的原理:

下图给出了蝶形变换的运算流图,可由两个N/2点的FFT (Y[k]和Z[k]得出N 点FFT X[k])。同理,每个N/2点的FFT 可以由两个N/4点的FFT 求得。按这种方法,该过程可延迟后推到2点的FFT 。 下图为N=8的分解过程。图中最右边的为8个时域采样点的8点FFTX[k],由偶编号采样点的4点FFT 和奇编号采样点的4点得到。这4点偶编号又由偶编号的偶采样点的2点FFT 和奇编号的偶采样点的2点FFT 产生。相同的4点奇编号也是如此。依次往左都可以用相同的方法算出,最后由偶编号的奇采样点和奇编号的偶采样点的2点FFT 算出。图中没2点FFT 成为蝶形,第一级需要每组一个蝶形的4组,第二级有每组两个蝶形的两组,最后一级需要一组4个蝶形。 四、实验内容 1.定义函数disbutterfly ,程序根据FFT 的定义:]2 [][][N n x n x n y + +=、n N W N n x n x n z -+ -=])2 [][(][,将序列x 分解为偶采样点y 和奇采样点z 。

MATLAB中FFT的使用方法

MATLAB FFT 的使用方法 2009-08-22 11:00 说明:以下资源来源于《数字信号处理的 MATLAB ;现》万永革主编 一.调用方法 X=FFT(X); X=FFTX, N); x=IFFT(X); x=IFFT(X,N) 用MATLAB!行谱分析时注意: (1) 函数FFT 返回值的数据结构具有对称性。 例: N=8; n=0:N-1; xn=[4 3 2 6 7 8 9 0]; Xk=fft(xn) Xk = -10.7782 + 6.2929i 7.7071i 5.0000 0 + 5.0000i -10.7782 - 6.2929i Xk 与xn 的维数相同,共有8个元素。Xk 的第一个数对应于直流分量,即频率值 为00 (2) 做FFT 分析时,幅值大小与FFT 选择的点数有关,但不影响分析结果。在 IFFT 时已经做了处理。要得到真实的振幅值的大小,只要将得到的变换后结果 乘以2除以N 即可。 二.FFT 应用举例 例 1: x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t) 。采样频率 fs=100Hz,分别绘 制N=128 1024点幅频图。 clf; fs=100;N=128; %采样频率和数据点数 n=0:N-1;t=n/fs; %时间序歹 U x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t); % 信号 39.0000 5.0000i 4.7782 + 7.7071i 4.7782

y=fft(x,N); %对信号进行快速Fourier变换 mag=abs(y); %求得Fourier变换后的振幅 f=n*fs/N; %? 率序列 subplot(2,2,1),plot(f,mag); %绘出随频率变化的振幅 xlabel('频率/Hz'); ylabel('振幅');title('N=128');grid on; subplot(2,2,2),plot(f(1:N⑵,mag(1:N⑵);%绘出Nyquist 频率之前随频率变化的振幅 xlabel('频率/Hz'); ylabel('振幅');title('N=128');grid on; 以对信号采样数据为1024点的处理 fs=100;N=1024;n=0:N-1;t=n/fs; x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t); % 信号 y=fft(x,N); %对信号进行快速Fourier变换 mag=abs(y); %求取Fourier变换的振幅 f=n*fs/N; subplot(2,2,3),plot(f,mag); % 绘出随频率变化的振幅 xlabel('频率/Hz'); ylabel('振幅');title('N=1024');grid on; subplot(2,2,4) plot(f(1:N/2),mag(1:N/2)); % 绘出Nyquist频率之前随频率变化的振幅xlabel('频率/Hz'); ylabel('振幅');title('N=1024');grid on; 运行结果:

用matlab进行fft谐波分析

FFT是离散傅立叶变换的快速算法,可以将一个信号变换到频域。有些信号在时域上是很难看出什么特征的,但是如果变换到频域之后,就很容易看出特征了。这就是很多信号分析采用FFT变换的原因。另外,FFT可以将一个信号的频谱提取出来,这在频谱分析方面也是经常用的。 虽然很多人都知道FFT是什么,可以用来做什么,怎么去做,但是却不知道FFT之后的结果是什意思、如何决定要使用多少点来做FFT。 现在就根据实际经验来说说FFT结果的具体物理意义。一个模拟信号,经过ADC采样之后,就变成了数字信号。采样定理告诉我们,采样频率要大于信号频率的两倍,这些我就不在此罗嗦了。 采样得到的数字信号,就可以做FFT变换了。N个采样点,经过FFT之后,就可以得到N个点的FFT结果。为了方便进行FFT运算,通常N取2的整数次方。 假设采样频率为Fs,信号频率F,采样点数为N。那么FFT之后结果就是一个为N点的复数。每一个点就对应着一个频率点。这个点的模值,就是该频率值下的幅度特性。具体跟原始信号的幅度有什么关系呢?假设原始信号的峰值为A,那么FFT的结果的每个点(除了第一个点直流分量之外)的模值就是A的N/2倍。而第一个点就是直流分量,它的模值就是直流分量的N倍。而每个点的相位呢,就是在该频率下的信号的相位。第一个点表示直流分量(即0Hz),而最后一个点N的再下一个点(实际上这个点是不存在的,这里是假设的第N+1个点,也可以看做是将第一个点分做两半分,另一半移到最后)则表示采样频率Fs,这中间被N-1个点平均分成N等份,每个点的频率依次增加。例如某点n所表示的频率为:Fn=(n-1)*Fs/N。由上面的公式可以看出,Fn所能分辨到频率为为Fs/N,如果采样频率Fs 为1024Hz,采样点数为1024点,则可以分辨到1Hz。1024Hz的采样率采样1024点,刚好是1秒,也就是说,采样1秒时间的信号并做FFT,则结果可以分析到1Hz,如果采样2秒时间的信号并做FFT,则结果可以分析到0.5Hz。如果要提高频率分辨力,则必须增加采样点数,也即采样时间。频率分辨率和采样时间是倒数关系。 假设FFT之后某点n用复数a+bi表示,那么这个复数的模就是An=根号a*a+b*b,相位就是Pn=atan2(b,a)。根据以上的结果,就可以计算出n点(n≠1,且n<=N/2)对应的信号的表达式为:An/(N/2)*cos(2*pi*Fn*t+Pn),即2*An/N*cos(2*pi*Fn*t+Pn)。对于n=1点的信号,是直流分量,幅度即为A1/N。 由于FFT结果的对称性,通常我们只使用前半部分的结果,即小于采样频率一半的结果。 好了,说了半天,看着公式也晕,下面以一个实际的信号来做说明。 假设我们有一个信号,它含有2V的直流分量,频率为50Hz、相位为-30度、幅度为3V 的交流信号,以及一个频率为75Hz、相位为90度、幅度为1.5V的交流信号。用数学表达式就是如下: S=2+3*cos(2*pi*50*t-pi*30/180)+1.5*cos(2*pi*75*t+pi*90/180) 式中cos参数为弧度,所以-30度和90度要分别换算成弧度。我们以256Hz的采样率对这个信号进行采样,总共采样256点。按照我们上面的分析,Fn=(n-1)*Fs/N,我们可以知道,每两个点之间的间距就是1Hz,第n个点的频率就是n-1。我们的信号有3个频率:0Hz、50Hz、75Hz,应该分别在第1个点、第51个点、第76个点上出现峰值,其它各点应该接近0。实际情况如何呢?我们来看看FFT的结果的模值如图所示。

相关文档
最新文档