快速傅里叶变换程序设计.doc

沈 阳 工 程 学 院课 程 设 计设计题目： 快速傅里叶变换程序设计 沈阳工程学院课程设计任务书课程设计题目： 快速傅里叶变换程序设计 教研室主任 年 月 日批准1.设计主要内容及要求；编写正弦信号发生器程序。要求：1）研究FFT原理以及利用DSP实现的方法。 2）编写FFT程序。 3）调试程序，观察结果。2.对设计论文撰写内容、格式、字数的要求；（1）.课程设计论文是体现和总结课程设计成果的载体，一般不应少于3000字。（2）.学生应撰写的内容为：中文摘要和关键词、目录、正文、参考文献等。课程设计论文的结构及各部分内容要求可参照沈阳工程学院毕业设计（论文）撰写规范执行。应做到文理通顺，内容正确完整，书写工整，装订整齐。（3）.论文要求打印，打印时按沈阳工程学院毕业设计（论文）撰写规范的要求进行打印。（4）. 课程设计论文装订顺序为：封面、任务书、成绩评审意见表、中文摘要和关键词、目录、正文、参考文献。3.时间进度安排；顺序阶段日期计 划 完 成 内 容备注17月12日教师讲解题目，学生查阅相关资料27月13日确定FFT算法以及程序流程37月14日编写程序47月15日调试程序57月16日撰写论文，程序验收DSP技术课程设计成绩评定表 指 导 教 师 评 审 意 见评价内容具 体 要 求权重评 分加权分调研论证能独立查阅文献,收集资料；能制定课程设计方案和日程安排。0.15432工作能力态度工作态度认真，遵守纪律，出勤情况是否良好，能够独立完成设计工作， 0.25432工作量按期圆满完成规定的设计任务，工作量饱满，难度适宜。0.25432说明书的质量说明书立论正确，论述充分，结论严谨合理，文字通顺，技术用语准确，符号统一，编号齐全，图表完备，书写工整规范。0.55432指导教师评审成绩（加权分合计乘以12） 分加权分合计指 导 教 师 签 名： 年 月 日评 阅 教 师 评 审 意 见评价内容具 体 要 求权重评 分加权分查阅文献查阅文献有一定广泛性；有综合归纳资料的能力0.25432工作量工作量饱满，难度适中。0.55432说明书的质量说明书立论正确，论述充分，结论严谨合理，文字通顺，技术用语准确，符号统一，编号齐全，图表完备，书写工整规范。0.35432评阅教师评审成绩（加权分合计乘以8）分加权分合计评 阅 教 师 签 名： 年 月 日课 程 设 计 总 评 成 绩分快速傅里叶变换程序设计中文摘要数字信号处理 (Digital Signal Processing，DSP)是一门涉及许多科学而又广泛应用于众多领域的新兴学科。步入21世纪以后，社会进入数字化时代，而DSP正是这场数字化的核心。简单的说，数字信号处理器就是把信号用数字符号表示成序列，通过计算机或专用信号处理设备，用数字的数值计算方法进行处理（如滤波、变换、压缩、增强、估计、识别等），以达到提取有用信息便于应用的目的。本次课程设计用的是TMS320C54x系列芯片，TMS320C54x系列芯片是TMS320C5000平台下的定点DSP芯片。54x系列芯提供了低成本、低功耗、高性能的处理能力，在各个领域应用日益广泛。本文就是利用它来实现快速傅里叶变换这种运算的。本次我课程设计的题目是快速傅里叶变换的DSP实现方法，作为数字信号处理的一种算法，快速傅里叶变换日益广泛的应用于实时控制和信号处理等各个领域。快速傅里叶变换(Fast Fourier Transform)是实现离散傅里叶变换(DFT)的一种快速高效的运算方法，是数字信号处理中最为重要的工具之一。它使DFT的运算效率提高12个数量级，为数字信号处理技术应用于各种高速信号的实时处理创造了良好的条件，从而大大推动了数字信号处理技术的发展。本次课程设计完成的是长度为256点的FFT运算，它的运算可以用一个流程图来描述，因为流程图的外形像一只蝴蝶，所以称之为蝶形图，一个蝶形图包含一次复乘，两次复加。通过计算蝶距和旋转因子就可画出每级的蝶形图。本次程序通过编写蝶距和旋转因子的子程序，每次都调用这两个子程序，就能计算出输入的数据所对应的输出。关键词 数字信号处理（DSP），快速傅里叶变换（FFT）目 录中文摘要I1 设计任务描述11.1 设计题目11.2 设计要求11.2.1 设计目的11.3 基本要求12 设计思路22.1 FFT算法的由来22.2 FFT算法原理22.3 基2FFT的蝶形运算流图32.4 时间抽取算法FFT的运算特点42.4.1 原位运算42.4.2 输入、输出序列的倒位序规律52.4.3 蝶距的计算52.4.4 旋转因子的计算62.5 FFT算法的DSP的实现方法62.6 FFT运算中应注意的问题63 软件流程图74 各部分程序设计及参数计算84.1 程序的初始化84.2 位反转子程序84.3 旋转因子的软件实现94.4 实现点复数FFT运算94.4.1 第一级蝶形运算104.4.2 第二级蝶形运算104.4.3 第三级第级蝶形运算114.4.4功率谱计算的实现134.5 参数计算145 程序的调试156 工作过程分析166.1 程序的初始化166.2 位倒序子程序166.3 FFT计算166.4 功率谱的计算17小结18致谢19参考文献20附录A1 程序清单21附录A2 程序图形30I快速傅里叶变换程序设计1 设计任务描述1.1 设计题目 快速傅里叶变换程序设计1.2 设计要求1.2.1 设计目的1）理解FFT的算法以及利用DSP实现的方法。2）能熟练的调试程序并能观察其结果。3）熟悉TMS320C54x系列DSP芯片的软件设计方法。1.3 基本要求1）研究FFT原理以及利用DSP实现的方法。2）编写FFT程序。3）调试程序，观察结果。2 设计思路2.1 FFT算法的由来傅里叶变换是数字信号处理领域中的一种分析工具，它可以将信号从时域变换到频域，称之为傅里叶正变换，也可以把信号从频域变换到时域，称之为傅里叶逆变换。傅里叶变换分为连续傅里叶变换和离散傅里叶变换。离散傅里叶变换简称DFT（Discrete Fourier Transform），是对离散信号进行傅里叶变换的方法，其运算量大、复杂度与变换点数的二次方成正比，因而不适用于进行实时信号处理。为了提高DFT的运算速度，在20世纪60年代由Cooley和Turkey提出了快速傅里叶变换的思想，简称FFT（Fast Fourier Transform），它是一种高效实现DFT的算法，能够明显降低DFT运算的复杂度，使DFT得到了广泛的应用。2.2 FFT算法原理若给定由个信号样本（0），（1），（-1）组成的信号序列（），DFT可用式2-1给出： =0,1，-1 （2-1）式2-1中，称为旋转因子或蝶形因子，=。从中可以看出：当信号样本为复数时，计算单个需经过次复数乘法和-1次复数加法运算，相当于4次实数乘法和2（2-1）次实数加法。完成全部点DFT共需次复数乘法和（-1）复数加法运算。可见，随着不断增加，整个DFT运算量是相当庞大的，而FFT算法通过对计算过程的深入分析，利用旋转因子具有的周期性与对称性，实现了降低运算复杂度的目的。当序列长度为偶数时，信号序列（）可被分解为奇、偶两个子序列，相应的点DFT被分解为两个/2点的DFT： =0,1, ，/2-1 （2-2） =0,1, ，/2-1 （2-3）式（2-2）和（2-3）中，和分别表示（）分解后得到的/2点偶序列点奇序列的DFT。式（2-2）和式（2-3）表明，只要求出和，（）前/2点和后/2点的DFT就得到了，整个序列的DFT也就得到了。这样做的好处是计算点DFT只需要约/2次复数乘法，总运算量约为直接DFT运算量的一半同理，当/2为偶数时，每个/2点的DFT又可被分解成两个/4点的DFT，进一步减少了DFT运算的复杂度。依次类推，直到不能继续分解为止。分解结束时，最小DFT的点数称为称为基数，当=（为正整数）时，经过-1次分解，点DFT最终可被分解为/2个两点的DFT，即得到基数为2的FFT运算，使得DFT所需复数乘法次数降至。2.3 基2FFT的蝶形运算流图基2FFT的蝶形运算过程可用图2-1所示，此时=8，=3。图2-1 8点基2FFT运算过程观察图2-1，根据DFT的基2FFT算法，可以总结出以下几条规律：（1）点FFT运算从输入端开始，逐级进行，共需经过级运算；在第（=1，2，）级中存在个相似的蝶形运算组（除输入数据不同外）；每个组内蝶形运算的个数为，参与每个蝶形运算的两个输入数据相距个点。（2）中间数据的存储，可采用原位存储法。即每次蝶形运算的结果存储在与原数据相同的内存单元内。（3）为了保证输出数据按自然数序排列，在进行FFT之前输入数据需要按照特定的顺序存放，通过位倒序寻址可以满足这种要求。2.4 时间抽取算法FFT的运算特点快速傅立叶变换(FFT)算法基本上分为两大类：按时间抽取的FFT运算和按频率抽取的FFT运算。两者在算法的时间和空间复杂度上是一致的，只是序列在计算前后的排列有所不同。在本论文里，采用的是按时间抽取的FFT(DITFFT)算法。2.4.1 原位运算当数据输入到存储器中以后，每一级运算的结果仍然存储在同一组存储器中，直到最后输出，中间无需其它存储器，这叫原位运算。DITFFT的运算就是原位计算，从图2-2可以看出这种运算是很有规律的，每级（每列）计算都是由N/2个蝶形运算构成的，每一个蝶形结构完成下述基本迭代运算 （2-4）式中：m表示第m列迭代； 则分别为该蝶形单元两个输入数据所在行数。式（2-4）的蝶形运算如图2-2 所示。-1图2-2 按时间抽取算法基本蝶形运算单元由前面介绍过的完整的DIT-FFT运算流程图2-2可以看出，第级蝶形运算的输出数据仅与该级蝶形运算的输入数据有关，与前级蝶形的输入数据无关，且第级蝶形运算的输出数据为第级蝶形运算的输入数据。某任何一个蝶形运算的两个输入节点和的节点变量进行蝶形运算后，得到的结果为该蝶形运算两个输出节点，两节点的节点变量，而和其他节点变量无关，因而可以采用原位运算。这种原位运算结构可节省存储单元，降低设备成本，还可节省找地址的时间。进行点的FFT运算时，只需要个寄存器存储节点变量及/2个寄存器存储/2系数，共个存储单元即可。2.4.2 输入、输出序列的倒位序规律由流程图2-1可以看出，当进行原位运算时，发现当运算完成后，FFT的输出按自然顺序排列在存储单元中，即按，的顺序排列；但是这是输入却不是按自然顺序存储的，而是按，的顺序存入存储单元，这种方式就称之为倒位序。当用二进制表示这个顺序时，它正好是“码位倒置”的顺序。例如，原来的自然顺序应是的地方，现在存放，用二进制码表示这一规律时，则是在处存放，出存放，即将自然顺序的二进制码位倒置过来，第一位码变成最末位码，这样倒置以后的顺序正是输入所需要的顺序。表2-1中列出的是=8时按码位倒置规律所得的顺序，其结果与按时间抽取算法FFT流程图中的输入顺序是一致的。同理，当=256时，亦可以采用同样的方法进行位倒序操作。表2-1 码位倒置顺序自然顺序二进码表示码位倒置倒位序00000000100110042010010230111106410000115101101561100113711111172.4.3 蝶距的计算设=，则整个运算流图中包含级蝶形运算，每一级则有个蝶形单元。蝶距即每个蝶形单元两个输入（出）节点的序号差。以为例，结合图2-1可知共包含3级蝶形运算，每一级蝶形单元的蝶距如下：第一级，蝶距为1，可以看作由得到：第二级，蝶距为2，可以看作由得到；第三级，蝶距为4，可以看作由得到。因此得：第级蝶形单元的蝶距为：。2.4.4 旋转因子的计算由FFT算法原理过程可知，若=，则共有级蝶形运算，各级蝶形运算中旋转因子分别如下：第级的旋转因子为（=0,1，）；第-1级的旋转因子为（=0,1，）；第一级的旋转因子为（=0,1，）。由此可见， 第级蝶形运算中旋转因子为，=0，1，。2.5 FFT算法的DSP的实现方法设FFT运算的输入数据为实数，则2点实数FFT算法的实现步骤为：第一步，把2实数输入序列组合成点的复数序列。然后把该复数序列进行位倒序操作后存储在输入区。第二步，进行的FFT运算。第三步，把点FFT输出拆成2的复数序列，这2的复数序列对应于2点时实数输入序列的DFT输出。第四步，结果输出及功率谱计算。2.6 FFT运算中应注意的问题为了避免可能的结果溢出，在编写程序时我们应该注意对每次蝶形运算的结果都右移一位，即除以2。为了减少FFT的运算时间和充分利用C54xDSP资源，编程时应尽可能多的采用并行指令。3 软件流程图结束程序初始化开始送入数据调入系数表输入数据位码倒置FFT的蝶形运算是否发生溢出？归一化输入数据结束？各图形输出图3-1 256点实序列FFT运算程序流程图4 各部分程序设计及参数计算4.1 程序的初始化输入数据的旋转因子表由文件输入 .dataDATA .space 1024 ;复输出数据的起始地址 .copy mdata.inc ;输入数据N .set 128 ;复数点数LOGN .set 7 ;蝶形级数 为输入数据和旋转因子表定义变量sav_grp .usect tempv,3 ;定义组变量值sav_sin .set sav_grp+1 ;定义旋转因子表索引值sav_idx .set sav_grp+2 ;定义输入数据索引值OUTPUT .usect OUTPUT,256 ;定义输出数据大小BOS .usect stack,0fh ;定义堆栈TOS .usect stack,1 .copy twiddle1.inc ;旋转因子sine表 .copy twiddle2.inc ;旋转因子cosine表在此段程序中设置了复数数据的个数为128，FFT运算的级数为7级，还设置了正弦表和余弦表。4.2 位反转子程序*输入数据位码倒置所使用寄存器定义如下：AR0-位翻转寻址索引AR2-以位翻转顺序指向已处理的数据AR3-指向原始输入数据AR7-数据的起始地址* STM #2*N, BK STM #INPUT,AR3 ;AR3指向第一个输入XR0 STM #DATA,AR7 ;AR7中存储数据的起始地址 MVMM AR7,AR2 ;AR2指向第一个被处理的数据R0 STM #N-1, BRC RPTBD p1end - 1 STM #N,AR0 ;AR0赋值为循环缓冲器的大小的一半 LDM AR3 , A READA *AR2+ ADD #1, A READA *AR2+ MAR *AR3+0B ;位翻转寻址 此段程序实现的是位倒序操作，使用C54xDSP提供的位倒序寻址方式，可以方便的实现位倒序操作。在该寻址方式下，用AR0存放数据，用另外的辅助寄存器ARx存放原始数据的首地址，位倒序操作把AR0中的数加到ARx中，且以向右进位的方式进行，则所得地址将以位倒序的方式产生。4.3 旋转因子的软件实现旋转因子是复数，可表示为： （4-1）由式（4-1）可以看出旋转因子的实部为余弦函数，虚部为正弦函数。为了获得FFT运算中需要的全部旋转因子，需要分别存储正弦表和余弦表，且每个表长度为，对应于0180，同时，采用循环寻址方式对正弦表和余弦表进行寻址。4.4 实现点复数FFT运算当2=256时，128点复数FFT运算过程分成七级、3个阶段实现，这三个阶段是：第一级蝶形运算，第二级蝶形运算，第三级第级蝶形运算，以下分别介绍这三个阶段的蝶形运算。*蝶形运算每一级的所有输出都除以2以防溢出，所使用寄存器定义如下：第一级和第二级蝶形运算 AR0-到下一个蝶形运算的偏移量 AR2-指向第一个蝶形的输入数据PR和PI AR3-指向第二个蝶形的输入数据QR和QI AR7-数据的起始地址剩余级的蝶形运算 AR0-旋转因子表索引 AR1-组计数器 AR2-指向WR（cosine表） AR3-指向WI(sine表) AR6-蝶形运算次数计数器 AR7-蝶形级数计数器*4.4.1 第一级蝶形运算第一级蝶形运算的程序如下： STM #0,BK ;循环缓冲器大小BK=0 LD #-1,ASM ;每一级的输出都除以2 MVMM AR7,AR2 ;AR2指向第一个蝶形的输入数据实部PR STM #DATA+2,AR3 ;AR3指向第二个蝶形的输入数据实部QR STM #N/2-1,BRC LD *AR2,16,A RPTBD s1end-1 ;A:=PR STM #3,AR0 ;块重复 SUB *AR3,16,A,B ;B：=PR-QR ADD *AR3,16,A ;A：=PR+QR STH A,ASM,*AR2+ ;PR:=(PR+QR)/2 ST B,*AR3+ ;QR:=(PR-QR)/2 |LD *AR2,A ;A:=PI SUB *AR3,16,A,B ;B:=PI-QI ADD *AR3,16,A ;A:=PI+QI STH A,ASM,*AR2+0 ;PI:=(PI+QI)/2 ST B,*AR3+0% ;QI:=(PI-QI)/2 |LD *AR2,A ;A:=PR在进行第一级蝶形运算时，仅用到旋转因子，而=1，故此时蝶形运算化简为： （4-2） （4-3）在式（4-2）和（4-3）中，和表示蝶形输入，和表示蝶形输出。由于该蝶形运算仅包含加减运算，容易实现，因此在编程时单独处理了。还有，在程序中用到了许多并行指令，其目的是使程序运行的更快即减少了FFT运算的时间。4.4.2 第二级蝶形运算第二级蝶形运算的程序如下：MVMM AR7,AR2 ;AR2指向PRSTM #DATA+4,AR3 ;AR3指向QRSTM #N/4-1,BRCLD *AR2,16,A ;A:=PR RPTBD s2end-1STM #5,AR0SUB *AR3,16,A,B ;B:=PR-QR ADD *AR3,16,A ;A:=PR+QR STH A,ASM,*AR2+ ;PR:=(PR+QR)/2 ST B,*AR3+ ;QR:=(PR-QR)/2 |LD *AR2,A ;A:=PI SUB *AR3,16,A,B ;B:=PI-QI ADD *AR3,16,A ;A:=PI+QI STH A,ASM,*AR2+ ;PI:=(PI+QI)/2 STH B,ASM,*AR3+ ;QI:=(PI-QI)/2 MAR *AR3+ ADD *AR2,*AR3,A ;A:=PR+QI SUB *AR2,*AR3-,B ;B:=PR-QI STH A,ASM,*AR2+ ;PR:=(PR+QI)/2 SUB *AR2,*AR3,A ;A:=PI-QR ST B,*AR3 ;QR:=(PR-QI)/2 |LD *AR3+,B ;B:=QR ST A,*AR2 ;PI:=(PI-QR)/2 |ADD *AR2+0%,A ;A:=PI+QR ST A,*AR3+0% ;QI:=(PI+QR)/2 |LD *AR2,A ;A:=PR 此时每组蝶形运算用到的旋转因子都是两个：和。参与的蝶形运算同第一级运算时相同，而=，则其参与的蝶形运算可用式4-4和4-5表示为： （4-4） （4-5）由式（4-4）和（4-5）可以看出在第二级运算时也可不涉及乘法运算，因此在编程时也单独处理了。4.4.3 第三级第级蝶形运算第三级到最后一级蝶形运算的程序如下： STM #512,BK ;循环缓冲器大小BK=512 ST #128,sav_sin ;初始化旋转因子表 ;如第3级128，.,第8级为4 STM #128,AR0 ;第3级的旋转因子 STM #TWI2,AR4 ;AR4指向WR STM #TWI1,AR5 ;AR5指向WI STM #-3+LOGN,AR7 ;初始化级计数器 ST #-1+N/8,sav_grp ;初始化组计数器 STM #3,AR6 ;初始化蝶形运算计数器 ST #8,sav_idx ;初始化输让入数据索引stage: STM #DATA,AR2 LD sav_idx,A ADD *(AR2),A STLM A,AR3 MVDK sav_grp,AR1group: MVMD AR6,BRC RPTBD bend-1 LD *AR4,T MPY *AR3+,A MACR *AR5+0%,*AR3-,A ADD *AR2,16,A,B ST B,*AR2 |SUB *AR2+,B ST B,*AR3 |MPY *AR3+,A MASR *AR3,*AR4+0%,A ADD *AR2,16,A,B ST B,*AR3+ |SUB *AR2,B LD *AR4,T ST B,*AR2+ |MPY *AR3+,Abend: PSHM AR0 MVDK sav_idx,AR0 MAR *AR2+0 MAR *AR3+0 BANZD group,*AR1- POPM AR0 MAR *AR3- LD sav_idx,A SUB #1,A,B STLM B,AR6 STL A,1,sav_idx LD sav_grp,A STL A,ASM,sav_grp LD sav_sin,A STL A,ASM,sav_sin BANZD stage,*AR7- MVDK sav_sin,AR0在第三级到最后一级蝶形运算中，组内蝶形运算的次数逐渐增多，而且多是乘法运算，单独处理是不可取的，因此这段运算采用相同的程序块来实现。最终所得的结果可用式（4-6）表示：0，1， （4-6）式（4-6）中，和分别表示复数的实部和虚部。4.4.4功率谱计算的实现 power: STM #OUTPUT,AR3 ;AR3指向输出缓冲地址 STM #255,BRC ;块循环计数器设置为255 RPTBD power_end-1 ;带延迟方式的重复指令执行 STM #DATA,AR2 ;AR2指向AR0 SQUR *AR2+,A ;A:=AR2 SQURA *AR2+,A ;A:=AR2+AI2 STH A,7,*AR3 ;将A中的数据存入输出缓冲中 ANDM #7FFFH,*AR3+ ;避免输出数据过大显示错误power_end: B power_end .end为了便于观察FFT的运算结果，需要求出信号的频谱。经过第三级到最后一级蝶形运算之后，已经得到式（4-7）所示： （4-7）故功率谱可以通过式（4-8）计算得到： （4-8）为了省去开方的复杂运算，此程序中仅求的实部和虚部的二次方和既可，它同样可以反映功率谱的特征。通过软件仿真，可以观察到输入信号的时间波形和频谱波形以及输出信号的功率谱波形。4.5 参数计算由于程序中一共分为512个点，所以每个点所占的度数为360/512。那么假设从0开始计算，第个点的度数就为度,即（360/512）*n。然后求每个点不同度数下的正弦值，即sin（360/512）*n的值。函数值sin（360/512）*n+1都为从0到2的正数。每个字中能存入四位16进制数即7FFFH，换算成十进制数为32767，幅值应为其的一半等于16383.5，所以公式前还应该乘以它的幅值，综上所述，最后的输出结果应该是Y= 16383.5*sin（360/512）*n+1。根据同样的原理可以计算出余弦函数值。5 程序的调试我们编写的每一个程序都需要进行调试，恰好本次课程设计的主要目的就是能正确的调试程序，使程序最终输出正确的结果。在程序的调试过程中我遇到了许多困难。首先，在编译时有两个错误是认为因素造成的，有一个是因为把指令的拼音打错了，还有一个是因为在程序的最开始多留出一行，由于DSP的汇编指令要求非常严格，所以在编写程序时应该注意书写的格式。其次，还有两个错误就是，正弦表和余弦表这两个文件打不开或者是不存在，这两个错误我们几个同学改了很长时间也没找出原因，后来在吕老师的指点下，终于把这两个错误给解决了。我在创建工程时，直接把正弦表和余弦表这两个表放在了ASM文件中了，其实应该在所建工程里再另外创建两个窗口，把正弦表存在一个文件夹中，把余弦表存在一个文件夹中。最后还有一个问题就是，当编译和链接都没有问题时，在load program时发现.out文件又加载不过去了，经过查找发现MEMORY和SECTION这两个地方的指令都是小写的，把小写的字母都改成大写的之后，问题不存在了解决了。这可真的是应了那句话啊：实践出真知啊。6 工作过程分析6.1 程序的初始化256点实数FFT的程序主要由4部分组成，分别是位反转子程序，FFT核心计算子程序，奇偶分离及结果产生子程序，还有功率谱计算子程序。首先在程序的开始应对程序进行初始化，程序的初始化包括为输入数据和旋转因子表定义变量，还有是在此段程序中设置复数数据的个数以及FFT运算的级数，还设置了正弦表和余弦表。6.2 位倒序子程序在编写位倒序子程序时，使用C54xDSP提供的位倒序寻址方式，可以方便的实现位倒序操作。在该寻址方式下，用AR0存放数据，用另外的辅助寄存器ARx存放原始数据的首地址，位倒序操作把AR0中的数加到ARx中，且以向右进位的方式进行，则所得地址将以位倒序的方式产生。6.3 FFT计算当2=256时，128点复数FFT运算过程分成七级、3个阶段实现，这三个阶段是：第一级蝶形运算，第二级蝶形运算，第三级第级蝶形运算，以下分别分析这三个阶段蝶形运算的工作过程。第一级和第二级蝶形运算仅包含加减运算，容易实现，因此在编程时单独处理了。在第三级到最后一级蝶形运算中，组内蝶形运算的次数逐渐增多，而且多是乘法运算，单独处理是不可取的，因此这段运算采用相同的程序块来实现。蝶形运算的主要思想就是，输入数据乘以旋转因子得到下一级的输出，所以在运算之前先要计算出这一级的旋转因子。而旋转因子又是一个复数运算，因此在程序中设置了正弦表和余弦表，通过查这两个表来计算旋转因子。最后来确定FFT的输出值。正弦系数表和余弦系数表由文件（twiddle1.inc和twiddle2.inc）给出，主程序通过.copy 汇编命令将正弦系数和余弦系数与程序代码汇编在一起（也可以用.include命令从twiddle1.inc和twiddle2.inc文件中读入系数。数据文件twiddle1.inc和twiddle2.inc分别给出FFT的正弦系数、余弦系数各512个。利用此系数表可以完成81024点的FFT运算。6.4 功率谱的计算对信号进行傅里叶变换，取sin部分为实部，cos部分为虚部，直接算实部和虚部的平方和，得到的就是功率谱分布。功率信号的功率谱反应了信号功率随频率分布的特点，功率谱是信号先自相关再作FFT变换。在CCS下进行程序调试和结果显示为输入信号时域波形图、输入信号频域波形图和输出信号功率谱。三个图见附录。参考文献1 俞一彪，孙兵. 数字信号处理理论与应用.南京：东南大学出版社，20052 李行一. 数字信号处理.重庆：重庆大学出版社，20023 郭森楙，颜允圣. 数字信号处理器体系结构、实现与应用.北京：清华大学出版社，20054姜沫岐，许涵，俞鹏，段国强 . DSP原理与应用从入门到提高. 北京: 机械工业出版社， 20075 俞卞章. 数字信号处理. 西安: 西北工业大学出版社， 20026 王安国. 数字信号处理基础.北京: 电子工业出版社， 2003附录A1 程序清单fft.asm程序清单： .title cscs.asm .mmregs .global reset,start,sav_sin,sav_idx,sav_grp .def start,_c_int00 .dataDATA .space 1024 .copy mdata.incN .set 128LOGN .set 7sav_grp .usect tempv,3sav_sin .set sav_grp+1sav_idx .set sav_grp+2OUTPUT .usect OUTPUT,256BOS .usect stack,0FHTOS .usect stack,1 .copy twiddle1.inc .copy twiddle2.inc .text_c_int00 B start NOP NOPstart: STM #TOS,SP LD #0,DP SSBX FRCT STM #2*N,BK STM #INPUT,AR3 STM #DATA,AR7 MVMM AR7,AR2 STM #N-1,BRC RPTBD plend-1 STM #N,AR0 LDM AR3,A READA *AR2+ ADD #1,A READA *AR2+ MAR *AR3+0Bp1end: STM #0,BK LD #-1,ASM MVMM AR7,AR2 STM #DATA+2,AR3 STM #N/2-1,BRC LD *AR2,16,A RPTBD slend-1 STM #3,AR0 SUB *AR3,16,A,B ADD *AR3,16,A STH A,ASM,*AR2+ ST B,*AR3+ |LD *AR2,A SUB *AR3,16,A,B ADD *AR3,16,A STH A,ASM,*AR2+0 ST B,*AR3+0% |LD *AR2,As1end: MVMM AR7,AR2 STM #DATA+4,AR3 STM #N/4-1,BRC LD *AR2,16,A RPTBD s2end-1 STM #5,AR0 SUB *AR3,16,A,B ADD *AR3,16,A STH A,ASM,*AR2+ ST B,*AR3+ |LD *AR2,A SUB *AR3,16,A,B ADD *AR3,16,A STH A,ASM,*AR2+ STH B,ASM,*AR3+ MAR *AR3+ ADD *AR2,*AR3,A SUB *AR2,*AR3-,B STH A,ASM,*AR2+ SUB *AR2,*AR3,A ST B,*AR3 |LD *AR3+,B ST A,*AR2 |ADD *AR2+0%,A ST A,*AR3+0% |LD *AR2,As2end: STM #512,BK ST #128,sav_sin STM #128,AR0 STM #TWI2,AR4 STM #TWI1,AR5 STM #-3+LOGN,AR7 ST #-1+N/8,sav_grp STM #3,AR6 ST #8,sav_idxstage: STM #DATA,AR2 LD sav_idx,A ADD *(AR2),A STLM A,AR3 MVDK sav_gr