杜哈梅尔matlab,毕业论文范文——FFT算法的研究与Matlab编程实现

学号：

毕业设计

题 目：

FFT算法的研究与Matlab编程实现

作 者

届 别

系 别

专 业

电子信息工程

指导老师

职 称

完成时间

1

摘 要

快速傅里叶变 (|Fas Fourier Tranformation，FFT)是将一个大点数N的DFT分解为若干小点的D F T的组合。将用运算工作量明显降低， 从而大大提高 离散傅里叶变换(D F T) 的计算速度。因各个科学技术领域广泛的使用了FFT 技术它大大推动了信号处理技术的进步，现已成为数字信号处理强有力的工具，本论文将比较全面的叙述各种快速傅里叶变换算法原理、特点，并完成了基于MATLAB的实现。

关键词：离散傅立叶变换；快速傅立叶变换；蝶形单元；MATLAB

目 录

第一章　绪 论4

1.1 FFT算法的意义4

1.2 研究目标、内容4

第二章 基本理论6

2.1 FFT算法基本概念6

2.1.1离散傅里叶变换(DFT)6

2.1.2 快速傅里叶变换(FFT)

2.2 FFT算法分类7

2.2.1基2、DIT-FFT(按时间抽取)

2.2.2基2、DIF-FFT(按频率抽取)

2.2．3 基4、DIF-FFT(按频率抽取)

2.2．4 分裂基FFT算法

2.2.5 N为组合数的FFT——混合基算法

2.2.6 Chirp-z变换

2.3 MATLAB的应用13

2.3.1MATLAB主要功能21

2.4基本概念

2.4.130

2.4.2

2.4.2

2.4.4

第三章 FFT的MATLAB设计与实现22

3.122

3.2.22

3.325

3.4

3.5

第四章 FFT的分析31

4.131

4.231

4.331

第五章 总结与展望33

参考文献34

致谢35

第一章绪论

1.1． 引言

1965年，库利(J.W.Cooley)和图基(J.W.Tukey)在《计算数学》杂志上发表了“机器计算傅立叶级数的一种算法”的文章，这是一篇关于计算DFT的一种快速有效的计算方法的文章。它的思路建立在对DFT运算内在规律的认识之上。这篇文章的发表使DFT的计算量大大减少，并导致了许多计算方法的发现。这些算法统称为快速傅立叶变换(Fast Fourier Transform)，简称FFT，1984年，法国的杜哈梅尔(P.Dohamel)和霍尔曼(H.Hollmann)提出的分裂基快速算法，使运算效率进一步提高。FFT即为快速傅氏变换，是离散傅氏变换的快速算法，它是根据离散傅氏变换的奇、偶、虚、实等特性，对离散傅立叶变换的算法进行改进获得的。它对傅氏变换的理论并没有新的发现，但是对于在计算机系统或者说数字系统中应用离散傅立叶变换，可以说是进了一大步。

随着科学的进步，FFT算法的重要意义已经远远超过傅里叶分析本身的应用。FFT算法之所以快速，其根本原因在于原始变化矩阵的多余行，此特性也适用于傅里叶变换外的其他一些正交变换，例如，快速沃尔什变换、数论变换等等。在FFT的影响下，人们对于广义的快速正交变换进行了深入研究，使各种快速变换在数字信号处理中占据了重要地位。因此说FFT对数字信号处理技术的发展起了重大推动作用。[2]

1.2 FFT算法的研究与发展

近十多年来数字信号处理技术同数字计算机、大规模集成电路等先进技术一样，有了突飞猛进的发展，日新月异，已经形成了一门具有强大生命力的技术科学。由于它本身具有一系列的优点，所以能有效地促进各工程技术领域的技术改造和学科发展，应用领域也更加广泛、深入，越来越受到人们的重视。

在数字信号处理中，离散傅里叶变换(Discrete Fourier Transform, DFT)是常用的变换方法，它在各种数字信号处理系统中扮演着重要的角色。快速傅里叶变换〔Fast Fourier Transfonn, FFT〕并不是与离散傅里叶变换不同的另一种变换，而是为了减少DFT 计算次数的一种快速有效的算法。

