DSP课程设计基于MATLAB的FFT算法实现

1、1 引言 12基于MATLAB勺FFT算法实现 21.1 系统总体流程图 21.2 FFT 运算规律及编程思想 31.2.1 语音信号的采集 31.2.2 DIT-FFT 算法的基本原理 31.2.3 DIT-FFT 算法的运算规律及编程思想 53 Matlab 程序实现 104 系统人机对话界面 134.1 GUI 简介 134.2 界面设计 134.3 运行调试 145 心得体会 16参考文献 17附录I 18附录H 211 引言MATLA是矩阵实验室(Matrix Laboratory )的简称，是美国 MathWorks公 司出品的商数学软件，用于算法开发、 数据可视化、 数据分析以及

2、数值计算的高 级技术计算语言和交互式环境，主要包括MATLAB口 Simulink两大部分。MATLAB勺应用范围非常广，包括信号和图像处理、通讯、控制系统设计、 测试和测量、 财务建模和分析以及计算生物学等众多应用领域。 附加的工具箱 (单 独提供的专用MATLAB函数集)扩展了 MATLAB环境，以解决这些应用领域内特 定类型的问题。它以矩阵运算为基础，把计算、可视化、程序设计融合在一个简 单易用的交互式工作环境中， 是一款数据分析和处理功能都非常强大的工程适用 软件。 它可以将声音文件变换为离散的数据文件， 然后利用其强大的矩阵运算能 力处理数据，如数据滤波、傅立叶变换、时域和频域分析、

