FFT算法的DSP实现

FFT算法移植到DSP的过程

一， 学习导向

    为什么要学习5000系列DSP？我想这就像为什么要学习《信号与系统》，《模拟电路》一样，他不仅是理论性强的东西，更是直接具有应用价值的东西。废话不多说，来看一个非常基础的5000系列DSP的应用案例--FFT算法移植到DSP芯片的过程。

    这个过程牵涉到比较多的关于CCS3.3的基础性操作问题，鉴于很基础，很有用，所以会配合很多图形来讲述。也算是留下记录以备日后需要。

二， 例程的实际意义

    这个例子的实践意义在于，通过利用移植过来的FFT算法对已知的输入信号进行频谱分析，产生FFT结果，将结果显示出来，并且验证这个结果的正确性。鉴于CCS集成开发环境自带有FFT工具并且具有很方便的显示能力，所以验证过程可以以它为参照。

三， 实现步骤

    1， 开发环境配置

    本例采用5502DSP做软件仿真，环境配置如下：

    选择好之后点击“Add”即可添加到左边，然后就可退出设置向导启动CCS。

    2， 工程的建立/编译/装载/断点/运行

    首先建立工程，Project\new，弹出的对话框填好工程文档名字，名字就写FFT即可。然后需要我们手动的把CCS目录下面Myprojects文件夹下面的FFT文件夹内部添加必要的工程文件。它们是本工程需要3个最起码的文件，FFT.c，FFT.cmd，rts55.lib三个。我将会在后面附录中给出相应文件代码，另外rts55.lib库文件可以在CCS安装目录CCS\C5500\cgtools\lib目录下面找到。做完这一步之后回到CCS开发环境，对着FFT工程下面的文档右键--Add files to project选项，将刚才复制过去的文档全部添加进来。至此，工程创建成功。

    编译程序，点击“Compile File”，没有错误，进入下一步。

    点击“Rebuild All”，没有错误，进入下一步。

    装载程序，File->Load Program->FFT->Debug->FFT.out。

    断点设置，需进入FFT.c主程序，本例程只需要执行一次就可得到全部数据，故将断点设置到main()函数的while(1)处，具体做法是先把光标放到while(1)这一行，点击“Toggle Breakpoint”。做完之后，while(1)这一行代码左边会有一个暗红色的圆圈，后面程序运行的时候到了这里就会停下来。

    运行程序，点击“Run”。

    3. 观测结果

    本例程采用两个不同频率的余弦信号直接相加，参数如下：

    信号1频率：60Hz

    信号2频率：180Hz

    采样频率：800Hz

    FFT点数：256

    以上参数均在源代码开始部分的宏部分定义和修改。合成之后的离散时间信号样本存放于SignalInput[256]数组内，FFT变换之后的256个样本存放于FFT_W[256]数组内。

    首先来观察SignalInput的时域波形，View->Graph->Time\Frequency，弹出的窗口按一下设置方式设置：

    单击“OK”，将会弹出如下显示窗口：

    上图便是输入信号的时域波形，下面来看它的频谱，我们先用CCS自带的FFT数据窗口工具直接观察结果，不需要再额外的打开显示窗口了，直接对着上图右键Properties更改参数即可，按如下方式更改：

    点击“OK”，会弹出频谱图如下：

    不去质疑这个结果的正确性，因为他必须是对的。

    我想下面才是本次例程的重点，接下来利用我们自己的FFT程序产生的结果来观察，这个结果存放于FFT_W[256]数组内部，依然如同上面一样，利用窗口显示，首先还是要设置好参数：

    点击“OK”，将会弹出结果如下：

    接下来，我将在下面部分来验证这个结果的可靠性。

四， FFT结果数组的物理意义分析

    FFT因为广泛应用与各个领域，所以不同领域的人会有不同的表示习惯，特别是结果的显示方式，我在之前的博客曾经写过一篇DFT/FFT的文章，详细论述了这一差异的根源以及频谱分辨率概念的具体意义。这里不再赘述，只给出一个具体的FFT横坐标转换关系。

    在本例中，沿用了计算机编程爱好者的习惯，将256点FFT的结果样本值存放在一个数组FFT_w[256]的数组中，这样做的好处是方便理解程序，方便后续转化。但是这种显示方式不符合我们通信人的习惯，我们需要的是具体的横坐标频率值。一言蔽之，0-256个FFT样本点涵盖了一个低频信号0-FS的所有频率信息，FS是采样频率。

    假定有如下参数：

    N点FFT结果数组X[N]，其中一个元素记为X[i]；

    信号采样频率FS；

    某一个点的信号频率记为f；

    以下的对应关系必定成立：

    上式便是N点FFT数组形式与实际物理频率之间的转化关系。下面我们来利用这个转化关系来看看，上面实验的最后一步的结果，通过肉眼大致观察到数组元素两个峰值点对应的横坐标是：19和58。

   带入上式，f1=800/256*19=59.375(Hz)

                    f2=800/256*58=181.25(Hz)

   很明显，这与实际的60Hz，180Hz是基本吻合的，这证明了，这套FFT程序是可行的。

