课程设计（论文）基于matlab的fft算法设计

1、目目 录录 1 引言.1 2 课程设计内容及要求.2 2.1 课程设计内容 .2 2.2 课程设计要求 .2 2.3 课程设计目的 .2 2.3 课程设计平台.2 3 基于 matlab 的 fft 算法设计原理.3 3.1 总体设计流程图 .3 3.2 语音信号的采集 .3 3.3 语音信号的时频分析 .3 3.4 快速傅里叶变换 .6 3.4.1 fft的运算规律.8 3.4.2 基于 matlab 的 fft 所编写程序的框图 .12 3.5 自编算法与机带算法仿真波形比较 .13 4 设计总结.16 参考文献.17 附录.18 1 引言 随着信息时代，数字时代的到来，数字信号处理已经成

2、为一门极其重要的学 科和技术领域。以 dsp 为核心芯片的处理系统日益变成了数字信号处理系统的主 流。它广泛用于电子信息、通信、图像处理、语音处理、生物医学、自动控制、 地质探测等领域，受到工程设计和使用人员的青睐。 matlab，它是美国 math works 公司推出的一种面向工程和科学计算的交互 式计算软件。它以矩阵运算为基础，把计算、可视化、程序设计融合在一个简单 易用的交互式工作环境中，是一款数据分析和处理功能都非常强大的工程适用软 件。通过本次实习我们学会了分析和处理音频信号，首先要对声音信号进行采集， matlab 的数据采集工具箱提供了一整套命令和函数，通过调用这些函数和命令，

3、 可直接控制声卡进行数据采集。window 自带的录音机程序也可驱动声卡来采集语 音信号，并能保存为 wav 格式文件，供 matlab 相关函数直接读取、写入或播 放。 matlab 语言是一种数据分析和处理功能十分强大的计算机应用软件，它可 以将声音文件变换位离散的数据文件，然后利用其强大的矩阵运算能力处理数据， 如数据滤波、傅立叶变换、时域和频域分析、声音回放以及各种图的呈现等，它 的信号处理与分析工具箱位语音信号分析提供了十分丰富的功能函数，利用这些 功能函数可以快捷而又方便的完成语音信号的处理和分析以及信号的可视化，是 人机交互更加便捷。信号处理是 matlab 重要应用的领域之一。

4、 语音信号处理是研究用数字信号处理技术和语音学知识对语音信号进行处理 的新兴的学科，是目前发展最为迅速的信息科学研究领域的核心技术之一。通过 语音传递信息是人类最重要、最有效、最常用和最方便的交换信息形式。 语音信号的处理与滤波的设计主要是用 matlab 作为工具平台，设计中涉 及到声音的录制、播放、存储和读取，语音信号的抽样、频谱分析，滤波器的设 计及语音信号的滤波，通过数字信号处理课程的理论知识的综合运用。从实践上 初步实现对数字信号的处理。 2 课程设计内容及要求 2.1 课程设计内容 录制一段个人自己的语音信号，并对录制的信号进行采样；画出采样后语音 信号的时域波形和频谱图；在 ma

5、tlab 环境下编写基 2 dit-fft 算法；利用自己编 写的算法对已采集的语音信号进行频谱分析，并画出语音信号的时域与频谱图， 并与 matlab 数字信号处理工具箱中的 fft 函数进行对比研究，验证自编算法的正确 性。 2.2 课程设计要求 1.完成语音信号的采集，利用 windows 自带的录音机或其他软件，录制一段 语音，时间在 1s 以内； 2.在 matlab 中编写程序，实现输入信号的倒序； 3.编写程序，实现蝶形运算； 4.画出语音信号的频谱图，与 matlab 数字信号处理工具箱中的 fft 函数进行对 比研究，并对设计结果进行独立思考和分析； 2.3 课程设计目的 1

6、.学会 matlab 的使用，掌握 matlab 的程序设计方法。 2.掌握在 windows 环境下语音信号采集的方法。 3.掌握数字信号处理的基本概念、基本理论和基本方法。 4.掌握 matlab 设计 fir 和 iir 数字滤波器的方法。 5.学会用 matlab 对信号进行分析和处理。 2.4 课程设计平台 matlab7.1 软件 matlab 是由美国 mathworks 公司发布的主要面对科学计算、可视化以及交 互式程序设计的高科技计算环境。它将数值分析、矩阵计算、科学数据可视化以 及非线性动态系统的建模和仿真等诸多强大功能集成在一个易于使用的视窗环境 中，为科学研究、工程设计

