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

1 引言 (1)
2 基于MATLAB的FFT算法实现 (2)
2.1系统总体流程图 (2)
2.2 FFT运算规律及编程思想 (3)
2.2.1语音信号的采集 (3)
2.2.2 DIT-FFT算法的基本原理 (3)
2.2.3 DIT-FFT算法的运算规律及编程思想 (5)
3 Matlab程序实现 (10)
4 系统人机对话界面 (13)
4.1 GUI简介 (13)
4.2 界面设计 (13)
4.3 运行调试 (14)
5 心得体会 (16)
参考文献 (17)
附录Ⅰ (18)
附录Ⅱ (21)
MATLAB是矩阵实验室（Matrix Laboratory）的简称，是美国MathWorks公司出品的商数学软件，用于算法开发、数据可视化、数据分析以及数值计算的高级技术计算语言和交互式环境，主要包括MATLAB和Simulink两大部分。

MATLAB 的应用范围非常广，包括信号和图像处理、通讯、控制系统设计、测试和测量、财务建模和分析以及计算生物学等众多应用领域。

附加的工具箱（单独提供的专用 MATLAB 函数集）扩展了 MATLAB 环境，以解决这些应用领域内特定类型的问题。

它以矩阵运算为基础，把计算、可视化、程序设计融合在一个简单易用的交互式工作环境中，是一款数据分析和处理功能都非常强大的工程适用软件。

它可以将声音文件变换为离散的数据文件，然后利用其强大的矩阵运算能力处理数据，如数据滤波、傅立叶变换、时域和频域分析、声音回放以及各种图的呈现等，它的信号处理与分析工具箱位语音信号分析提供了十分丰富的功能函数，利用这些功能函数可以快捷而又方便的完成语音信号的处理和分析以及信号的可视化。

数字信号处理是MATLAB重要应用的领域之一。

对于有限长序列x(n),若要求其N点的傅里叶变换（DFT）需要经过2N次复数乘法运算和N*(N-1)次复数加法运算。

随着N的增加，运算量将急剧增加，而在实际问题中，N往往是较大的，如当N=1024时，完成复数乘法和复数加法的次数分别为百万以上，无论是用通用计算机还是用DSP芯片，都需要消耗大量的时间和机器内存，不能满足实时的要求。

因此，DFT的这种运算只能进行理论上的计算，不适合对实时处理要求高的场合。

因此，研究作为DSP的快速算法的FFT是相当必要的，快速傅里叶变换（FFT）是为提高DFT运算速度而采用的一种算法，快速算法的种类很多，而且目前仍在改进和提高，它是根据离散傅里叶变换的奇、偶、虚、实等特性，对离散傅立叶变换的算法进行改进获得的。

基于本学期所学的DIT-FFT的运算规律和编程思想以及Matlab的学习和使用，本课设要求在Matlab环境下编写基2 DIT-FFT算法实现对离散信号的快速傅里叶变换，再与Matlab软件自带的FFT函数实现对离散信号的傅里叶变换进行比较，如果得到的频谱相同，那么我们编写的程序就是正确的。

其中离散信号是通过PC自带的录音机录制一段wav语音信号，用Matlab采样得到离散序列x1。

如果有能力可以选做系统人机对话界面。

用GUI界面完成人机交互方便使用的。

本课程设计主要是对数字信号的分析。

2 基于MATLAB的FFT算法实现
2.1系统总体流程图
本设计要求录制一段个人自己的语音信号，并对录制的信号进行采样；画出采样后语音信号的时域波形和频谱图；在Matlab环境下编写基2 DIT-FFT算法；利用自己编写的算法对已采集的语音信号进行频谱分析，并画出语音信号的时域与频谱图，并与Matlab数字信号处理工具箱中的fft函数进行对比研究，验证自编算法的正确性。

所以得到系统总体流程图如图1所示。

图1 系统总体流程图
2.2 FFT 运算规律及编程思想
2.2.1语音信号的采集
利用PC 机自带的录音机，录制一段语音信号，保存格式为wave 的文件，并将其保存在电脑中。

