实验5 fft实验.doc

基础实验5：快速傅里叶变换FFT 5.1 实验目的 (1) 了解FFT 的原理； (2) 了解使用Matlab 语言实现FFT 的方法； (3) 了解在DSP 中FFT 的设计及编程方法； (4) 熟悉对FFT 的调试方法； 5.2 实验内容 本试验要求使用FFT 变换求一个时域信号的频域特性，并从这个频域特性求出该信号的频率值。使用Matlab 语言实现对FFT 算法的仿真，然后使用DSP C语言实现对FFT 的DSP 编程。 5.3 实验原理对于有限长离散数字信号{x[n]}，0≤x≤N-1，它的频谱离散数学值{XK)}可由离散傅氏变换（DFT）求得。DFT 定义为： 也可以方便的把它改写成如下形式： 式中WN（有时简写为W）代表。不难看出，是周期性的，且周 期为N，即 的周期性是DFT 的关键之一。 由DFT 的定义可以看出，在x[n]为复数序列的情况下，完全可以直接运算N点DFT 需要N*N次复数乘法和N*(N-1）次复数加法。因此，对于一些相当大的N 值（如1024 点）来说，直接计算它的DFT 所需要的计算量是很大的。一个优化的实数FFT 算法是一个组合以后的算法。原始的2N 个点的实输入序列组合成一个N 点的复序列，然后对复序列进行N 点的FFT 运算，最后再由N 点复数输出拆散成2N 点的复数序列，这2N 点的复数序列与原始的2N 点的实数输入序列的DFT 输出一致。 FFT 的基本思想：将原来的N 点序列分成两个较短的序列，这些序列的DFT 可很简单地组合起来得到原序列的DFT。 例如，若N为偶数，将原有的N 点序列分成两个（N/2）点序列，那么计算N 点DFT 将只需要约[（N/2）*(N/2)*2]=N*N/2次复数乘法。即比直接计算少做一半乘法。 因此（N/2）*(N/2) 表示直接计算（N/2）点DFT 所需要的乘法次数，而乘数2 代表必须完成两个DFT。上述处理方法可以反复使用，即（N/2）点的DFT 计算也可以化成两个（N/4）点DFT（假定N/2 为偶数），从而又少做一半乘法。使用这种方法，在组合输入和拆散输出的操作中，FFT 的运算量减半。这样，利用实数FFT 算法来计算实输入序列的DFT 的速度几乎是一般复FFT 算法的两倍。 假定序列{x[n]}的点数N是2的幂，定义两个分别为x[n]的偶数项和基数项的（N/2）点序列x1[n]和x2[n]，即： x1[n]=x[2n] n=0，1，…，（N/2）-1 x2[n]=x[2n+1] n=0，1，…，（N/2）-1 {x[n]}的N 点DFT 可写成： 因考虑到可写成： 故X（k）可写为： X1 (K) 和X 2 (K) 是x1 (n ) 和x2 (n) 的（N/2）点DFT。上式表明，N点DFT X（k）可分解为按上式的规则加以组合的两个（N/2）点DFT。 依此类推，当N为2的整数次幂时()，由于每分解一次降低一阶幂次，所以通过M次的分解，最后全部成为一系列2点DFT运算。以上就是按时间抽取的快速傅立叶变换(FFT)算法。当需要进行变换的序列的长度不是2的整数次方的时候，为了使用以2为基的FFT，可以用末尾补零的方法，使其长度延长至2的整数次方。 2．利用FFT进行频谱分析 若信号本身是有限长的序列，计算序列的频谱就是直接对序列进行FFT运算求得，就代表了序列在之间的频谱值。 幅度谱 相位谱 5.4 FFT的matlab编程 FFT 算法的Matlab 语言编程十分简单，直接使用Matlab 下的FFT 函数就可 以了，可以参考FFT 函数的帮助文件，使得Matlab 编程简单方便。 N=128; Len=128; n=1:1:Len; fs=10000; f1=100; f2=160; x = cos(2*pi*f1*n/fs)+cos(2*pi*f2*n/fs); plot(n,x); title([ Signal sample data length= ,int2str(Len) ]); xlabel(time (milliseconds)); Y = fft(x,N); Pyy = Y.* conj(Y) / N; f =fs*(0:N/2)/N; plot(f,Pyy(1:N/2+1)) title([ Frequency content of y ,N=,int2str(N), Len=,int2str(Len) ]); xlabel(frequency (Hz)); 5.5 FFT的高级编程 FFT的DSP实现步骤主要有以下四步： （1） 产生信号序列 （2） 将输入序列位倒序。 （3） N 点的复数FFT。 （4） 产生信号的幅度频谱特性。 #include ma