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

《课程设计（论文）_基于matlab的fft算法设计.doc》由会员分享，可在线阅读，更多相关《课程设计（论文）_基于matlab的fft算法设计.doc（22页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、目 录1 引言12课程设计内容及要求2课程设计内容2课程设计要求22.3课程设计目的2课程设计平台23 基于MATLAB的FFT算法设计原理3总体设计流程图3语音信号的采集3语音信号的时频分析3快速傅里叶变换63 fft的运算规律8基于MATLAB的FFT所编写程序的框图12自编算法与机带算法仿真波形比拟134设计总结16参考文献17附录181 引言随着信息时代，数字时代的到来，数字信号处理已经成为一门极其重要的学科和技术领域。以DSP为核心芯片的处理系统日益变成了数字信号处理系统的主流。它广泛用于电子信息、通信、图像处理、语音处理、生物医学、自动控制、地质探测等领域，受到工程设计和使用人员的

2、青睐。MATLAB，它是美国Math Works公司推出的一种面向工程和科学计算的交互式计算软件。它以矩阵运算为根底，把计算、可视化、程序设计融合在一个简单易用的交互式工作环境中，是一款数据分析和处理功能都非常强大的工程适用软件。通过本次实习我们学会了分析和处理音频信号，首先要对声音信号进行采集，MATLAB的数据采集工具箱提供了一整套命令和函数，通过调用这些函数和命令，可直接控制声卡进行数据采集。Window自带的录音机程序也可驱动声卡来采集语音信号，并能保存为WAV格式文件，供MATLAB相关函数直接读取、写入或播放。MATLAB语言是一种数据分析和处理功能十分强大的计算机应用软件，它可以

3、将声音文件变换位离散的数据文件，然后利用其强大的矩阵运算能力处理数据，如数据滤波、傅立叶变换、时域和频域分析、声音回放以及各种图的呈现等，它的信号处理与分析工具箱位语音信号分析提供了十分丰富的功能函数，利用这些功能函数可以快捷而又方便的完成语音信号的处理和分析以及信号的可视化，是人机交互更加便捷。信号处理是MATLAB重要应用的领域之一。语音信号处理是研究用数字信号处理技术和语音学知识对语音信号进行处理的新兴的学科，是目前开展最为迅速的信息科学研究领域的核心技术之一。通过语音传递信息是人类最重要、最有效、最常用和最方便的交换信息形式。语音信号的处理与滤波的设计主要是用MATLAB作为工具平台，

4、设计中涉及到声音的录制、播放、存储和读取，语音信号的抽样、频谱分析，滤波器的设计及语音信号的滤波，通过数字信号处理课程的理论知识的综合运用。从实践上初步实现对数字信号的处理。 2 课程设计内容及要求录制一段个人自己的语音信号，并对录制的信号进行采样；画出采样后语音信号的时域波形和频谱图；在Matlab环境下编写基2 DIT-FFT算法；利用自己编写的算法对已采集的语音信号进行频谱分析，并画出语音信号的时域与频谱图，并与Matlab数字信号处理工具箱中的fft函数进行比照研究，验证自编算法的正确性。1.完成语音信号的采集，利用windows自带的录音机或其他软件，录制一段语音，时间在1s以内；2

5、.在Matlab中编写程序，实现输入信号的倒序；3.编写程序，实现蝶形运算；4.画出语音信号的频谱图，与Matlab数字信号处理工具箱中的fft函数进行比照研究，并对设计结果进行独立思考和分析；1.学会 MATLAB 的使用，掌握 MATLAB 的程序设计方法。 2.掌握在 Windows 环境下语音信号采集的方法。3.掌握数字信号处理的根本概念、根本理论和根本方法。 4.掌握 MATLAB 设计 FIR 和IIR 数字滤波器的方法。5.学会用 MATLAB 对信号进行分析和处理。MATLAB软件MATLAB是由美国mathworks公司发布的主要面对科学计算、可视化以及交互式程序设计的高科技

6、计算环境。它将数值分析、矩阵计算、科学数据可视化以及非线性动态系统的建模和仿真等诸多强大功能集成在一个易于使用的视窗环境中，为科学研究、工程设计以及必须进行有效数值计算的众多科学领域提供了一种全面的解决方案，并在很大程度上摆脱了传统非交互式程序设计语言如C、Fortran的编辑模式，代表了当今国际科学计算软件的先进水平。3 基于MATLAB的FFT算法设计原理3.1 总体设计流程图 在一个相对较安静的环境下，录下1s左右的wav声音信号，然后对声音进行采样，画出其时域波形和频谱图，其流程图如图1所示：开始输入声音信号对声音信号采样蝶形运算原信号fft运算结束 图1设计流程图在实际工作中，我们可

