《dsp基于matlab的fft算法实现精品资料.doc》由会员分享，可在线阅读，更多相关《dsp基于matlab的fft算法实现精品资料.doc（33页珍藏版）》请在taowenge.com淘文阁网|工程机械CAD图纸|机械工程制图|CAD装配图下载|SolidWorks_CaTia_CAD_UG_PROE_设计图分享下载上搜索。
1、课程设计说明书目 录1 摘要12 设计目的和内容23 基2 DIT-FFT算法33.1 DIT-FFT算法的基本原理43.2 DIT-FFT算法的运算规律及编程思想43.3 原位计算53.4 倒序计算53.5 蝶形运算94 MATLAB运行界面图124.1 fs=1000;n=2000时的原始的语音信号时域图124.1 fs=1000;n=2000时的原始的语音信号频域图124.3 原始语音信号FFT频谱图与原始语音信号自编FFT频谱图比较144.4 原始语音信号FFT频谱图与原始语音信号自编FFT频谱图比较155 设计总结16参考文献19附录201 引言 傅里叶变换在信号处理中具有十分重要的
2、作用，但是基于离散时间的傅里叶变换具有很大的时间复杂度，根据傅里叶变换理论，对一个有限长度且长度为的离散信号，做傅里叶变换的时间复杂度为，当很大时，其实现的时间是相当惊人的（比如当为时，其完成时间为（为计算机的时钟周期），故其实现难度是相当大的，同时也严重制约了DFT 在信号分析中的应用，故需要提出一种快速的且有效的算法来实现。 tt 正是鉴于DFT 极其复杂的时间复杂度，1965 年和巧妙地利用因子的周期性和对称性，提出了一个DFT 的快速算法，即快速傅里叶变换（FFT），从而使得DFT 在信号处理中才得到真正的广泛应用。 本文基于时间抽选奇偶分解，利用Matlab 软件实现快速傅里叶变换。
3、基于所编的FFT 源程序应用的一个实例，本文对有限长度离散时间和连续时间信号进行频谱分析。 DFT是一种应用广泛的数学变换工具，MATLAB是一款功能强大的科学计算语言。MATLAB提供的fft函数解决了DFT的快速计算问题，但由于它是内建函数而不能了解到软件实现的过程。文章以按时间抽取的基2FFT算法为例，根据快速傅里叶变换的原理和规律，绘出了算法实现的程序框图，列出了MATLAB环境下软件实现的程序，建立了从算法理论到程序实现的完整概念。 在信号处理中，DFT（离散傅里叶变换）的计算具有举足轻重的地位。但是基于其复杂的计算，直接应用起来十分麻烦，基于此，本文利用Matlab 软件对有限长度
4、信号的DFT 进行改进，提出FFT（快速傅里叶变换），并利用FFT 对所给连续时间和离散时间信号做了频谱分析。 语音信号的处理与滤波的设计主要是用MATLAB作为工具平台，设计中涉及到声音的录制、播放、存储和读取，语音信号的抽样、频谱分析，滤波器的设计及语音信号的滤波，通过数字信号处理课程的理论知识的综合运用。从实践上初步实现对数字信号的处理。 2 设计目的和内容 MATLAB全称是Matrix Laboratory，是一种功能强大、效率高、交互性好的数值和可视化计算机高级语言，它将数值分析、矩阵运算、信号处理和图形显示有机地融合为一体，形成了一个极其方便、用户界面友好的操作环境。经过多年的发
5、展，已经发展成为一种功能全面的软件，几乎可以解决科学计算中所有问题。MATLAB软件还提供了非常广泛和灵活的用于处理数据集的数组运算功能。综合运用本课程的理论知识进行频谱分析以及滤波器设计，通过理论推导得出相应结论，并利用MATLAB作为工具进行实现，从而复习巩固课堂所学的理这次课程设计的主要目的是综合运用本课程的理论知识进行频谱分析以及滤波器设计，通过理论推导得出相应结论，并利用MATLAB或者开发系统作为工具进行实现，从而复习巩固课堂所学的理论知识，提高对所学知识的综合应用能力，并从实践上初步实现对数字信号的处理。通过对声音的采样，将声音采样后的频谱与滤波。MATLAB全称是Matrix
6、Laboratory，是一种功能强大、效率高、交互性好的数值和可视化计算机高级语言，它将数值分析、矩阵运算、信号处理和图形显示有机地融合为一体，形成了一个极其方便、用户界面友好的操作环境。经过多年的发展，已经发展成为一种功能全面的软件，几乎可以解决科学计算中所有问题。MATLAB软件还提供了非常广泛和灵活的用于处理数据集的数组运算功能。在本次课程设计中，主要通过MATLAB来编程对语音信号处理与滤波，设计滤波器来处理数字信号并对其进行分析。通过理论推导得出相应结论，并利用MATLAB或者开发系统作为工具进行实现，从而复习巩固课堂所学的理论知识，提高对所学知识的综合应用能力，并从实践上初步实现对
7、数字信号的处理。通过对声音的采样，将声音采样后的频谱与滤波。 录制一段个人自己的语音信号，并对录制的信号进行采样；画出采样后语音信号的时域波形和频谱图；在Matlab环境下编写基2 DIT-FFT算法；利用自己编写的算法对已采集的语音信号进行频谱分析，并画出语音信号的时域与频谱图，并与Matlab数字信号处理工具箱中的fft函数进行对比研究，验证自编算法的正确性。 分析和处理音频信号，首先要对声音信号进行采集，MATLAB的数据采集工具箱提供了一整套命令和函数，通过调用这些函数和命令，可以直接控制声卡进行数据采集。windows自带的录音机程序也可驱动声卡来采集语音信号，并能保存为wav格式文
8、件，提供MATLAB相关函数直接读取、写入和播放。课设以wav格式音频信号作为分析处理的输入数据。用MATLAB处理音频信号的基本流程是：现将wav格式音频信号经wavread函数转化成MATLAB列数组变量；再用MATLAB强大的运算能力进行数据分析和处理，如时域分析、频域分析、数字滤波、信号合成、信号变换等等；处理后的数据如是音频数据，可用sound、wavplay等函数直接回放。我们可以利用前面设计好的各种滤波器来实现对某些声音信号的滤波，从而改变输入信号所包含的频率成分的相对比例或滤出某些频率成分。下面介绍一下具体的操作过程。(1)利用wavread函数把指定路径下的WAV格式的声音文
9、件读入矢量中，采样率是22050Hz。x,fs,bits=wavread(D:3.wav);(2)利用sound函数以同样的码率对原声进行再现。sound(x,22050);(3)利用filter函数对样本信号进行滤波处理，以设计好的FIR滤波器为例，z=filter(bd,ad,x); (4)利用sound函数对滤波之后的信号进行再现。sound(z);(5)利用plot函数分别作出信号在滤波前后的时域波形图，以作对比。plot(x(:,1); plot(z(:,1);(6)利用fft函数对滤波前后的信号进行1024点的快速傅里叶变换FFT，采用一次基2按时间抽取。X=fft(x);Z=ff
10、t(z);(7)利用plot函数分别作出信号在滤波前后的频域波形图，以作对比。plot(f,abs(X(1:n/2); plot(f,abs(Z(1:n/2);MATLAB全称是Matrix Laboratory，是一种功能强大、效率高、交互性好的数值和可视化计算机高级语言，它将数值分析、矩阵运算、信号处理和图形显示有机地融合为一体，形成了一个极其方便、用户界面友好的操作环境。经过多年的发展，已经发展成为一种功能全面的软件，几乎可以解决科学计算中所有问题。MATLAB软件还提供了非常广泛和灵活的用于处理数据集的数组运算功能。在本次课程设计中，主要通过MATLAB来编程对语音信号处理与滤波，设计
11、滤波器来处理数字信号并对其进行分析。通过理论推导得出相应结论，并利用MATLAB或者开发系统作为工具进行实现，从而复习巩固课堂所学的理论知识，提高对所学知识的综合应用能力，并从实践上初步实现对数字信号的处理。通过对声音的采样，将声音采样后的频谱与滤波。 录制一段个人自己的语音信号，并对录制的信号进行采样；画出采样后语音信号的时域波形和频谱图；在Matlab环境下编写基2 DIT-FFT算法；利用自己编写的算法对已采集的语音信号进行频谱分析，并画出语音信号的时域与频谱图，并与Matlab数字信号处理工具箱中的fft函数进行对比研究，验证自编算法的正确性。 3 基2 DIT-FFT算法3.1 FF
12、T算法 有限长序列x(n)的N点DFT定义为：，式中，其整数次幂简称为旋转因子。直接进行DFT运算大约需要次三角函数计算、次实数乘法计算和次实数加法计算，且需许多索引和寻址操作2。文3列出了直接DFT的MATLAB程序，这种直接DFT运算概念清楚、编程简单，但占用内存大、运算速度低，在实际工作中并不实用。基2FFT算法的基本思想是把原始的N点序列依次分解成一系列短序列，充分利用旋转因子的周期性和对称性，分别求出这些短序列对应的DFT，再进行适当的组合，得到原N点序列的DFT，最终达到减少运算次数，提高运算速度的目的。按时间抽取的基2FFT算法，先是将N点输入序列x(n)在时域按奇偶次序分解成2
13、个N/2点序列x1(n)和x2(n)，再分别进行DFT运算，求出与之对应的X1(k)和X2(k)，然后利用图1所示的运算流程进行蝶形运算，得到原N点序列的DFT。只要N是2的整数次幂，这种分解就可一直进行下去，直到其DFT就是本身的1点时域序列。一个完整的8点DIT-FFT运算流程如图2所示4。图中的输入序列不再是顺序排列但有规律可循，数组A(存储地址)用于存放输入数据和每级运算的结果。3.2 DIT-FFT算法的运算规律及编程思想 并绘出程序框图。由图2可知，DIT-FFT算法的运算过程很有规律。直接DFT的MATLAB程序，这种直接DFT运算概念清楚、编程简单，但占用内存大、运算速度低，在
14、实际工作中并不实用。基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点时域序列。一个完整的8点DIT-FFT运算流程如图
15、2所示4。图中的输入序列不再是顺序排列但有规律可循，数组A(存储地址)用于存放输入数据和每级运算的结果。直接DFT的MATLAB程序，这种直接DFT运算概念清楚、编程简单，但占用内存大、运算速度低，在实际工作中并不实用。基2FFT算法的基本思想是把原始的N点序列依次分解成一系列短序列，充分利用旋转因子的周期性和对称性，分别求出这些短序列对应的DFT，再进行适当的组合，得到原N点序列的DFT，最终达到减少运算次数，提高运算速度的目的。按时间抽取的基2FFT算法，先是将N点输入序列x(n)在时域按奇偶次序分解成2个N/2点序列x1(n)和x2(n)，再分别进行DFT运算，求出与之对应的X1(k)和
16、X2(k)，然后利用图1所示的运算流程进行蝶形运算，得到原N点序列的DFT。只要N是2的整数次幂，这种分解就可一直进行下去，直到其DFT就是本身的1点时域序列。一个完整的8点DIT-FFT运算流程如图2所示4。图中的输入序列不再是顺序排列但有规律可循，数组A(存储地址)用于存放输入数据和每级运算的结果。直接DFT的MATLAB程序，这种直接DFT运算概念清楚、编程简单，但占用内存大、运算速度低，在实际工作中并不实用。基2FFT算法的基本思想是把原始的N点序列依次分解成一系列短序列，充分利用旋转因子的周期性和对称性，分别求出这些短序列对应的DFT，再进行适当的组合，得到原N点序列的DFT，最终达
17、到减少运算次数，提高运算速度的目的。按时间抽取的基2FFT算法，先是将N点输入序列x(n)在时域按奇偶次序分解成2个N/2点序列x1(n)和x2(n)，再分别进行DFT运算，求出与之对应的X1(k)和X2(k)，然后利用图1所示的运算流程进行蝶形运算，得到原N点序列的DFT。只要N是2的整数次幂，这种分解就可一直进行下去，直到其DFT就是本身的1点时域序列。一个完整的8点DIT-FFT运算流程如图2所示4。图中的输入序列不再是顺序排列但有规律可循，数组A(存储地址)用于存放输入数据和每级运算的结果。3.3 原位计算对点的FFT共进行M级运级由N/2个蝶形运算组成。在同一级中，每个蝶的输入数据只
18、对本蝶有用，且输出节点与输入节点在同一水平线上，得到原N点序列的DFT，最终达到减少运算次数，提高运算速度的目的。按时间抽取的基2FFT算法，先是将N点输入序列x(n)在时域按奇偶次序分解成2个N/2点序列x1(n)和x2(n)，再分别进行DFT运算，求出与之对应的X1(k)的基本思想是把原始的N点序列依次分解成一系列短序列，充分利用旋转因子的周期性和对称性，分别求出这些短序列对应的DFT，再进行适当的组合，得到原N点序列的DFT，最终达到减少运算次数，提高运算速度的目的。按时间抽取的基2FFT算法，先是将N点输入序列x(n)在时域按奇偶次序分解成2个N/2点序列x1(n)和x2(n)，再分别
19、进行DFT运算，求出与之对应的X1(k)元素(存储单元)，这种原位(址)计算的方法可节省大量内存。先是将N点输入序列x(n)在时域按奇偶次序分解成2个N/2点序列x1(n)和x2(n)，再分别进行DFT运算，求出与之对应的X1(k)元素(存储单元)，这种原位(址)计算的方法可节省大量内存。这样数据倒序后的运算可用三重循环程序实现。图1运算流程图3.4 倒序运算快速傅里叶变换（FFT）是为提高DFT运算速度而采用的一种算法。对一个有限长度序列x(n)的N点的DFT为：X（k）=x（n）WknN (k=0,1,，N-1；n=0,1,，N-1;W=e-j2/N)当N=4时，X（k）可展开为：从上式可
20、以看出,要求4点的DFT,需要16次的复数乘法运算,12次复数乘法运算算。由此类推，要求出N点的DFT，需要N2次复数乘法运算，N*(N-1)次复数加法运算。当N值较大时，要完成的复数乘法运算和复数加法运算得次数都非常多，无论是用通用计算机还是用DSP芯片，都需要消耗大量的时间，不适合于对实时处理要求高的场合。为了能实时处理DFT,要想减少DFT的运算量可以有两个途径：第一是降N，N的值减小了，运算量就减少了；第二是利用旋转因子的周期性和对称性，可约性。利用这两个途径实现DFT的快速傅里叶变换（FFT），FFT算法基本上可分为时域抽取法和频域抽取法。当N值较大时，要完成的复数乘法运算和复数加法
21、运算得次数都非常多，无论是用通用计算机还是用DSP芯片，都需要消耗大量的时间，不适合于对实时处理要求高的场合。为了能实时处理DFT,要想减少DFT的运算量可以有两个途径：第一是降N，N的值减小了，运算量就减少了；第二是利用旋转因子的周期性和对称性，可约性。输入序列先按自然顺序存入存储单元,然后经变址运算来实现倒位序排列,用J表示倒序的十进制数，对N=2M，M位的二进制数从左到右各位数权值位N/2，N/4，N/82，1。因此，最高位加1相当于J+N/2。.如果最高位为0，则直接得到下一个倒序值，J+N/2；.如果最高位为1，则最高位为0(J-N/2)，次高位加1(J+N/4)。.以此类推，直到最
22、后一位二进制数字。LH=N/2N1=N-1J=0T=x(J+1)x(J+1)=x(I+1) x(I+1)=TI=0：N1I=JK=LHJKJ=J+KJ=J-KK=K/2YNYN倒序的流程图 图2 倒序的流程图3.5 蝶形运算 实现FFT运算的核心是蝶形运算，找出蝶形运算的规律是编程的基础。蝶形运算是分级进行的；每级的蝶形运算可以按旋转因子的指数大小排序进行；如果指数大小一样则可从上往下依次蝶算。对点的FFT共有M级运算，用L表示从左到右的运算级数(L=1，2，M )。第L级共有个不同指数的旋转因子，用R表示这些不同指数旋转因子从上到下的顺序(R=0，1，B-1)。图3 8点DIT-FFT运算流
23、程第一节点标号k从R开始，由于本级中旋转因子指数相同的蝶共有个，且这些蝶的相邻间距为，故旋转因子指数为P的最后一个蝶的第一节点标号k为：，本级中各蝶的第二个节点与第一个节点都相距B点。应用原位计算，蝶形运算可表示成如下形式。 总结上述运算规律，可采用如下运算方法进行DIT-FFT运算。首先读入数据，根据数据长度确定运算级数M，运算总点数，不足补0处理。然后对读入数据进行数据倒序操作。数据倒序后从第1级开始逐级进行，共进行M级运算。在进行第L级运算时，先算出该级不同旋转因子的个数(也是该级中各个蝶形运算两输入数据的间距)，再从R=0开始按序计算，直到R=B-1结束。每个R对应的旋转因子指数，旋转
24、因子指数相同的蝶从上往下依次逐个运算，各个蝶的第一节点标号k都是从R开始，以为步长，到(可简取极值N-2)结束。考虑到蝶形运算有两个输出，且都要用到本级的两个输入数据，故第一个输出计算完毕后，输出数据不能立即存入输入地址，要等到第二个输出计算调用输入数据完毕后才能覆盖。这样数据倒序后的运算可用三重循环程序实现。整个运算流程如图3所示。蝶形运算的流程图如图4：开始开始开始J=0：B-1N=2MB=2L/2L=1:M输入x(n)K=J+1:2L:NP=J*2(M-L)T=x(K)+x(K+B)* WNp x(K+B)=x(K)-x(K+B)* WNp x(K)=T倒序结束输出 1 2 3图4 蝶形
25、运算的流程图4 MATLAB运行界面图 4.1 fs=1000;n=2000时的原始的语音信号时域图图5 语音信号时域图 4.2 fs=1000;n=2000时的原始语音信号频率特性图图6 频率特性图4.3原始语音信号和自编FFT频谱图比较图7 FFT频谱图比较 4.4 原始语音信号和自编FFT频谱图比较图 8 FFT频谱图比较 5 设计总结 在之前数字信号处理的学习以及完成实验的过程中，已经使用过MATLAB，对其有了一些基础的了解和认识，通过这次的课程设计使我进一步了解了信号的产生，采样及频谱分析的方法，以及其中产生信号和绘制信号的基本命令和一些基础编程语言。让我感受到只有在了解课本知识的
26、前提下，才能更好的应用这个工具，并且熟练的应用MATLAB也可以很好的加深我对课程的理解，方便我的思维。这次课程设计使我了解了MATLAB的使用方法，提高了自己的分析和动手实践能力。同时我相信，进一步加强对MATLAB的学习与研究对我今后的学习将会起到很大的帮助。总的来说，这次程序设计让我获益匪浅，对MATLAB也有了进一步的认识：MATLAB功能强、使用灵活方便等。MATLAB是在国内外广泛使用的一种数字信号处理系统，相信除了以上优点，还有许多我还未发现，希望能在以后的学习中有更深入的认识。通过本次数字信号处理课程设计，我觉得对自己提高很大：克服了的偷懒的毛病，这在我以后的学习和工作中的心理
27、定位与调节有很大的帮助。我感受到了编程是一项非常烦琐周密的活动，不但需要一个人周密的思考问题的能力，处理问题的能力，还需要有足够的耐心和严谨治学的作风，来不得半点马虎。DFT是一种应用广泛的数学变换工具，MATLAB是一款功能强大的科学计算语言。MATLAB提供的fft函数解决了DFT的快速计算问题，但由于它是内建函数而不能了解到软件实现的过程。文章以按时间抽取的基2FFT算法为例，根据快速傅里叶变换的原理和规律，绘出了算法实现的程序框图，列出了MATLAB环境下软件实现的程序，建立了从算法理论到程序实现的完整概念。基于数字滤波和快速傅立叶变换的数字信号处理在上世纪70年代有了飞速的发展，到了
28、80年代，数字信号处理已经开始应用到了各个工程技术领域，不管在军用还是民用系统中都发挥了积极的作用。进入90年代，微电子工业飞速发展，DSP技术已经成为工程实用技术，DSP技术的成本大幅度降低，DSP处理器由通用的芯片向专用的DSP芯片发展，价格也由80年代的数千美元单片降到数十美元每片，甚至更低。以DSP为核心芯片的处理系统日益变成了数字信号处理系统的主流。它广泛用于电子信息、通信、图像处理、语音处理、生物医学、自动控制、地质探测等领域，受到工程设计和使用人员的青睐。MATLAB，它是美国Math Works公司推出的一种面向工程和科学计算的交互式计算软件。它以矩阵运算为基础，把计算、可视化
29、、程序设计融合在一个简单易用的交互式工作环境中，是一款数据分析和处理功能都非常强大的工程适用软件。通过本次实习我们学会了分析和处理音频信号，首先要对声音信号进行采集，MATLAB的数据采集工具箱提供了一整套命令和函数，通过调用这些函数和命令，可直接控制声卡进行数据采集。Window自带的录音机程序也可驱动声卡来采集语音信号，并能保存为WAV格式文件，供MATLAB相关函数直接读取、写入或播放。MATLAB语言是一种数据分析和处理功能十分强大的计算机应用软件，它可以将声音文件变换位离散的数据文件，然后利用其强大的矩阵运算能力处理数据，如数据滤波、傅立叶变换、时域和频域分析、声音回放以及各种图的呈
30、现等，它的信号处理与分析工具箱位语音信号分析提供了十分丰富的功能函数，利用这些功能函数可以快捷而又方便的完成语音信号的处理和分析以及信号的可视化，是人机交互更加便捷。信号处理是MATLAB重要应用的领域之一。语音信号处理是研究用数字信号处理技术和语音学知识对语音信号进行处理的新兴的学科，是目前发展最为迅速的信息科学研究领域的核心技术之一。通过语音传递信息是人类最重要、最有效、最常用和最方便的交换信息形式。语音信号的处理与滤波的设计主要是用MATLAB作为工具平台，设计中涉及到声音的录制、播放、存储和读取，语音信号的抽样、频谱分析，滤波器的设计及语音信号的滤波，通过数字信号处理课程的理论知识的综
31、合运用。从实践上初步实现对数字信号的处理。 在信号处理中，DFT（离散傅里叶变换）的计算具有举足轻重的地位。但是基于其复杂的计算，直接应用起来十分麻烦，基于此，本文利用Matlab 软件对有限长度信号的DFT 进行改进，提出FFT（快速傅里叶变换），并利用FFT 对所给连续时间和离散时间信号做了频谱分析。 傅里叶变换在信号处理中具有十分重要的作用，但是基于离散时间的傅里叶变换具有很大的时间复杂度，根据傅里叶变换理论，对一个有限长度且长度为的离散信号，做傅里叶变换的时间复杂度为，当很大时，其实现的时间是相当惊人的（比如当为时，其完成时间为（为计算机的时钟周期），故其实现难度是相当大的，同时也严重
32、制约了DFT 在信号分析中的应用，故需要提出一种快速的且有效的算法来实现。 t t 正是鉴于DFT 极其复杂的时间复杂度，1965 年和巧妙地利用因子的周期性和对称性，提出了一个DFT 的快速算法，即快速傅里叶变换（FFT），从而使得DFT 在信号处理中才得到真正的广泛应用。 本文基于时间抽选奇偶分解，利用Matlab 软件实现快速傅里叶变换。基于所编的FFT 源程序应用的一个实例，本文对有限长度离散时间和连续时间信号进行频谱分析。 有限长序列x(n)的N点DFT定义为：，式中，其整数次幂简称为旋转因子。直接进行DFT运算大约需要次三角函数计算、次实数乘法计算和次实数加法计算，且需许多索引和寻
33、址操作。文3列出了直接DFT的MATLAB程序，这种直接DFT运算概念清楚、编程简单，但占用内存大、运算速度低，在实际工作中并不实用。基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
34、的整数次幂，这种分解就可一直进行下去，直到其DFT就是本身的1点时域序列。一个完整的8点DIT-FFT运算流程如图2所示4。图中的输入序列不再是顺序排列但有规律可循，数组A(存储地址)用于存放输入数据和每级运算的结果。参考文献1范寿康 DSP技术与DSP芯片.北京：电子工业出版社2程佩青.数字信号处理教程.北京：清华大学出版社出版，20013高西全, 丁玉美等 .数字信号处理. 北京：电子工业出版社，20094 余成波,陶红艳.数字信号处理及MATLAB实现.北京：清华大学出版社，2008 5（美）Edward W.Kamen, Bonnie S.Heck 著,高强译. 信号与系统基础教程，北
35、京：电子工业出版社，2007 6 曹弋,赵阳.MATLAB. 实用教程，北京：电子工业出版社，2007 信号与系统课程设计。 附录x1,fs,bits=wavread(song.wav,N); %读取语音信号的数据，赋给变量x1sound(x1,2*N); %播放语音信号x1=reshape(x1,1,2*n);y1=fft(x1); figure(1)plot(x1) %做原始语音信号的时域图形title(原始语音信号)xlabel(n);ylabel(幅值);M=nextpow2(x1); % 求x的长度对应的2的最低幂次m N=2M;if length(x1)N x1=x1,zeros(
36、1,N-length(x1); % 若x的长度不是2的幂，补零到2的整数幂 end NV2=N/2;NM1=N-1;I=0;J=0;while INM1 if IJ T=x1(J+1); x1(J+1)=x1(I+1); x1(I+1)=T; end K=NV2; while K=J J=J-K; K=K/2; end J=J+K; I=I+1;end %x1; y=x1; % 将x倒序排列作为y的初始值 WN=exp(-i*2*pi/N); for L=1:M B=2L/2;%第L级中，每个蝶形的两个输入数据相距B个点，每级有B个不同的旋转因子 for J=0:B-1 % J代表了不同的旋转
37、因子 p=J*2(M-L); WNp=WNp; for k=J+1:2L:N % 本次蝶形运算的跨越间隔为2L kp=k+B; % 蝶形运算的两个因子对应单元下标的关系 t=y(kp)*WNp; % 蝶形运算的乘积项 y(kp)=y(k)-t; % 蝶形运算， 注意必须先进行减法运算，然后进行加法运算，否则要使用中间变量来传递y(k) y(k)=y(k)+t; % 蝶形运算 end end end %y figure(2)x1,w1=freqz(x1,1) ; %绘制原始语音信号的频率图plot(w1/pi,20*log10(abs(x1);title(频率特性图)xlabel(归一化频率);
38、ylabel(幅度/DB);figure(3)subplot(2,1,1);plot(abs(y1) %做原始语音信号的FFT频谱图title(原始语音信号FFT频谱)xlabel(K);ylabel(Y1(k);subplot(2,1,2);plot(abs(y);title(语音信号FFT频谱) xlabel(K);ylabel(Y(k); 附录资料：matlab画二次曲面一、螺旋线1.静态螺旋线a=0:0.1:20*pi;h=plot3(a.*cos(a),a.*sin(a),2.*a,b,linewidth,2);axis(-50,50,-50,50,0,150);grid onset
39、(h,erasemode,none,markersize,22);xlabel(x轴);ylabel(y轴);zlabel(z轴);title(静态螺旋线); 2.动态螺旋线t=0:0.1:10*pi;i=1;h=plot3(sin(t(i),cos(t(i),t(i),*,erasemode,none);grid onaxis(-2 2 -2 2 0 35)for i=2:length(t) set(h,xdata,sin(t(i),ydata,cos(t(i),zdata,t(i); drawnow pause(0.01)endtitle(动态螺旋线);(图略) 3.圆柱螺旋线t=0:0.
40、1:10*pi;x=r.*cos(t);y=r.*sin(t);z=t;plot3(x,y,z,h,linewidth,2);grid onaxis(square)xlabel(x轴);ylabel(y轴);zlabel(z轴);title(圆柱螺旋线) 二、旋转抛物面b=0:0.2:2*pi;X,Y=meshgrid(-6:0.1:6);Z=(X.2+Y.2)./4;meshc(X,Y,Z);axis(square)xlabel(x轴);ylabel(y轴);zlabel(z轴);shading flat;title(旋转抛物面)或直接用：ezsurfc(X.2+Y.2)./4) 三、椭圆柱
41、面load clownezsurf(2*cos(u),4*sin(u),v,0,2*pi,0,2*pi)view(-105,40) %视角处理shading interp %灯光处理colormap(map) %颜色处理grid on %添加网格线axis equal %使x,y轴比例一致xlabel(x轴);ylabel(y轴);zlabel(z轴);shading flat;title(椭圆柱面) %添加标题四、椭圆抛物面b=0:0.2:2*pi;X,Y=meshgrid(-6:0.1:6);Z=X.2./9+Y.2./4;meshc(X,Y,Z);axis(square)xlabel(x
42、轴);ylabel(y轴);zlabel(z轴);shading flat;title(椭圆抛物面)或直接用：ezsurfc(X.2./9+Y.2./4)五、双叶双曲面ezsurf(8*tan(u)*cos(v),8.*tan(u)*sin(v),2.*sec(u),-pi./2,3*pi./2,0,2*pi)axis equalgrid onaxis squarexlabel(x轴);ylabel(y轴);zlabel(z轴);shading flat;title(双叶双曲面)六、双曲柱面load clownezsurf(2*sec(u),2*tan(u),v,-pi/2,pi/2,-3*p
43、i,3*pi)hold on %在原来的图上继续作图ezsurf(2*sec(u),2*tan(u),v,pi/2,3*pi/2,-3*pi,3*pi)colormap(map)shading interpview(-15,30)axis equalgrid onaxis equalxlabel(x轴);ylabel(y轴);zlabel(z轴);shading flat;title(双曲柱面)七、双曲抛物面（马鞍面）X,Y=meshgrid(-7:0.1:7);Z=X.2./8-Y.2./6;meshc(X,Y,Z);view(85,20)axis(square)xlabel(x轴);ylabel(y轴);zlabel(z轴);shading flat;title(双曲抛物面)或直接用：ezsurfc(X.2./8-Y.2./6) 八、抛物柱面X,Y=meshgrid(-7:0.1:7);Z=Y.2./8;h=mesh(Z);rotate(h,1 0 1,180) %旋转处理%axis(-8,8,-8,8,-2,6);axis(square)xlabel(x轴);ylabel(y轴);zlabel(z轴);shading flat;title(抛物柱面)或直接用：ezsurfc(Y.2./8) 九、环面ezmesh(5+2*co