dsp5000指导书第八章.doc

FFT的C语言实现第8章 FFT的C语言实现TI公司的TMS320C55x DSP函数库DSPLIB是专门为TMS320C55x DSP芯片使用的C语言优化程序函数库。它包含了50多种函数数字信号处理程序，为可以用C语言调用的汇编优化目标代码。它们是用于实时计算应用的典型的程序，经过优化处理后其执行效率是非常高的。使用TI公司的TMS320C55x DSP函数库DSPLIB，不但使编程容易得多及程序容易被人理解，而且实现同样的功能代码变得非常短小，因而缩短了应用开发的时间。用DSP实现快速傅立叶变换（FFT），可以使用标准的C语言ANSI或使用汇编语言编程，但是使用DSPLIB库函数实现FFT，可以取得比标准的C语言快得多的执行速度。以下给出的实验是分别用两种不同方法来实现FFT，目的是为了让读者初步认识和使用DSPLIB。8.1 利用标准C语言实现FFT实验8.1.1实验目的n 了解FFT基本原理；n 学习55x DSP中FFT的设计和编程思想；8.1.2 实验设备n PC兼容机一台；操作系统为Windows2000 (或WindowsNT、Windows98、WindowsXP）；n 计算机安装CCS 5000或CCS 3.1。8.1.3实验要求使用CCS集成仿真环境，完成建立工程、源文件、命令文件，保存和添加文件到工程，进行编译、链接、运行和调试等操作。使用属性窗口设置项目及观察输出波形等。8.1.4实验原理1. FFT基本原理快速傅立叶变换（FFT）是离散傅立叶变换（DFT）的一种高效运算方法。对于有限长离散数字信号x(n),0nN-1,其离散傅立叶变换为 k=0,1,2,N-1 (8.1)式中称为旋转因子或蝶形因子。由上式可以看出，计算N点DFT的所有X(k)值，需要N2次复数乘法和N(N-1)次复数加法。当N较大时, N(N-1)N2,计算X(k)的运算量几乎与N2成正比，因此直接计算DFT的运算量很大，即使采用计算机也很难做到实时处理，必须加以改进。FFT主要利用DFT旋转因子的周期性与对称性来减少运算量。 周期性 (8.2) 对称性 (8.3)利用周期性与对称性，一方面可以在DFT的运算过程中合并某些项，另一方面可以把长序列的DFT分解成若干个短序列的DFT,由于DFT运算量与变换长度的平方成正比，因此分解成短序列的DFT可以大大减少运算量。常用的FFT算法有两大类，一类是按时间抽取的FFT算法（DIT-FFT），另一类是按频率抽取的FFT算法（DIF-FFT）。这里只介绍DIT-FFT，至于DIF-FFT原理类似。设有限长序列x(n)的长度为N=2L,L为整数，若不满足该条件，加零补足，显然N为偶数。此时定义两个x(n)的偶数项和奇数项序列,长度均为N/2, 即 (8.4) (8.5)则x(n)的N点序列DFT可写成 因为 所以 (8.6) 式中，X1(k)与X2(k)分别是x1(n)与x2(n)的N/2点DFT。上式表明，一个N点DFT可以分解为两个N/2点DFT。此时需要2(N/2)2+N/2N2/2次复数乘法，N(N/2-1)+N=N2/2次复数加法,可见，经过一次分解运算量就减少了接近一半。由于N/2依然是偶数，故可将两个N/2点的DFT按同样方法分解成四个N/4点的DFT，四个N/4点的DFT继续分解为八个N/8点的DFT，如此进行下去，经过L-1次分解后，就把一个N点DFT分解成N/2个2点的DFT,至此分解结束。一个完整的N=8点DIT-FFT运算蝶形图如图8-1所示。图8-1 N=8 DIT-FFT运算蝶形图DIT-FFT运算量为：当N=2L时，经过L-1次分解，整个运算过程有L级蝶形，每一级蝶形有N/2个蝶形运算，每一个蝶形运算有一次复数乘法和两次复数加法，所以整个运算过程的运算量如下 复数乘法 复数加法 直接计算DFT需要N2次复数乘法，N(N-1)次复数加法，直接计算DFT 与DIT-FFT复数乘法的运算量之比为： N越大，DIT-FFT运算量就减少得越多，FFT的优越性就更加突出。例如，当N=256时，直接计算中复数乘法次数为65536，FFT算法中复数乘法次数为1024，速度提高倍数为64。2. FFT流程图 FFT流程图如图8-2所示。按照编码逆序排列输入序列开 始波形发生FFT初始化工作变量调用波形发生子程序产生波形（2个正弦波）调用FFT子程序计算功率谱计算步长用标准C的sin函数计算当前波形值（256点）结 束用蝶形算法计算计算功率谱返回计算结果图8-2 FFT流程图8.1.5实验步骤1. 新建立工程、源文件和命令文件 分别在源文件和命令文件中,输入8.1.6程序清单中FFT的C语言源程序命令文件，并把这两个文件添加文件到工程。2. 编译工程并运行程序编译工程过程中,如有错误，修改错误，编译成功后，装载输出文件。3. 设置断点观察窗口（1）设置断点 在程序中的break point处设置断点，如图8-2所示。图8-2 在程序中的break point处设置断点4. 设置波形观察窗口（1）查看输入信号波形从主菜单中选择ViewGraphTime/Frequency命令，出现的如图8-3所示图形属性对话框，按照图8-3所示进行相应属性修改，修改好后，单击OK按钮。将程序运行到第二个断点处，得到输入信号波形如图8-4所示。 图 8-3 图形属性对话框图 8-4 输入信号波形 （2）使用CCS提供的工具观察输入信号FFT变换结果从主菜单中选择ViewGraphTime/Frequency命令，出现的如图8-5所示图形属性对话框，按照图8-5所示进行相应属性修改，修改好后，单击OK按钮。可以看到使用CCS提供的工具对输入信号形进行FFT变换结果，如8-6所示。图 8-5 图形属性对话框 图 8-6 CCS工具FFT变换结果（3）使用C语言程序计算FFT变换结果图 8-7 图形属性对话框 8-8 C语言程序计算FFT变换结果从主菜单中选择ViewGraphTime/Frequency命令，出现的如图8-7所示图形属性对话框，按照图8-5所示进行相应属性修改，修改好后，单击OK按钮。将程序运行到第三个断点处，可以看到如图-8所示的C语言程序计算FFT变换结果。对照以上波形情况，说明参数变化时，波形基本无变化，程序及其FFT结果正确。对CCS工具FFT变换结果和C语言程序计算FFT变换结果进行比较分析可以看到两个结果是一致的。读者可以修改两个信号的频率，重新编译程序，运行程序，比较CCS工具FFT变换结果和C语言程序计算FFT变换结果。8.1.6 FFT的C语言程序清单参考1. FFT.c文件#include #define pi 3.1415926#define sample_l 256#define signal_1_f 60;#define signal_2_f 200;#define signal_sample_f 512;int inputsample_l;float fWaveRsample_l,fWaveIsample_l,wsample_l;float sin_tabsample_l;float cos_tabsample_l;void init_fft_tab();void input_data();void fft(float dataRsample_l,float dataIsample_l);void main()int i;init_fft_tab();input_data();/break point for ( i=0;isample_l;i+ )fWaveRi=inputi;fWaveIi=0.0f;wi=0.0f;fft(fWaveR,fWaveI);/break pointwhile(1);/break pointvoid init_fft_tab()float wt1;float wt2; int i; for (i=0;isample_l;i+) wt1=2*pi*i*signal_1_f; wt1=wt1/signal_sample_f; wt2=2*pi*i*signal_2_f; wt2=wt2/signal_sample_f; inputi=(cos(wt1)+cos(wt2)/2*32768; void input_data()int i;for(i=0;isample_l;i+)sin_tabi=sin(2*pi*i/sample_l);cos_tabi=cos(2*pi*i/sample_l);void fft(float dataRsample_l,float dataIsample_l)int x0,x1,x2,x3,x4,x5,x6,x7,xx;int i,j,k,b,p,L;float TR,TI,temp;/* following code invert sequence */for ( i=0;isample_l;i+ )x0=x1=x2=x3=x4=x5=x6=0;x0=i&0x01; x1=(i/2)&0x01; x2=(i/4)&0x01; x3=(i/8)&0x01;x4=(i/16)&0x01; x5=(i/32)&0x01; x6=(i/64)&0x01;x7=(i/128)&0x01;xx=x0*128+x1*64+x2*32+x3*16+x4*8+x5*4+x6*2+x7;dataIxx=dataRi;for ( i=0;isample_l;i+ )dataRi=dataIi; dataIi=0; /* following code FFT */for ( L=1;L0 ) b=b*2; i-; /* b= 2(L-1) */for ( j=0;j0 ) /* p=pow(2,7-L)*j; */p=p*2; i-;p=p*j;for ( k=j;k256;k=k+2*b ) /* for (3) */TR=dataRk; TI=dataIk; temp=dataRk+b;dataRk=dataRk+dataRk+b*cos_tabp+dataIk+b*sin_tabp;dataIk=dataIk-dataRk+b*sin_tabp+dataIk+b*cos_tabp;dataRk+b=TR-dataRk+b*cos_tabp-dataIk+b*sin_tabp;dataIk+b=TI+temp*sin_tabp-dataIk+b*cos_tabp; /* END for (3) */ /* END for (2) */ /* END for (1) */for ( i=0;i DARAM .vectors: VECT .trcinit: DARAM .gblinit: DARAM frt: DARAM .cinit: DARAM .pinit: DARAM .sysinit: DARAM .bss: DARAM2 .far: DARAM2 .const: DARAM2 .switch: DARAM2 .sysmem: DARAM2 .cio: DARAM2 .MEM$obj: DARAM2 .sysheap: DARAM2 .sysstack DARAM2 .stack: DARAM2 .input: DARAM2 .fftcode: DARAM2 8.2 利用DSP库函数的C语言实现FFT实验8.2.1实验目的n 认识和了解TI公司的TMS320C55x DSP库函数DSPLIB；n 学习调用TMS320C55x DSP函数库实现实数FFT“rfft”的编程。8.2.2 实验设备n PC兼容机一台；操作系统为Windows2000 (或WindowsNT、Windows98、WindowsXP）；n 计算机安装CCS 5000或CCS 3.1。8.2.3实验要求使用软件编程产生两个周期性信号，调用DSP函数库对信号实现实数FFT计算。8.2.4实验原理1. 函数库DSPLIB介绍TI公司TMS320C55x DSP函数库DSPLIB的内部函数库一般包括以下八种类型：1) 快速傅立叶变换FFT； 2) 滤波和卷积；3) 自适应滤波；4) 相关运算；5) 数学运算；6) 三角运算；7) 混合运算；8) 矩阵运算。DSPLIB函数库由四个部分组成：1) dsplib.h由C程序组成的头文件；2) 55xdsp.lib由目标代码组成的库文件；3) 55xdsp.src允许用户增添和修改的库源文件；4) 在“55x_test”子目录下的程序范例和链接命令文件。调用DSPLIB函数库方法分三个步骤：1) 在用户应用程序中添加dsplib.h包含文件；语法：#include 2) 将用户应用程序与DSPLIB函数库的目标文件相连接，即添加55xdsp.lib或55xdspx.lib；语法：#include 3) 使用正确的链接命令文件，描述用户C55x目标板上可用的存储器空间配置。8.2.5实验步骤1. 新建工程、源文件和命令文件新建工程、C源文件和命令文件，把后两个文件添加到工程，在C源文件中输入C程序，在命令文件中输入命令文件的内容。2.编译器的选项设置1) 选择ProjectBuild Options，选择Compiler选项卡, 进行编译器选项设置。在Category栏中选处理器Preprocessor （图8-9左边框中倒数第二个高亮处），在Preprocessor栏的Include Search Path中输入库文件所在的路径：C:CCStudio_v3.1C5500dsplibinclude。一般编译器不能自动搜索库文件所在的路径，如果不给编译器正确的路径，则在调试中编译器会报出找不到库文件的错误信息。 2) 再选择Linker选项卡，进行连接器选项设置。在Category栏中选Basic（图8-10左边框中第一个高亮处），在Library Search Path中输入：C:CCStudio_v3.1C5500dsplib键入rts55.lib;55xdsp.lib，这里的两个库文件用“；”间隔。单击OK确认，如图8-10所示3. 编译工程，装载输出文件4. 设置观察时域图和FFT结果属性对话框从主菜单中选择ViewGraphTime/Frequency命令，如图8-11所示图形属性对话框进行修改，单击OK按钮，观察时域图；从主菜单中选择ViewGraphTime/Frequency命令，如图8-12所示图形属性对话框进行修改，单击OK按钮，观察频域图。5.观察程序运行结果在程序注释断点处设断点，全速运行程序，运行到第一个断点可以看到如图8-13所示图形，上半图为输入为1024点16-bit 实数矩形波，下半图为的FFT结果。运行到第二个断点可以看到如图8-14所示图形，上半图为输入为1024点16-bit 实数三角波，下半图为的FFT结果。图8-9 连接器选项设置（一） 图8-10 连接器选项设置（二）图8-11 时域图属性对话框图 8-12 频域图属性对话框图8-13 输入为矩形波（上半图）与FFT结果（下半图）图8-14 输入为三角波（上半图）与FFT结果（下半图）8.2.6程序清单和运行结果1.1024点16-bit 实数FFT程序#include #include #include #define pi 3.1415926#define NX 1024#define PEROID 32 /Period & Amplitude of the waveform1#define HPERIOD 64 /Period & Amplitude of the waveform2#define AMPLITUDE 0.5#define signal_1_f 60;#define signal_2_f 200;#define signal_sample_f 512;#define Q15c 32768 /231DATA xNX; /Waveform Data#pragma DATA_SECTION (X,.input)DATA XNX; /Twiddle DataLDATA PXNX/2; /Power Frequency SpectrumLDATA k;int i,j;void Square_Wave();void Triangle_wave();void sin_tab();void fft();void main() Square_Wave();fft(); Triangle_wave();/Break Pointfft(); exit(1); /Break Point /Terminate the programvoid Square_Wave() k=(LDATA)(AMPLITUDE*Q15c);/Waveform Generation: (Square Wave)i=0;while(1)j=0;while (1)xi=k;if (j+=PEROID) break;if (+i=NX) break; k=-k;if (i=NX) break;void Triangle_wave()float triangle_step;/Waveform Generation: Triangle wavei=0;triangle_step=2*AMPLITUDE*Q15c/HPERIOD;k=1;while(1)j=0;while (1)xi=(triangle_step*j-AMPLITUDE*Q15c)*k;if (j+=