实验2_FFT的实现.doc

1、电 子 科 技 大 学实 验 报 告学生姓名： 学 号： 指导教师： 一、实验室名称：数字信号处理实验室二、实验项目名称：FFT 的实现三、实验原理：一FFT 算法思想：1DFT 的定义：对于有限长离散数字信号xn，0 n N-1，其离散谱xk可以由离散付氏变换（DFT）求得。DFT 的定义为：，k=0,1,N-1210NjkNnXkxe通常令 ，称为旋转因子。2jNeW2直接计算 DFT 的问题及 FFT 的基本思想：由 DFT 的定义可以看出，在 xn为复数序列的情况下，完全直接运算 N 点DFT 需要（N-1） 2次复数乘法和 N（N-1）次加法。因此，对于一些相当大的 N值（如 102

2、4）来说，直接计算它的 DFT 所作的计算量是很大的。FFT 的基本思想在于，将原有的 N 点序列分成两个较短的序列，这些序列的 DFT 可以很简单的组合起来得到原序列的 DFT。例如，若 N 为偶数，将原有的 N 点序列分成两个（N/2）点序列，那么计算 N 点 DFT 将只需要约(N/2) 2 2=N2/2 次复数乘法。即比直接计算少作一半乘法。因子（N/2） 2表示直接计算（N/2）点 DFT 所需要的乘法次数，而乘数 2 代表必须完成两个 DFT。上述处理方法可以反复使用，即（N/2）点的 DFT 计算也可以化成两个（N/4）点的DFT（假定 N/2 为偶数） ，从而又少作一半的乘法。

3、这样一级一级的划分下去一直到最后就划分成两点的 FFT 运算的情况。3基 2 按时间抽取（DIT）的 FFT 算法思想：设序列长度为 ，L 为整数（如果序列长度不满足此条件，通过在后2N面补零让其满足） 。将长度为 的序列 ，先按 n 的奇偶分成两组：2LN(0,1.)xnN，r=0,1,N/2-12rDFT 化为： 1/21/21(21)000/21/212200/21/21/2/200NNNnkrk rkNnrNrk kNNr rrkkrkr rXkDFTxWxxWxx 上式中利用了旋转因子的可约性，即： 。又令/2rkrkNW，则上式可以写成：/21/211/2/200,NNrk rkN

4、r rXkxWXx(k=0,1,N/2-1)12kNX可以看出， 分别为从 中取出的 N/2 点偶数点和奇数点序列12,k的 N/2 点 DFT 值，所以，一个 N 点序列的 DFT 可以用两个 N/2 点序列的 DFT 组合而成。但是，从上式可以看出，这样的组合仅表示出了 前 N/2 点的 DFTXk值，还需要继续利用 表示 的后半段本算法推导才完整。利用旋12,XkXk转因子的周期性，有： ，则后半段的 DFT 值表达式：(/2)/rrNNW，同样， /21/21()21/ /2100Nrkrkr rXkxx22NXk（k=0,1,N/2-1） ，所以后半段（k=N/2,N-1）的 DFT

5、 值可以用前半段 k 值表达式获得，中间还利用到 ，得到后半段的 值表达()22NkkkW式为： （k=0,1,N/2-1） 。1kXkWX这样，通过计算两个 N/2 点序列 的 N/2 点 DFT ，可以12,xn12,Xk组合得到 N 点序列的 DFT 值 ，其组合过程如下图所示：k1Xk 12kNXW-1 2XknkNW12kNX比如，一个 N = 8 点的 FFT 运算按照这种方法来计算 FFT 可以用下面的流程图来表示：W0W0W2W0W2W0W1W2W3x(0)x(1)x(2)x(3)x(4)x(5)x(6)x(7) X(7)X(6)X(5)X(4)X(3)X(2)X(1)X(0)

6、W0W0W04基 2 按频率抽取（DIF）的 FFT 算法思想：设序列长度为 ，L 为整数（如果序列长度不满足此条件，通过在后2N面补零让其满足） 。在把 按 k 的奇偶分组之前，把输入按 n 的顺序分成前后两半：X1/21100/2/21/21()200/2120,0,1.NNNnknknkNnN nknkNNn knnDFTxxWxxWx A因为 ，则有 ，所以：2NW2(1)Nkk/210(),01,.2NknkNnXkxxWA按 k 的奇偶来讨论，k 为偶数时： /2120,.NrnNnrxkk 为奇数时：/21(21)0,0,.1rnnXWN A前面已经推导过 ，所以上面的两个等式可

7、以写为：2/2rkrkNW/21 /20,01,./2rnNnrxA/21 /20,./1NnrnXWN 通过上面的推导， 的偶数点值 和奇数点值 分别可以XkX2Xr由组合而成的 N/2 点的序列来求得，其中偶数点值 为输入 xn的前半段和后半段之和序列的 N/2 点 DFT 值，奇数点值 为输入 xn的前半段和21r后半段之差再与 相乘序列的 N/2 点 DFT 值。nNW令 ， ，则有：12xx22nNxnxWA/2 /11/22/20 0,0,1.2N Nrn rnNn nXrX A这样，也可以用两个 N/2 点 DFT 来组合成一个 N 点 DFT，组合过程如下图所示：x 2xn-1

