快速傅立叶变换FFT的DSP实现.docx
《快速傅立叶变换FFT的DSP实现.docx》由会员分享，可在线阅读，更多相关《快速傅立叶变换FFT的DSP实现.docx（10页珍藏版）》请在冰豆网上搜索。
快速傅立叶变换FFT的DSP实现
快速傅立叶变换（FFT）的DSP实现
一、摘要
在数字信号处理系统中，FFT作为一个非常重要的工具经常使用，甚至成为DSP运算能力的一个考核因素。
FFT是一种高效实现离散付氏变换的算法。
离散付氏变换的目的是把信号由时域变换到频域，从而可以在频域分析处理信息，得到的结果再由付氏逆变换到时域。
本实验的目的在于学习FFT算法，及其在TMS320C54X上的实现，并通过编程掌握C54X的存储器管理、辅助寄存器的使用、位倒序寻址方式等技巧，同时练习使用CCS的探针和图形工具。
另外在BIOS子目录下是一个使用DSP/BIOS工具实现FFT的程序。
通过该程序，你可以使用DSP/BIOS提供的分析工具评估FFT代码执行情况。
二、实验原理
1基—2按时间抽取FFT算法
对于有限长离散数字信号{x[n]}，0nN-1，其离散谱{x[k]}可以由离散付氏变换（DFT）求得。
DFT的定义为
可以方便的把它改写为如下形式：
不难看出，WN是周期性的，且周期为N，即
WN的周期性是DFT的关键性质之一。
为了强调起见，常用表达式WN取代W以便明确其周期是N。
由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运算的情况。
比如，一个N=8点的FFT运算按照这种方法来计算FFT可以用下面的流程图来表示：
关于蝶形结运算的具体原理及其推导可以参照讲义，在此就不再赘述。
2实数FFT运算
对于离散傅立叶变换（DFT）的数字计算，FFT是一种有效的方法。
一般假定输入序列是复数。
当实际输入是实数时，利用对称性质可以使计算DFT非常有效。
一个优化的实数FFT算法是一个组合以后的算法。
原始的2N个点的实输入序列组合成一个N点的复序列，之后对复序列进行N点的FFT运算，最后再由N点的复数输出拆散成2N点的复数序列，这2N点的复数序列与原始的2N点的实数输入序列的DFT输出一致。
使用这种方法，在组合输入和拆散输出的操作中，FFT运算量减半。
这样利用实数FFT算法来计算实输入序列的DFT的速度几乎是一般复FFT算法的两倍。
本实验就用这种方法实现了一个256点实数FFT（2N=256）运算。
三、系统设计框架
⒈实数FFT运算序列的存储分配
如何利用有限的DSP系统资源，合理的安排好算法使用的存储器是一个比较重要的问题。
参见FFT实验程序的CMD文件：
MEMORY
{
PAGE0:
IPROG:
origin=0x3080,len=0x1F80
VECT:
origin=0x3000,len=0x80
EPROG:
origin=0x38000,len=0x8000
PAGE1:
USERREGS:
origin=0x60,len=0x1c
BIOSREGS:
origin=0x7c,len=0x4
IDATA:
origin=0x80,len=0xB80
EDATA:
origin=0xC00,len=0x1400
}
SECTIONS
{
.vectors:
{}>VECTPAGE0
.sysregs:
{}>BIOSREGSPAGE1
.trcinit:
{}>IPROGPAGE0
.gblinit:
{}>IPROGPAGE0
.bios:
{}>IPROGPAGE0
frt:
{}>IPROGPAGE0
.text:
{}>IPROGPAGE0
.cinit:
{}>IPROGPAGE0
.pinit:
{}>IPROGPAGE0
.sysinit:
{}>IPROGPAGE0
.data{}>EDATAPAGE1
.bss:
{}>IDATAPAGE1
.far:
{}>IDATAPAGE1
.const:
{}>IDATAPAGE1
.switch:
{}>IDATAPAGE1
.sysmem:
{}>IDATAPAGE1
.cio:
{}>IDATAPAGE1
.MEM$obj:
{}>IDATAPAGE1
.sysheap:
{}>IDATAPAGE1
}
从上面的连接定位CMD文件可以了解到，程序代码安排在0x3000开始的存储器中。
其中0x3000-0x3080存放中断向量表。
FFT程序使用的正弦表、余弦表数据（.data段）安排在0xc00开始的地方。
变量（.bss段定义）存放在0x80开始的地址中。
另外，本256点实数FFT程序的输入数据缓冲为0x2300-0x23ff，FFT后功率谱的计算结果存放在0x2200-0x22ff中。
⒉基二实数FFT运算的算法
FFT的计算可以分为三步：
首先将1个N点的时域信号分成N个1点的时域信号，然后计算这N个1点时域信号的频域，得到N个频域的点，然后将这个N个频域的点按照一定的顺序加起来，就得到了我们需要的频谱。
这里每个点的意思是复数，都有实部和虚部。
第一步的信号分解按照下面的规律执行：
可以看出它是按照比特反转顺序来分解的。
第二步是计算每个点的频谱：
这一步很简单，因为一个时域的点的频谱的数值就是它自己，所以这一步什么也不需做，但需明白这时候N个点不是时域信号了，而是频域信号。
第三步是将这N个频域信号结合起来
这一步是最麻烦的一步。
就是和前面时域分解的顺序相反，将2个1点的频域信号变成1个2点的频域信号，再将2个2点的频域信号变成1个4点的频域信号，一直到结束。
这里看下如何将2个4点的频域信号变成1个8点的频域信号。
首先对1个4点的频域信号进行复制，这样能稀释时域信号，也对另1个4点的频域信号进行复制不过复制之前需要乘上正弦函数，这样得到的稀释时域信号时经过了平移的，然后将这两个频域信号加起来，如下图所示。
之所以这么做的目的是在时域分解的时候就是用这种交织的分解方式的。
以下是基本的运算，称为蝶形运算，它将2个1点的复数变成1个2点的复数。
三、程序
#include"math.h"
#definesample_1256
#definesignal_1_f128
#definesignal_2_f64
#definesignal_3_f32
#definesignal_sample_f1024
#definepi3.1415926
intinput[sample_1];
intoutput[sample_1];
floatfwaver[sample_1],fwavei[sample_1],w[sample_1];
floatsin_tab[sample_1];
floatcos_tab[sample_1];
voidinit_fft_tab（）;
voidinput_data（）;
voidfft（floatdatar[sample_1],floatdatai[sample_1]）;
voidmain（）
{
inti;
init_fft_tab（）;
input_data（）;
for（i=0;i{
fwaver[i]=input[i];
fwavei[i]=0.0f;
w[i]=0.0f;
}
fft（fwaver,fwavei）;
//while
（1）;
}
voidinit_fft_tab（）
{
floatwt1;
floatwt2;
floatwt3;
inti;
for（i=0;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;
wt3=2*pi*i*signal_3_f;
wt3=wt3/signal_sample_f;
input[i]=（cos（wt1）+cos（wt2）+cos（wt3））/3*1024;//两个信号叠加
}
}
voidinput_data（）
{
inti;
for（i=0;i{
sin_tab[i]=sin（2*pi*i/sample_1）;
cos_tab[i]=cos（2*pi*i/sample_1）;
}
}
voidfft（floatdatar[sample_1],floatdatai[sample_1]）
{
intx0,x1,x2,x3,x4,x5,x6,x7,xx;
inti,j,k,b,p,L;
floatTR,TI,temp;
for（i=0;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{
datar[i]=datai[i];datai[i]=0;
}
for（L=1;L<=8;L++）
{
b=1;i=L-1;
while（i>0）
{
b=b*2;i--;
}
for（j=0;j<=b-1;j++）
{
p=1;i=8-L;
while（i>0）
{
p=p*2;i--;
}
p=p*j;
for（k=j;k<256;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{
w[i]=sqrt（datar[i]*datar[i]+datai[i]*datai[i]）;
}
for（j=0;j{
if（j<128）output[j]=datar[i++];
else{i=0;}
if（j>127）output[j]=datai[i++];
}
}
原理图设计
四、实验结果
1.下图为输入波形对应频率为128Hz、64Hz和32Hz的三个波形相叠加的结果。
2、下图是经过FFT变换后的波形，表现为三个尖峰，因为是三个波形的叠加，出现了三个尖峰
升级会员 