FFT算法的DSP实现
合集下载
相关主题
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
.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
2）N点复数FFT运算
在数据处理器里进行N点复数FFT运算。由于在FFT运算中要用到旋转因子 ，它是一个复数。我们把它分为正弦和余弦部分，用Q15格式将它们存储在两个分离的表中。每个表中有128项，对应从 ~ 。因为采用循环寻址地址来对表寻址，128= < ，因此每张表排队的开始地址就必须是8个LSB位为0的地址。按照系统的存储器配置，把正弦表第一项“sine_table”放在0x0D00的位置，余弦表第一项“cos-table”放在0x0E00的位置。
使用这种方法，在组合输入和拆散输出的操作中，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文件程序如下：
如图一所示，输入实数序列为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]表示虚部：
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时，通过指令：
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的编程中为了节约程序运行时间，提高代码的效率，往往要用到并行程序指令。比如：
这一步中，实现FFT计算的具体程序如下：
Fft:
；计算FFT的第一步，两点的FFT
.asg AR1, GROUP_COUNTER；定义FFT计算的组指针
.asg AR2, PX；AR2为指向参加蝶形运算第一个
；数据的指针
.asg AR3, QX；AR2为指向参加蝶形运算第二个
；数据的指针
.asg AR4, WR；定义AR4为指向余弦表的指针
FFT算法的DSP实现
对于离散傅里叶变换(DFT)的数字计算，FFT是一种有效的方法。一般假定输入序列是复数。当实际输入是实数时，利用对称性质可以使计算DFT非常有效。
一个优化的实数FFT算法是一个组合以后的算法。原始的2N个点的实输入序列组合成一个N点的复序列，之后对复序列进行N点的FFT运算，最后再由N点的复数输出拆散成2N点的复数序列，这2 N点的复数序列与原始的2N点的实数输入序列的DFT输出一致。
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)}。
：放入输入地址
STM#fft_data, DATA_PROC_BUF：在AR7（DATA_PROC_BUF）中
：放入处理后输出的地址
MVMMDATA_PROC_BUF，REORDERED_DA
：AR2（REORDERED_DATA）中
：装入第一个位倒序数据指针
STM#K_FFT_SINE-1, BRC
ST B, *AR3
‖LD *AR3+, B
并行指令的执行效果是，使原本分开要两个指令周期才能执行完的两条指令在一个指令周期中就能执行完。上述指令时将B移位(ASM-16)所决定的位数，存于AR3所指定的存储单元中，同时并行执行，将AR3所指的单元中的值装入到累加器B的高位中。由于指令的src和dst都是Acc、B，所以存入*AR3中的值是这条指令执行以前的值。
.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”。
IDATA: origin=0x80,len=0xB80
EDATA: origin=0xC00,len=0x1400
}
SECTIONS
{
.vectors: { } > VECT PAGE 0
.sysregs: { } > BIOSREGS PAGE 1
.trcinit: { } > IPROG PAGE 0
d[n]=r[n]+j i[n]
按位倒序的方式存储d[n]到数据处理缓冲中，如图二所示。
2200h
2200hຫໍສະໝຸດ Baidu
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
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）
LD *PX, 16, A；AH: =PR
STM #fft_data+K_DATA_IDX_1, QX；指向参加蝶形运算的第二个数的实
；部(QR)
STM #K_FFT_SIZE/2-1, BRC；设置块循环计数器
RPTBD STAGELEND-1；语句重复执行的范围到地址
；STAGELEND-1处
STM #K_DATA_IDX_1+1, AR0；延迟执行的两字节的指令（该指令
根据公式
k=0,1,…,N-1
利用蝶形结对d[n]进行N=128点复数FFT运算，其中
所需的正弦值和余弦值分别以Q15的格式存储于内存区以0x0D00开始的正弦表和以0x0E00开始的余弦表中。把128点的复数FFT分为七级来算，第一级是计算2点的FFT蝶形结，第二级是计算4点的FFT蝶形结，然后是8点、16点、32点、64点、128点的蝶形结计算。最后所有的结果表示为
：将原始输入缓冲中的数据放入到位倒序
：缓冲中去之后，输入缓冲（AR3）指针
：减1，位倒序缓冲（AR2）指针加1，
：以保证位倒序寻址正确
MAR *ORIGINAL_INPUT+0B：按位倒序寻址方式修改AR3
bit_rev_end
注意，在上面的程序中，输入缓冲指针AR3（即ORIGINAL_INPUT）在操作时先加1再减1，是因为我们把输入数据相邻的两个字看成一个复数，在用寄存器间接寻址移动了一个复数（两个字的数据）之后，对AR3进行位倒序寻址之前要把AR3的值恢复到这个复数的首字的地址，这样才能保证位倒序寻址的正确。
.gblinit: { } > IPROG PAGE 0
.bios: { } > IPROG PAGE 0
frt: { } > IPROG PAGE 0
.text: { } > IPROG PAGE 0
.cinit: { } > IPROG PAGE 0
.pinit: { } > IPROG PAGE 0
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]
STM#K_FFT_SINE, AR0：AR0的值是输入数据数目的一半=128
RPTBbit_rev_end
MVDD*ORIGINAL_INPUT+, *REORDERED_DATA+
：将原始输入缓冲中的数据放入到位倒序
：缓冲中去之后，输入缓冲（AR3）指针
：加1，位倒序缓冲（AR2）指针也加1
MVDD*ORIGINAL_INPUT+, *REORDERED_DATA+
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]
2305h
A[5]
2305h
A[5]
2306h
A[6]
.asg AR5, WT；定义AR5为指向正弦表的指针
.asg AR6, BUTTERFLY_COUNTER；定义AR6为指向蝶形结的指针
.asg AR7, DATA_PROC_BUF；定义在第一步中的数据处理缓冲指针
.asg AR7, STAGE_COUNTER；定义剩下几步中的数据处理缓冲指针
pshm st0
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
MAR AR2+0B
就可以得到AR2位倒序寻址后的值为0x0010。
下面是0x0060(1100000)与0x0040(1000000)以位倒序方式相加的过程：
1100000
+ 1000000
0010000
实现256点数据位倒序存储的具体程序段如下：
bit_rev:
STM#d_input_addr, ORIGINAL_INPUT：在AR3(ORIGINAL_INPUT)中
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]
.data:{ } > EDATA PAGE 1
.bss:{ } > IDATA PAGE 1
.far:{ } > IDATA PAGE 1
.const:{ } > IDATA PAGE 1
.switch:{ } > IDATA PAGE 1
.sysmem:{ } > IDATA PAGE 1
2）N点复数FFT运算
在数据处理器里进行N点复数FFT运算。由于在FFT运算中要用到旋转因子 ，它是一个复数。我们把它分为正弦和余弦部分，用Q15格式将它们存储在两个分离的表中。每个表中有128项，对应从 ~ 。因为采用循环寻址地址来对表寻址，128= < ，因此每张表排队的开始地址就必须是8个LSB位为0的地址。按照系统的存储器配置，把正弦表第一项“sine_table”放在0x0D00的位置，余弦表第一项“cos-table”放在0x0E00的位置。
使用这种方法，在组合输入和拆散输出的操作中，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文件程序如下：
如图一所示，输入实数序列为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]表示虚部：
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时，通过指令：
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的编程中为了节约程序运行时间，提高代码的效率，往往要用到并行程序指令。比如：
这一步中，实现FFT计算的具体程序如下：
Fft:
；计算FFT的第一步，两点的FFT
.asg AR1, GROUP_COUNTER；定义FFT计算的组指针
.asg AR2, PX；AR2为指向参加蝶形运算第一个
；数据的指针
.asg AR3, QX；AR2为指向参加蝶形运算第二个
；数据的指针
.asg AR4, WR；定义AR4为指向余弦表的指针
FFT算法的DSP实现
对于离散傅里叶变换(DFT)的数字计算，FFT是一种有效的方法。一般假定输入序列是复数。当实际输入是实数时，利用对称性质可以使计算DFT非常有效。
一个优化的实数FFT算法是一个组合以后的算法。原始的2N个点的实输入序列组合成一个N点的复序列，之后对复序列进行N点的FFT运算，最后再由N点的复数输出拆散成2N点的复数序列，这2 N点的复数序列与原始的2N点的实数输入序列的DFT输出一致。
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)}。
：放入输入地址
STM#fft_data, DATA_PROC_BUF：在AR7（DATA_PROC_BUF）中
：放入处理后输出的地址
MVMMDATA_PROC_BUF，REORDERED_DA
：AR2（REORDERED_DATA）中
：装入第一个位倒序数据指针
STM#K_FFT_SINE-1, BRC
ST B, *AR3
‖LD *AR3+, B
并行指令的执行效果是，使原本分开要两个指令周期才能执行完的两条指令在一个指令周期中就能执行完。上述指令时将B移位(ASM-16)所决定的位数，存于AR3所指定的存储单元中，同时并行执行，将AR3所指的单元中的值装入到累加器B的高位中。由于指令的src和dst都是Acc、B，所以存入*AR3中的值是这条指令执行以前的值。
.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”。
IDATA: origin=0x80,len=0xB80
EDATA: origin=0xC00,len=0x1400
}
SECTIONS
{
.vectors: { } > VECT PAGE 0
.sysregs: { } > BIOSREGS PAGE 1
.trcinit: { } > IPROG PAGE 0
d[n]=r[n]+j i[n]
按位倒序的方式存储d[n]到数据处理缓冲中，如图二所示。
2200h
2200hຫໍສະໝຸດ Baidu
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
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）
LD *PX, 16, A；AH: =PR
STM #fft_data+K_DATA_IDX_1, QX；指向参加蝶形运算的第二个数的实
；部(QR)
STM #K_FFT_SIZE/2-1, BRC；设置块循环计数器
RPTBD STAGELEND-1；语句重复执行的范围到地址
；STAGELEND-1处
STM #K_DATA_IDX_1+1, AR0；延迟执行的两字节的指令（该指令
根据公式
k=0,1,…,N-1
利用蝶形结对d[n]进行N=128点复数FFT运算，其中
所需的正弦值和余弦值分别以Q15的格式存储于内存区以0x0D00开始的正弦表和以0x0E00开始的余弦表中。把128点的复数FFT分为七级来算，第一级是计算2点的FFT蝶形结，第二级是计算4点的FFT蝶形结，然后是8点、16点、32点、64点、128点的蝶形结计算。最后所有的结果表示为
：将原始输入缓冲中的数据放入到位倒序
：缓冲中去之后，输入缓冲（AR3）指针
：减1，位倒序缓冲（AR2）指针加1，
：以保证位倒序寻址正确
MAR *ORIGINAL_INPUT+0B：按位倒序寻址方式修改AR3
bit_rev_end
注意，在上面的程序中，输入缓冲指针AR3（即ORIGINAL_INPUT）在操作时先加1再减1，是因为我们把输入数据相邻的两个字看成一个复数，在用寄存器间接寻址移动了一个复数（两个字的数据）之后，对AR3进行位倒序寻址之前要把AR3的值恢复到这个复数的首字的地址，这样才能保证位倒序寻址的正确。
.gblinit: { } > IPROG PAGE 0
.bios: { } > IPROG PAGE 0
frt: { } > IPROG PAGE 0
.text: { } > IPROG PAGE 0
.cinit: { } > IPROG PAGE 0
.pinit: { } > IPROG PAGE 0
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]
STM#K_FFT_SINE, AR0：AR0的值是输入数据数目的一半=128
RPTBbit_rev_end
MVDD*ORIGINAL_INPUT+, *REORDERED_DATA+
：将原始输入缓冲中的数据放入到位倒序
：缓冲中去之后，输入缓冲（AR3）指针
：加1，位倒序缓冲（AR2）指针也加1
MVDD*ORIGINAL_INPUT+, *REORDERED_DATA+
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]
2305h
A[5]
2305h
A[5]
2306h
A[6]
.asg AR5, WT；定义AR5为指向正弦表的指针
.asg AR6, BUTTERFLY_COUNTER；定义AR6为指向蝶形结的指针
.asg AR7, DATA_PROC_BUF；定义在第一步中的数据处理缓冲指针
.asg AR7, STAGE_COUNTER；定义剩下几步中的数据处理缓冲指针
pshm st0
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
MAR AR2+0B
就可以得到AR2位倒序寻址后的值为0x0010。
下面是0x0060(1100000)与0x0040(1000000)以位倒序方式相加的过程：
1100000
+ 1000000
0010000
实现256点数据位倒序存储的具体程序段如下：
bit_rev:
STM#d_input_addr, ORIGINAL_INPUT：在AR3(ORIGINAL_INPUT)中
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]