傅里叶变换已有一百多年的历史了，我们知道频域分析常常比时域分析更优越，不仅简单，且易于分析复杂信号。但用较精确的数字方法，即 DFT 进行谱分析，在FFT 出现以前是不切实际的。这是因为 DFT 计算量太大。直到 1965 年出现了FFT。

应该指出，当时电子数字计算机的条件也促成了这个算法的提出。1967 年至 1968年间 FFT 的数字硬件就制成了。至此 DFT 的运算大为简化，运算时间一般可降低 1-2个数量级。因而各个科学技术领域广泛地采用了 FFT 技术，它大大推动了近 30 年来信号处理技术的发展，成为数字信号处理应用领域强有力的工具，广泛应用于雷达、声纳、通信、地质劫探、图像处理、生物医学等领域中。

Sande 提出了按照频率抽取的 FFT 算法，Bergland 提出了采用高基数结构的算法,Winograd 博士提出的，可以称为 WFTA 算法，Rader 和 Brenner 提出的余割因子算法，王中德提出的对称分解法，Vetlerli 和 Nussbaumer 提出的 DFT DCT 算法，其中最具代表性的是法国的 Duhamel 和 Hollman 提出的分裂基算法。各种 DFT 的快速算法，都利用了 的周期性和对称性，通过将一个大点数N的DFT分解为若干小点数的DFT的组合,来减少运算量

第二章 基本理论

2.1 FFT算法的基本概念

离散傅立叶变换(DFT)是信号分析与处理中的一种重要的变换。DFT 的计算，需要作 N2次复数乘法运算和 N(N－1)次复数加法运算，当N较大时，计算量太大，为了快速计算 DFT，近半个世纪以来，人们对离散傅立叶变换的计算进行了大量的研究，提出了很多有效的快速计算 DFT 的方法。这些算法，称之为快速傅立叶变换(Fast Fourier Transform，FFT)。快速博里叶变换并不是与DF T 不同的另外一种变换，而是为减少 DFT 计算次数的一种快速有效的算法。其突出的优点在于能快速高效地和比较精确地完成 DFT 的计算。

2.1.1离散傅里叶变换(DFT)

随着科学技术的进步，人们越来越正视频率分析技术的发展与应用。以前要进行一次频率分析，其计算的工作量大的惊人。现在，电子计算机的飞速发展使得计算的速度加快。但是，电子计算机不可能对连续的信号进行处理，它只能处理有限的离散数据。为了便于利用电子计算机进行傅立叶级数与傅立叶积分交换的运算，需要对连续信号进行采样，从而得到一系列离散数据。这种对离散量的傅立叶变换，称为离散傅立叶变换，记作 DFT。

DFT 的定义：

设x(n) 是一个长度为 N 的有限长序列，定义x(n) 的 N 点离散傅立叶变换为

(2．2.1)

其中k=0，1,…，N－1

X (k) 的傅立叶反变换为

(2．2.2)

其中n= 0,1,…,N-1

在(2-2-1)式和(2-2-2)式中， x(n) 和 X(k)均可以是复数。因为在(2-2-1)式和(2-2-2)式的右边仅在WN 指数上差一个符号，并相差一个比例因子 1／N，所以有关(2-2-1)式计算步骤的讨论稍加修改可以直接用于(2-2-2)式。

DFT隐含周期性，在DFT变换对中，x(n) 和 X(k)均为有限长序列，设

=，由的周期性，使得x(n) 和 X(k)隐含周期性，且周期为 N。

2.1.2 快速傅里叶变换(FFT)

从前一节的讨论中知道，DFT 的计算，需要作 N2次复数乘法运算和 N(N－1)次复数加法运算，运算量大。为了快速计算 DFT，近半个世纪以来，人们对离散傅立叶变换的计算进行了大量的研究，提出了很多有效的快速计算 DFT 的方法。这些算法，称之为快速傅立叶变换(Fast Fourier Transform，FFT)。快速博里叶变换并不是与DF T 不同的另外一种变换，而是为减少 DFT 计算次数的一种快速有效的算法。其突出的优点在于能快速高效地和比较精确地完成 DFT 的计算。

2.2.FFT算法的分类

