DFT的快速算法分析及FFT的DSP实现.docx

DFT的快速算法分析及FFT的DSP实现.docx

《DFT的快速算法分析及FFT的DSP实现.docx》由会员分享，可在线阅读，更多相关《DFT的快速算法分析及FFT的DSP实现.docx（20页珍藏版）》请在冰豆网上搜索。

DFT的快速算法分析及FFT的DSP实现

《现代信号处理课程设计》

课程设计报告

设计题目DFT的快速算法分析及FFT的DSP实现

第1章绪论..............................................3

1.1设计背景.............................................3

1.2数字信号处理技术现状.................................3

1.3设计思路简介.........................................5

第2章系统开发平台与环境.................................7

2.1CCS开发环境.........................................7

2.2ICETEKDSP教学实验箱................................7

第3章DFT的快速算法分析与研究.........................8

3.1快速傅立叶变换原理...................................8

3.2FFT快速算法分析.................................9

3.2.1 按时间抽取的基2DIT-FFT…............................9

3.2.2 按频率抽取的基2DIT-FFT..............................10

3.2.3 按频率抽取的基2DIT-FFT.............................10

第4章系统软件设计..................................12

4.1DFT的编程思想......................................12

4.2程序源代码.........................................12

第5章系统仿真..........................................15

参考文献.................................................19

第1章绪论

1.1设计背景

快速傅氏变换（FFT）是离散傅氏变换的快速算法，它是根据离散傅氏变换的奇、偶、虚、实等特性，对离散傅立叶变换的算法进行改进获得的。

它对傅氏变换的理论并没有新的发现，但是对于在计算机系统或者说数字系统中应用离散傅立叶变换，可以说是进了一大步。

设x（n）为N项的复数序列，由DFT变换，任一X（m）的计算都需要N次复数乘法和N-1次复数加法，而一次复数乘法等于四次实数乘法和两次实数加法，一次复数加法等于两次实数加法，即使把一次复数乘法和一次复数加法定义成一次“运算”（四次实数乘法和四次实数加法），那么求出N项复数序列的X（m）,即N点DFT变换大约就需要N2次运算。

当N=1024点甚至更多的时候，需要N2=1048576次运算，在FFT中，利用WN的周期性和对称性，把一个N项序列（设N=2k,k为正整数），分为两个N/2项的子序列，每个N/2点DFT变换需要（N/2）2次运算，再用N次运算把两个N/2点的DFT变换组合成一个N点的DFT变换。

这样变换以后，总的运算次数就变成N+2（N/2）2=N+N2/2。

继续上面的例子，N=1024时，总的运算次数就变成了525312次，节省了大约50%的运算量。

而如果我们将这种“一分为二”的思想不断进行下去，直到分成两两一组的DFT运算单元，那么N点的DFT变换就只需要Nlog2N次的运算，N在1024点时，运算量仅有10240次，是先前的直接算法的1%，点数越多，运算量的节约就越大，这就是FFT的优越性。

数字信号处理器（DSP）是一种可编程的高性能处理器，近年来发展很快．它不仅适用于数字信号处理，而且在图像处理、语音处理、通信等领域得到了广泛的应用．通用的微处理器在运算速度上很难适应信号实时处理的要求．联沪处理器中集成有高速的乘法器硬件，能快速地进行大量数据的乘法和加法运算。

快速傅里叶变换（FFT）的出现使得DFr在实际应用中得到了广泛的应用．

1．2数字信号处理技术现状

信号处理（slgnalProcessing）的本质是信息的安排和提取，将信息通过模拟、数字或光学的方法从各种噪声、干扰的环境中提取出来，并变换为一种便于人或计算机所使用的形式。

通常的待测信号大多是模拟信号，数字信号处理技术先将模拟信号数字化后（包括时域的离散化和幅值的量化处理），再进行数字处理，最后还原成模拟信号，相对于传统的模拟处理技术，数字信号处理的突出优点是精度高、灵活性强、可靠性高、易于大规模集成等等。

在当前信息、控制等很多应用领域都涌现出了所谓的“实时系统”

