FFT算法的DSP实现.docx
《FFT算法的DSP实现.docx》由会员分享，可在线阅读，更多相关《FFT算法的DSP实现.docx（23页珍藏版）》请在冰点文库上搜索。
FFT算法的DSP实现
FFT算法的DSP实现
对于离散傅里叶变换（DFT）的数字计算，FFT是一种有效的方法。
一般假定输入序列是复数。
当实际输入是实数时，利用对称性质可以使计算DFT非常有效。
一个优化的实数FFT算法是一个组合以后的算法。
原始的2N个点的实输入序列组合成一个N点的复序列，之后对复序列进行N点的FFT运算，最后再由N点的复数输出拆散成2N点的复数序列，这2N点的复数序列与原始的2N点的实数输入序列的DFT输出一致。
使用这种方法，在组合输入和拆散输出的操作中，FFT运算量减半。
这样利用实数FFT算法来计算实输入序列的DFT的速度几乎是一般FFT算法的两倍。
下面用这种方法来实现一个256点实数FFT（2N=256）运算。
1.实数FFT运算序列的存储分配
如何利用有限的DSP系统资源，合理的安排好算法使用的存储器是一个比较重要的问题。
本文中，程序代码安排在0x3000开始的存储器中，其中0x3000~0x3080存放中断向量表，FFT程序使用的正弦表和余弦表数据（.data段）安排在OxcOO开始的地方，变量（.bss段定
义）存放在0x80开始的地址中。
另外，本文中256点实数FFT程序的数据缓冲位Ox23OO~Ox23ff,FFT变换后功率谱的计算结果存放在Ox22OO~Ox22ff中。
连续定位.cmd文件程序如下：
MEMORY{
PAGEO:
IPROG:
origin=Ox3O8O,len=Ox1F8O
VECT:
lorigin=Ox3OOO,len=Ox8O
EPROG:
origin=Ox38OOO,len=Ox8OOO
PAGE1:
USERREGS:
origin=Ox6O,len=Ox1c
BIOSREGS:
origin=Ox7c,len=Ox4
IDATA:
origin=Ox8O,len=OxB8O
}
SECTIONS
{
EDATA:
origin=OxCOO,len=Ox14OO
{
.vectors:
{}>VECTPAGEO
.sysregs:
.trcinit:
.gblinit:
{}>BIOSREGSPAGE1{}>IPROGPAGEO{}>IPROGPAGEO
.bios:
frt:
{}>IPROGPAGEO{}>IPROGPAGEO
.text:
{}>IPROGPAGEO
.cinit:
{}>IPROGPAGEO
.pinit:
{}>IPROGPAGEO
.sysinit:
{}>IPROGPAGEO
.data:
.bss:
.far:
.const:
{}>EDATAPAGE1
{}>IDATAPAGE1
{}>IDATAPAGE1
{}>IDATAPAGE1
.switch:
{}>IDATAPAGE1
.sysmem:
{}>IDATAPAGE1
•cio:
{}>IDATAPAGE1
.MEM$obj:
{}>IDATAPAGE1
.sysheap:
{}>IDATAPAGE1
}
2.基2实数FFT运算的算法
该算法主要分为以下四步进行：
1）输入数据的组合和位排序
首先，原始输入的2N=256个点的实数序列复制放到标记有“d_input_addr"的相邻单元，
当成N=128点的复数序列d[n]，其中奇数地址是d[n]实部，偶数地址是d[n]的虚部，这个过程叫做组合（n为序列变量，N为常量）。
然后，把输入序列作位倒序，是为了在整个运算最后的输出中得到的序列是自然顺序，复数序列经过位倒序，存储在数据处理缓冲其中，标记为"”_data”
如图一所示，输入实数序列为a[n],n=0,1,2,3,,,255。
分离a[n]成两个序列，如图二所示，原始的输入序列是从地址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]到数据处理缓冲中，如图二所示。
2200h
2200h
R[0]=a[0]
2201h
2201h
I[0]=a[1]
2202h
2202h
R[64]=a[128]
2203h
2203h
l[64]=a[129]
2204h
2204h
R[32]=a[64]
2205h
2205h
I[32]=a[65]
2206h
2206h
R[96]=a[192]
2207h
2207h
I[96]=a[193]
2208h
2208h
R[16]=a[32]
2209h
2209h
I[16]=a[33]
220ah
220ah
R[80]=a[160]
220bh.
220bh.
I[80]=a[161]
22feh
R[127]=a[254]
22ffh
22ffh
l[127]=a[255]
2300h
A[0]
2300h
A[0]
2301h
A[1]
2301h
A[1]
2302h
A[2]
2302h
A[2]
2303h
A[3]
2303h
A[3]
2304h
A[4]
2304h
A[4]
2305h
A[5]
2305h
A[5]
2306h
A[6]
2306h
A[6]
2307h
A[7]
2307h
A[7]
2308h
:
A[8]
2308h
A[8]
2309h
A[9]
2309h
A[9]
230ah
A[10]
230ah
A[10]
230bh
A[11]
230bh
A[11]
A[254]
23ffh
A[255]
23ffh
A[255]
程序设计中，在用C54x进行位倒序组合时，使用位倒序寻址方式可以大大提高程序执行的速度和使用存储器的效率。
在这种寻址方式中，AR0存放的整数N是FFT点数的一半，
一个辅助寄存器指向一个数据存放的单元。
当使用位倒序寻址把AR0加到辅助寄存器中时，
地址以位倒序的方式产生，即进位是从左到右，而不是从右到左。
例如，当AR0=0x0060，
AR2=0x0040时，通过指令：
MARAR2+0B
就可以得到AR2位倒序寻址后的值为0x0010。
下面是0x0060（1100000）与0x0040（1000000）以位倒序方式相加的过程：
1100000
+1000000
0010000
实现256点数据位倒序存储的具体程序段如下：
:
AR2（REORDERED_DATA）
:
装入第一个位倒序数据指针
注意，在上面的程序中，输入缓冲指针AR3（即ORIGINAL_INPUT）在操作时先加1再
减1，是因为我们把输入数据相邻的两个字看成一个复数，在用寄存器间接寻址移动了一个复数（两个字的数据）之后，对AR3进行位倒序寻址之前要把AR3的值恢复到这个复数的首字的地址，这样才能保证位倒序寻址的正确。
2）N点复数FFT运算
在数据处理器里进行N点复数FFT运算。
由于在FFT运算中要用到旋转因子WN，它
是一个复数。
我们把它分为正弦和余弦部分，用Q15格式将它们存储在两个分离的表中。
每个表中有128项，对应从0°〜1800。
因为采用循环寻址地址来对表寻址，128=27<28，
因此每张表排队的开始地址就必须是8个LSB位为0的地址。
按照系统的存储器配置，把
正弦表第一项"sine_table”放在0x0D00的位置，余弦表第一项"cos-table”放在OxOEOO的位置。
根据公式
N-1
D（k）八d[n]W「k=0,1,,,N-1
n=0
利用蝶形结对d[n]进行N=128点复数FFT运算，其中
nk」（2-/N）nk2■■
WN=ecos（nk）—jsin（nk）
NN
所需的正弦值和余弦值分别以Q15的格式存储于内存区以0x0E00开始的余弦表中。
把128点的复数FFT分为七级来算，形结，第二级是计算4点的FFT蝶形结，然后是8点、16点、结计算。
最后所有的结果表示为
0x0D00开始的正弦表和以
第一级是计算2点的FFT蝶
32点、64点、128点的蝶形
D[K]=F{d[n]}=R[k]+jl[k]
其中，R[k]、I[k]分别是D[K]的实部和虚部。
FFT完成以后，结果序列D[K]就存储到数据处理缓冲器的上半部分，如图三所示，下
半部分仍然保留原始的输入序列a[n],这半部分将在第三步中被改写。
这样原始的a[n]序列
的所有DFT的信息都在D（k）中了，第三步中需要做的就是把D（k）变为最终的2N=256点复
合序列，A[k]=F{a（n）}。
2200h
R[0]
2201h
I[0]
2202h
R[1]
2203h
I[1]
2204h
R[2]
2205h
I[2]
2206h
R[3]:
2207h
I[3]
2208h
R[4]
2209h
I[4]
220ah
R[5]
220bh
I[5]
22feh
R[127]
22ffh
l[127]
2300h
A[0]
2301h
A[1]
2302h
A[2]
2303h
A[3]
2304h
A[4]
2305h
A[5]
2306h
A[6]
2307h
A[7]
2308h
A[8]
2309h
A[9]「
230ah
A[10]
230bh
A[11]
23feh
A[254]
23ffh
A[255]
注意，在实际的DSP的编程中为了节约程序运行时间，提高代码的效率，往往要用到并行程序指令。
比如：
STB,*AR3
IILD*AR3+,B
并行指令的执行效果是，使原本分开要两个指令周期才能执行完的两条指令在一个指令周期中就能执行完。
上述指令时将B移位（ASM-16）所决定的位数，存于AR3所指定的存
储单元中，同时并行执行，将AR3所指的单元中的值装入到累加器B的高位中。
由于指令
的src和dst都是Acc、B，所以存入*AR3中的值是这条指令执行以前的值。
这一步中，实现FFT计算的具体程序如下：
Fft:
;计算FFT的第一步，两点的FFT
•asg
AR1,
GROUP_COUNTER
;定义FFT计算的组指针
•asg
AR2,
PX
；AR2为指向参加蝶形运算第一个
；数据的指针
•asg
AR3,
QX
；AR2为指向参加蝶形运算第二个
；数据的指针
•asg
AR4,
WR
;定义AR4为指向余弦表的指针
•asg
AR5,
WT
;定义AR5为指向正弦表的指针
•asg
AR6,
BUTTERFLY_COUNTER
;定义AR6为指向蝶形结的指针
•asg
AR7,
DATA_PROC_BUF
；定义在第一步中的数据处理缓冲指针
•asg
AR7,
STAGE_COUNTER
;定义剩下几步中的数据处理缓冲指针
pshm
st0
pshm
ar0
pshm
bk
;保存坏境变量
SSBX
SXM
；开启符合扩展模式
STM
#K_ZERO_BK,BK
;让BK=0使*ARn+O%=*ARn+o
LD
#-1,ASM
；为避免溢出在每一步输出时右移一位
MVMMDATA_PROC_BUF,PX
；PX指向参加蝶形结运算的第一个数
；的实部（PR）
LD
*PX,16,A
；AH:
=PR
STM
#fft_data+K_DATA_IDX_1,
QX；指向参加蝶形运算的第二个数的实
;部（QR）
STM
#K_FFT_SIZE/2-1,BRC
；设置块循环计数器
RPTBDSTAGELEND-1
;语句重复执行的范围到地址
;STAGELEND-1处
STM
#K_DATA_IDX_1+1,AR0
;延迟执行的两字节的指令（该指令
;不重复执行）
SUB
*QX,16,A,B
;BH:
=PR-QR
ADD
*QX,16,A
;AH:
=PR+QR
STH
A,ASM,*PX+
;PR':
=（PR+QR）/2
STB,
*QX+
;QR':
=（PR-QR）/2
IILD
*PX,A
;AH:
=PI
SUB
*QX,16,A,B
;BH：
=PI-QI
ADD
*QX,16,A
;AH：
=PI+QI
STH
A,ASM,*PX+0%
;PI'：
=（PI+QI）/2
STB,*QX+0%
;QI'：
=（PI-QI）/2
ILD
*PX,A
;AH:
=nextPR
stagelend:
;计算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
RPTBDSTAGEL2END-1
;语句重复执行的范围到地址
;STAGEL1END-1处
STM
#K_DATA_IDX_2+1,AR0
;初始化AR0以被循环寻址以下是第
;一个蝶形
;运算过程
SUB
*QX,16,A,B
;BH:
=PR-QR
ADD
*QX,16,A
;AH:
=PR+QR
STH
A,ASM,*PX+
;PR':
=（PR+QR）/2
STB,*QX+
;QR':
=（PR-QR）/2
ILD
*PX,A
;AH:
=PI
SUB
*QX,16,A,B
;BH:
=PI-QI
ADD
*QX,16,A
;AH:
=PI+QI
STH
A,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
STH
A,ASM,*PX+
；QR':
=（PR+QI）/2
SUB
*PX,*QX,A
；AH:
=PI+QR
ST
B,*QX
；QR':
=（PR-QI）/2
IILD*QX,B
；BH:
=OR
STA,*PX
；PI':
=（PI-QR）/2
IADD*PX+0%,A
；AH:
=PI+QR
ST
A,*QX+0%
；QIR':
=（PI+QR）/2
ILD*PX,A
；AH:
=PR
Stage2end:
；stage3thrustagelogN-1:
从第三个蝶
；形到第六个蝶形的过程如下
STM
#K_TWID_TBL_SIZE,BK
；BK=旋转因子表格的大小值
ST
#K_TWID_IDX_3,d_twid_idx
；初始化旋转表格索引值
STW
#K_TWID_IDX_3,AR0
；AR0=旋转表格初始索引值
STW
#COS_TABLE,WR
；初始化WR指针
STW
#SINE_TABLE,WI
；初始化WI指针
STW
#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）
LD
d_data_idx,A
ADD
*（PX）,A
STLM
A,QX
；QX指向参加蝶形运算第二个数据；的实部（QR）
MVDK
d_grps_cnt,GROUP_COUNTER
；AR1是组个数计数器
group:
；以下是每一组的运算过程
MVMD
BUTTERFLY_COUNTER,BRC
；将每一组中的蝶形的个数装入BRC
RPTBD
butterflyend-1
；重复执行至butterflyend-1处
LD
*WR,T
MPY
*QX+,A
；A:
=QR*WRIQX*QI
MAC
*WI+0%,*QX-,A
；A:
=QR*WR+QI*WI
ADD
*PX,16,A,B
；B:
=（QR*WR+QI*WI）+PR
；IQX指向QR
ST
B,*PX
；PR':
=（（QR*WR+QI*WI）+PR）/2
ISUB
*PX+B
；B:
=PR-（QR*WR+QI*WI）
STB,*QX
；QR':
=（PR-（QR*WR+QI*WI））/2
IMPY
*QX+,A
；A:
=QR*WI[T=WT]
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
IISUB*PX,B
LD*WR,T
STB,*PX+
IMPY*QX+,A
butterflyend:
PSHMAR0
MVDKd_data_idx,AR0
MAR*PX+0
MAR*QX+0
BANZSgroup,*GROUP_COUNTER-
POPMAR0
MAR*QX-
;IQX指向QR
;B:
=PI-（QR*WI-QI*WR）
;T:
=WR
;PI':
=（PI-（QR*WI-QI*WR））/2
;IPX指向PR
;A:
=QR*WRIQX指向QI
;更新指针以准备下一组蝶形的运算
;保存AR0
;ARO中装入在该步运算中每一组所
；用的蝶形的数目
;增加PX准备进行下一组的运算
;增加QX准备进行下一组的运算
；当组计数器减一后不等于零时，延
;迟跳转至group处
;恢复ARO（一字节）
；修改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
POPMbk
POPMar0
POPMst0
Fft_end:
RET
3）
；修改蝶形个数计算器
；下一步计算的数据索引翻倍
；下一步计算的组数目减少一半
；下一步计算的旋转因索引减少一半
;人只0=旋转因子索引（两字节）
；恢复环境变量
分离复数FFT的输出为奇部分和偶部分
分离FFT输出为相关的四个序列：
RPRMIP和IM,即偶实数、奇实数、偶虚数和奇
虚数四部分，以便第四部形成最终结果。
利用信号分析的理论我们把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*（l[K]+l[N-K]）
IM[K]=-IM[N-K]=0.5*（I[K]-I[N-K]）
RP[O]=R[O]
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]
RM[K]
图四显示了第三步完成以后存储器中的数据情况，RP[K]和IP[K]存储在上半部分,
和IM[K]存储在下半部分。
2200h
Rp[O]=R[O]
2201h
lp[O]=l[O]
2202h
RP[1]
2203h
IP[1]
2204h
RP[2]
2205h
IP[2]
2206h
RP[3]
2207h
IP[3]
2208h
RP[4]
2209h
IP[4]
220ah
RP[5]
220bh
IP[5]
22feh
RP[127]
22ffh
IP[127]
2300h
A[0]
2301h
A[1]
2302h
IM[127]
2303h
RM[127]
2304h
IM[126]
2305h
RM[126]
2306h
IM[125]
2307h
RM[125]
2308h
IM[124]
2309h
RM[124]
230ah
IM[123]
230bh
RM[123]
23feh
IM[1]
23ffh
RM[1]
这一过程的程序代码如下删除两个字
unpack
•asgAR2,XP_K
.asgAR3,XP_Nminusk
.asgAR6,XM_K
.asgAR7,XM_Nminusk
STM#fft_data+2,XP_K
；AR2指向R[K]（tempRP[K]）
ST
升级会员 