3、声音回放以及各种图 的呈现等， 它的信号处理与分析工具箱位语音信号分析提供了十分丰富的功能函 数， 利用这些功能函数可以快捷而又方便的完成语音信号的处理和分析以及信号 的可视化。数字信号处理是MATLABt要应用的领域之一。对于有限长序列x(n),若要求其N点的傅里叶变换(DFT需要经过M次复 数乘法运算和N*(N-1)次复数加法运算。随着N的增加，运算量将急剧增加，而 在实际问题中，N往往是较大的，如当N=1024时，完成复数乘法和复数加法的 次数分别为百万以上，无论是用通用计算机还是用 DSPE片，都需要消耗大量的 时间和机器内存，不能满足实时的要求。因此， DFT的这种运算只能进行理论上

4、 的计算，不适合对实时处理要求高的场合。因此，研究作为 DSP 的快速算法的 FFT是相当必要的，快速傅里叶变换(FFT)是为提高DFT运算速度而采用的一 种算法， 快速算法的种类很多，而且目前仍在改进和提高， 它是根据离散傅里叶 变换的奇、偶、虚、实等特性，对离散傅立叶变换的算法进行改进获得的。基于 本学期所学的 DIT-FFT 的运算规律和编程思想以及Matlab 的学习和使用，本课设要求在Matlab 环境下编写基2 DIT-FFT 算法实现对离散信号的快速傅里叶变换，再与Matlab 软件自带的FFT 函数实现对离散信号的傅里叶变换进行比较，如果得到的频谱相同，那么我们编写的程序就是正

5、确的。其中离散信号是通过PC自带的录音机录制一段 wav语音信号，用Matlab采样得到离散序列x1。如果 有能力可以选做系统人机对话界面。 用GUI界面完成人机交互方便使用的。 本课 程设计主要是对数字信号的分析。2基于MATLAB勺FFT算法实现2.1 系统总体流程图本设计要求录制一段个人自己的语音信号，并对录制的信号进行采样；画出采样后语音信号的时域波形和频谱图；在Matlab环境下编写基2 DIT-FFT算法;利用自己编写的算法对已采集的语音信号进行频谱分析，并画出语音信号的时域与频谱图，并与Matlab数字信号处理工具箱中的fft函数进行对比研究，验证 自编算法的正确性。所以得到系统

6、总体流程图如图1所示。语音信号采集编写fft程序，画出信号频谱图完成信号时域图兀 成 信 号 频 率 响 应实 现 输 入 信 号 的 倒 序实现一级中不同种蝶形算运实 现级 中 相 同 种 蝶 形 运Matlab自带的FFT比较图1系统总体流程图2.2 FFT运算规律及编程思想2.2.1 语音信号的采集利用PM自带的录音机，录制一段语音信号，保存格式为 wave的文件，并 将其保存在电脑中。在 MATLAB 中， fn=input(' Enter WAV filename:','s'); x,fs,nb=wavread(fn,n1 n2);用于读取语音，采样值

7、放在向量 x中，fs表示采样频率(Hz) , nb表示采样位数。n1 n2表示读取从n1点到n2点的值(若 只有一个n的点则表示读取前n点的采样值)。sound(x,fs,nb); 用于对声音的 回放。向量x则就代表了一个信号(也即一个复杂的“函数表达式”)也就是说可以像处理一个信号表达式一样处理这个声音信号。采集到语音信号之后，需要对语音信号进行分析，如语音信号的时域分析、 频谱分析、谱图分析。2.2.2 DIT-FFT算法的基本原理快速傅里叶变换(FFT)是为提高DFT运算速度而采用的一种算法。 对一个有限长度序列x(n)的N点的DFT为：N 1X k xn 01 Nx n 一WNnk ,

8、 k 0, 1, 2,., N 1k Wn nk,n 0, 1, 2,., N 1所以，要求N点的DFT,需要2次的复数乘法运算，N*(N-1)次复数乘法运算算。随着N的增加，运算量将急剧增加，而在实际问题中，N往往是较大的，如当N=1024时，完成复数乘法和复数加法的次数分别为百万以上，无论是用通用 计算机还是用DSP芯片，都需要消耗大量的时间，不能满足实时的要求不适合于对实时处理要求高的场合。为了能实时处理DFT,要想减少DFT的运算量可以有两个途径：第一是降N, N的值减小了，运算量就减少了；第二是利用旋转 因子的周期性，对称性和可约性。利用这两个途径实现DFT的快速傅里叶变换(FFT)

9、, FFT算法基本上可分为按时间抽取的 FFT算法(DIT-FFT)和按频率抽 设序列x(n)的长度为N,且N满足N=2,M为正整数。若N不能满足上述关系， 可以将序列x(n)补零实现。按时间抽取基 2-FFT算法的基本思路是将 N点序列 按时间下标的奇偶分为两个 N/2点序列，计算这两个N/2点序列的N/2点DFI 计算量可减小约一半；每一个N/2点序列按照同样的划分原则，可以划分为两个 N/4点序列，最后，将原序列划分为多个 2点序列，将计算量大大降低。取的FFT算法(DIF-FFT)旋转因子的性质：(1)周期性(2)共腕对称性(3)可约性本次课设要求用用基WW;n WWN啊*WNWN)n

10、)* W；n wN k)n* W；(n)*Wkn WCWkn WNk肝2的按时间抽取的FFT算法(DIT-FFT)实现FFT功能,按时间下标的奇偶将N点x(n)分别抽取组成两个N/2点序列，分别记为x1(n) 和x2(n),将x(n)的DFTW化为x1(n)和x2(n)的DFT的计算。r 0,1,x(2r) xi(r)x(2r 1)x2(r)n 0,2,4kx nx n WNnkn 1,3,5N12x 2rW2rkr 0,1N12x 2r 1 Wn2r 1 kr 0,12rkxr Wnr 0,1N 122r 1 kx2r Wnr 0,1利用旋转因子的可约性，即:W；*X kj2_2rk e N