（Realtimesystem），其特点在于系统的输入、处理和输出等各个阶段都具有绝对的时间限制，若超过所限定的时间，系统将被认为不起作用而失败。

因此处理单元的速度必须要足够快，以满足操作的实时性限制。

如电话网上的一个用于回音消除的语音信号处理系统，处理数据的吞吐率必须要跟得上语音信号的输入和输出速率，否则通信信道将无法工作。

实时数字信号处理问题的结构层次划分如图卜3所示，总的可以分为算法的开发和硬件体系结构的开发。

50年代末期在时域上发展了离散傅立叶变换（DFT）、相关等处理算法，这类算法的复杂度为O（Ⅳ2）。

根据当时计算机的处理速度是做不到实时性的。

到了70年代，大规模集成电路（LsI）技术的发展，产生了微处理器（Micro—processorunit，MPu）和各种基于快速傅立叶变换（FFT）算法的实时信号处理硬件，如高速并行乘法器等；同时，1965年J．w．Colley和J．w．Tukey提出的FFT，将DFT算法的复杂度降为0（Ⅳlog：

Ⅳ）。

以FFT为标志的使用快速算法的提出，促进了专门用来实现数字信号处理的一类MPu的分化和发展，并在此基础上出现了面向实时性、重复性、数值运算密集型的特殊MPu，即所谓的DsP。

DsP为实时数字信号处理提供了廉价、高效、灵活的工具，促进了器件实用化和推广。

DsP是实时数字信号处理技术的核心和标志，从功能上DSP可分为专用和通用两大类，专用DSP一般采用定点的数据格式，结构简单，速度很快；通用DsP的灵活性好、功能强，但处理速度有所降低。

通用DsP片内取指、译码、执行等都采用流水线（Pipeline）技术，缩短了执行时间；针对滤波、矩阵运算、FFT等需要大量乘加法运算的特点，DsP内部配有独立的乘法器和加法器；许多DsP带有直接数据传输（DMA）通道控制器和多种外设接口，提高了数据块传输的速度；DSP还普遍带有中断控制器和定时器，有的甚至还有A／D、D／A、PwM、调制解调器（Modem）等，方便了实时处理系统的构建。

随着近20年来超大规模集成电路（VLsI）技术的发展，DsP的处理速度越来越块，内存越来越大，功能越来越复杂，价格也越来越低，出现了面向通讯、图像处理、多媒体、仪器仪表等各个领域、各个应用层次的产品，成为了电子技术的主要热点。

数字信号处理的应用领域是多种多样的，各种应用的抽样率变化范围很大，但数字信号处理基本上从两个方面来解决信号的处理问题：

一个是时域方法，即数字滤波；另一个是频域方法，即频谱分析。

数字信号处理的任务大致分为三类：

1卷积：

用于各种滤波器，对给定频率范围的原始信号进行加工（通过或滤除）来提高信噪比。

2相关：

用于信号比较，去除随机事件而放大重复的信号。

3变换：

用于分析信号的频率内容，寻找特征“指纹”。

这些信号处理算法的主要特点是：

1包含大量数据。

2同样数据多次使用。

3存在反馈路径。

4都有快速乘加这样的基本运算。

5大量中间数据需要存储和管理。

6涉及大量地址产生。

7算法基本上针对并行运算。

1.3设计思路简介

快速傅立叶变换（FFT）是实现离散傅立叶变换的快速算法，它在信号分析与处理领域有着十分广泛的应用，是数字信号处理中最为重要的工具之一。

高速F”处理技术是FFT技术应用的一个主要方向。

发展高速FFT处理技术是十分必要，因为：

1FFT作为时域和频域转换的基本运算，是数字谱分析的必要条件，在通信等数字处理领域有着极为广泛的应用，部分应用更有高速实时运算的要求。

2处理场合要求实时快速进行FFT运算，如电力系统中电压电流信号的谐波分析、数字滤波、振动分析和故障诊断和图像处理等方面，这些都需要使用专用的FFT计算设备进行实时处理。

通常，提高FFT的运算速度有两种途径：

一种是改进运算工具；另一种是改进FFT算法本身。