7、以及必须进行有效数值计算的众多科学领域提供了一 种全面的解决方案，并在很大程度上摆脱了传统非交互式程序设计语言（如 c、fortran）的编辑模式，代表了当今国际科学计算软件的先进水平。 3 基于 matlab 的 fft 算法设计原理 3.1 总体设计流程图 在一个相对较安静的环境下，录下 1s 左右的 wav 声音信号，然后对声音进行 采样，画出其时域波形和频谱图，其流程图如图 1 所示： 图 1 设计流程图 3.2 语音信号的采集 在实际工作中，我们可以利用 windows 自带的录音机录制语音文件。采集到 语音信号之后，需要对语音信号进行分析，如语音信号的时域分析、频谱分析、 语谱图分

8、析。在 matlab 中，我们可以通过y,fs,bits=wavread(语音信号路径 ,n1 n2)语句。用于读取语音，采样值放在向量 y 中，fs 表示采样频率(hz)，bits 表示采样位数。n1 n2表示读取从 n1 点到 n2 点的值（若只有一个 n 的点则表 示读取前 n 点的采样值）。向量 y 则就代表了一个信号（也即一个复杂的“函数表 达式”）也就是说可以像处理一个信号表达式一样处理这个声音信号。 3.3 语音信号的时频分析 利用matlab中的“wavread”命令来读入（采集）语音信号，将它赋值给某一 开始 输入声音信号 对声音信号采样 蝶形运算 原信号 fft 运算 结束

9、 向量。再对其进行采样，记住采样频率和采样点数。下面介绍wavread 函数几种 调用格式。 1.y=wavread（file） 功能说明：读取 file 所规定的 wav 文件，返回采样值放在向量 y 中。 2.y,fs,nbits=wavread(file) 功能说明：采样值放在向量 y 中，fs 表示采样频率（hz），nbits 表示采样位 数。 3.y=wavread（file，n） 功能说明：读取钱 n 点的采样值放在向量 y 中。 4.y=wavread（file，n1,n2） 功能说明：读取从 n1 到 n2 点的采样值放在向量 y 中。 接下来，对语音信号 speech off

10、.wav 进行采样。其程序如下： y,fs,nbits=wavered (speech off.wav); 功能说明：把语音信号加载入 matlab 仿真软件平台中 然后，画出语音信号的时域波形，再对语音信号进行频谱分析。matlab 提供了 快速傅里叶变换算法 fft 计算 dft 的函数 fft,其调用格式如下： xk=fft(xn,n) 参数 xn 为被变换的时域序列向量，n 是 dft 变换区间长度，当 n 大于 xn 的长度 时，fft 函数自动在 xn 后面补零。，当 n 小于 xn 的长度时，fft 函数计算 xn 的前 n 个元素，忽略其后面的元素。 原始信号的时域波形图如图

11、3 所示： 图 3 原始信号的时域波形图 原始信号的频域特性图如图 4 所示： 图 4 原始信号的频域特性图 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）可展开为： x（0）= x（0）w0*4+ x（1）w0*4 +x（2）w0*4+ x（3）w0*4 x（1）= x（0）w0*4+ x（1）w1*4 +x（2）w2*4+ x（3）w3*4 x（2）= x（0）w0*4+

12、x（1）w2*4 +x（2）w4*4+ x（3）w6*4 x（3）= x（0）w0*4+ x（1）w3*4 +x（2）w6*4+ x（3）w9*4 从上式可以看出,要求 4 点的 dft,需要 16 次的复数乘法运算,12 次复数乘法运 算算。由此类推，要求出 n 点的 dft，需要 n2 次复数乘法运算，n*(n-1)次复 数加法运算。当 n 值较大时，要完成的复数乘法运算和复数加法运算得次数都非 常多，无论是用通用计算机还是用 dsp 芯片，都需要消耗大量的时间，不适合于 对实时处理要求高的场合。为了能实时处理 dft,要想减少 dft 的运算量可以有两 个途径：第一是降 n，n 的值减小

13、了，运算量就减少了；第二是利用旋转因子的 周期性和对称性，可约性。利用这两个途径实现 dft 的快速傅里叶变换（fft） ， fft 算法基本上可分为时域抽取法和频域抽取法。 w=e-j2/n 的性质： （1）周期性 （2）共轭对称性 （3）可约性 本程序是用基 2 的按时间抽取的 fft 算法（dit-fft） ，设序列 x(n)的长度为 n，且 n 满足 n=2m,m 为正整数。若 n 不能满足上述关系，可以将序列 x(n)补零 实现，则 x(n)的 n 点 dft 为： x（k）=x（n）wknn (k=0,1,，n-1；n=0,1,，n-1;w=e- j2/n)将 n 分为奇数与偶数两