DFT 分解法基本上分为两类：一类是将时间序列x(n) (n 为时间标号)进行逐次分解，由此得到的 FFT 算法称为按时间抽取(Decimation-in-time)算法，另一类是 将 傅 立 叶 变 换 序 列 X(k) (k为 频率标号)进行分 解 ，叫做按频率抽取(Decimation-in-frequency)算法。对这两种算法，库利—图基和桑德－图基进行了理论的推导，故又称为库利—图基〔Cooley－Tukey〕算法和桑德—图基(Sande－Tukey)算法。对每一算法，按基本的蝶形运算的构成又可分为基 2、基 4、基 8 以及任意因子等的 FFT 算法。不同基的 FFT 算法所需的计算量略有差异。之所以说略有差异是指并无数量级的差别，甚至无成倍的差别。只是某种基的算法比另一种省几分之几而已。表 2.2.1 列出了 N=4096 时各种基的 FFT 算法所需运算量的比较。

算法

实数乘法次数

实数加法次数

基2

81 924

139 266

基4

57 348

126 978

基8

49 156

126 978

基16

48 132

125 422

(表2.2.1 算法运算量的比较)

表 2.2.1 给出了N是2的任意次乘方时所需的运算量。这里假定每一种算法都以最有效的方法来完成尽可能多的计算。

从表格可以看出，一般地说，基数越高，总计算量愈少。但是我们不能光从运算量的角度来简单地判断说基数越高的算法越好。判断一个算法不仅要从计算量而且还应从算法的复杂性方面来加以考察。一般来说，基数愈高，算法本身也愈复杂。所以，选择算法要依据具体情况进行考虑。如考虑用何种语言设计程序，在什么机器上运行程序，或者采用何种逻辑器件和系统结构来构成FFT硬件等等。

2.2.1 基2、DIT-FFT(按时间抽取)

令，则有

蝶形运算单元如下所示：

2.2.2　基2、DIF-FFT(按频率抽取)

则有

蝶形运算单元如下所示：

由前面的分析可知，DIT(按时间抽取)算法与DIF(按频率抽取)算法没有本质上的区别，只是复数加减法与旋转因子乘法的次序有区别，两种方法的运算量是一样的。

在基2算法中，每个蝶形运算单元都包括1次复数乘法、2次复数加法。N(N= )点序列的运算流图应有M级蝶形，每一级都由N/2个蝶形运算组成，所以N点序列的基2FFT算法，总的运算量为次复数乘法，次复数加法。直接DFT运算量为次复数乘法、次复数加法。可见，FFT算法大大减少了运算量，当N越大时，FFT算法的优越性越明显。

2.2．3 基4、DIF-FFT(按频率抽取)

令：

则有

其蝶形运算单元如下所示：

2.2.4 分裂基FFT算法

1984年提出的分裂基(split-radix)算法同时使用基2和基4算法，被认为是目前对于N=2L各类算法中最为理想的一种。

当n=pq，且p=N/4，q=4时，n可表示为

并有

再将上式中的k表示为

可得

对k0=0，1,2,3,并用k表示k1，用n表示n0，可以写出

(1-1)

(1.2)

令

(1.3)

则(1.2)式可写成如下更简明的形式

(1.4)

得以上分裂基第一次分解L形流图

例如，N=16，第一次抽选分解时，由式(4.4.3)得

(n)=x(n)+x(n+8), 0≤n≤7

(n)={［x(n)-x(n+8)］-j［x(n+4)-x(n+12)］},

0≤n≤3

(n)={［x(n)-x(n+8)］+j［x(n+4)-x(n+12)］}

0≤n≤3

把上式代入式(4.4.4)，可得

X(2k)=DFT［ (n)］, 0≤k≤7

X(4k+1)=DFT［(n)］, 0≤k≤3

X(4k+3)=DFT［(n)］, 0≤k≤3

图(1.2) 分裂基FFT算法L形排列示意图与结构示意图如下

(a)分裂基FFT算法L形排列示意图；

(b)分裂基FFT算法运算流图结构示意图

。

(图1.2)

【16点分裂基第一次分解L形流图(图中省去箭头)】

第二次分解：

先对图4.4.3中N/2点DFT进行分解。

令 (l)=X(2l)，则有

(2l)=DFT［y2(n)］, 0≤l≤3

(4l+1)=DFT［(n)］, 0≤l≤1