五，小结

    本例程是5000系列DSP的一个最基本最经典的应用案例，虽然整个过程并未涉及到硬件，但是要明白，5000系列DSP是为通信系统众多复杂算法服务的，所以并非纸上谈兵。如果实际运用需要，那就是输入信号采集这个过程需要单独考虑硬件设计，以及接口的考虑而已。

    后面程序附件部分列出了256点FFT程序和配套的CMD文件源代码，如果需要使用1024点FFT程序的话，需要做如下更改。

    首先要改CMD文件，因为1024点的话，数据量比较大，存储空间分配的不够，所以相关代码需要改成下面这样：

MEMORY{
DATA:       origin = 0x6000, len = 0x8000
		PROG : origin = 0x200, len = 0x5e00
		   VECT : origin = 0xF000, len = 0x100
}

    其次，FFT.c源代码的void fft(float datar[Sample_Numb],float datai[Sample_Numb])函数体内，需要新增便面x8，x9，相关代码要替换成下面这样：

            int x0,x1,x2,x3,x4,x5,x6,x7,x8,x9,xx;
            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;
            x8=(i/256)&0x01;
            x9=(i/512)&0x01;
            xx=x0*512+x1*256+x2*128+x3*64+x4*32+x5*16+x6*8+x7*4+x8*2+x9*1;

六，程序附件

//1. 主程序FFT.c
#include "math.h"
#define Sample_Numb 256   // FFT点数
#define S1_Freq 60       //信号1频率
#define S2_Freq 180      // 信号2频率
#define SampleFreq 800   //采样率
#define pi 3.1415926
#define M log(Sample_Numb)/log(2) //变更FFT点数修改fft子程序方便的需要
int SignalInput[Sample_Numb];//输入信号，S=S1+S2
float FFT_Re[Sample_Numb];// FFT实部
float FFT_Im[Sample_Numb];//FFT虚部
float FFT_W[Sample_Numb];    //功率谱
float sin_tab[Sample_Numb];
float cos_tab[Sample_Numb];
void init_fft_tab();
void input_data();
void fft(float datar[Sample_Numb],float datai[Sample_Numb]);


 void init_fft_tab(void) //输入波形的初始化
  {
   float wt1;
   float wt2;
   int i;
   for (i=0;i<Sample_Numb;i++)
    {
wt1=2*pi*i*S1_Freq;
wt1=wt1/SampleFreq;


wt2=2*pi*i*S2_Freq;
wt2=wt2/SampleFreq;
SignalInput[i]=(cos(wt1)+cos(wt2))/2*100;
}
  }


  void FFT_WNnk(void)//蝶形运算系数表计算
   {
    int i;
for(i=0;i<Sample_Numb;i++)
{
 sin_tab[i]=sin(2*pi*i/Sample_Numb);
 cos_tab[i]=cos(2*pi*i/Sample_Numb);
}
}


void fft(float datar[Sample_Numb],float datai[Sample_Numb])
{
 int x0,x1,x2,x3,x4,x5,x6,x7,xx;
 int i,j,k,b,p,L;
 float TR,TI,temp;
 
 for(i=0;i<Sample_Numb;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;
   datai[xx]=datar[i];
 }
 for(i=0;i<Sample_Numb;i++)
 {
    datar[i]=datai[i];datai[i]=0;
 }


 for(L=1;L<=M;L++)
  {
    b=1;i=L-1;
while(i>0)
{
b=b*2;i--;
}
for(j=0;j<=b-1;j++)
{
p=1;i=M-L;
while(i>0) { p=p*2;i--;}  
p=p*j;


        for(k=j;k<Sample_Numb;k=k+2*b)
{
  TR=datar[k];TI=datai[k];temp=datar[k+b];
  datar[k]=datar[k]+datar[k+b]*cos_tab[p]+datai[k+b]*sin_tab[p];
  datai[k]=datai[k]-datar[k+b]*sin_tab[p]+datai[k+b]*cos_tab[p];
  datar[k+b]=TR-datar[k+b]*cos_tab[p]-datai[k+b]*sin_tab[p];
  datai[k+b]=TI+temp*sin_tab[p]-datai[k+b]*cos_tab[p];
    }
  }
}


for(i=0;i<Sample_Numb/2;i++)
{
  FFT_W[i]=sqrt(datar[i]*datar[i]+datai[i]*datai[i]);
}
}


void main()
{
  int i;
  init_fft_tab();
  FFT_WNnk();
  for (i=0;i<Sample_Numb;i++)
   {
     FFT_Re[i]=SignalInput[i];
FFT_Im[i]=0.0f;
FFT_W[i]=0.0f;
}
fft(FFT_Re,FFT_Im);
while(1);
 }

//2. CMD文件
MEMORY {
   DATA:       origin = 0x6000,        len = 0x4000
   PROG:       origin = 0x200,         len = 0x5e00
   VECT:       origin = 0xd000,        len = 0x100
}


SECTIONS
{
    .vectors: {} > VECT
    .trcinit: {} > PROG
    .gblinit: {} > PROG
     frt:     {} > PROG
    .text:    {} > PROG
    .cinit:   {} > PROG
    .pinit:   {} > PROG
    .sysinit: {} > PROG
    .bss:     {} > DATA
    .far:     {} > DATA
    .const:   {} > DATA
    .switch:  {} > DATA
    .sysmem:  {} > DATA
    .cio:     {} > DATA
    .MEM$obj: {} > DATA
    .sysheap: {} > DATA
    .stack:   {} > DATA
    .sysstack {} > DATA
}