14、部分。 按时间抽取基 2-fft 算法的基本思路是将 n 点序列按时间下标的奇偶分为两 个 n/2 点序列，计算这两个 n/2 点序列的 n/2 点 dft，计算量可减小约一半；每 一个 n/2 点序列按照同样的划分原则，可以划分为两个 n/4 点序列，最后，将原 序列划分为多个 2 点序列，将计算量大大降低。 1. 按时间下标的奇偶将 n 点 x(n)分别抽取组成两个 n/2 点序列，分别记为 x1(n)和 x2(n)，将 x(n)的 dft 转化为 x1(n)和 x2(n)的 dft 的计算。 )()(nnk n nnk n kn n www *)(*)( nk n nk n kn n w

15、ww mkn mn kn n mkn mn kn n wwww / / , 1 2 , 1 , 0, )() 12( )()2( 2 1 n r rxrx rxrx 1 n 0 21 nn 0,2,4.1,3,5. 11 22 212 nn 0,10,1 11 22 212 1n2n 0,10,1 221 n nk n nn nknk nn nn rkrk rr nn rkrk rr x kx n w x n wx n w xr wxrw x r wxr w 2 j 2 j2 2 2 2 ee rk n rk rkrk n nn ww 用蝶形运算可表式为： 以此类推，还可以把 x1(n)和 x

16、2(n)按 n 值得奇偶分为两个序列，这样就达到了 降 n 得目的，从而减少了运算量。 fft 对 dft 的数学运算量改进： 直接采用 dft 进行计算，运算量为 n2 次复数乘法和 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 减少运算量的特点表现的越明显。 3.4.1 fft 的运算规律 （1）原位运算 n=2m的f

17、ft共m级运算，每级有n/2蝶形原位计算，当数据输入到存储器 以后，每一组蝶形运算后，结果仍然存放在这同一组存储器的同一位置，不需要 另辟存储空间，直接最后输出。 同一级的蝶形运算每个蝶形运算的输入数据对其他级输入没有影响。 （2）倒序运算的规律 输入序列先按自然顺序存入存储单元,然后经变址运算来实现倒位序排列,用 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)

18、。 .以此类推，直到最后一位二进制数字。 例如 ，n=8 时如下图 5 所示： 11 22 12 22 00 12 , 01 nn rkkrk nnn rr k n x kx r wwxr w x kw xkn （）（k） 顺序倒序 十进制二进制二进制十进制 00000000 10011004 20100102 30111106 41000011 51011015 61100113 71111117 图 5 码位倒序（n=8） 倒序的流程图如图 3 所示： 图 3 倒序的流程图 （3）蝶形运算的运算规律 设序列x（n）倒序后存放在数组a中，如果蝶形运算两个输入相距b=2(l-1). lh=n/

19、2 n1=n-1 j=0 t=x(j+1) x(j+1)=x(i+1) x(i+1)=t i=0：n1 i=j k=lh jk j=j+k j=j-k k=k/2 y n y n 采用原位运算：蝶形运算可表示为： x(k)=x(k)+x(k+b)*wnp; x(k+b)=x(k)- x(k+b)*wnp.其中p=j*2（m-1）,(j=0,1,2,3,2(l-1)-1) 旋转因子的确定： 第l级共有2(l-1)个旋转因子，以(n=23=8为例)： l=1时，w=w j=0； p n j l2 l=2时，w=w j=0，1； p n j l2 l=3时，w=w j=0，1，2，3； p n j