在
MATLAB
中
，
fn=input('
Enter
WAV
filename:','s');
[x,fs,nb]=wavread(fn,[n1 n2]); 用于读取语音，采样值放在向量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 点的DFT,需要N 2次的复数乘法运算,N*(N-1)次复数乘法运算算。

随着N 的增加，运算量将急剧增加，而在实际问题中，N 往往是较大的，如当N=1024时，完成复数乘法和复数加法的次数分别为百万以上，无论是用通用计算机还是用DSP 芯片，都需要消耗大量的时间，不能满足实时的要求,，不适合于对实时处理要求高的场合。

为了能实时处理DFT,要想减少DFT 的运算量可以有两个途径：第一是降N ，N 的值减小了，运算量就减少了；第二是利用旋转因子的周期性，对称性和可约性。

利用这两个途径实现DFT 的快速傅里叶变换（FFT ），FFT 算法基本上可分为按时间抽取的FFT 算法（DIT-FFT ）和按频率抽取的FFT 算法（DIF-FFT ）。

旋转因子的性质： （1）周期性 （2）共轭对称性 （3）可约性
本次课设要求用用基2的按时间抽取的FFT 算法（DIT-FFT ）实现FFT 功能，
)
()(N n k N n N k N kn N W W W ++==*)(*)(][][n k N n k N kn N W W W --==*)(*)(][][n k N n k N kn N W W W --==m
kn m
N kn N mkn mN kn N W W W W //,
==
设序列x(n)的长度为N ，且N 满足N=2M ,M 为正整数。

若N 不能满足上述关系，可以将序列x(n)补零实现。

按时间抽取基2-FFT 算法的基本思路是将N 点序列按时间下标的奇偶分为两个N/2点序列，计算这两个N/2点序列的N/2点DFT ，计算量可减小约一半；每一个N/2点序列按照同样的划分原则，可以划分为两个N/4点序列，最后，将原序列划分为多个2点序列，将计算量大大降低。

按时间下标的奇偶将N 点x(n)分别抽取组成两个N/2点序列，分别记为x1(n)和x2(n)，将x(n)的DFT 转化为x1(n)和x2(n)的DFT 的计算。

利用旋转因子的可约性，即：
用蝶形运算可表示为如图2所示：
以此类推，还可以把x1(n)和x2(n)按n 值得奇偶分为两个序列，这样就达到了降N 得目的，从而减少了运算量。

FFT 对DFT 的数学运算量改进：
()()()()()()(
)()()(
)1N 0
2
1
N
N
0,2,4...1,3,5 (112)
2
212N
N
0,10,1112
2
2121N 2N
0,10,1
221N nk
n N N nk nk n n N N r k rk r r N N r k
rk
r r X k x n W x n W
x n W
x r W x r W x r W x r W -=--==--+==--+====
+
=
+
+=
+
∑∑∑∑∑∑
∑
2j 2j 2222
e
e rk N
rk
rk rk N
N N W W π
π--===()()()1122
122
2
12,01N N rk k rk
N N N r r k N X k x r W W x r W X k W X k N --===+=+≤≤-∑∑（）（k）图2 DIT-FFT 蝶形运算流图符号
12
,
,1,0,)()12()()2(21-=⎭
⎬⎫=+=N
r r x r x r x r x
直接采用DFT 进行计算，运算量为N 2次复数乘法和N*(N-1)次复数乘法。

当采用M 次FFT 时，由N=2M 求得M=logN ，运算流图有M 级蝶形，每一级都由N/2个蝶形运算构成，这样每一级蝶形运算都需要N/2次复数乘法和N 次复数加法。

M 级运算共需要复数乘法次数为C=N/2*M,复数加法次数为C=N*M 。

当N 值较大时，FFT 减少运算量的特点表现的越明显。

2.2.3 DIT-FFT 算法的运算规律及编程思想
为了编写DIT-FFT 算法的运算程序，首先要分析其运算规律，总结编程思想并绘出程序框图。

1. 原位计算
对M
N 2=点的FFT 共进行M 级运算，每级由N/2个蝶形运算组成。

在同一级中，每个蝶的输入数据只对本蝶有用，且输出节点与输入节点在同一水平线上，这就意味着每算完一个蝶后，所得数据可立即存入原输入数据所占用的数组元素(存储单元)，这种原位(址)计算的方法可节省大量内存。

2. 蝶形运算
实现FFT 运算的核心是蝶形运算，找出蝶形运算的规律是编程的基础。

蝶形运算是分级进行的；每级的蝶形运算可以按旋转因子的指数大小排序进行；如果指数大小一样则可从上往下依次蝶算。

对M
N 2=点的FFT 共有M 级运算，用L 表示从左到右的运算级数(L=1，2，…，M )。

第L 级共有12-=L B 个不同指数的旋转因子，用R 表示这些不同指数旋转因子从上到下的顺序(R=0，1，…，B-1)。

第R 个旋转因子的指数R P L M -=2，旋转因子指数为P 的第一个蝶的第一节点标号k 从R 开始，由于本级中旋转因子指数相同的蝶共有L M -2个，且这些蝶的相邻间距为L 2，故旋转因子指数为P 的最后一个蝶的第一节点标号k 为：
R N R L L L M +-=+⋅--22)12(，本级中各蝶的第二个节点与第一个节点都相距
B 点。

应用原位计算，蝶形运算可表示成如下形式：
L A (J)= 1-L A (J)+ 1-L A (J+B)* P
N W
L A (J+B)= 1-L A (J)-1-L A (J+B)* P N W
总结上述运算规律，可采用如下运算方法进行DIT-FFT 运算。

首先读入数据，根据数据长度确定运算级数M ，运算总点数M
N 2=，不足补0处理。

然后对读入数据进行数据倒序操作。

数据倒序后从第1级开始逐级进行，共进行M 级运算。

在进行第L 级运算时，先算出该级不同旋转因子的个数12-=L B (也是该级中各个蝶形运算两输入数据的间距)，再从R=0开始按序计算，直到R=B-1结束。

每个R 对应的旋转因子指数R P L M -=2，旋转因子指数相同的蝶从上往下依次逐个运算，
各个蝶的第一节点标号k 都是从R 开始，以L 2为步长，到R N L
+-2(可简取极值
N-2)结束。

考虑到蝶形运算有两个输出，且都要用到本级的两个输入数据，故第一个输出计算完毕后，输出数据不能立即存入输入地址，要等到第二个输出计算调用输入数据完毕后才能覆盖。

这样数据倒序后的运算可用三重循环程序实现。

整个蝶形运算流程图如图3所示。

3.序列倒序
为了保证运算输出的X(k)按顺序排列，要求序列x(n)倒序输入，即在运算前要先对输入的序列进行位序颠倒。

如果总点数为M
=的x(n)的顺序数是用M
N2
位二进制数表示，则倒序数只需将顺序数的二进制位倒置即可，按照这一规律用硬件电路和汇编语言很容易产生倒序数。

但用MATLAB等高级语言实现倒序时，直接倒置二进制数位的方法不可取，还须找出产生倒序的十进制规律。

将十进制顺序数用I表示，与之对应的二进制数用IB表示。

十进制倒序数用J表示，与之对应的二进制数用JB表示。

JB是IB的位倒置结果，十进制顺序数I增加1，相当于IB最低位加1且逢2向高位进1，即相当于JB最高位加1且逢2向低位进1。

JB的变化规律反映到J的变化分二种情况：如果JB的最高位是0)2/
J<，
(N
则直接由加1)2/
J≥，则⇐得到下一个倒序值；如果JB的最高位是1)2/
(N
(N
J
J+
要先将最高位变0)2/
J
J+
(N
⇐。

但次高位加1时，
⇐，再在次高位加1)4/
J
(N
J-
同样要判断0、1值，如果是0 )4/
J+
J
⇐，否则要先将
J<，则直接加1)4/
(N
(N
次高位变0)4/
J-
⇐，再判断下一位。

依此类推，直到完成最高位加1，逢2 J
(N
向右进位的运算。

利用这一算法可按顺序数I的递增顺序，依次求得与之对应的倒序数J。

为了节省内存，数据倒序可原址进行，当I = J时不需要交换，当I ≠J时需要交换数据。

另外，为了避免再次调换前面已经调换过的一对数据，只对I<J的情况进行数据交换即可实现数据倒序操作。

图3中数据倒序的程序流程图如图4所示。

例如，N=8时，序列倒序结果如表1所示。

表1 码位倒序（N=8）
自然顺序二进制倒位序二进制倒位顺序
0 0 0 0 0 0 0 0
1 0 0 1 1 0 0 4
2 0 1 0 0 1 0 2
3 0 1 1 1 1 0 6
4 1 0 0 0 0 1 1
5 1 0 1 1 0 1 5
6 1 1 0 0 1 1 3
7 1 1 1 1 1 1 7
图4 倒序流程图
3 Matlab程序实现
用ringin.wav作为例子，运行调试程序。

1.程序运行开始时，要求输入采样频率fs（fs=1024），所需要变换的起点N1（N1=1000）和终点N2（N2=5095）以及要采样的语音文件（ringin.wav）；其中N1和N2的差值必须在语音信号长度范围内，但不能太小，否则听不到较为清晰的语音。

其输入窗口显示如图5所示。

图5 程序运行开始的输入窗口
将fs=1024，N1=1000，N2=5095和语音信号ringin.wav输入后得到采样后的语音信号x1的时域波形，如图6所示，其频率响应特性如图7所示，Matlab 自带的FFT函数实现的x1的频谱与编写的FFT程序实现的x1的频谱的比较如图8所示。

图6 采样后的语音信号的时域波形 050010001500200025003000350040004500-0.8-0.6
-0.4
-0.2
0.20.40.6
0.8语音信号时域波形n 幅值图7 x1的频率响应特性
00.51 1.52 2.53 3.5010203040
50
60x1的频率响应特性频率幅度
取不同的点数进行FFT 变换，经观察，编写FFT 程序得到的语音信号的频谱图与Matlab 中自带的FFT 函数得到的语音信号频谱图总是基本一致，但是如果输入的N1和N2差值加1不是2的整数次幂就会有细微差别，这是因为编写的快速傅里叶变换计算信号的N 点傅里叶变换要求N 为2 的整数次幂，不够的话信号将会被补零后运算，即参与运算的信号已经不同（差别很小），所以FFT 运算后的结果也不尽相同，所示频谱图自然有细微的差别。

由用 MATLAB 自带FFT 函数实现的频谱图与用MATLAB 编写的FFT 程序实现的频谱图相比较，可知，两个算法计算后的结果几乎相同，验证了自编算法的正确性。

如果改变输入的值和语音信号，那么得到不一样的时域波形，频率响应，通过FFT 得到频谱都会不同。

需要注意的是自由输入的语音信号不同，其长度也不同，所以根据要输入的语音信号输入N1，N2的值。

图8 FFT 函数实现与编写的FFT 程序实现的比较
050010001500200025003000350040004500050
100
150
200
Matlab 自带的FFT 函数实现的x1的频谱K Y 1(k )050010001500200025003000350040004500050100
150
200编写的FFT 程序实现的的x1的频谱
K Y (k )
4 系统人机对话界面
4.1 GUI简介
图形用户界面（GUI），是一种提供人机交互的工具和方法。

GUI是包含图形对象，如窗口、图标、菜单和文本等图文并茂的用户界面。

4.2 界面设计
用MATLAB图形用户界面开发环境设计GUI点的一般步骤是：
1.进行界面设计。

2.设计控件属性。

3.进行M语言编程。

以本设计要求为例介绍。

第一步，该选择本图形用户界面需要的控件：
三个文本框(Edit Text)用来输入N1，N2和fs。

三个静态文本(static text)用来提示我们需要人为输入的内容。

两个推按钮(Push button)，用来运行和退出。

四个轴对象(axes)用来显示信号时域波形，频谱图，和两种不同FFT实现的频谱图。

完成人机设计界面如图9所示：
图9 完成人机设计界面
第二步，设置控件属性：
双击组件可以设置文本框，推按钮的属性，如显示大小，名称和默认值等。

第三步，编写回调函数。

组件事件的发生是通过回调函数进行工作的。

控件设置完成后保存，然后运行GUI（操作为ctrl+T），就会进入editor窗口，加入各个控件功能的函数代码。

完成后保存即可。

第四步，运行GUI。

运行editor窗口的程序后，会弹出已经激活的人机对话界面。

系统人机对话界面如图10所示。

图10 系统人机对话界面
4.3 运行调试
“Input 运行GUI，已经弹出图10的系统人机对话界面，根据提示“Input N1：”，
N2：”和“Input fs：”分别输入N1，N2和fs，然后点击“运行”推按钮就会得到本设计的要求。

运行结果如图11所示。

如果运行结束，可以通过点击“退出”推按钮退出该人机对话界面，返回Matlab。

图11 运行结果
5 心得体会
本次实习的主要内容是通过用Matlab实现FFT的设计，可以实现对一段自己录制的语音信号进行分析，并画出采样信号的时域与频域图。

把自己编写的FFT算法与Matlab自带FFT算法进行比较。

程序运行调试时，自己选择输入要采样的语音信号，采样频率以及要变换的范围，可以实现对不同信号的信号采样和进行不同点的FFT运算。

在之前数字信号处理的学习以及完成实验的过程中，已经使用过Matlab，对其有了一些基础的了解和认识，通过这次的课程设计使我进一步了解了信号的产生，采样及频谱分析的方法，以及其中产生信号和绘制信号的基本命令和一些基础编程语言。

让我感受到只有在了解课本知识的前提下，才能更好的应用这个工具，并且熟练的应用Matlab也可以很好的加深我对课程的理解，方便我的思维。

这次课程设计使我了解了Matlab的使用方法，提高了自己的分析和动手实践能力。

同时我相信，进一步加强对MATLAB的学习与研究对我今后的学习将会起到很大的帮助。

这次的课程设计是对本学期所学知识的一次重要巩固，使得在课堂上掌握的知识得到了真正的运用。

在学习的过程中和同学讨论，更明白了理论知识与实践的联系。

书到用时方恨少，有些知识学会是一回事，掌握是一回事，但应用起来，确实不是那么简单的，需要很多知识的融会贯通。

程序运行调试初期，曾经多次出现错误、不能产生图形等问题，但在我翻阅资料认真改正及老师同学的帮助下基本功能还是完成了，经过1个星期的上机实习，程序已得到一些完善，能完成基本的要求的功能。

最后经过努力，又深入学习了图形用户界面（GUI），完成了选做要求的人机对话界面。

学习就是一个了解，疑惑，进而解惑的过程，这次实习就是提供了这样一个发现自己知识漏洞，与同学老师探讨进行解惑的的机会。

通过这次课程设计实习，我更深刻的了解了Matlab的运用，重新复习了FFT 中的重要的序列倒序和蝶形变换的程序，对课本上的知识有了更深的理解，使我对数字信号处理有了系统的认知。

在这里特别感谢董老师和李老师，他们给了我们很大的发挥空间，让我们真正自己动手真正掌握了知识，感谢他们细心指导。

也非常感谢我的同学，他们解开了我在实习中出现的诸多知识死角，谢谢大家！
参考文献
[1]范寿康 DSP 技术与DSP芯片.北京：电子工业出版社
[2]程佩青.数字信号处理教程.北京：清华大学出版社出版，2001
[3]高西全丁玉美等.数字信号处理.北京：电子工业出版社，2009
[4] 李勇徐震.MATLAB辅助现代工程数字信号处理.西安电子科技大学出版社
[5] 陈杰.Matlab宝典.电子工业出版社
[6] 苏金明张莲花刘波.MATLAB工具箱应用，电子工业出版社
附录Ⅰ
fs=input('输入采样频率fs='); %语音信号采样频率为fs
N1=input('输入所需变换的起点N1=');
N2=input('输入所需变换的终点N2=');
fn=input(' Enter WAV filename:','s'); %获取一个*.wav的文件
[x,fs,nb]=wavread(fn,[N1 N2]); %读取语音信号的数据
sound(x,fs,nb); %播放语音信号
%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的最低幂次m
N=2^M;
if length(x1)<N
x1=[x1,zeros(1,N-length(x1))]; % 若x的长度不是2的幂，补零到2的整数幂
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;
end
J=J+K;
end
%x1;
y=x1; % 将x倒序排列作为y的初始值
WN=exp(-i*2*pi/N);
for L=1:M
B=2^L/2;%第L级中，每个蝶形的两个输入数据相距B个点，每级有B个不同的旋转因子
for J=0:B-1 % J代表了不同的旋转因子
p=J*2^(M-L);
WNp=WN^p;
for k=J+1:2^L:N % 本次蝶形运算的跨越间隔为2^L
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,abs(x1));
title('x1的频率响应特性')
xlabel('频率');
ylabel('幅度');
figure(3)
subplot(2,1,1);
plot(abs(y1)) %Matlab自带的FFT函数实现的语音信号的FFT频谱图
title('Matlab自带的FFT函数实现的x1的频谱')
xlabel('K');
ylabel('Y1(k)');
subplot(2,1,2);
plot(abs(y)); %编写的FFT程序实现的语音信号的FFT频谱图title('编写的FFT程序实现的的x1的频谱')
xlabel('K');
附录Ⅱ
function varargout = gaoli208(varargin)
% GAOLI208 M-file for gaoli208.fig
% GAOLI208, by itself, creates a new GAOLI208 or raises the existing % singleton*.
%
% H = GAOLI208 returns the handle to a new GAOLI208 or the handle to
% the existing singleton*.
%
% GAOLI208('CALLBACK',hObject,eventData,handles,...) calls the local
% function named CALLBACK in GAOLI208.M with the given input arguments.
%
% GAOLI208('Property','Value',...) creates a new GAOLI208 or raises the
% existing singleton*. Starting from the left, property value pairs are
% applied to the GUI before gaoli208_OpeningFunction gets called. An
% unrecognized property name or invalid 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 response to help gaoli208
% Last Modified by GUIDE v2.5 11-Jan-2010 14:14:19
% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name', mfilename, ...
'gui_Singleton', gui_Singleton, ...
'gui_OpeningFcn', @gaoli208_OpeningFcn, ...
'gui_OutputFcn', @gaoli208_OutputFcn, ...
'gui_LayoutFcn', [] , ...
'gui_Callback', []);
if nargin && ischar(varargin{1})
gui_State.gui_Callback = str2func(varargin{1});
end
if nargout
[varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:}); else
gui_mainfcn(gui_State, varargin{:});
end
% End initialization code - DO NOT EDIT
% --- 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 data (see GUIDATA)
% varargin command line arguments to gaoli208 (see VARARGIN)
% Choose default command line output for gaoli208
handles.output = hObject;
% Update handles structure
guidata(hObject, handles);
% UIWAIT makes gaoli208 wait for user response (see UIRESUME)
% uiwait(handles.figure1);
% --- Outputs 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 structure with handles and user data (see GUIDATA)
% Get default command line output from handles structure
varargout{1} = 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 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)
% 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.
if ispc && isequal(get(hObject,'BackgroundColor'), 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 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)
% hObject 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.
if ispc && isequal(get(hObject,'BackgroundColor'), 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 handles 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)
% hObject 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.
if ispc && isequal(get(hObject,'BackgroundColor'), 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 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:\Matlab\work\ringin.wav',[N1 N2]);
sound(x,fs,bits); %播放语音信号
%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的最低幂次m
N=2^M;
if length(x1)<N
x1=[x1,zeros(1,N-length(x1))]; % 若x的长度不是2的幂，补零到2的整数幂
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;
end
J=J+K;
end
%x1;
y=x1; % 将x倒序排列作为y的初始值
WN=exp(-i*2*pi/N);
for L=1:M
B=2^L/2;%第L级中，每个蝶形的两个输入数据相距B个点，每级有B个不同的旋转因子
for J=0:B-1 % J代表了不同的旋转因子
p=J*2^(M-L);
WNp=WN^p;
for k=J+1:2^L:N % 本次蝶形运算的跨越间隔为2^L
kp=k+B; % 蝶形运算的两个因子对应单元下标的关系
t=y(kp)*WNp; % 蝶形运算的乘积项
y(kp)=y(k)-t; % 蝶形运算，注意必须先进行减法运算，然后进行加法运算，否则要使用中间变量来传递y(k)
y(k)=y(k)+t; % 蝶形运算
end
end
end
%y
axes(handles.tu2);
[x1,w1]=freqz(x1,1); %绘制原始语音信号的频率图
plot(w1,abs(x1));
title('x1的频率响应特性')
xlabel('频率');
ylabel('幅度');
axes(handles.tu3);
plot(abs(y1)) %Matlab自带的FFT函数实现的语音信号的FFT频谱图
title('Matlab自带的FFT函数实现的x1的频谱')
xlabel('K');
ylabel('Y1(k)');
axes(handles.tu4);
plot(abs(y)); %编写的FFT程序实现的语音信号的FFT频谱图
title('编写的FFT程序实现的的x1的频谱')
xlabel('K');
ylabel('Y(k)');
% --- Executes on button press in quit.
function quit_Callback(hObject, eventdata, handles)
% hObject handle to quit (see GCBO)
% eventdata reserved - to be defined in a future version of MATLAB % handles structure with handles and user data (see GUIDATA) close all;
% --- Executes on selection change in ringin.
function ringin_Callback(hObject, eventdata, handles)
% hObject handle to ringin (see GCBO)
% eventdata reserved - to be defined in a future version of MATLAB % handles structure with handles and user data (see GUIDATA)
% Hints: contents = get(hObject,'String') returns ringin contents as cell array
% contents{get(hObject,'Value')} returns selected item from ringin
% --- Executes during object creation, after setting all properties. function ringin_CreateFcn(hObject, eventdata, handles)
% hObject handle to ringin (see GCBO)
% eventdata reserved - to be defined in a future version of MATLAB % handles empty - handles not created until after all CreateFcns called
% Hint: popupmenu controls usually have a white background on Windows. % See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
set(hObject,'BackgroundColor','white');
end。