8、 2NxnnNW2nNxW二在 FFT 计算中使用到的 MATLAB 命令：函数 fft(x)可以计算 R 点序列的 R 点 DFT 值；而 fft(x,N)则计算 R 点序列的 N 点 DFT，若 RN，则直接截取 R 点 DFT 的前 N 点，若 RN，则 x 先进行补零扩展为 N 点序列再求 N 点 DFT。函数 ifft(X)可以计算 R 点的谱序列的 R 点IDFT 值；而 ifft(X,N)同 fft(x,N)的情况。四、实验目的：离散傅氏变换（DFT）的目的是把信号由时域变换到频域，从而可以在频域分析处理信息，得到的结果再由逆 DFT 变换到时域。FFT 是 DFT 的一种快速算

9、法。在数字信号处理系统中，FFT 作为一个非常重要的工具经常使用，甚至成为 DSP 运算能力的一个考核因素。本实验通过直接计算 DFT，利用 FFT 算法思想计算 DFT，以及使用 MATLAB函数中的 FFT 命令计算离散时间信号的频谱，以加深对离散信号的 DFT 变换及FFT 算法的理解。五、实验内容：a) 计算实数序列 的 256 点 DFT。5()cos,0261xnnb) 计算周期为 1kHz 的方波序列（占空比为 50，幅度取为/-512，采样频率为 25kHz，取 256 点长度） 256 点 DFT。六、实验器材（设备、元器件）：安装 MATLAB 软件的 PC 机一台，DSP

10、 实验演示系统一套。七、实验步骤：（1） 先利用 DFT 定义式，编程直接计算 2 个要求序列的 DFT 值。（2） 利用 MATLAB 中提供的 FFT 函数，计算 2 个要求序列的 DFT 值。（3） (拓展要求)不改变序列的点数，仅改变 DFT 计算点数（如变为计算 1024 点 DFT 值） ，观察画出来的频谱与前面频谱的差别，并解释这种差别。通过这一步骤的分析，理解频谱分辨力的概念，解释如何提高频谱分辨力。（4） 利用 FFT 的基本思想（基 2DIT 或基 2DIF） ，自己编写 FFT计算函数，并用该函数计算要求序列的 DFT 值。并对前面 3 个结果进行对比。（5） (拓展要求

11、)尝试对其他快速傅立叶变换算法（如 Goertzel 算法）进行 MATLAB 编程实现，并用它来计算要求的序列的 DFT 值。并与前面的结果进行对比。（6） （拓展要求）在提供的 DSP 实验板上演示要求的 2 种序列的 FFT算法（基 2DIT） ，用示波器观察实际计算出来的频谱结果，并与理论结果对比。八、实验数据及结果分析：注：本次实验在寝室电脑上完成，所用 MATLAB 版本为 MATLAB R2010b程序：（1） 对要求的 2 种序列直接进行 DFT 计算的程序clc,clearclose all%正弦序列Fs = 1; %设采样频率为 1HzTs = 1/Fs; %采样周期N =