所以在实际的应用系统中，一方面我们要通过优化算法来提高运算速度，另一方面选取高效、快速、实用的FFT计算平台也十分关键。

目前，实现FFT算法的计算平台主要有以下三种：

1）利用单片机通过软件编程来实现

这种方法实现起来具有很大的灵活性，但是受到单片机本身性能的限制，例如没有专门的乘法运算部件，所以难以实现高速、大规模的FFT计算。

2）应用专用FFT芯片或用户定制的大规模集成电路来实现

这种方法可以实现很高的运行速度，但是使用中灵活性较差，特别是在使用用户定制的大规模集成电路时，需要很高的开发、研制费用。

3）利用专门的DSP芯片通过软件编程实现

这种方法作为FFT计算的平台具有十分突出的优势。

上个世纪八十年代，专用数字信号处理器（DsP）诞生，用于实现数字信号处理的计算机结构发生了重大变化。

它们大都具备了硬件乘法器，执行一次乘法运算所需的时间与加法相同：

采用多级流水线操作，取指、译码、读数、执行等步骤并行进行；DSP的程序总线和数据总线分开，访问数据和读取程序能够同时进行：

并可在一个指令周期内进行多个操作数的读取等等。

总之，DsP技术的迅猛发展使其以优良的性能价格比和优化的结构为FFT的实现提供了一个有效的平台。

正是由于采用DsP实现快速傅立叶变换（FFT）的上述优越性，本论文采用TI公司的32位浮点DsP—TMs320Vc33为核心处理单元来实现了FFT算法程序。

该芯片具有快速的指令周期，高精度的浮点运算能力，具有专用的硬件乘法器可以在一个指令周期内完成一次乘法和一次加法，满足了FFr算法中有大量的乘法运算的需求，这样大大提高了程序的执行效率，同时，vC33为浮点DsP，在FFl运算中不需要考虑溢出问题，因而具有比定点DsP具有更高的精确度，此外，它还有反序间接寻址指令，它是专门为n叮的变换寻址而设计的，这样可以大大提高程序的执行效率。

本课程设计主要是实现了对各种FFT算法的研究、分析和比较。

编写各种算法实现程序，这些程序可以方便地作为子函数调用，以满足各种应用需要。

第2章系统开发平台与环境

2.1CCS开发环境

CCS提供了配置、建立、调试、跟踪和分析程序的工具，它便于实时、嵌入式信号处理程序的编制和测试，它能够加速开发进程，提高工作效率。

CCS提供了基本的代码生成工具，它们具有一系列的调试、分析能力。

CCS支持如下图1.1所示的开发周期的所有阶段。

图1.1

2.2ICETEKDSP教学实验箱

实验箱部分：

􀁺一个数字信号发生器，可同时提供四种波形、三路输出；信号的波形、频率、幅度可调。

-频率调整：

在每个频率段范围内进行频率调整。

-波形切换：

提供4种波形（方波，三角波，正弦波，白噪声），可通过拨动开关进行选

择。

-幅值微调：

0—3.3V平滑调整。

-信号接插孔：

4路A/D输入（ADCIN0-3），4路D/A输出（DACOUT1-4），每路均提供

信号和地。

􀁺多种直流电源输出。

+5V（5A），+12V（1A），+9V（0.5A），地。

􀁺底板提供插座，可使用插座完成DSP评估板上的A/D信号输入和D/A输出。

􀁺测试模块：

提供14个测试点，可以测量PWM输出、AD输入和DA输出波形。

􀁺双信号发生器设计，更加贴近DSP的实际应用，许多实际的情况都是需要对两个信号

进行相关分析。

通用DSP开发系统部分：

􀁺USB2.0接口开发系统，支持C2000/VC33/C5000/C6000的开发应用。

􀁺支持CCS。

􀁺通用开发系统和DSP控制板分离，有利于将来DSP的升级。

同时，也可以脱离实验箱

单独从事科研开发使用。

通用控制模块部分：

􀁺显示输出：

-液晶显示（LCD）：

128×64点阵图形显示屏，可调整显示对比度。