7、以利用windows自带的录音机录制语音文件。采集到语音信号之后，需要对语音信号进行分析，如语音信号的时域分析、频谱分析、语谱图分析。在MATLAB中，我们可以通过y,fs,bits=wavread(语音信号路径,N1 N2)语句。用于读取语音，采样值放在向量y中，fs表示采样频率(Hz)，bits表示采样位数。N1 N2表示读取从N1点到N2点的值假设只有一个N的点那么表示读取前N点的采样值。向量y那么就代表了一个信号也即一个复杂的“函数表达式也就是说可以像处理一个信号表达式一样处理这个声音信号。利用MATLAB中的“wavread命令来读入采集语音信号，将它赋值给某一向量。再对其进行采样，

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

9、功能说明：把语音信号加载入Matlab 仿真软件平台中然后，画出语音信号的时域波形，再对语音信号进行频谱分析。MATLAB提供了快速傅里叶变换算法FFT计算DFT的函数fft,其调用格式如下：Xk=fft(xn,N)参数xn为被变换的时域序列向量，N是DFT变换区间长度，当N大于xn的长度时，fft函数自动在xn后面补零。，当N小于xn的长度时，fft函数计算xn的前N个元素，忽略其后面的元素。原始信号的时域波形图如图3所示：图 3 原始信号的时域波形图 原始信号的频域特性图如图4所示：图 4 原始信号的频域特性图3.4 快速傅里叶变换快速傅里叶变换FFT是为提高DFT运算速度而采用的一种算法

10、。对一个有限长度序列x(n)的N点的DFT为：Xk=xnWknN (k=0,1,，N-1；n=0,1,，N-1;W=e-j2/N)当N=4时，Xk可展开为：X0= x0W0*4+ x1W0*4 +x2W0*4+ x3W0*4X1= x0W0*4+ x1W1*4 +x2W2*4+ x3W3*4X2= x0W0*4+ x1W2*4 +x2W4*4+ x3W6*4X3= x0W0*4+ x1W3*4 +x2W6*4+ x3W9*4从上式可以看出,要求4点的DFT,需要16次的复数乘法运算,12次复数乘法运算算。由此类推，要求出N点的DFT，需要N2次复数乘法运算，N*(N-1)次复数加法运算。当N值

11、较大时，要完成的复数乘法运算和复数加法运算得次数都非常多，无论是用通用计算机还是用DSP芯片，都需要消耗大量的时间，不适合于对实时处理要求高的场合。为了能实时处理DFT,要想减少DFT的运算量可以有两个途径：第一是降N，N的值减小了，运算量就减少了；第二是利用旋转因子的周期性和对称性，可约性。利用这两个途径实现DFT的快速傅里叶变换FFT，FFT算法根本上可分为时域抽取法和频域抽取法。W=e-j2/N的性质：1周期性 2共轭对称性 3可约性 本程序是用基2的按时间抽取的FFT算法DIT-FFT，设序列x(n)的长度为N，且N满足N=2M,M为正整数。假设N不能满足上述关系，可以将序列x(n)补

12、零实现，那么x(n)的N点DFT为：Xk=xnWknN (k=0,1,，N-1；n=0,1,，N-1;W=e-j2/N)将n分为奇数与偶数两局部。按时间抽取基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的计算。用蝶形运算可表式为：以此类推

13、，还可以把x1(n)和x2(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减少运算量的特点表现的越明显。 FFT的运算规律1原位运算N=2M的FFT共M级运算，每级有N/2蝶形原位计算，当数据输入到存储器以后，每一组蝶形运

14、算后，结果仍然存放在这同一组存储器的同一位置，不需要另辟存储空间，直接最后输出。同一级的蝶形运算每个蝶形运算的输入数据对其他级输入没有影响。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)。.以此类推，直到最后一位二进制数字。例如 ，N=8时如下列图5所示：顺序倒序十进制二进制二进制十进制0000000

15、010011004201001023011110641000011510110156110011371111117图5 码位倒序N=8倒序的流程图如图3所示：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图3 倒序的流程图3蝶形运算的运算规律设序列xn倒序后存放在数组A中，如果蝶形运算两个输入相距B=2(L-1).采用原位运算：蝶形运算可表示为： x(K)=x(K)+x(K+B)*WNp; x(K+B)=x(K)-x(K+B)*WNp.其中P=J*2M-1,(J=0,1,2,3,2(

16、L-1)-1) 旋转因子确实定：第L级共有2(L-1)个旋转因子，以(N=23=8为例)：L=1时，W=W J=0；L=2时，W=W J=0，1；L=3时，W=W J=0，1，2，3；W=W=W P=J*2(M-L)，用来确定第L级旋转因子从输入端开始，共进行M级运算，在进行第L级运算时，依次求出2(L-1)个旋转因子，然后计算每个旋转因子所对应的2(M-L)个蝶形元素。第L级的蝶形运算中，每个蝶形运算的两个输入相距B=2L-1,同一旋转因子对应的蝶形运算相隔2L个，同一旋转因子对应的蝶形运算有2M-L个。注：1：控制第L级顺序运算2：控制不同种的旋转因子3：控制同种旋转因子所对应的蝶形运算蝶

