快速傅立叶变换FFT的实现.docx
《快速傅立叶变换FFT的实现.docx》由会员分享，可在线阅读，更多相关《快速傅立叶变换FFT的实现.docx（24页珍藏版）》请在冰豆网上搜索。
快速傅立叶变换FFT的实现
快速傅立叶变换（FFT）的实现
一、实验目的
在数字信号处理系统中，FFT作为一个非常重要的工具经常使用，甚至成为DSP运算能力的一个考核因素。
FFT是一种高效实现离散付氏变换的算法。
离散付氏变换的目的是把信号由时域变换到频域，从而可以在频域分析处理信息，得到的结果再由付氏逆变换到时域。
本实验的目的在于学习FFT算法，及其在TMS320C54X上的实现，并通过编程掌握C54X的存储器管理、辅助寄存器的使用、位倒序寻址方式等技巧，同时练习使用CCS的探针和图形工具。
另外在BIOS子目录下是一个使用DSP/BIOS工具实现FFT的程序。
通过该程序，你可以使用DSP/BIOS提供的分析工具评估FFT代码执行情况。
二、实验原理
1基—2按时间抽取FFT算法
对于有限长离散数字信号{x[n]}，0£n£N-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运算的算法
该算法主要分为四步：
第一步，输入数据的组合和位倒序
把输入序列作位倒序，是为了在整个运算最后的输出中得到的序列是自然顺序。
首先，原始的输入的2N=256个点的实数序列复制放到标记有“d_input_addr”的相邻单元，当成N=128点的复数序列d[n]。
奇数地址是d[n]的实部，偶数地址是d[n]的虚部。
这个过程叫做组合（n是从0到无穷，指示时间的变量，N是常量）。
然后，复数序列经过位倒序，存储在数据处理缓冲器中，标记为“fft-data”。
①如图2，输入实数序列为a[n]，n=0，1，2，3，…，255。
分离a[n]成两个序列，如图3所示。
原始的输入序列是从地址0x2300到0x23FF，其余的从0x2200到0x22FF的是经过位倒序之后的组合序列：
n=0，1，2，3，…，127。
②d[n]表示复合FFT的输入，r[n]表示实部，i[n]表示虚部：
d[n]=r[n]+ji[n]
按位倒序的方式存储d[n]到数据处理缓冲中，如图2。
图2
图2
*编程技巧：
在用C54x进行位倒序组合时，使用位倒序寻址方式可以大大提高程序执行的速度和使用存储器的效率。
在这种寻址方式中，AR0存放的整数N是FFT点数的一半，一个辅助寄存器指向一数据存放的单元。
当使用位倒序寻址把AR0加到辅助寄存器中时，地址以位倒序的方式产生，即进位是从左到右，而不是从右到左。
例如，当AR0=0x0060，AR2=0x0040时，通过指令：
MARAR2+0B
我们就可以得到AR2位倒序寻址后的值为0x0010。
下面是0x0060（1100000）与0x0040（1000000）以位倒序方式相加的过程：
实现256点数据位倒序存储的具体程序段如下：
bit_rev:
STM#d_input_addr,ORIGINAL_INPUT；在AR3（ORIGINAL_INPUT）中
；放入输入地址
STM#fft_data,DATA_PROC_BUF；在AR7（DATA_PROC_BUF）中
；放入处理后输出的地址
MVMMDATA_PROC_BUF,REORDERED_DATA;AR2（REORDERED_DATA）
;中装入第一个位倒序数据指针
STM#K_FFT_SIZE-1,BRC
STM#K_FFT_SIZE,AR0;AR0的值是输入数据数目的一半=128
RPTBbit_rev_end
MVDD*ORIGINAL_INPUT+,*REORDERED_DATA+;将原始输入缓冲中的数据
;放入到位倒序缓冲中去之
;后输入缓冲（AR3）指针加1
;位倒序缓冲（AR2）指针也加
;一
MVDD*ORIGINAL_INPUT-,*REORDERED_DATA+;将原始输入缓冲中的数据
;放入到位倒序缓冲中去之
;后输入缓冲（AR3）指针减一
;位倒序缓冲（AR2）指针加一
;以保证位倒序寻址正确
MAR*ORIGINAL_INPUT+0B；按位倒序寻址方式修改AR3
bit_rev_end:
注意，在上面的程序中。
输入缓冲指针AR3（即ORIGINAL_INPUT）在操作时先加一再减一，是因为我们把输入数据相邻的两个字看成一个复数，在用寄存器间接寻址移动了一个复数（两个字的数据）之后，对AR3进行位倒序寻址之前要把AR3的值恢复到这个复数的首字的地址，这样才能保证位倒序寻址的正确。
第二步，N点复数FFT
在数据处理缓冲器里进行N点复数FFT运算。
由于在FFT运算中要用到旋转因子WN，它是一个复数。
我们把它分为正弦和余弦部分，用Q15格式将它们存储在两个分离的表中。
每个表中有128项，对应从0度到180度。
因为采用循环寻址来对表寻址，128=27<28，因此每张表排队的开始地址就必须是8个LSB位为0的地址。
参照前面图1DES系统的存储区分配，所以我们必须把正弦表第一项“sine_table”放在0x0D00的位置,余弦表第一项“cos_table”放在0x0E00的位置。
1根据公式
利用蝶形结对d[n]进行N=128点复数FFT运算，其中
所需的正弦值和余弦值分别以Q15的格式存储于内存区以0x0D00开始的正弦表和以0x0E00开始的余弦表中。
我们把128点的复数FFT分为七级来算，第一级是计算两点的FFT蝶形结，第二级是计算四点的FFT蝶形结，然后是八点、十六点、三十二点六十四点、一百二十八点的蝶形结计算。
最后所得的结果表示为：
D[k]=F{d[n]}=R[k]+jI[k]
其中，R[k]、I[k]分别是D[k]的实部和虚部。
图3
2FFT完成以后，结果序列D[k]就存储到数据处理缓冲器的上半部分，如图3所示，下半部分任然保留原始的输入序列a[n]，这半部分将在第三步中被改写。
这样原始的a[n]序列的所有DFT的信息都在D[k]中了，第三步中需要做的就是就是把D[k]变为最终的2N=256点复合序列，A[k]=F{a（n）}。
*编程技巧：
在实际的DSP的编程中为了节约程序运行时间，提高代码的效率，往往要用到并行程序指令。
比如：
STB,*AR3
||LD*AR3+,B
并行指令的执行效果是，使原本分开要两个指令周期才能执行完的两条指令在一个指令周期中中就能执行完。
上述指令是将B移位（ASM-16）所决定的位数，存于AR3所指定的存储单元中，同时并行执行，将AR3所指的单元中的值装入到累加器B的高位中。
由于指令的src和dst都是accB，所以存入*AR3中的值是这条指令执行以前的值。
这一步中，实现FFT计算的具体程序如下：
fft:
;---------stage1：
计算FFT的第一步，两点的FFT
.asgAR1,GROUP_COUNTER；定义FFT计算的组指针
.asgAR2,PX；定义AR2为指向参加蝶形运算第一个数据的指针
.asgAR3,QX；定义AR3为指向参加蝶形运算第二个数据的指针
.asgAR4,WR；定义AR4为指向余弦表的指针
.asgAR5,WI；定义AR5为指向正弦表的指针
.asgAR6,BUTTERFLY_COUNTER；定义AR6为指向蝶形结的指针
.asgAR7,DATA_PROC_BUF；定义在第一步中的数据处理缓冲指针
.asgAR7,STAGE_COUNTER；定义剩下几步中的数据处理缓冲指针
pshmst0
pshmar0
pshmbk；保存环境变量
SSBXSXM；开启符号扩展模式
STM#K_ZERO_BK,BK；让BK=0使*ARn+0%=*ARn+0
LD#-1,ASM；为避免溢出在每一步输出时右移一位
MVMMDATA_PROC_BUF,PX；PX指向参加蝶形结运算的第一个数的实部（PR）
LD*PX,16,A；AH:
=PR
STM#fft_data+K_DATA_IDX_1,QX；QX指向参加蝶形结运算的第二个数的实部（QR）
STM#K_FFT_SIZE/2-1,BRC；设置块循环计数器
RPTBDstage1end-1；语句重复执行的范围到地址stage1end-1处
STM#K_DATA_IDX_1+1,AR0；延迟执行的两字节的指令（该指令不重复执行）
SUB*QX,16,A,B；BH:
=PR-QR
ADD*QX,16,A；AH:
=PR+QR
STHA,ASM,*PX+；PR':
=（PR+QR）/2
STB,*QX+；QR':
=（PR-QR）/2
||LD*PX,A；AH:
=PI
SUB*QX,16,A,B；BH:
=PI-QI
ADD*QX,16,A；AH:
=PI+QI
STHA,ASM,*PX+0%；PI':
=（PI+QI）/2
STB,*QX+0%；QI':
=（PI-QI）/2
||LD*PX,A；AH:
=nextPR
stage1end:
;-------Stage2：
计算FFT的第二步，四点的FFT
MVMMDATA_PROC_BUF,PX；PX指向参加蝶形结运算第一个数据的实部（PR）
STM#fft_data+K_DATA_IDX_2,QX；QX指向参加蝶形结运算第二个数据的实部（QR）
STM#K_FFT_SIZE/4-1,BRC；设置块循环计数器
LD*PX,16,A；AH:
=PR
RPTBDstage2end-1；语句重复执行的范围到地址stage1end-1处
STM#K_DATA_IDX_2+1,AR0；初始化AR0以被循环寻址
;以下是第二步运算的第一个蝶形结运算过程
SUB*QX,16,A,B；BH:
=PR-QR
ADD*QX,16,A；AH:
=PR+QR
STHA,ASM,*PX+；PR':
=（PR+QR）/2
STB,*QX+；QR':
=（PR-QR）/2
||LD*PX,A；AH:
=PI
SUB*QX,16,A,B；BH:
=PI-QI
ADD*QX,16,A；AH:
=PI+QI
STHA,ASM,*PX+；PI':
=（PI+QI）/2
STHB,ASM,*QX+；QI':
=（PI-QI）/2
;以下是第二步运算的第二个蝶形结运算过程
MAR*QX+；QX中的地址加一
ADD*PX,*QX,A；AH:
=PR+QI
SUB*PX,*QX-,B；BH:
=PR-QI
STHA,ASM,*PX+；PR':
=（PR+QI）/2
SUB*PX,*QX,A；AH:
=PI-QR
STB,*QX；QR':
=（PR-QI）/2
||LD*QX+,B；BH:
=QR
STA,*PX；PI':
=（PI-QR）/2
||ADD*PX+0%,A；AH:
=PI+QR
STA,*QX+0%；QI':
=（PI+QR）/2
||LD*PX,A；AH:
=PR
stage2end:
;Stage3thruStagelogN-1：
从第三步到第六步的过程如下
STM#K_TWID_TBL_SIZE,BK；BK=旋转因子表格的大小值
ST#K_TWID_IDX_3,d_twid_idx；初始化旋转表格索引值
STM#K_TWID_IDX_3,AR0；AR0=旋转表格初始索引值
STM#cos_table,WR；初始化WR指针
STM#sine_table,WI；初始化WI指针
STM#K_LOGN-2-1,STAGE_COUNTER；初始化步骤指针
ST#K_FFT_SIZE/8-1,d_grps_cnt；初始化组指针
STM#K_FLY_COUNT_3-1,BUTTERFLY_COUNTER；初始化蝶形结指针
ST#K_DATA_IDX_3,d_data_idx；初始化输入数据的索引
stage:
；以下是每一步的运算过程
STM#fft_data,PX；PX指向参加蝶形结运算第一个数据的实部（PR）
LDd_data_idx,A
ADD*（PX）,A
STLMA,QX；QX指向参加蝶形结运算第二个数据的实部（QR）
MVDKd_grps_cnt,GROUP_COUNTER；AR1是组个数计数器:
group:
；以下是每一组的运算过程
MVMDBUTTERFLY_COUNTER,BRC；将每一组中的蝶形结的个数装入BRC
RPTBDbutterflyend-1；重复执行至butterflyend-1处
LD*WR,T
MPY*QX+,A；A:
=QR*WR||QX*QI
MAC*WI+0%,*QX-,A；A:
=QR*WR+QI*WI
ADD*PX,16,A,B；B:
=（QR*WR+QI*WI）+PR
;||QX指向QR
STB,*PX；PR':
=（（QR*WR+QI*WI）+PR）/2
||SUB*PX+,B；B:
=PR-（QR*WR+QI*WI）
STB,*QX；QR':
=（PR-（QR*WR+QI*WI））/2
||MPY*QX+,A；A:
=QR*WI[T=WI]
;||QX指向QI
MAS*QX,*WR+0%,A；A:
=QR*WI-QI*WR
ADD*PX,16,A,B；B:
=（QR*WI-QI*WR）+PI
STB,*QX+；QI':
=（（QR*WI-QI*WR）+PI）/2
;||QX指向QR
||SUB*PX,B；B:
=PI-（QR*WI-QI*WR）
LD*WR,T；T:
=WR
STB,*PX+；PI':
=（PI-（QR*WI-QI*WR））/2
;||PX指向PR
||MPY*QX+,A；A:
=QR*WR||QX指向QI
butterflyend:
；更新指针以准备下一组蝶形结的运算
PSHMAR0；保存AR0
MVDKd_data_idx,AR0；AR0中装入在该步运算中每一组所用的蝶形结的数目
MAR*PX+0；增加PX准备进行下一组的运算
MAR*QX+0；增加QX准备进行下一组的运算
BANZDgroup,*GROUP_COUNTER-；当组计数器减一后不等于零时，延迟跳转至group处
POPMAR0；恢复AR0（一字节）
MAR*QX-；修改QX以适应下一组的运算
;
；更新指针和其他索引数据以变进入下一个步骤的运算
LDd_data_idx,A
SUB#1,A,B；B=A-1
STLMB,BUTTERFLY_COUNTER；修改蝶形结个数计数器
STLA,1,d_data_idx；下一步计算的数据索引翻倍
LDd_grps_cnt,A
STLA,-1,d_grps_cnt；下一步计算的组数目减少一半
LDd_twid_idx,A
STLA,-1,d_twid_idx；下一步计算的旋转因子索引减少一半
BANZDstage,*STAGE_COUNTER-
MVDKd_twid_idx,AR0；AR0=旋转因子索引（两字节）
popmbk
popmar0
popmst0；恢复环境变量
fft_end:
RET
第三步，分离复数FFT的输出为奇部分和偶部分
分离FFT输出为相关的四个序列：
RP、RM、IP和IM，即偶实数，奇实数、偶虚数和奇虚数四部分，以便第四步形成最终结果。
1利用信号分析的理论我们把D[k]通过下面的公式分为偶实数RP[k]、奇实数RM[k]、偶虚数IP[k]和奇虚数IM[k]：
RP[k]=RP[N-k]=0.5*（R[k]+R[N-k]）
RM[k]=-RM[N-k]=0.5*（R[k]-R[N-k]）
IP[k]=IP[N-k]=0.5*（I[k]+I[N-k]）
IM[k]=-IM[N-k]=0.5*（I[k]-I[N-k]）
RP[0]=R[0]
IP[0]=I[0]
RM[0]=IM[0]=RM[N/2]=IM[N/2]=0
RP[N/2]=R[N/2]
IP[N/2]=I[N/2]
2
下面的图4显示了第三步完成以后存储器中的数据情况，RP[k]和IP[k]存储在上半部分，RM[k] 和IM[k]存储在下半部分。
图4
这一过程的程序代码如下所示：
unpack
.asgAR2,XP_k
.asgAR3,XP_Nminusk
.asgAR6,XM_k
.asgAR7,XM_Nminusk
STM#fft_data+2,XP_k;AR2指向R[k]（tempRP[k]）
STM#fft_data+2*K_FFT_SIZE-2,XP_Nminusk;AR3指向R[N-K]（tempRP[N-K]）
STM#fft_data+2*K_FFT_SIZE+3,XM_Nminusk;AR7指向tempRM[N-K]
STM#fft_data+4*K_FF
升级会员 