-发光二极管显示阵列：

8×8点阵。

-发光二极管。

􀁺音频输出：

可由DSPI/O脚控制的蜂鸣器；D/A输出提供音频插座，可直接接插耳机。

􀁺键盘：

17键数字键盘（标准PS2接口）。

􀁺步进电机：

四相步进电机，步距角5.625，起动频率≥300PPS，运行频率≥900PPS，可由DSPI/O端口控制旋转和方向、速度。

􀁺直流电机：

空载转速3050转/分，输出功率1.35W，启动力矩21.3N，可以接收DSP输

出的PWM控制信号，实现电机的转速和方向控制。

􀁺拨动开关（DIP）：

4路，可实现复位和设置DSP应用板参数。

DSP主处理板部分：

􀁺（见第一部分）支持：

ICETEK-VC5416-A板、ICETEK-VC33-A板、ICETEK-VC5416-A

板、ICETEK-VC33-A板、ICETEK-C6713-A板、ICETEK-LF2407-A板、ICETEK-F2407-A

板等。

第3章DFT的快速算法分析与研究

FFT算法的基本原理是把长序列的DFT逐次分解为较短序列的DFT。

按照抽取方式的不同可分为DIT-FFT（按时间抽取）和DIF-FFT（按频率抽取）算法。

按照蝶形运算的构成不同可分为基2、基4、基8以及任意因子（2n,n为大于1的整数），基2、基4算法较为常用。

3.1快速傅立叶变换原理

从定义式（2—2一卜1）计算N点DFT需要Ⅳ2次乘法和N（N一1）次加法运算，当N较大时，运算量是非常大的，其中还不包括计算缈。

所需的运算量，为此寻求其快速算法是非常必要的。

c001ey和Tukey将DFT运算量从O（Ⅳ2）降为0（Ⅳlog：

Ⅳ）。

目前已有各种各样的计算DFT的快速算法。

这些方法大抵可以分为两类，一类是将DFT转变为卷积，利用计算卷积的方法计算，其代表是winograd傅立叶变换算法（wFTA）和素因子（FPA），另一类是递归型算法，是将一维DFT转化为容易计算的二维或多维DFT，且这个过程可以重复，其代表算法是cooley—Tukey算法、Rader—Brenner算法和分裂基算法等。

FFT的基本思想在于，将原有的N点序列分解成两个或更多的较短序列，这些较短序列的DFT可重新组合成元序列的DFT，而总的运算次数却比直接的DFT运算少的多，从而达到提高速度的目的。

这种分解基本上可分为两类：

一类是时间序列x（n）进行逐次分解，称为按时间抽取算法（DecimationInTime）；另一类是将傅立叶变换序列x（k）进行分解，称为按频率抽取算法（DecimationInFrequency）。

3.2FFT算法分析

FFT算法的基本原理是把长序列的DFT逐次分解为较短序列的DFT。

按照抽取方式的不同可分为DIT-FFT（按时间抽取）和DIF-FFT（按频率抽取）算法。

按照蝶形运算的构成不同可分为基2、基4、基8以及任意因子（2n,n为大于1的整数），基2、基4算法较为常用。

3.2.1基2DIT-FFT（按时间抽取）：

令

，

，则有：

蝶形运算单元如下所示：

3.2.2基2DIF-FFT（按频率抽取）：

则有：

蝶形运算单元如下所示：

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

在基2算法中，每个蝶形运算单元都包括1次复数乘法、2次复数加法。

N（N=

）点序列的运算流图应有M级蝶形，每一级都由N/2个蝶形运算组成，所以N点序列的基2FFT算法，总的运算量为

次复数乘法，

次复数加法。

直接DFT运算量为

次复数乘法、

次复数加法。

可见，FFT算法大大减少了运算量，当N越大时，FFT算法的优越性越明显。

3.2.3基4DIF-FFT（按频率抽取）

令：

则有：

蝶形运算单元如下所示：

由上图可知每个基4蝶形运算单元包括3次复数乘法、8次复数加法。

N（N=

