实验2_FFT的实现.doc

电 子 科 技 大 学
实 验 报 告
学生姓名： 学 号： 指导教师：
一、实验室名称：数字信号处理实验室
二、实验项目名称：FFT的实现
三、实验原理：
一． FFT算法思想：
1． DFT的定义：
对于有限长离散数字信号{x[n]}，0 £ n £ N-1，其离散谱{x[k]}可以由离散付氏变换（DFT）求得。DFT的定义为：
，k=0,1,…N-1
通常令，称为旋转因子。
2． 直接计算DFT的问题及FFT的基本思想：
由DFT的定义可以看出，在x[n]为复数序列的情况下，完全直接运算N点DFT需要（N-1）2次复数乘法和N（N-1）次加法。因此，对于一些相当大的N值（如1024）来说，直接计算它的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为偶数），从而又少作一半的乘法。这样一级一级的划分下去一直到最后就划分成两点的FFT运算的情况。
3． 基2按时间抽取（DIT）的FFT算法思想：
设序列长度为，L为整数（如果序列长度不满足此条件，通过在后面补零让其满足）。
将长度为的序列，先按n的奇偶分成两组：
，r=0,1,…,N/2-1
DFT化为：

上式中利用了旋转因子的可约性，即：。又令，则上式可以写成：
(k=0,1,…,N/2-1)
可以看出，分别为从中取出的N/2点偶数点和奇数点序列的N/2点DFT值，所以，一个N点序列的DFT可以用两个N/2点序列的DFT组合而成。但是，从上式可以看出，这样的组合仅表示出了前N/2点的DFT值，还需要继续利用表示的后半段本算法推导才完整。利用旋转因子的周期性，有：，则后半段的DFT值表达式：
，同样， （k=0,1,…,N/2-1），所以后半段（k=N/2,…,N-1）的DFT值可以用前半段k值表达式获得，中间还利用到，得到后半段的值表达式为：（k=0,1,…,N/2-1）。
这样，通过计算两个N/2点序列的N/2点DFT，可以组合得到N点序列的DFT值，其组合过程如下图所示：

-1

比如，一个N = 8点的FFT运算按照这种方法来计算FFT可以用下面的流程图来表示：

4． 基2按频率抽取（DIF）的FFT算法思想：
设序列长度为，L为整数（如果序列长度不满足此条件，通过在后面补零让其满足）。
在把按k的奇偶分组之前，把输入按n的顺序分成前后两半：

因为，则有，所以：

按k的奇偶来讨论，k为偶数时：

k为奇数时：
前面已经推导过，所以上面的两个等式可以写为：

通过上面的推导，的偶数点值和奇数点值分别可以由组合而成的N/2点的序列来求得，其中偶数点值为输入x[n]的前半段和后半段之和序列的N/2点DFT值，奇数点值为输入x[n]的前半段和后半段之差再与相乘序列的N/2点DFT值。
令，，则有：

这样，也可以用两个N/2点DFT来组合成一个N点DFT，组合过程如下图所示：

-1

二． 在FFT计算中使用到的MATLAB命令：
函数fft(x)可以计算R点序列的R点DFT值；而fft(x,N)则计算R点序列的N点DFT，若R>N，则直接截取R点DFT的前N点，若R<N，则x先进行补零扩展为N点序列再求N点DFT。函数ifft(X)可以计算R点的谱序列的R点IDFT值；而ifft(X,N)同fft(x,N)的情况。
四、实验目的：
离散傅氏变换（DFT）的目的是把信号由时域变换到频域，从而可以在频域分析处理信息，得到的结果再由逆DFT变换到时域。FFT是DFT的一种快速算法。在数字信号处理系统中，FFT作为一个非常重要的工具经常使用，甚至成为DSP运算能力的一个考核因素。
本实验通过直接计算DFT，利用FFT算法思想计算DFT，以及使用MATLAB函数中的FFT命令计算离散时间信号的频谱，以加深对离散信号的DFT变换及FFT算法的理解。
五、实验内容：
a) 计算实数序列的256点DFT。
b) 计算周期为1kHz的方波序列（占空比为50％，幅度取为＋/-512，采样频率为25kHz，取256点长度） 256点DFT。
六、实验器材（设备、元器件）：
安装MATLAB软件的PC机一台，DSP实验演示系统一套。
七、实验步骤：
（1） 先利用DFT定义式，编程直接计算2个要求序列的DFT值。
（2） 利用MATLAB中提供的FFT函数，计算2个要求序列的DFT值。
（3） (拓展要求)不改变序列的点数，仅改变DFT计算点数（如变为计算1024点DFT值），观察画出来的频谱与前面频谱的差别，并解释这种差别。通过这一步骤的分析，理解频谱分辨力的概念，解释如何提高频谱分辨力。
（4） 利用FFT的基本思想（基2－DIT或基2－DIF），自己编写FFT计算函数，并用该函数计算要求序列的DFT值。并对前面3个结果进行对比。
（5） (拓展要求)尝试对其他快速傅立叶变换算法（如Goertzel算法）进行MATLAB编程实现，并用它来计算要求的序列的DFT值。并与前面的结果进行对比。
（6） （拓展要求）在提供的DSP实验板上演示要求的2种序列的FFT算法（基2－DIT），用示波器观察实际计算出来的频谱结果，并与理论结果对比。

