FFT算法的DSP实现
FFT 算法的DSP 实现
对于离散傅里叶变换（DFT）的数字计算，FFT是一种有效的方法。一般假定输入序列是复数。当实际输入是实数时，利用对称性质可以使计算DFT 非常有效。
一个优化的实数FFT算法是一个组合以后的算法。原始的2N个点的实输入序列组合成一个N 点的复序列，之后对复序列进行N 点的FFT 运算，最后再由N 点的复数输出拆散成2N点的复数序列，这 2 N点的复数序列与原始的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 {
PAGE O: IPROG: origin = Ox3O8O,len=Ox1F8O
VECT: lorigin=Ox3OOO,len=Ox8O
EPROG: origin=Ox38OOO,len=Ox8OOO
PAGE 1:
USERREGS: origin=Ox6O,len=Ox1c
BIOSREGS: origin=Ox7c,len=Ox4
IDATA: origin=Ox8O,len=OxB8O
}
SECTIONS
{
EDATA: origin=OxCOO,len=Ox14OO
{
.vectors: { } > VECT PAGE O
.sysregs:
.trcinit:
.gblinit: { } > BIOSREGS PAGE 1 { } > IPROG PAGE O { } > IPROG PAGE O
.bios:
frt:
{ } > IPROG PAGE O { } > IPROG PAGE O
.text: { } > IPROG PAGE O
.cinit: { } > IPROG PAGE O
.pinit: { } > IPROG PAGE O
.sysinit: { } > IPROG PAGE O
.data: .bss: .far:
.const: { } > EDATA PAGE 1 { } > IDATA PAGE 1 { } > IDATA PAGE 1 { } > IDATA PAGE 1
.switch: { } > IDATA PAGE 1 .sysmem: { } > IDATA PAGE1
?cio:{ } > IDATA PAGE1
.MEM$obj: { } > IDATA PAGE1
.sysheap: { } > IDATA PAGE1
}
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]+j i[n]
程序设计中，在用 C54x 进行位倒序组合时，使用位倒序寻址方式可以大大提高程序执 行的速度和使
用存储器的效率。
在这种寻址方式中，AR0存放的整数N 是FFT 点数的一半，
一个辅助寄存器指向一个数据存放的单元。 当使用位倒序寻址把 AR0加到辅助寄存器中时，
地址以位倒序的方式产生，即进位是从左到右，而不是从右到左。例如，当
AR0=0x0060 ，
AR2=0x0040时，通过指令： MAR AR2+0B
就可以得到 AR2位倒序寻址后的值为 0x0010。
下面是0x0060(1100000)与0x0040(1000000)以位倒序方式相加的过程：
1100000 + 1000000 0010000
实现256点数据位倒序存储的具体程序段如下： :AR2( REORDERED_DATA )
:装入第一个位倒序数据指针
bit_rev:
STM #d 」n put_addr, ORIGINAL_INPUT STM #fft_data, DATA_PROC_BUF
MVMM
DATA_PROC_BUF ， :在 AR3(ORIGINAL_INPUT) :放入输入地址
:在 AR7( DATA_PROC_BUF ) :放入处理后输出的地址
REORDERED DA
STM STM RPTB MVDD
MVDD
#K_FFT_SINE-1, BRC
#K_FFT_SINE, AR0 bit_rev_e nd *ORIGINAL_INPUT+, *REORDERED_DATA+
:将原始输入缓冲中的数据放入到位倒序 :缓冲中去之后，输入缓冲( AR3 )指针 :力口 1,位倒序缓冲(AR2 )指针也加1
*REORDERED_DATA+
:将原始输入缓冲中的数据放入到位倒序 :缓冲中去之后，输入缓冲( AR3 )指针 :减1，位倒序缓冲(AR2 )指针加1， :以保证位倒序寻址正确
:按位倒序寻址方式修改
AR3
*ORIGINAL_INPUT+,
MAR *ORIGINAL_INPUT+0B bit rev end
:AR0的值是输入数据数目的一半 =128
注意，在上面的程序中，输入缓冲指针 AR3 (即ORIGINAL_INPUT )在操作时先加1再
减1，是因为我们把输入数据相邻的两个字看成一个复数，在用寄存器间接寻址移动了一 个复数(两个字的数据)之后，对AR3进行位倒序寻址之前要把 AR3的值恢复到这个复数 的首字的地址，这样才能保证位倒序寻址的正确。
2) N 点复数FFT 运算
在数据处理器里进行 N 点复数FFT 运算。由于在FFT 运算中要用到旋转因子 W N ，它 是一个复数。我们把它分为正弦和余弦部分，用
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)nk
2
■■
W N =e
cos( n k) — js in(
nk)
N N
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)}。
所需的正弦值和余弦值分别以 Q15的格式存储于内存区以 0x0E00开始的余弦表中。把128点的复数FFT 分为七级来算， 形结，第二级是计算 4点的FFT 蝶形结，然后是8点、16点、 结计算。最后所有的结果表示为
注意，在实际的DSP的编程中为了节约程序运行时间，提高代码的效率，往往要用到并行程序指令。比如：
ST B, *AR3
II LD *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, BUTTERFL Y_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