17、形运算的流程图如图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 蝶形运算的流程图基于MATLAB的FFT所编写程序的框图 基于MATLAB的FFT所编写程序的框图如图5所示：编写fft程序，画出语音信号频谱图语音信号采集完成语音信号时域图完成语音信号频率特性图实现输入信号的倒序实现一级中不同种蝶形算运实现一级中相同种蝶形运算与matlab自带的fft比拟图 5 程序框图3.5 自编算法与机带算法仿真

18、波形比拟 我们知道MATLAB软件自带FFT算法，我们可以通过比拟自编算法仿真结果与机带算法仿真结果来检验自编算法的正确性。自编算法与FFT算法幅值比拟图如图6所示：图 6 幅值比拟图自编算法FFT算法频谱分贝计较图如图7所示：图 7 分贝比拟图由以上两图可以看出，经过蝶形运算得出的频谱图和信号直接FFT得出的频谱图一致。该程序严格按程序框图编写，思路清晰、容易理解，程序的运行过程在命令窗中一目了然。通过与FFT函数运算的结果比对，程序编写正确。4 设计总结为期四天的课设很快接近尾声了，基于MATLAB的FFT算法实现设计已按方案如期全部完成，通过这次DSP课程设计，我对课堂上所学到的理论知识

19、的理解加深了许多， 自己动脑、动手设计的能力也得到了较大提高。在这次课程设计的过程中，我对 MATLAB 语言有了更深的认识。现在仔细想想，这次课程设计使得我对 MATLAB 语言的理解与应用能力得到了较大的提升，也让我认识到只要深入学习，提升的空间永远是存在的。在设计的过程中我遇到了一些问题，如：编写源程序中出现了语法错误等。通过查阅书本和以前设计的程序我发现了产生错误的原因并解决了问题完成了设计。经过反思我发现较大一局部错误是因为操作的不熟练造成的，这也让我明白了要保持设计的高效率 必须经常练习。另一方面我也发现了动手实践的重要性。动手实践是理论知识得以灵活运用的必要前提，也是今后走上工作

20、岗位之后能够很好的完成设计工作的技术保证。只有遇到实际问题并根据自己对课堂上获得的专业知识的理解来解决才能真正的提高自己的能力。这次程序设计让我获益匪浅，对MATLAB也有了进一步的认识：MATLAB功能强、使用灵活方便等。MATLAB是在国内外广泛使用的一种数字信号处理系统，相信除了以上优点，还有许多我还未发现，希望能在以后的学习中有更深入的认识。在此，我也要感谢一下董翠英老师对我这次课设的大力支持和指导，如果没有老师的帮助我很难在这么短的时间内完成这次课设。参考文献1范寿康 DSP技术与DSP芯片.北京：电子工业出版社2程佩青.数字信号处理教程.北京：清华大学出版社出版，20013高西全

21、丁玉美等 数字信号处理. 北京：电子工业出版社，20214 余成波，陶红艳。数字信号处理及MATLAB 实现，北京：清华大学出版社，2021 5 曹弋，赵阳。MATLAB 实用教程，北京：电子工业出版社，2007 附录源程序：Clear all;x,fs,bits=wavread(gl.wav,2048); %读取声音x1=reshape(x,1,4096);sound(x1,fs,bits); %播放语音信号y1=fft(x1)N=length(x1); %计数读取信号的点数t=(1:N)/fs; %信号的时域采样点figure(1)plot(t, x1); %画出声音采样后的时域波形tit

22、le(原声音信号的时域波形); %给图形加注标签说明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

23、=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代表了不同的旋转因子 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; % 蝶形运算， 注意必须先进行减法运算，然后进行加法运算，否那么要使用中间变量

24、来传递y(k) y(k)=y(k)+t; % 蝶形运算 end end end %yfigure(2)x2,w1=freqz(x1,1) ; %绘制原始语音信号的频率图plot(w1/pi,20*log10(abs(x2);title(频率特性图)xlabel(归一化频率);ylabel(幅度/DB);grid; figure(3)subplot(2,1,1); %把画图区域划分为2行1列，指定第1 个图plot(abs(y1); %绘制原始语音信号的幅频响应图title(直接运算 语音信号FFT频谱特性); %给图形加注标签说明xlabel(n);ylabel(幅值/A);grid; %添加网格subplot(2,1,2); %把画图区域划分为2行1列，指定第2个图plot(abs(y);title(蝶形运算 语音信号FFT频谱特性); xlabel(n);ylabel(幅值/A);grid; %添加网格figure(4)subplot(2,1,1)if y=0 %判断指数是否为0 plot(20*log10(abs(y1); %画信号频谱的分贝图endxlabel(n);ylabel(振幅/db);title(直接运算 语音信号FFT频谱特性);grid ; %添加网格