(4l+3)=DFT［(n)］, 0≤l≤1

其中

y2(n)=x2(n)+x2(n+4), 0≤n≤3

(n)={［x2(n)-x2(n+4)］-［x2(n+2)x(n+6)］}Wn8,n=0,1

(n)={［x2(n)-x2(n+4)］+j［x2(n+2)x2(n+6)］}W3n8,n=0,1

(N/2点DFT的分解L形流图)

(4点分裂基L形运算流图)

(16点分裂基FFT运算流图 )

2.2．4 分裂基FFT算法的运算量

设第j级有lj个L形，j=1,2,…,M-1,M=log2N，则有l1=N/4。由图4.4.2(b)可见，第j-1列中的L形包含了第j列中的一部分结点的计算，即空白部分，所占结点数刚好等于第j-1列中所有L形对应结点的一半，所以第j列L形个数就减少l j-1/2个，即

由于每个L形有两次复(数)乘运算，所以全部复乘次数为

2.2．5 N为组合数的FFT——混合基算法

以2为基数的FFT算法，即N=2M，这种情况实际上使用得最多。其优点：程序简单，效率高，使用方便。

实际应用时，有限长序列的长度N很大程度上由人为因素确定，因此多数场合可取 N=2M，从而直接使用以2为基数的FFT算法。

如N不能人为确定 , N的数值也不是以2为基数的整数次方,处理方法有两种:

(1)补零:将x(n)补零 ,使N=2M.

例如N=30,补上x(30)=x(31)=0两点,使N=32=25,这样可直接采用以2为基数M=5的FFT程序。

有限长度序列补零后并不影响其频谱 X(ejw) ，只是频谱的采样点数增加了,上例中由30点增加到32点，所以在许多场合这种处理是可接受的。

(2)如要求准确的N点DFT值,可采用任意数为基数的 FFT 算法 , 其 计算效率低于以2为基数FFT算法。

如 N 为复合数，可分解为两个整数p与q的乘积,像前面以2为基数时一样,FFT的基本思想是将DFT的运算尽量分小,因此,在N=pq情况下,也希望将N点的DFT分解为p个q点DFT或q个p点DFT,以减少计算量。步骤如下：

，分别为0, 1,…，Q-1，，分别为0, 1,…，P-1

N点DFT可以重新写成为

考虑到

令

再令

以P=3，Q=4, N=12为例

(1) 先将 x(n)通过 x(n1Q+n0)改写成 x(n1,n0)。因为 Q=4, n1=0,1,2, n0=0,1,2,3,故输入是按自然顺序的，即

x(0,0)=x(0) x(0,1)=x(1) x(0,2)=x(2) x(0,3)=x(3)

x(1,0)=x(4) x(1,1)=x(5) x(1,2)=x(6) x(1,3)=x(7)

x(2,0)=x(8) x(2,1)=x(9) x(2,2)=x(10) x(2,3)=x(11)

(2) 求Ｑ个Ｐ点的DFT

(３)X1(k0,n0)乘以得到X1′(k0,n0)。

(４)求P个Ｑ点的DFT，参变量是k0

(５)将X2(k0,k1)通过X(k0+k1P)恢复为X(k)

(N=12为组合数时的FFT)

(１)求Ｑ个P点DFT需要ＱP2次复数乘法和Ｑ·P·(P-1)次复数加法；

(２)乘Ｎ个W因子需要Ｎ次复数乘法；

(３)求P个Ｑ点DFT需要PQ2 次复数乘法和P·Ｑ(Ｑ-1)次复数加法。

总的复数乘法量: QP2+N+PQ2=N(P+Q+1)；

总的复数加法量: Q·P(P-1)+P·Q·(Q-1)=N(P+Q-2)

上述分解原则可推广至任意基数的更加复杂的情况。

例如, 如果N可分解为m个质数因子p1,p2,…，pm，即

N=p1p2p3…pm

则

第一步：可把N先分解为两个因子N=p1q1，其中q1=p2p3…pm ，并用上述讨论的方法将DFT分解为p1个q1点DFT；

第二步：将q1分解为q1=p2q2，q2=p3p4…pm，然后将每个q1 点DFT再分解为p2个q2点DFT；

依此类推，