11、N 12 .rkx r Wn r 02.2jNrkTN2 k WnrWn：x20r WN"X 1(k) W；X2(k), 0 k N 1用蝶形运算可表示为如图2所示:(k) 人 一_ Y X (#)=国 + 际X" (k)M+幻=X") -际与图2 DIT-FFT蝶形运算流图符号以此类推，还可以把x1(n)和x2(n)按n值得奇偶分为两个序列，这样就达到了降N得目的，从而减少了运算量FFT对DFT的数学运算量改进:直接采用DFT进行计算，运算量为 M次复数乘法和N*(N-1)次复数乘法。当采用M次FFT时，由N=皆求得M=logN,运算流图有M级蝶形，每一级都 由

12、N/2个蝶形运算构成，这样每一级蝶形运算都需要 N/2次复数乘法和N次复数 加法。M级运算共需要复数乘法次数为 C=N/2*M,复数加法次数为C=N*M当N值较大时，FFT减少运算量的特点表现的越明显。2.2.3 DIT-FFT 算法的运算规律及编程思想为了编写 DIT-FFT 算法的运算程序， 首先要分析其运算规律， 总结编程思想 并绘出程序框图。1. 原位计算对N 2M点的FFT共进行M级运算，每级由N/2个蝶形运算组成。在同一级 中，每个蝶的输入数据只对本蝶有用，且输出节点与输入节点在同一水平线上，这就意味着每算完一个蝶后， 所得数据可立即存入原输入数据所占用的数组元素 ( 存储单元 )

13、 ，这种原位( 址 ) 计算的方法可节省大量内存。2. 蝶形运算实现FFT运算的核心是蝶形运算，找出蝶形运算的规律是编程的基础。蝶形 运算是分级进行的； 每级的蝶形运算可以按旋转因子的指数大小排序进行； 如果 指数大小一样则可从上往下依次蝶算。对 N 2M点的FFT共有M级运算，用L表 示从左到右的运算级数(L=1, 2,M)。第L级共有B 2L1个不同指数的旋转 因子，用R表示这些不同指数旋转因子从上到下的顺序(R=0, 1,，B-1) 0第R个旋转因子的指数P 2MLR,旋转因子指数为P的第一个蝶的第一节点标号k 从R开始，由于本级中旋转因子指数相同的蝶共有2ML个，且这些蝶的相邻间距为

14、2L ， 故旋 转 因子指数为 P 的最后 一个蝶 的 第 一节 点标号 k 为 ：MLLL(21) 2 R N 2 R ，本级中各蝶的第二个节点与第一个节点都相距B 点。应用原位计算，蝶形运算可表示成如下形式：AL(J)= AL-1(J)+ AL-1 (J+B)* WNPAL(J+B)= AL-1(J)- AL-1(J+B)* WNP总结上述运算规律， 可采用如下运算方法进行DIT-FFT 运算。 首先读入数据，根据数据长度确定运算级数 M,运算总点数N 2M ,不足补0处理。然后对读入 数据进行数据倒序操作。数据倒序后从第1级开始逐级进行，共进行 M级运算。在进行第 L 级运算时，先算出该

