实验5fft实验

基础实验5：快速傅里叶变换FFT5.1 实验目的(1) 了解FFT 的原理；(2) 了解使用Matlab 语言实现FFT 的方法；(3) 了解在DSP 中FFT 的设计及编程方法；(4) 熟悉对FFT 的调试方法；5.2 实验内容本试验要求使用FFT 变换求一个时域信号的频域特性，并从这个频域特性求出该信号的频率值。使用Matlab 语言实现对FFT 算法的仿真，然后使用DSP C语言实现对FFT 的DSP 编程。5.3 实验原理1快速傅立叶变换(FFT)算法对于有限长离散数字信号xn，0xN-1，它的频谱离散数学值X(K)可由离散傅氏变换（DFT）求得。DFT 定义为：也可以方便的把它改写成如下形式：式中WN（有时简写为W）代表。不难看出，是周期性的，且周期为N，即的周期性是DFT 的关键之一。由DFT 的定义可以看出，在xn为复数序列的情况下，完全可以直接运算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 算法的两倍。假定序列xn的点数N是2的幂，定义两个分别为xn的偶数项和基数项的（N/2）点序列x1n和x2n，即：x1n=x2n n=0，1，（N/2）-1x2n=x2n+1 n=0，1，（N/2）-1xn的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 <math.h>#define pai 3.14159265358979#define Len 128#define N 1024#define L 10typedef struct int real; int imag;complex;complex dataLen;complex data_zeropadN;int outputN;static short f1=100,f2=160;static short fs=10000;void fft(complex *a,int l);main() unsigned short h; /*generate data signal*/ for(h=0x0;h<Len;h+) datah.real=( cos(2*pai*f1*h/fs) + cos(2*pai*f2*h/fs) )*0x80; datah.imag=0; for(h=0x0;h<N;h+) if(h<Len) data_zeropadh.real=datah.real; data_zeropadh.imag=0; else data_zeropadh.real=0; data_zeropadh.imag=0; fft(data_zeropad,L); for(h=0x0;h<N;h+) outputh=abs(sqrt(data_zeropadh.real*data_zeropadh.real+data_zeropadh.imag*data_zeropadh.imag); while(1) ;/*end of main*/void fft(complex *a,int l)complex w,t;unsigned short n=1,nv2,nm1,k,le,lei,ip;unsigned short i,j,m;double tmp; n<<=l; nv2=n>>1;/N/2nm1=n-1; /N-1j=0;for(i=0;i<nm1;i+)/*bit reverse*/ if(i<j) t=aj; aj=ai; ai=t; k=nv2; while(k<=j)/*if k>N/2*/ j-=k; k>>=1; j+=k;le=1;for(m=1;m<=l;m+)lei=le;le<<=1; for(j=0;j<lei;j+) tmp=j*pai/lei; w.real= 32767*cos(tmp);/*w=complex(cos(tmp),-sin(tmp);*/ w.imag=-32767*sin(tmp); for(i=j;i<n;i+=le)/*如何考虑数据溢出*/ ip=i+lei; t.real=(aip.real*w.real - aip.imag*w.imag)>>15;/*t=aip*u;*/t.imag=(aip.real*w.imag + aip.imag*w.real)>>15; aip.real=ai.real- t.real;/*aip=ai-t;*/ aip.imag=ai.imag- t.imag; ai.real=ai.real + t.real;/*ai+=t;*/ ai.imag=ai.imag + t.imag; /u*=w; /return 0;5.6 实验步骤1参照基础实验1的步骤正确连接实验设备2设置Code Composer Studio2.2 在硬件仿真(Emulator)方式下运行3. 启动Code Composer Studio2.24执行Projects->workspaceLoad Workspace 装载工作台文件F:lab5fft20100605.wks，打开lab5fft工程。5 编译工程，并按照链接命令文件lnk.cmd进行存储器配置，选择File->Load Program装载程序6、运行程序至fft函数开始处，选择菜单栏的view->Graph->Time/Frequency进入信号的时域和频域观察窗口。7. 在弹出窗口中，按照以下参数进行配置，在时域观察产生的余弦信号序列8. 在显示类型选择FFT magnitude，借助于CCS软件的DFT功能，可从频域观察信号序列9. 而实验程序用1024个点的FFT观察信号序列所得的结果如下，此时通过补零，1024个点的FFT能够分辨的出原始的两个信号，但是信号的形状没有发生变化。10. 将FFT的点数N改为128点，级数L改为7级时，通过频谱观察窗口，我们已不能分辨两个原有的信号。11. 将信号序列的长度增加到1024点时，可以得到其时域和频域的波形，此时序列的内容也更为丰富。12.通过运用ccs的Profiler功能，可以测出fft程序的运行时间5.7思考题1 对于不同的N，幅频特性会相同吗？为什么？2 FFT进行谱分析，可以应用的什么方面？5.8实验报告要求3 简述实验原理及目的；4 结合实验中所给定典型序列幅频特性曲线，与理论结果比较，并分析说明误差产生的原因以及用FFT作谱分析时有关参数的选择方法。5 总结实验所得主要结论。