八、实验数据及结果分析：
注：本次实验在寝室电脑上完成，所用MATLAB版本为MATLAB R2010b
程序：
（1） 对要求的2种序列直接进行DFT计算的程序
clc,clear
close all

%正弦序列
Fs = 1; %设采样频率为1Hz
Ts = 1/Fs; %采样周期
N = 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-domain 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 k'), ylabel('Magnitude')
title('Magnitude of DFT samples(direct method)')

函数文件：myDFT.m
function X = myDFT(x,N)

%利用定义式计算序列的DFT
%输入参数：
% x-离散时间信号
% N-计算的DFT点数
%输出参数：
% X-N点DFT序列

WN = exp( -1i*(2*pi/N) ); %旋转因子

for t1 = 1:N
k = t1-1;
X(t1) = 0;
for t2 = 1:N
n = t2-1;
X(t1) = X(t1) + x(t2)*WN^( n*k );
end
end

（2） 对要求的2种序列进行基2－DIT和基2－DIF FFT算法程序

序列生成的代码同（1）。这里给出FFT算法的程序：
基2-DIT FFT算法程序：
函数文件：myfitfft.m
function X = myditfft(x)

%按时间抽选的基2-FFT算法
%输入参数：
% x-离散时间信号
%
%输出参数：
% X-序列x的N点DFT（N是序列长度，必须是2的整数次幂）

M = nextpow2(length(x)); % x的长度对应的2的最低次幂
N = 2^M; % 序列长度
if length(x) < N
x = [x,zeros(1,N-length(x))]; % 若x的长度不是2的整数次幂，则补零直到长度为N
end
nxd = bin2dec(fliplr(dec2bin([1:N]-1,M))) + 1; %求1：N序列序号的倒位序
X = x(nxd); %调整x输入顺序后的序列，并作为X的初始化

WN = exp(-1i*2*pi/N); %旋转因子

for L = 1:M
B = 2^(L-1); %第L级中，每个蝶形的两个输入数据相距B个点，共有B个不同的旋转因子
for J = 0 : B-1 %第L级中不同的旋转因子
p = J*2^(M-L); %旋转因子的指数
WNp = WN^p; %旋转因子的值
for k = J+1 : 2^L : N %蝶形运算
t = X(k+B)*WNp;
X(k+B) = X(k)-t;
X(k) = X(k)+t;
end
end
end

（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(A2_fft)),
xlabel('Frequency index k'), ylabel('Magnitude')
title('Magnitude of DFT samples-Periodic square wave(fft method)')

结果：
（1） 对2种要求的序列直接进行DFT计算的频域波形

（2） 对2种要求的序列进行基2－DIT和基2－DIF FFT算法频域波形

（3） 对2种要求的序列用MATLAB中提供的FFT函数计算的频域波形。

（4）（拓展要求）分析利用上面的方法画出的信号频谱与理论计算出来的频谱之间的差异，并解释这种差异。
（5）（拓展要求）保持序列点数不变，改变DFT计算点数（变为1024点），观察频谱的变化，并分析这种变化，由此讨论如何提高频谱分辨力的问题。

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

报告评分：
指导教师签字：