通过m次分解，一直分到最少点数的DFT运算，从而获得最高的运算效率。其运算量近似为N(p1+p2+…+pm)次复数乘法和复数加法。

但计算效率的提高是要以编程的复杂性为代价的，一般较少应用。

p1=p2 = …=pm =2，为基2 FFT 算法。

当组合数 N=P1P2P3…Pm中所有的Pi均为4时，就是基四FFT算法

以N=43为例,第一级运算的一般形式为：

2.2．6 Chirp-z变换

采用FFT可以算出全部N点DFT值，即z变换X(z)在z平面单位圆上的等间隔取样值，但要求 N 为复合数。

问题的提出：

(1)不需要计算整个单位圆上z变换的取样，如对于窄带信号，只需要对信号所在的一段频带进行分析，这时,希望频谱的采样集中在这一频带内，以获得较高的分辨率，而频带以外的部分可不考虑。

(2)对其它围线上的z变换取样感兴趣，例如语音信号处理中，需要知道z变换的极点所在频率，如极点位置离单位圆较远，则其单位圆上的频谱就很平滑，如果采样不是沿单位圆而是沿一条接近这些极点的弧线进行，则在极点所在频率上的将出现明显的尖峰，由此可较准确地测定极点频率。

(3)要求能有效地计算当N是素数时序列的DFT

算法原理

螺旋线采样是一种适合于这种需要的变换，且可以采用FFT来快速计算，这种变换也称作Chirp-z变换。

已知x(n)，0≤n≤N-1

令zk=AW-k ，k=0，…，M-1 ，M ：采样点数，A、W：任意复数

其中：

A0表示起始取样点的半径长度，通常A0≤1

θ0表示起始取样点z0的相角

φ0表两相邻点之间的等分角

W0螺旋线的伸展率，W0〈1则线外伸，W0〉1则线内缩(反时针)，W0=1则表示半径为A0的一段圆弧，若A0=1，这段圆弧则是单位圆的一部分。 计算z变换在采样点的值

k=0,1,… ,M-1

显然，按照以上公式计算出全部M点采样值需要 NM 次复乘和(N-1)M次复加，当N及M较大时，计算量迅速增加，以上运算可转换为卷积形式，从而可采用FFT进行，这样可大大提高计算速度。

利用zk的表示式代入

nk可以用以下表示式来替换

则

定义

则有

上式说明，如对信号x(n)先进行一次加权处理，加权系数为 ，然后通过一个单位脉冲响应为h(n)的线性系统，最后，对该系统的前M点输出再作一次 的加权，就可得到全部M点螺旋线采样值。

X(zk)

x(n)

X

X

g(n)

系统的单位脉冲响应 ,与频率随时间成线性增加的线性调频信号相似，因此称为Chirp -z变换。

2.3 MATLAB在FFT算法中的应用

MATLAB是一套用于科学计算的可视化高性能语言与软件环境。它集数值分析、矩阵运算、信号处理和图形显示于一体，构成了一个界面友好的用户环境。它的信号处理工具箱包含了各种经典的和现代的数字信号处理技术，是一个非常优秀的算法研究与辅助设计的工具。在进行FFT算法时，通常采用MATLAB来进行辅助设计和

1MATLAB主要功能

高级语言可用于技术计算开发环境可对代码、文件和数据进行管理

交互式工具可以按迭代的方式探查、设计及求解问题数学函数可用于线性代数、统计、傅立叶分析、筛选、优化以及数值积分等二维和三维图形函数可用于可视化数据各种工具可用于构建自定义的图形用户界面各种函数可将基于 MATLAB 的算法与外部应用程序和语言

第三章 FFT具体应用及其Matlab仿真

源程序如下

t=0:0.001:1.023;

x=sin(2*pi*50*t)+sin(2*pi*25*t); %周期信号的有限化和离散化

%subplot(2,2,1)

figure(1)

plot(t,x)

xlabel('时间轴t')

ylabel('信号值f(t)')

title('正弦波','FontSize',10)

num=length(x); cc=num;

lnum=log(num)/log(2); %获取蝶形运算的级数

pow=1;

a=[0,1];

for m=1:(lnum-1)

pow=pow*2;

for n=1:pow