15、级不同旋转因子的个数B 2L1 (也是该级中各个蝶形运算两输入数据的间距)，再从R=0开始按序计算，直到R=B-1结束。每个 R对应的旋转因子指数P 2MLR，旋转因子指数相同的蝶从上往下依次逐个运算， 各个蝶的第一节点标号k都是从R开始，以2L为步长，到N 2L R(可简取极值N-2）结束。考虑到蝶形运算有两个输出，且都要用到本级的两个输入数据，故第 一个输出计算完毕后，输出数据不能立即存入输入地址，要等到第二个输出计算 调用输入数据完毕后才能覆盖。这样数据倒序后的运算可用三重循环程序实现 整个蝶形运算流程图如图3所示。3. 序列倒序为了保证运算输出的X(k)按顺序排列，要求序列x(n)倒序

16、输入，即在运算 前要先对输入的序列进行位序颠倒。如果总点数为n 2M的x(n)的顺序数是用M位二进制数表示，则倒序数只需将顺序数的二进制位倒置即可，按照这一规律用硬件电路和汇编语言很容易产生倒序数。但用MATLA舒高级语言实现倒序时，直接倒置二进制数位的方法不可取，还须找出产生倒序的十进制规律。 将十进制顺序数用I表示，与之对应的二进制数用 旧表示。十进制倒序数用J表示，与 之对应的二进制数用JB表示。JB是旧的位倒置结果，十进制顺序数I增加1, 相当于 旧最低位加1且逢2向高位进1,即相当于JB最高位加1且逢2向低位 进10JB的变化规律反映到J的变化分二种情况：如果JB的最高位是0(J n

17、/2), 则直接由加1(J J n/2)得到下一个倒序值；如果JB的最高位是1(J n/2),则 要先将最高位变0(J J n/2),再在次高位加1(J J n/4)。但次高位加1时， 同样要判断0、1值，如果是0 (J N/4),则直接加1(J J N/4),否则要先将 次高位变0(J J n/4),再判断下一位。依此类推，直到完成最高位加1,逢2向右进位的运算。利用这一算法可按顺序数I的递增顺序，依次求得与之对应的 倒序数Jo为了节省内存，数据倒序可原址进行，当I = J时不需要交换，当I* J时需要交换数据。另外，为了避免再次调换前面已经调换过的一对数据，只对 I <J的情况进行数

18、据交换即可实现数据倒序操作。图 3中数据倒序的程序流程图 如图4所示。例如，N=8时，序列倒序结果如表1所示。表1码位倒序(N=8)00 0 00 0 0010 0 11 0 0420 1 00 1 0230 1 11 1 0641 0 00 0 1151 0 11 0 1561 1 00 1 1371 1 11 1 17图4倒序流程图3 Matlab程序实现用ringin.wav作为例子，运行调试程序。1.程序运行开始时，要求输入采样频率fs (fs=1024),所需要变换的起点N1 (N1=100O和终点N2 (N2=509§以及要采样白语音文件(ringin.wav );其 中

19、N1和N2的差值必须在语音信号长度范围内， 但不能太小，否则听不到较为清 晰的语音。具输入窗口显示如图 5所示。图5程序运行开始的输入窗口将fs=1024, N1=100Q N2=5095和语音信号ringin.wav 输入后得到采样后 的语音信号x1的时域波形，如图6所示，其频率响应特性如图7所示，Matlab 自带的FFT函数实现的x1的频谱与编写的FFT程序实现的x1的频谱的比较如图 8所示。语音信号时域波形0.80.60.40.20-0.2-0.4-0.6-0.8 匚1LJ050010001500200025003000350040004500图6采样后的语音信号的时域波形60x1的频

20、率响应特性50403020100.511.522.5频率3.5图7 x1的频率响应特性Matlab自带的FFT函数实现的x1的频谱200150 -k 100 一-Y|50 -I ,-JhLI11 JH 10仁；. N'一 2、:-中, 一050010001500200025003000350040004500KK图8 FFT函数实现与编写的 FFT程序实现的比较取不同的点数进行FFT变换，经观察，编写FFT程序得到的语音信号的频谱 图与Matlab中自带的FFT函数得到的语音信号频谱图总是基本一致，但是如果 输入的N1和N2差值加1不是2的整数次幕就会有细微差别，这是因为编写的快 速傅

21、里叶变换计算信号的 N点傅里叶变换要求N为2的整数次幕，不够的话信 号将会被补零后运算，即参与运算的信号已经不同（差别很小），所以FFT运算后的结果也不尽相同，所示频谱图自然有细微的差别。由用MATLAB自带FFT函数实现的频谱图与用 MATLA踹写的FFT程序实现 的频谱图相比较，可知，两个算法计算后的结果几乎相同，验证了自编算法的正 确性。如果改变输入的值和语音信号，那么得到不一样的时域波形，频率响应，通 过FFT得到频谱都会不同。需要注意的是自由输入的语音信号不同，其长度也不同，所以根据要输入的语音信号输入 N1, N2的值。4系统人机对话界面4.1 GUI简介图形用户界面(GUI),是

22、一种提供人机交互的工具和方法。GUI是包含图形对象，如窗口、图标、菜单和文本等图文并茂的用户界面。4.2 界面设计用MATLA图形用户界面开发环境设计 GUI点的一般步骤是：1 .进行界面设计。2 .设计控件属性。3 .进行M语言编程。以本设计要求为例介绍。第一步，该选择本图形用户界面需要的控件：三个文本框(Edit Text) 用来输入N1, N2和fs。三个静态文本(static text)用来提示我们需要人为输入的内容。两个推按钮(Push button),用来运行和退出。四个轴对象(axes)用来显示信号时域波形，频谱图，和两种不同FFT实现的频谱 图。完成人机设计界面如图9所示：图9

23、完成人机设计界面第二步，设置控件属性：双击组件可以设置文本框，推按钮的属性，如显示大小，名称和默认值等。第三步，编写回调函数。组件事件的发生是通过回调函数进行工作的。控件设置完成后保存，然后运行GUI （操彳乍为ctrl+T ）,就会进入editor窗口，加入各个控件功能的函数代码。 完成后保存即可。第四步，运行GUI。运行editor窗口的程序后，会弹出已经激活的人机对话界面。系统人机对话界面如图10所示。图10系统人机对话界面4.3运行调试运行GUI,已经弹出图10的系统人机对话界面，根据提示“Input N1:"，“Input N2:"和“Input fs :&quo

24、t;分别输入N1, N2和fs ,然后点击“运行”推按钮就会得 到本设计的要求。运行结果如图11所示。如果运行结束，可以通过点击“退出” 推按钮退出该人机对话界面，返回 Matlab。匚庆图11运行结果5 心得体会本次实习的主要内容是通过用 Matlab 实现FFT 的设计，可以实现对一段自己录制的语音信号进行分析，并画出采样信号的时域与频域图。把自己编写的FFT 算法与 Matlab 自带 FFT 算法进行比较。程序运行调试时，自己选择输入要采样的语音信号， 采样频率以及要变换的范围， 可以实现对不同信号的信号采样和进行不同点的FFT运算。在之前数字信号处理的学习以及完成实验的过程中，已经使

25、用过Matlab ，对其有了一些基础的了解和认识， 通过这次的课程设计使我进一步了解了信号的产生， 采样及频谱分析的方法， 以及其中产生信号和绘制信号的基本命令和一些基础编程语言。 让我感受到只有在了解课本知识的前提下， 才能更好的应用这个工具，并且熟练的应用 Matlab 也可以很好的加深我对课程的理解，方便我的思维。这次课程设计使我了解了 Matlab 的使用方法，提高了自己的分析和动手实践能力。同时我相信，进一步加强对MATLAB勺学习与研究对我今后的学习将会起到很大的帮助。这次的课程设计是对本学期所学知识的一次重要巩固， 使得在课堂上掌握的知识得到了真正的运用。 在学习的过程中和同学讨

26、论， 更明白了理论知识与实践的联系。 书到用时方恨少， 有些知识学会是一回事， 掌握是一回事， 但应用起来，确实不是那么简单的，需要很多知识的融会贯通。程序运行调试初期， 曾经多次出现错误、不能产生图形等问题，但在我翻阅资料认真改正及老师同学的帮助下基本功能还是完成了， 经过 1 个星期的上机实习，程序已得到一些完善，能完成基本的要求的功能。最后经过努力，又深入学习了图形用户界面（GUI）,完成了选做要求的人机对话界面。学习就是一个了解， 疑惑，进而解惑的过程，这次实习就是提供了这样一个发现自己知识漏洞，与同学老师探讨进行解惑的的机会。通过这次课程设计实习， 我更深刻的了解了 Matlab 的

27、运用， 重新复习了 FFT中的重要的序列倒序和蝶形变换的程序， 对课本上的知识有了更深的理解， 使我对数字信号处理有了系统的认知。在这里特别感谢董老师和李老师， 他们给了我们很大的发挥空间， 让我们真正自己动手真正掌握了知识，感谢他们细心指导。也非常感谢我的同学， 他们解开了我在实习中出现的诸多知识死角，谢谢大家！参考文献1范寿康DSP技术与DSP芯片.北京：电子工业出版社2 程佩青 . 数字信号处理教程. 北京：清华大学出版社出版， 20013 高西全 丁玉美等 . 数字信号处理. 北京：电子工业出版社， 20094李勇徐震.MATLAB辅助现代工程数字信号处理.西安电子科技大学出版社5 陈

28、杰 .Matlab 宝典 . 电子工业出版社6苏金明张莲花刘波.MATLAB工具箱应用，电子工业出版社附录Ifs=input(' 输入采样频率fs='); % 语音信号采样频率为 fsN1=input(' 输入所需变换的起点N1=');N2=input(' 输入所需变换的终点N2=');fn=input(' Enter WAV filename:','s'); %获取一个 *.wav 的文件x,fs,nb=wavread(fn,N1 N2); % 读取语音信号的数据sound(x,fs,nb); % 播放语音信号%

29、n=N2-N1+1;脸语音信号文件较大时用这两条%x1=reshape(x,1,2*n);% 语句替换 x1=x'x1=x'y1=fft(x1);figure(1)plot(x1) % 做原始语音信号的时域图形title(' 语音信号时域波形 ')xlabel('n');ylabel(' 幅值 ');M=nextpow2(x1); % 求x的长度对应的2的最低幕次mN=2AM;if length(x1)<Nx1=x1,zeros(1,N-length(x1);%若x的长度不是2的幕，补零到2的整数幂end%数据倒序操作J=0

30、;%给倒序数赋初值for I=0:N-1;% 按序交换数据和算倒序数if I<J;% 条件判断及数据交换T=x1(I+1);x1(I+1)=x1(J+1);x1(J+1)=T;end% 算下一个倒序数K=N/2;while J>=K;J=J-K;K=K/2;endJ=J+K;end %x1;y=x1; %将x倒序排列作为y的初始值WN=exp(-i*2*pi/N);for L=1:MB=2AL/2;% 第L级中，每个蝶形的两个输入数据相距B个点，每级有Bt不同的旋转因子p=J*2A(M-L);WNp=WNAp;for k=J+1:2AL:N %kp=k+B; %t=y(kp)*WN

31、p; % y(kp)=y(k)-t; %for J=0:B-1 % J 代表了不同的旋转因子本次蝶形运算的跨越间隔为 2AL蝶形运算的两个因子对应单元下标的关系蝶形运算的乘积项蝶形运算， 注意必须先进行减法运算，然后进行加法运算，否则要使用中间变量来传递y(k)y(k)=y(k)+t; %蝶形运算endendend%yfigure(2)x1,w1=freqz(x1,1); % 绘制原始语音信号的频率图plot(w1,abs(x1);title('x1 的频率响应特性')xlabel(' 频率');ylabel(' 幅度');figure(3)su

32、bplot(2,1,1);plot(abs(y1) %Matlab 自带的FFT®数实现的语音信号的FF砌谱图title('Matlab 自带的FFT®数实现的x1的频谱')xlabel('K');ylabel('Y1(k)');subplot(2,1,2);plot(abs(y); % 编写的FFTS序实现的语音信号的FFT1谱图 title(' 编写的FFTS序实现的的x1的频谱') xlabel('K');附录nfunction varargout = gaoli208(varargin)

33、% GAOLI208 M-file for gaoli208.fig% GAOLI208, by itself, creates a new GAOLI208or raises the existing% singleton*.% H = GAOLI208 returns the handle to a new GAOLI208 or the handleto% the existing singleton*.% GAOLI208('CALLBACK',hObject,eventData,handles,.) calls the localfunction named CALL

34、BACKin GAOLI208.M with the given inputarguments.%GAOLI208('Property','Value',.)creates a newGAOLI208or raisesthe% existing singleton*. Starting from the left, property value pairs are% applied to the GUI before gaoli208_OpeningFunction gets called.An% unrecognized property name or in

35、valid value makes property application% stop. All inputs are passed to gaoli208_OpeningFcn via varargin.% *See GUI Options on GUIDE's Tools menu. Choose "GUI allows only one% instance to run (singleton)".% See also: GUIDE, GUIDATA, GUIHANDLES % Edit the above text to modify the respons

36、e to help gaoli208% Last Modified by GUIDE v2.5 11-Jan-2010 14:14:19 % Begin initialization code - DO NOT EDITgui_Singleton = 1;gui_State = struct('gui_Name', mfilename, .'gui_Singleton', gui_Singleton, .'gui_OpeningFcn', gaoli208_OpeningFcn, .'gui_OutputFcn', gaoli20

37、8_OutputFcn, .'gui_LayoutFcn', , .'gui_Callback', );if nargin && ischar(varargin1)gui_State.gui_Callback = str2func(varargin1);endif nargoutvarargout1:nargout = gui_mainfcn(gui_State, varargin:);elsegui_mainfcn(gui_State, varargin:);end% End initialization code - DO NOT EDIT%

38、 - Executes just before gaoli208 is made visible.function gaoli208_OpeningFcn(hObject, eventdata, handles, varargin)% This function has no output args, see OutputFcn.% hObject handle to figure% eventdata reserved - to be defined in a future version of MATLAB% handles structure with handles and user

39、data (see GUIDATA)% varargin command line arguments to gaoli208 (see VARARGIN)% Choose default command line output for gaoli208handles.output = hObject;% Update handles structureguidata(hObject, handles);% UIWAIT makes gaoli208 wait for user response (see UIRESUME) % uiwait(handles.figure1);% - Outp

40、uts from this function are returned to the command line. function varargout = gaoli208_OutputFcn(hObject, eventdata, handles) % varargout cell array for returning output args (see VARARGOUT);% hObject handle to figure% eventdata reserved - to be defined in a future version of MATLAB% handles structu

41、re with handles and user data (see GUIDATA)% Get default command line output from handles structure varargout1 = handles.output;function N1_Callback(hObject, eventdata, handles)% hObject handle to N1 (see GCBO)% eventdata reserved - to be defined in a future version of MATLAB% handles structure with

42、 handles and user data (see GUIDATA)% Hints: get(hObject,'String') returns contents of N1 as text% str2double(get(hObject,'String') returns contents of N1 as a double% - Executes during object creation, after setting all properties.function N1_CreateFcn(hObject, eventdata, handles)%

43、hObject handle to N1 (see GCBO)% eventdata reserved - to be defined in a future version of MATLAB%handles empty - handles not created until after all CreateFcns called% Hint: edit controls usually have a white background on Windows.% See ISPC and COMPUTER.ifispc&& isequal(get(hObject,'Ba

44、ckgroundColor'),get(0,'defaultUicontrolBackgroundColor')set(hObject,'BackgroundColor','white');end function N2_Callback(hObject, eventdata, handles)% hObject handle to N2 (see GCBO)% eventdata reserved - to be defined in a future version of MATLAB% handles structure with

45、handles and user data (see GUIDATA)% Hints: get(hObject,'String') returns contents of N2 as text% str2double(get(hObject,'String') returns contents of N2 as a double% - Executes during object creation, after setting all properties.function N2_CreateFcn(hObject, eventdata, handles)% h

46、Object handle to N2 (see GCBO)% eventdata reserved - to be defined in a future version of MATLAB%handles empty - handles not created until after all CreateFcns called% Hint: edit controls usually have a white background on Windows.% See ISPC and COMPUTER.ifispc&& isequal(get(hObject,'Bac

47、kgroundColor'),get(0,'defaultUicontrolBackgroundColor')set(hObject,'BackgroundColor','white');end function fs_Callback(hObject, eventdata, handles)% hObject handle to fs (see GCBO)% eventdata reserved - to be defined in a future version of MATLAB% handles structure with h

48、andles and user data (see GUIDATA)% Hints: get(hObject,'String') returns contents of fs as text% str2double(get(hObject,'String') returns contents of fs as a double% - Executes during object creation, after setting all properties. function fs_CreateFcn(hObject, eventdata, handles)% h

49、Object handle to fs (see GCBO)% eventdata reserved - to be defined in a future version of MATLAB%handles empty - handles not created until after all CreateFcns called% Hint: edit controls usually have a white background on Windows.% See ISPC and COMPUTER.ifispc&& isequal(get(hObject,'Bac

50、kgroundColor'),get(0,'defaultUicontrolBackgroundColor') set(hObject,'BackgroundColor','white');end% - Executes on button press in begin.function begin_Callback(hObject, eventdata, handles)% hObject handle to begin (see GCBO)% eventdata reserved - to be defined in a future

51、 version of MATLAB% handles structure with handles and user data (see GUIDATA)N1=str2double(get(handles.N1,'String');N2=str2double(get(handles.N2,'String');fs=str2double(get(handles.fs,'String');x,fs,bits=wavread('D:Matlabworkringin.wav',N1 N2);sound(x,fs,bits); % 播放语

52、音信号%n=N2-N1+1;脸语音信号文件较大时用这两条%x1=reshape(x,1,2*n);% 语句替换 x1=x'x1=x'y1=fft(x1);axes(handles.tu1);plot(x1) % 做原始语音信号的时域图形title(' 语音信号时域波形 ')xlabel('n');ylabel(' 幅值 ');M=nextpow2(x1); % 求x的长度对应的2的最低幕次mN=2AM;if length(x1)<Nx1=x1,zeros(1,N-length(x1);%若x的长度不是2的幕，补零到2的整数幂

53、end%数据倒序操作J=0;%给倒序数赋初值for I=0:N-1;% 按序交换数据和算倒序数if I<J;% 条件判断及数据交换T=x1(I+1);x1(I+1)=x1(J+1);x1(J+1)=T;end% 算下一个倒序数K=N/2;while J>=K;J=J-K;K=K/2;endJ=J+K;end%x1;y=x1; %将x倒序排列作为y的初始值WN=exp(-i*2*pi/N);for L=1:MB=2AL/2;% 第L级中，每个蝶形的两个输入数据相距B个点，每级有Bt不同的旋转因子for J=0:B-1 % J 代表了不同的旋转因子WNp=WNAp;for k=J+1:

54、2AL:N %kp=k+B; %t=y(kp)*WNp; % y(kp)=y(k)-t; %p=J*2A(M-L);本次蝶形运算的跨越间隔为 2AL蝶形运算的两个因子对应单元下标的关系蝶形运算的乘积项蝶形运算， 注意必须先进行减法运算，然后进行加法运算，否则要使用中间变量来传递y(k)y(k)=y(k)+t; %蝶形运算endendend%yaxes(handles.tu2);x1,w1=freqz(x1,1); % 绘制原始语音信号的频率图plot(w1,abs(x1);title('x1 的频率响应特性')xlabel(' 频率');ylabel(' 幅度');axes(handles.tu3);plot(abs(y1) %Matlab自带的FFT®数实现的语音信号的FF砌谱图title('Matlab自带的FFT®数实现的x1的频谱')xlabel('K');ylabel('Y1(k)');axes(handles.tu4);plot(abs(y); % 编写的FFTS序实现的语音信号的FF砌谱