12、 256; %采样点数n = 0:N-1;x = cos(5*pi/16*n);A = myDFT(x,256); %直接法计算DFT%周期方波序列N2 = 256;Fs2 = 25*1e3;Ts2 = 1/Fs2;f2 = linspace(0,Fs2,N2);n2 = 1:N2;x2 = 512*square(2*pi*1000/Fs2*n2, 50);A2 = myDFT(x2,256);%-作图-figure,subplot(211), plot(n,x)xlabel(Time index n), ylabel(Amplitude)title(Sinusoidal, time-doma

13、in sequence)subplot(212), stem(n,abs(A)xlabel(Frequency index k), ylabel(Magnitude)title(Magnitude of DFT samples(direct method)figure, subplot(211), plot(n2,x2)xlabel(Time index n), ylabel(Amplitude)title(Periodic square wave, time-domain sequence)subplot(212), stem(n2,abs(A2)xlabel(Frequency index

14、 k), ylabel(Magnitude)title(Magnitude of DFT samples(direct method)函数文件：myDFT.mfunction X = myDFT(x,N)%利用定义式计算序列的DFT%输入参数：% x-离散时间信号% N-计算的DFT点数%输出参数：% X-N点DFT序列WN = exp( -1i*(2*pi/N) ); %旋转因子for t1 = 1:Nk = t1-1;X(t1) = 0;for t2 = 1:Nn = t2-1;X(t1) = X(t1) + x(t2)*WN( n*k );endend（2） 对要求的 2 种序列进行基

15、2DIT 和基 2DIF FFT 算法程序序列生成的代码同（1） 。这里给出 FFT 算法的程序：基 2-DIT FFT 算法程序：函数文件：myfitfft.mfunction X = myditfft(x)%按时间抽选的基2-FFT算法%输入参数：% x-离散时间信号% %输出参数：% X-序列x的N点DFT （N是序列长度，必须是2 的整数次幂）M = nextpow2(length(x); % x的长度对应的2的最低次幂N = 2M; % 序列长度if length(x) Nx = x,zeros(1,N-length(x); % 若x的长度不是2的整数次幂，则补零直到长度为Nendn

16、xd = bin2dec(fliplr(dec2bin(1:N-1,M) + 1; %求1：N序列序号的倒位序X = x(nxd); %调整x输入顺序后的序列，并作为X 的初始化WN = exp(-1i*2*pi/N); %旋转因子for L = 1:MB = 2(L-1); %第L级中，每个蝶形的两个输入数据相距B 个点，共有B个不同的旋转因子for J = 0 : B-1 %第L 级中不同的旋转因子p = J*2(M-L); %旋转因子的指数WNp = WNp; %旋转因子的值for k = J+1 : 2L : N %蝶形运算t = X(k+B)*WNp;X(k+B) = X(k)-t;

17、X(k) = X(k)+t;endendend（3） 对要求的 2 种序列用 MATLAB 中提供的 FFT 函数进行计算的程序只要在程序（1）中的代码后加上以下这段即可：%-用 MATLAB提供的 fft函数计算DFT-A_fft = fft(x,256);A2_fft = fft(x2,256);figure, stem(n,abs(A_fft), xlabel(Frequency index k), ylabel(Magnitude)title(Magnitude of DFT samples-Sinusoidal wave(fft method)figure, stem(n2,abs(

18、A2_fft),xlabel(Frequency index k), ylabel(Magnitude)title(Magnitude of DFT samples-Periodic square wave(fft method)结果：（1） 对 2 种要求的序列直接进行 DFT 计算的频域波形0 50 100 150 200 250 300-1-0.500.51Time index nAmplitudeSinusoidal, time-domain sequence0 50 100 150 200 250 300050100150Frequency index kMagnitudeMagni

19、tude of DFT samples(direct method)0 50 100 150 200 250 300-1000-50005001000Time index nAmplitudePeriodic square wave, time-domain sequence0 50 100 150 200 250 30002468x 104Frequency index kMagnitudeMagnitude of DFT samples(direct method)（2） 对 2 种要求的序列进行基 2DIT 和基 2DIF FFT 算法频域波形（3） 对 2 种要求的序列用 MATLAB

20、 中提供的 FFT 函数计算的频域波形。0 50 100 150 200 250 300020406080100120140Frequency index kMagnitudeMagnitude of DFT samples-Sinusoidal wave(fft method)0 50 100 150 200 250 300012345678x 104Frequency index kMagnitudeMagnitude of DFT samples-Periodic square wave(fft method)（4） （拓展要求）分析利用上面的方法画出的信号频谱与理论计算出来的频谱之间的

21、差异，并解释这种差异。（5） （拓展要求）保持序列点数不变，改变 DFT 计算点数（变为 1024 点） ，观察频谱的变化，并分析这种变化，由此讨论如何提高频谱分辨力的问题。0 200 400 600 800 1000 1200020406080100120140MagnitudeFrequency index k1024 DFT samples of sinusoidal sequence0 200 400 600 800 1000 12000123456789x 104MagnitudeFrequency index k1024 DFT samples of square wave seq

22、uence可见，将 DFT 计算点数由 256 点增加到 1024 点后，频谱变得更密了。这是由于“栅栏效应 ”通过栅栏观察频谱，可能造成一些频谱丢失。增加 DFT计算点数相当于时域补零，使得谱线更密。但是补零并不能提高频谱分辨率，提高频谱分辨率的有效方法是增加采样信号的实际长度。九、实验结论：1. 直接用 DFT 定义和用 FFT 算法计算的 DFT 结果一致，但是 FFT 算法的效率更高。2. 时域补零并不能提高频谱分辨率，有效的方法是增加采样信号的实际长度。十、总结及心得体会：通过本实验，加强了 MATLAB 使用技能，练习了用定义和 FFT 算法计算序列的 DFT；对于 FFT 算法，用基 2-DIT 实现，由于时间原因，没有做基 2-DIF，但基本思想一致，只是变成频域抽样；实验加深了我对 FFT 算法思想的理解，学到了不少新知识、新技巧，将许多知识运用到了实践中。十一、对本实验过程及方法、手段的改进建议：1.建议在实验指导书中增加一些内容，如编程思想、码位倒序等，因为具体的编程实现和仅了解算法思想之间还是有一定差别的。2.建议适当调整实验时间段，不要和考试时间段重合，让同学们有机会认真研究实验内容，并有时间和老师交流。报告评分：指导教师签字：