20、l2 w=w=w p=j*2(m-l)，用来确定第l级旋转因子 p n j l2 )(2*lmj n 从输入端开始，共进行 m 级运算，在进行第 l 级运算时，依次求出 2(l-1)个 旋转因子，然后计算每个旋转因子所对应的 2(m-l)个蝶形元素。第 l 级的蝶形运 算中，每个蝶形运算的两个输入相距 b=2l-1,同一旋转因子对应的蝶形运算相隔 2l 个，同一旋转因子对应的蝶形运算有 2m-l 个。 注：1：控制第 l 级顺序运算 2：控制不同种的旋转因子 3：控制同种旋转因子所对应的蝶形运算 蝶形运算的流程图如图 4 所示： ： 1 2 3 图 4 蝶形运算的流程图 3.4.2 基于 ma

21、tlab 的 fft 所编写程序的框图 基于 matlab 的 fft 所编写程序的框图如图 5 所示： 开始开始开始 j=0：b-1 n=2m b=2l/2 l=1:m 输入 x(n) k=j+1:2l:n p=j*2(m-l) t=x(k)+x(k+b)* wnp x(k+b)=x(k)-x(k+b)* wnp x(k)=t 倒序 结束 输出 图 5 程序框图 3.5 自编算法与机带算法仿真波形比较 我们知道 matlab 软件自带 fft 算法，我们可以通过比较自编算法仿真结果与 机带算法仿真结果来检验自编算法的正确性。 自编算法与 fft 算法幅值比较图如图 6 所示： 语音信号采集

22、完 成 语 音 信 号 时 域 图 完 成 语 音 信 号 频 率 特 性 图 编写 fft 程序，画出语音信号 频谱图 实 现 输 入 信 号 的 倒 序 实 现 一 级 中 不 同 种 蝶 形 算 运 实 现 一 级 中 相 同 种 蝶 形 运 算 与 m atl ab 自 带 的 fft 比 较 图 6 幅值比较图 自编算法 fft 算法频谱分贝计较图如图 7 所示： 图 7 分贝比较图 由以上两图可以看出，经过蝶形运算得出的频谱图和信号直接 fft 得出的频 谱图一致。该程序严格按程序框图编写，思路清晰、容易理解，程序的运行过程 在命令窗中一目了然。通过与 fft 函数运算的结果比对，

23、程序编写正确。 4 设计总结 为期四天的课设很快接近尾声了，基于 matlab 的 fft 算法实现设计已按计划 如期全部完成，通过这次 dsp 课程设计，我对课堂上所学到的理论知识的理解加 深了许多， 自己动脑、动手设计的能力也得到了较大提高。在这次课程设计的过 程中，我对 matlab 语言有了更深的认识。现在仔细想想，这次课程设计使得我 对 matlab 语言的理解与应用能力得到了较大的提升，也让我认识到只要深入学 习，提升的空间永远是存在的。在设计的过程中我遇到了一些问题，如：编写源 程序中出现了语法错误等。通过查阅书本和以前设计的程序我发现了产生错误的 原因并解决了问题完成了设计。经

24、过反思我发现较大一部分错误是因为操作的不 熟练造成的，这也让我明白了要保持设计的高效率 必须经常练习。另一方面我也 发现了动手实践的重要性。动手实践是理论知识得以灵活运用的必要前提，也是 今后走上工作岗位之后能够很好的完成设计工作的技术保证。只有遇到实际问题 并根据自己对课堂上获得的专业知识的理解来解决才能真正的提高自己的能力。 这次程序设计让我获益匪浅，对 matlab 也有了进一步的认识：matlab 功能强、 使用灵活方便等。matlab 是在国内外广泛使用的一种数字信号处理系统，相信除 了以上优点，还有许多我还未发现，希望能在以后的学习中有更深入的认识。 在此，我也要感谢一下董翠英老师

25、对我这次课设的大力支持和指导，如果没 有老师的帮助我很难在这么短的时间内完成这次课设。 参考文献 1范寿康 dsp 技术与 dsp 芯片.北京：电子工业出版社 2程佩青.数字信号处理教程.北京：清华大学出版社出版，2001 3高西全 丁玉美等 数字信号处理. 北京：电子工业出版社，2009 4 余成波，陶红艳。数字信号处理及 matlab 实现，北京：清华大学出版社， 2008 5 曹弋，赵阳。matlab 实用教程，北京：电子工业出版社，2007 附录 源程序： clear all; x,fs,bits=wavread(gl.wav,2048); %读取声音 x1=reshape(x,1,4

26、096); sound(x1,fs,bits); %播放语音信号 y1=fft(x1) n=length(x1); %计数读取信号的点数 t=(1:n)/fs; %信号的时域采样点 figure(1) plot(t, x1); %画出声音采样后的时域波形 title(原声音信号的时域波形); %给图形加注标签说明 xlabel(时间/t); ylabel(振幅/a); grid ; %添加网格 m=nextpow2(x1); % 求 x 的长度对应的 2 的最低幂次 m n=2m; if length(x1)n x1=x1,zeros(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;