FFT算法的DSP实现

.asg AR5, WT；定义AR5为指向正弦表的指针
.asg AR6, BUTTERFLY_COUNTER；定义AR6为指向蝶形结的指针
.asg AR7, DATA_PROC_BUF；定义在第一步中的数据处理缓冲指针
.asg AR7, STAGE_COUNTER；定义剩下几步中的数据处理缓冲指针
pshm st0
2）N点复数FFT运算
在数据处理器里进行N点复数FFT运算。由于在FFT运算中要用到旋转因子 ，它是一个复数。我们把它分为正弦和余弦部分，用Q15格式将它们存储在两个分离的表中。每个表中有128项，对应从 ~ 。因为采用循环寻址地址来对表寻址，128= < ，因此每张表排队的开始地址就必须是8个LSB位为0的地址。按照系统的存储器配置，把正弦表第一项“sine_table”放在0x0D00的位置，余弦表第一项“cos-table”放在0x0E00的位置。
根据公式
k=0,1,…,N-1
利用蝶形结对d[n]进行N=128点复数FFT运算，其中
所需的正弦值和余弦值分别以Q15的格式存储于内存区以0x0D00开始的正弦表和以0x0E00开始的余弦表中。把128点的复数FFT分为七级来算，第一级是计算2点的FFT蝶形结，第二级是计算4点的FFT蝶形结，然后是8点、16点、32点、64点、128点的蝶形结计算。最后所有的结果表示为
MAR AR2+0B
就可以得到AR2位倒序寻址后的值为0x0010。
下面是0x0060(1100000)与0x0040(1000000)以位倒序方式相加的过程：
1100000
+ 1000000
0010000
实现256点数据位倒序存储的具体程序段如下：
bit_rev:
STM#d_input_addr, ORIGINAL_INPUT：在AR3(ORIGINAL_INPUT)中
pshm ar0
pshm bk；保存环境变量
SSBX SXM；开启符合扩展模式
STM #K_ZERO_BK, BK；让BK=0使*ARn+0%=*ARn+o
LD #-1, ASM；为避免溢出在每一步输出时右移一位
MVMMDATA_PROC_BUF, PX；PX指向参加蝶形结运算的第一个数
；的实部（PR）
MEMORY
{
PAGE 0: IPROG: origin = 0x3080,len=0x1F80
VECT: lorigin=0x3000,len=0x80
EPROG: origin=0x38000,len=0x8000
PAGE 1: USERREGS: origin=0x60,len=0x1c
BIOSREGS: origin=0x7c,len=0x4
：将原始输入缓冲中的数据放入到位倒序
：缓冲中去之后，输入缓冲（AR3）指针
：减1，位倒序缓冲（AR2）指针加1，
：以保证位倒序寻址正确
MAR *ORIGINAL_INPUT+0B：按位倒序寻址方式修改AR3
bit_rev_end
注意，在上面的程序中，输入缓冲指针AR3（即ORIGINAL_INPUT）在操作时先加1再减1，是因为我们把输入数据相邻的两个字看成一个复数，在用寄存器间接寻址移动了一个复数（两个字的数据）之后，对AR3进行位倒序寻址之前要把AR3的值恢复到这个复数的首字的地址，这样才能保证位倒序寻址的正确。
R[127]=a[254]
22ffh
22ffh
I[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]
23]
2306h
A[6]
.cio:{ } > IDATA PAGE 1
.MEM$obj:{ } > IDATA PAGE 1
.sysheap:{ } > IDATA PAGE 1
}
2.基2实数FFT运算的算法
该算法主要分为以下四步进行：
1）输入数据的组合和位排序
首先，原始输入的2N=256个点的实数序列复制放到标记有“d_input_addr”的相邻单元，当成N=128点的复数序列d[n]，其中奇数地址是d[n]实部，偶数地址是d[n]的虚部，这个过程叫做组合（n为序列变量，N为常量）。然后，把输入序列作位倒序，是为了在整个运算最后的输出中得到的序列是自然顺序，复数序列经过位倒序，存储在数据处理缓冲其中，标记为”fft_data”。
D[K]=F{d[n]}=R[k]+jI[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
I[127]
2300h
A[0]
2301h
A[1]
IDATA: origin=0x80,len=0xB80
EDATA: origin=0xC00,len=0x1400
}
SECTIONS
{
.vectors: { } > VECT PAGE 0
.sysregs: { } > BIOSREGS PAGE 1
.trcinit: { } > IPROG PAGE 0
这一步中，实现FFT计算的具体程序如下：
Fft:
；计算FFT的第一步，两点的FFT
.asg AR1, GROUP_COUNTER；定义FFT计算的组指针
.asg AR2, PX；AR2为指向参加蝶形运算第一个
；数据的指针
.asg AR3, QX；AR2为指向参加蝶形运算第二个
；数据的指针
.asg AR4, WR；定义AR4为指向余弦表的指针
ST B, *AR3
‖LD *AR3+, B
并行指令的执行效果是，使原本分开要两个指令周期才能执行完的两条指令在一个指令周期中就能执行完。上述指令时将B移位(ASM-16)所决定的位数，存于AR3所指定的存储单元中，同时并行执行，将AR3所指的单元中的值装入到累加器B的高位中。由于指令的src和dst都是Acc、B，所以存入*AR3中的值是这条指令执行以前的值。
使用这种方法，在组合输入和拆散输出的操作中，FFT运算量减半。这样利用实数FFT算法来计算实输入序列的DFT的速度几乎是一般FFT算法的两倍。下面用这种方法来实现一个256点实数FFT（2N=256）运算。
1.实数FFT运算序列的存储分配
如何利用有限的DSP系统资源，合理的安排好算法使用的存储器是一个比较重要的问题。本文中，程序代码安排在0x3000开始的存储器中，其中0x3000~0x3080存放中断向量表，FFT程序使用的正弦表和余弦表数据(.data段)安排在0xc00开始的地方，变量(.bss段定义)存放在0x80开始的地址中。另外，本文中256点实数FFT程序的数据缓冲位0x2300~0x23ff , FFT变换后功率谱的计算结果存放在0x2200~0x22ff中。连续定位.cmd文件程序如下：
d[n]=r[n]+j i[n]
按位倒序的方式存储d[n]到数据处理缓冲中，如图二所示。
2200h
2200h
R[0]=a[0]
2201h
2201h
I[0]=a[1]
2202h
2202h
R[64]=a[128]
2203h
2203h
I[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.
.
22feh
I[80]=a[161]
.sysinit: { } > IPROG PAGE 0
.data:{ } > EDATA PAGE 1
.bss:{ } > IDATA PAGE 1
.far:{ } > IDATA PAGE 1
.const:{ } > IDATA PAGE 1
.switch:{ } > IDATA PAGE 1
.sysmem:{ } > IDATA PAGE 1
LD *PX, 16, A；AH: =PR
STM #fft_data+K_DATA_IDX_1, QX；指向参加蝶形运算的第二个数的实