，M为偶数）点序列的FFT运算若采用基4算法则有M/2级蝶形，每级由N/4个蝶形运算构成。

采用基4算法计算N点序列的FFT共需要

次复数乘法、

次复数加法。

由于主要的运算时间集中在乘法上面，可见基4算法的运算量较基2算法减少了25%，但运算量的减少是以硬件的复杂性及使用更多资源为代价的。

第4章系统软件设计

4.1DIF-FFT的编程思想

DIF-FFT的运算过程与DIT-FFT有类似规律。

N=点的FFT共进行M级运算，且每一级都由N/2个蝶形运算组成，后一级的节点数据由前一级同处一条水平线位置的节点数据产生，所以我们同样可以将后一级的节点数据储存到前一级的节点中，这样的方法就叫做原位计算。

运用原位计算不仅可以降低成本，简化运算流程，而且更重要的是它节省了大量的内存资源。

如要进行DIF-FFT运算得到正确，必须先将输入序列按照一定的方式进行重新排列。

在DIF-FFT运算中，每一级都由若干个蝶形运算构成，每次蝶形运算都要乘以一个旋转因子，为了方便编写程序，我们要找到其变化的规律。

通过对DIF-FFT运算规律的计算可知：

程序可由三层循环运算构成，他们分别是级数L，序数R，次数K。

在第L级中，每个蝶形运算的两个输入量相距B=2^（L-1）个点，每一级有J个旋转因子，可表示为P=（K-1）*B，且每计算一级，相应的参与运算的点数减少一半,最终得到蝶形运算后的输出序列，进行倒序，则可输出正确的FFT结果。

4.2程序源代码

/*****************fftprograme*********************/

#include"typedef.h"

#include"math.h"

structcompxEE（structcompxb1,structcompxb2）

{

   structcompxb3;

   b3.real=b1.real*b2.real-b1.imag*b2.imag;

   b3.imag=b1.real*b2.imag+b1.imag*b2.real;

   return（b3）;

}

voidFFT（structcompx*xin,intN）

{

   intf,m,nv2,nm1,i,k,j=1,l;

   /*intf,m,nv2,nm1,i,k,j=N/2,l;*/

   structcompxv,w,t;

   nv2=N/2;

   f=N;

   for（m=1;（f=f/2）!

=1;m++）

   {

       ;

   }

   nm1=N-1;

      /*变址运算*/

   for（i=1;i<=nm1;i++）

   {

       if（i

       {

           t=xin[j];

           xin[j]=xin[i];

           xin[i]=t;

       }

       k=nv2;

       while（k

       {

           j=j-k;

           k=k/2;

       }

       j=j+k;

   }

   {

       intle,lei,ip;

       floatpi;

       for（l=1;l<=m;l++）

       {

           le=pow（2,l）;

           //这里用的是L而不是1 

           lei=le/2;

           pi=3.14159;

           v.real=1.0;

           v.imag=0.0;

           w.real=cos（pi/lei）;

           w.imag=-sin（pi/lei）;

           for（j=1;j<=lei;j++）

           {

               /*doublep=pow（2,m-l）*j;

                 doubleps=2*pi/N*p;

                 w.real=cos（ps）;

                 w.imag=-sin（ps）;*/

               for（i=j;i<=N;i=i+le）

               {

                   /* w.real=cos（ps）;

                     w.imag=-sin（ps）;*/

                   ip=i+lei;

                   t=EE（xin[ip],v）;

                   xin[ip].real=xin[i].real-t.real;

                   xin[ip].imag=xin[i].imag-t.imag;

                   xin[i].real=xin[i].real+t.real;

                   xin[i].imag=xin[i].imag+t.imag;

               }

               v=EE（v,w）;

           }

       }

   }

   return;

}

/*****************mainprograme********************/

#include

#include

#include

#include"typedef.h"

floatresult[257];

structcompxs[257];

intNum=256;

constfloatpp=3.14159;

main（）

{

   inti=1;

   for（;i<0x101;i++）

   {

       s[i].real=sin（pp*i/32）;

       s[i].imag=0;

   }

   FFT（s,Num）;

   for（i=1;i<0x101;i++）