fft算法的应用研究大学论文.doc

数字信号处理课程设计报告FFT算法的应用研究专 业： 通信工程 班 级： 通信11级 组 次： 第20组 姓名及学号： 吴杨生 2011013847 姓名及学号： 朱泽章 2011013864 组 员承 担 任 务吴杨生共同完成任务朱泽章共同完成任务指导教师评价意见FFT算法的应用研究一、设计目的MATLAB是一种功能强大、效率高、交互性好的数值和可视化计算机高级语言，它将数值分析、矩阵运算、信号处理和图形显示有机地融合为一体，形成了一个极其方便、用户界面友好的操作环境。经过多年的发展，已经发展成为一种功能全面的软件，几乎可以解决科学计算中所有问题。MATLAB软件还提供了非常广泛和灵活的用于处理数据集的数组运算功能。FFT算法的应用研究很广泛，数字信号也有很多，本次课程设计采取对语音信号进行FFT算法的的应用研究：录制一段个人自己的语音信号，并对录制的信号进行采样；画出采样后语音信号的时域波形和频谱图；在Matlab环境下编写基2 DIT-FFT算法；利用自己编写的算法对已采集的语音信号进行频谱分析，并画出语音信号的时域与频谱图，并与Matlab数字信号处理工具箱中的fft函数进行对比研究，验证自编算法的正确性二、设计任务对语音信号进行FFT算法的的应用研究三、设计原理1系统总体流程图本设计要求录制一段个人自己的语音信号，并对录制的信号进行采样；画出采样后语音信号的时域波形和频谱图；在Matlab环境下编写基2 DIT-FFT算法；利用自己编写的算法对已采集的语音信号进行频谱分析，并画出语音信号的时域与频谱图，并与Matlab数字信号处理工具箱中的fft函数进行对比研究，验证自编算法的正确性。所以得到系统总体流程图如图1所示。语音信号采集完成信号时域图完成信号频率响应编写fft程序，画出信号频谱图实现输入信号的倒序实现一级中不同种蝶形算运实现一级中相同种蝶形运算与Matlab自带的FFT比较图1 系统总体流程图2 FFT运算规律及编程思想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 DIT-FFT算法的基本原理快速傅里叶变换（FFT）是为提高DFT运算速度而采用的一种算法。对一个有限长度序列x(n)的N点的DFT为：所以，要求N点的DFT,需要N2次的复数乘法运算,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功能，设序列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所示：图2 DIT-FFT蝶形运算流图符号以此类推，还可以把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减少运算量的特点表现的越明显。2.3 DIT-FFT算法的运算规律及编程思想为了编写DIT-FFT算法的运算程序，首先要分析其运算规律，总结编程思想并绘出程序框图。1. 原位计算对点的FFT共进行M级运算，每级由N/2个蝶形运算组成。在同一级中，每个蝶的输入数据只对本蝶有用，且输出节点与输入节点在同一水平线上，这就意味着每算完一个蝶后，所得数据可立即存入原输入数据所占用的数组元素(存储单元)，这种原位(址)计算的方法可节省大量内存。2. 蝶形运算实现FFT运算的核心是蝶形运算，找出蝶形运算的规律是编程的基础。蝶形运算是分级进行的；每级的蝶形运算可以按旋转因子的指数大小排序进行；如果指数大小一样则可从上往下依次蝶算。对点的FFT共有M级运算，用L表示从左到右的运算级数(L=1，2，M )。第L级共有个不同指数的旋转因子，用R表示这些不同指数旋转因子从上到下的顺序(R=0，1，B-1)。第R个旋转因子的指数，旋转因子指数为P的第一个蝶的第一节点标号k从R开始，由于本级中旋转因子指数相同的蝶共有个，且这些蝶的相邻间距为，故旋转因子指数为P的最后一个蝶的第一节点标号k为：，本级中各蝶的第二个节点与第一个节点都相距B点。应用原位计算，蝶形运算可表示成如下形式： (J)= (J)+ (J+B)* (J+B)= (J)-(J+B)* 总结上述运算规律，可采用如下运算方法进行DIT-FFT运算。首先读入数据，根据数据长度确定运算级数M，运算总点数，不足补0处理。然后对读入数据进行数据倒序操作。数据倒序后从第1级开始逐级进行，共进行M级运算。在进行第L级运算时，先算出该级不同旋转因子的个数(也是该级中各个蝶形运算两输入数据的间距)，再从R=0开始按序计算，直到R=B-1结束。每个R对应的旋转因子指数，旋转因子指数相同的蝶从上往下依次逐个运算，各个蝶的第一节点标号k都是从R开始，以为步长，到(可简取极值N-2)结束。考虑到蝶形运算有两个输出，且都要用到本级的两个输入数据，故第一个输出计算完毕后，输出数据不能立即存入输入地址，要等到第二个输出计算调用输入数据完毕后才能覆盖。这样数据倒序后的运算可用三重循环程序实现。整个蝶形运算流程图如图3所示。图3 整个蝶形运算流程图送入x(n),MN=2M倒序L=1:MB=2(L-1)J=0:B-1P=J*2(M-L) K=J:2L:N-1T=A(K)+A(K+B)*WNPA(K+B)=A(K)-A(K+B)*WNPA(K)=T输出开始 结束束束3.序列倒序为了保证运算输出的X(k)按顺序排列，要求序列x(n)倒序输入，即在运算前要先对输入的序列进行位序颠倒。如果总点数为的x(n)的顺序数是用M位二进制数表示，则倒序数只需将顺序数的二进制位倒置即可，按照这一规律用硬件电路和汇编语言很容易产生倒序数。但用MATLAB等高级语言实现倒序时，直接倒置二进制数位的方法不可取，还须找出产生倒序的十进制规律。将十进制顺序数用I表示，与之对应的二进制数用IB表示。十进制倒序数用J表示，与之对应的二进制数用JB表示。JB是IB的位倒置结果，十进制顺序数I增加1，相当于IB最低位加1且逢2向高位进1，即相当于JB最高位加1且逢2向低位进1。JB的变化规律反映到J的变化分二种情况：如果JB的最高位是0，则直接由加1得到下一个倒序值；如果JB的最高位是1，则要先将最高位变0，再在次高位加1。但次高位加1时，同样要判断0、1值，如果是0 ，则直接加1，否则要先将次高位变0，再判断下一位。依此类推，直到完成最高位加1，逢2向右进位的运算。利用这一算法可按顺序数I的递增顺序，依次求得与之对应的倒序数J。为了节省内存，数据倒序可原址进行，当I = J时不需要交换，当I J时需要交换数据。另外，为了避免再次调换前面已经调换过的一对数据，只对IJ的情况进行数据交换即可实现数据倒序操作。图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 倒序流程图LH=N/2 J=1 N1=N-1I=0，N-1数据交换算倒序数YNIJ?T=A(J)A(J)=A（I）A(I)=TK=N/2JK?J=J+KNJ=J-K K=K/2Y四、设计过程3 Matlab程序实现3.1源程序fs=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); %播放语音信号%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=2M;if length(x1)N x1=x1,zeros(1,N-length(x1); % 若x的长度不是2的幂，补零到2的整数幂 end %数据倒序操作J=0;%给倒序数赋初值for I=0:N-1;%按序交换数据和算倒序数 if I=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=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; % 蝶形运算， 注意必须先进行减法运算，然后进行加法运算，否则要使用中间变量来传递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);3.1实验结果图用1.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 采样后的语音信号的时域波形图7 X1的频率响应特性图8 FFT函数实现的x1的频谱与编写的FFT程序实现的x1的频谱的比较数据分析取不同的点数进行FFT变换，经观察，编写FFT程序得到的语音信号的频谱图与Matlab中自带的FFT函数得到的语音信号频谱图总是基本一致，但是如果输入的N1和N2差值加1不是2的整数次幂就会有细微差别，这是因为编写的快速傅里叶变换计算信号的N点傅里叶变换要求N为2 的整数次幂，不够的话信号将会被补零后运算，即参与运算的信号已经不同（差别很小），所以FFT运算后的结果也不尽相同，所示频谱图自然有细微的差别。由用 MATLAB自带FFT函数实现的频谱图与用MATLAB编写的FFT程序实现的频谱图相比较，可知，两个算法计算后的结果几乎相同，验证了自编算法的正确性。如果改变输入的值和语音信号，那么得到不一样的时域波形，频率响应，通过FFT得到频谱都会不同。需要注意的是自由输入的语音信号不同，其长度也不同，所以根据要输入的语音信号输入N1，N2的值。5 心得体会本次实习的主要内容是通过用Matlab实现FFT的设计，可以实现对一段自己录制的语音信号进行分析，并画出采样信号的时域与频域图。把自己编写的FFT算法与Matlab自带FFT算法进行比较。程序运行调试时，自己选择输入要采样的语音信号，采样频率以及要变换的范围，可以实现对不同信号的信号采样和进行不同点的FFT运算。在之前数字信号处理的学习以及完成实验的过程中，已经使用过Matlab，对其有了一些基础的了解和认识，通过这次的课程设计使我进一步了解了信号的产生，采样及频谱分析的方法，以及其中产生信号和绘制信号的基本命令和一些基础编程语言。让我感受到只有在了解课本知识的前提下，才能更好的应用这个工具，并且熟练的应用Matlab也可以很好的加深我对课程的理解，方便我的思维。这次课程设计使我了解了Matlab的使用方法，提高了自己的分析和动手实践能力。同时我相信，进一步加强对MATLAB的学习与研究对我今后的学习将会起到很大的帮助。这次的课程设计是对本学期所学知识的一次重要巩固，使得在课堂上掌握的知识得到了真正的运用。在学习的过程中和同学讨论，更明白了理论知识与实践的联系。书到用时方恨少，有些知识学会是一回事，掌握是一回事，但应用起来，确实不是那么简单的，需要很多知识的融会贯通。程序运行调试初期，曾经多次出现错误、不能产生图形等问题，但在我翻阅资料认真改正及老师的帮助下基本功能还是完成了，程序已得到一些完善，能完成基本的要求的功能。学习就是一个了解，疑惑，进而解惑的过程，这次课程设计就是提供了这样一个发现自己知识漏洞，与老师同学探讨进行解惑的的机会。通过这次课程设计实习，我更深刻的了解了Matlab的运用，重新复习了FFT 中的重要的序列倒序和蝶形变换的程序，对课本上的知识有了更深的理解，使我对数字信号处理有了系统的认知。在这里特别感谢金星老师，他给了我们很大的帮助与指导，让我们真正自己动手真正掌握了知识，感谢老师细心指导。