b(2*n-1)=a(n);

b(2*n)=a(n)+pow;

end

a=b;

end

for n=1:1024

a(n)=a(n)+1;

end

for m=1:num

y(m)=x(a(m));

end

x=y %码位倒置算法

w0=exp(-i*2*pi/num);

pow=1;

for m=1:lnum;

pow=pow*2;

num=num/2;

e=num;

for p=0:(num-1)

w=1;

for n=(1+p*pow):(pow/2+p*pow)

b(n)=x(n)+w*x(n+pow/2);

w=w*(w0.^e);

end

w=1;

for n=(1+pow/2+p*pow):(pow+p*pow)

b(n)=x(n-pow/2)-w*x(n);

w=w*(w0.^e);

end

end

x=b;

end

y=b

%FFT算法

Y=abs(y) %幅度值

f=1000*(0:256)/1024; %频率值

%subplot(2,2,2)

figure(2)

plot(f,Y(1:257)) %幅度特性曲线

xlabel('频率轴omega')

ylabel('频率幅值F(omega)')

title('信号频谱','FontSize',10)

%IFFT算法

num=length(y);

lnum=log(num)/log(2);

pow=1;

a=[0,1];

for m=1:(lnum-1)

pow=pow*2;

for n=1:pow

b(2*n-1)=a(n);

b(2*n)=a(n)+pow;

end

a=b;

end

for n=1:1024

a(n)=a(n)+1;

end

for m=1:num

z(m)=y(a(m));

end

y=z

w0=exp(i*2*pi/num);

pow=1;

for m=1:lnum;

pow=pow*2;

num=num/2;

e=num;

for p=0:(num-1)

w=1;

for n=(1+p*pow):(pow/2+p*pow)

b(n)=y(n)+w*y(n+pow/2);

w=w*(w0.^e);

end

w=1;

for n=(1+pow/2+p*pow):(pow+p*pow)

b(n)=y(n-pow/2)-w*y(n);

w=w*(w0.^e);

end

end

y=b;

end

z=b/cc

t=0:0.001:1.023;

%subplot(2,2,3)

figure(3)

plot(t,z)

xlabel('时间轴t')

ylabel('信号值f(t)')

title('正弦波','FontSize',10)

第四章FFT算法的分析

第五章 总结与展望

总结：在本文中我们具体研究了FFT算法的基本理论、实现方法。并通过学习数原理和特性。FFT算法的实现打下了坚实的理论基础。并在MATLAB环境下进行了FFT算法的实现，

参考文献

[1] YUAN Li-guo,FU Xin-chu,YU Rong-zhong.Admissibility Conditions for Symbolic

Sequences in Dynamics of Digital Filter with Two s complement Arithmetic[J].

Journal of Shanghai University,2005年05期.

[2] 郑君里等，离散傅里叶变换以及其他离散正交变换.[J]. 信号与系统(下)，2007年重印.

[4] 赵伶俐.浅谈MATLAB在数字信号处理课程中的应用[J].科技信息，2006年S3期.

[7] 廖宇，郭黎. 基于MATLAB的数字滤波器的设计[J].高科技与产业化，2008年9期.

致谢

经过半年的忙碌和工作，本次毕业设计已经接近尾声，作为一个本科生的毕业设计，由于经验的匮乏，难免有许多考虑不周全的地方，如果没有导师的督促指导，以及一起工作的同学们的支持，想要完成这个设计是难以想象的。

在这里首先要感谢我的导师李雯老师。李老师平日里工作繁多，但在我做毕业设计的每个阶段，从外出实习到查阅资料，设计草案的确定和修改，中期检查，后期详细设计，装配草图等整个过程中都给予了我悉心的指导。我的设计较为复杂烦琐，但是李老师仍然细心地纠正图纸中的错误。除了敬佩李老师的专业水平外，她的治学严谨和科学研究的精神也是我永远学习的榜样，并将积极影响我今后的学习和工作。

然后还要感谢大学四年来所有的老师，为我们打下电子专业知识的基础；同时还要感谢所有的同学们，正是因为有了你们的支持和鼓励。此次毕业设计才会顺利完成。

最后感谢我的母校—湖南理工学院南湖学院大学四年来对我的大力栽培。

30

展开阅读全文