FFT算法的DSP实现.doc

发布时间：2022-06-14 发布人：admin 分类：说明书资料大小：0.20M 资料格式：doc 举报版权申诉

c1705969-2e74-489c-9aab-dff82ae7443f.doc.pdf-第1页.png

第1页 / 共13页

c1705969-2e74-489c-9aab-dff82ae7443f.doc.pdf-第2页.png

第2页 / 共13页

c1705969-2e74-489c-9aab-dff82ae7443f.doc.pdf-第3页.png

第3页 / 共13页

c1705969-2e74-489c-9aab-dff82ae7443f.doc.pdf-第4页.png

第4页 / 共13页

c1705969-2e74-489c-9aab-dff82ae7443f.doc.pdf-第5页.png

第5页 / 共13页

c1705969-2e74-489c-9aab-dff82ae7443f.doc.pdf-第6页.png

第6页 / 共13页

c1705969-2e74-489c-9aab-dff82ae7443f.doc.pdf-第7页.png

第7页 / 共13页

c1705969-2e74-489c-9aab-dff82ae7443f.doc.pdf-第8页.png

第8页 / 共13页

文本预览

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 段)安排在 0xc00 开始的地方，变量(.bss 段定义 ) 存放在 0x80 开始的地址中。另外，本文中 256 点实数 FFT 程序的数据缓冲位 0x2300~0x23ff , FFT 变换后功率谱的计算结果存放在 0x2200~0x22ff 中。连续定位.cmd 文件程序如下： MEMORY { PAGE 0: PAGE 1: IPROG: VECT: EPROG: USERREGS: BIOSREGS: IDATA: EDATA: origin = 0x3080,len=0x1F80 lorigin=0x3000,len=0x80 origin=0x38000,len=0x8000 origin=0x60,len=0x1c origin=0x7c,len=0x4 origin=0x80,len=0xB80 origin=0xC00,len=0x1400 } SECTIONS { .vectors: .sysregs: .trcinit: .gblinit: .bios: frt: .text: .cinit: .pinit: .sysinit: .data: .bss: .far: .const: { } > VECT PAGE 0 { } > BIOSREGS PAGE 1 { } > IPROG PAGE 0 { } > IPROG PAGE 0 { } > IPROG PAGE 0 { } > IPROG PAGE 0 { } > IPROG PAGE 0 { } > IPROG PAGE 0 { } > IPROG PAGE 0 { } > IPROG PAGE 0 { } > EDATA PAGE 1 { } > IDATA PAGE 1 { } > IDATA PAGE 1 { } > IDATA PAGE 1

.switch: .sysmem: .cio: .MEM$obj: .sysheap: { } > IDATA PAGE 1 { } > IDATA PAGE 1 { } > IDATA PAGE 1 { } > IDATA PAGE 1 { } > IDATA PAGE 1 } 2. 基 2 实数 FFT 运算的算法该算法主要分为以下四步进行： 1）输入数据的组合和位排序首先，原始输入的 2N=256 个点的实数序列复制放到标记有“d_input_addr”的相邻单元，当成 N=128 点的复数序列 d[n]，其中奇数地址是 d[n]实部，偶数地址是 d[n]的虚部，这个过程叫做组合（n 为序列变量，N 为常量）。然后，把输入序列作位倒序，是为了在整个运算最后的输出中得到的序列是自然顺序，复数序列经过位倒序，存储在数据处理缓冲其中，标记为”fft_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] 按位倒序的方式存储 d[n]到数据处理缓冲中，如图二所示。 2200h 2201h 2202h 2203h 2204h 2205h 2206h 2207h 2208h 2209h 220ah 220bh. . . 22ffh 2300h 2301h 2302h 2303h 2304h 2305h 2306h . . . A[0] A[1] A[2] A[3] A[4] A[5] A[6] 2200h 2201h 2202h 2203h 2204h 2205h 2206h 2207h 2208h 2209h 220ah 220bh. . 22feh 22ffh 2300h 2301h 2302h 2303h 2304h 2305h 2306h R[0]=a[0] I[0]=a[1] R[64]=a[128] I[64]=a[129] R[32]=a[64] I[32]=a[65] R[96]=a[192] I[96]=a[193] R[16]=a[32] I[16]=a[33] R[80]=a[160] I[80]=a[161] R[127]=a[254] I[127]=a[255] A[0] A[1] A[2] A[3] A[4] A[5] A[6]

2307h 2308h 2309h 230ah 230bh . . . 23ffh A[7] A[8] A[9] A[10] A[11] . . . A[255] 2307h 2308h 2309h 230ah 230bh . . 23ffh A[7] A[8] A[9] A[10] A[11] . A[254] A[255] 程序设计中，在用 C54x 进行位倒序组合时，使用位倒序寻址方式可以大大提高程序执行的速度和使用存储器的效率。在这种寻址方式中，AR0 存放的整数 N 是 FFT 点数的一半，一个辅助寄存器指向一个数据存放的单元。当使用位倒序寻址把 AR0 加到辅助寄存器中时，地址以位倒序的方式产生，即进位是从左到右，而不是从右到左。例如，当 AR0=0x0060， AR2=0x0040 时，通过指令： MAR AR2+0B 就可以得到 AR2 位倒序寻址后的值为 0x0010。下面是 0x0060(1100000)与 0x0040(1000000)以位倒序方式相加的过程： 1100000 + 1000000 0010000 实现 256 点数据位倒序存储的具体程序段如下： bit_rev: STM #d_input_addr, ORIGINAL_INPUT ：在 AR3(ORIGINAL_INPUT)中 STM #fft_data, DATA_PROC_BUF ：放入输入地址：在 AR7（DATA_PROC_BUF）中：放入处理后输出的地址 MVMM DATA_PROC_BUF，REORDERED_DA ：AR2（REORDERED_DATA）中：装入第一个位倒序数据指针 STM STM RPTB MVDD #K_FFT_SINE-1, BRC #K_FFT_SINE, AR0 bit_rev_end *ORIGINAL_INPUT+, *REORDERED_DATA+ ：AR0 的值是输入数据数目的一半=128 ：将原始输入缓冲中的数据放入到位倒序：缓冲中去之后，输入缓冲（AR3）指针：加 1，位倒序缓冲（AR2）指针也加 1 MVDD *ORIGINAL_INPUT+, *REORDERED_DATA+ MAR *ORIGINAL_INPUT+0B ：将原始输入缓冲中的数据放入到位倒序：缓冲中去之后，输入缓冲（AR3）指针：减 1，位倒序缓冲（AR2）指针加 1，：以保证位倒序寻址正确：按位倒序寻址方式修改 AR3

bit_rev_end 注意，在上面的程序中，输入缓冲指针 AR3（即 ORIGINAL_INPUT）在操作时先加 1 再减 1，是因为我们把输入数据相邻的两个字看成一个复数，在用寄存器间接寻址移动了一个复数（两个字的数据）之后，对 AR3 进行位倒序寻址之前要把 AR3 的值恢复到这个复数的首字的地址，这样才能保证位倒序寻址的正确。 2）N 点复数 FFT 运算在数据处理器里进行 N 点复数 FFT 运算。由于在 FFT 运算中要用到旋转因子 NW ，它是一个复数。我们把它分为正弦和余弦部分，用 Q15 格式将它们存储在两个分离的表中。每个表中有 128 项，对应从 0o ~ 180o 。因为采用循环寻址地址来对表寻址，128= 72 < 82 ，因此每张表排队的开始地址就必须是 8 个 LSB 位为 0 的地址。按照系统的存储器配置，把正弦表第一项“sine_table”放在 0x0D00 的位置，余弦表第一项“cos-table”放在 0x0E00 的位置。根据公式 D(k) = NWN-1 d[n] nk n=0 利用蝶形结对 d[n]进行 N=128 点复数 FFT 运算，其中 k=0,1,…,N-1 2  N 2  N nk NW  e j (2 /  N nk )  cos( nk )  j sin( nk ) 所需的正弦值和余弦值分别以 Q15 的格式存储于内存区以 0x0D00 开始的正弦表和以 0x0E00 开始的余弦表中。把 128 点的复数 FFT 分为七级来算，第一级是计算 2 点的 FFT 蝶形结，第二级是计算 4 点的 FFT 蝶形结，然后是 8 点、16 点、32 点、64 点、128 点的蝶形结计算。最后所有的结果表示为 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 2201h 2202h 2203h 2204h 2205h 2206h 2207h 2208h 2209h 220ah 220bh …. R[0] I[0] R[1] I[1] R[2] I[2] R[3] I[3] R[4] I[4] R[5] I[5] ….

22feh 22ffh 2300h 2301h 2302h 2303h 2304h 2305h 2306h 2307h 2308h 2309h 230ah 230bh ….. 23feh 23ffh R[127] I[127] A[0] A[1] A[2] A[3] A[4] A[5] A[6] A[7] A[8] A[9] A[10] A[11] … A[254] A[255] 注意，在实际的 DSP 的编程中为了节约程序运行时间，提高代码的效率，往往要用到并行程序指令。比如： *AR3 *AR3+, B B, ST ‖LD 并行指令的执行效果是，使原本分开要两个指令周期才能执行完的两条指令在一个指令周期中就能执行完。上述指令时将 B 移位(ASM-16)所决定的位数，存于 AR3 所指定的存储单元中，同时并行执行，将 AR3 所指的单元中的值装入到累加器 B 的高位中。由于指令的 src 和 dst 都是 Acc、B，所以存入*AR3 中的值是这条指令执行以前的值。这一步中，实现 FFT 计算的具体程序如下： Fft: AR3, QX .asg .asg .asg .asg .asg .asg .asg .asg pshm pshm pshm SSBX AR1, GROUP_COUNTER AR2, PX ；计算 FFT 的第一步，两点的 FFT ；定义 FFT 计算的组指针；AR2 为指向参加蝶形运算第一个；数据的指针；AR2 为指向参加蝶形运算第二个；数据的指针；定义 AR4 为指向余弦表的指针 AR4, WR AR5, WT ；定义 AR5 为指向正弦表的指针 AR6, BUTTERFLY_COUNTER ；定义 AR6 为指向蝶形结的指针 AR7, DATA_PROC_BUF AR7, STAGE_COUNTER st0 ar0 bk SXM ；定义在第一步中的数据处理缓冲指针；定义剩下几步中的数据处理缓冲指针；保存环境变量；开启符合扩展模式

STM LD MVMM #K_ZERO_BK, BK #-1, ASM DATA_PROC_BUF, PX ；让 BK=0 使*ARn+0%=*ARn+o ；为避免溢出在每一步输出时右移一位；PX 指向参加蝶形结运算的第一个数；的实部（PR）；AH: =PR LD STM 16, A *PX, #fft_data+K_DATA_IDX_1, QX ；指向参加蝶形运算的第二个数的实 STM RPTBD #K_FFT_SIZE/2-1, BRC STAGELEND-1 ；部(QR) ；设置块循环计数器；语句重复执行的范围到地址；STAGELEND-1 处 STM #K_DATA_IDX_1+1, AR0 ；延迟执行的两字节的指令（该指令 SUB ADD STH STB, ‖LD SUB ADD STH ST ‖LD stagelend: 16, A, B 16, A *PX+ B *QX, *QX, A, ASM, *QX+ A 16, A, A *PX, *QX, * QX, 16, A, B, *PX, ASM, *QX+0% A *PX+0% ；不重复执行）；BH: =PR-QR ；AH: =PR+QR ；PR’: =(PR+QR)/2 ；QR’: =(PR-QR)/2 ；AH: ；BH：=PI-QI ；AH：=PI+QI ；PI’：=（PI+QI）/2 ；QI’：=（PI-QI）/2 ；AH: =next =PI PR MVMM DATA_PROC_BUF, PX STM #fft_data+K_DATA_IDX_2, QX STM LD RPTBD #K_FFT_SIZE/4-1, BRC *PX, STAGEL2END-1 16, A STM #K_DATA_IDX_2+1, AR0 ；计算 FFT 的第二步，四点的 FFT ；PX 指向参加蝶形运算第一个数据的；实部(PR) ；QX 指向参加蝶形运算第二个数据的；实部(QR) ；设置块循环计数器；AH: =PR ；语句重复执行的范围到地址；STAGEL1END-1 处；初始化 AR0 以被循环寻址以下是第 SUB ADD STH ST ‖LD SUB ADD STH STH 16, A, B 16, A *PX+ *QX, *QX, A, B, ASM, *QX+ *PX, A *QX, *QX, A, B, 16, A, B 16, A ASM, ASM, *PX+ *QX+ ；一个蝶形；运算过程；BH: =PR-QR ；AH:= PR+QR ；PR’:= (PR+QR)/2 ；QR’:= (PR-QR)/2 ；AH:=PI ；BH:=PI-QI ；AH:=PI+QI ；PI’:= (PI+QI)/2 ；QI’:= (PI-QI)/2

MAR ADD SUB STH SUB ST ‖LD ST ‖ADD ST ‖LD Stage2end: STM ST STW STW STW STW ST STM A *QX, *QX-, B *PX+ ASM, *QX, A *QX+ *PX, *PX, A, *PX, B, *QX *QX, B A, *PX *PX+0%, A, A *QX+0% *PX, A ；以下是第二个蝶形运算过程；QX 中的地址加一；AH:= PR+QI ；BH: =PR-QI ；QR’:= (PR+QI)/2 ；AH:=PI+QR ；QR’:= (PR-QI)/2 ；BH:=OR ；PI’:= (PI-QR)/2 ；AH:= PI+QR ；QIR’:= (PI+QR)/2 ；AH:=PR ；stage 3 thru stage logN-1: 从第三个蝶；形到第六个蝶形的过程如下；BK=旋转因子表格的大小值 #K_TWID_TBL_SIZE, BK #K_TWID_IDX_3, ；初始化旋转表格索引值；AR0=旋转表格初始索引值 #K_TWID_IDX_3, ；初始化 WR 指针 #COS_TABLE, WR ；初始化 WI 指针 #SINE_TABLE, WI #K_LOGN-2-1, STAGE_COUNTER ；初始化步骤指针 #K_FFT_SIZE/8-1, #K_FLY_COUNT_3-1, BUTTERFLY_COUNTER d_twid_idx AR0 ；初始化组指针 d_grps_cnt ST #K_DATA_IDX_3, d_data_idx Stage: STM #fft_data, PX LD ADD STLM d_data_idx, A *(PX), A A, QX MVDK group: MVMD RPTBD LD MPY MAC ADD d_grps_cnt, GROUP_COUNTER BUTTERFLY_COUNTER, BRC butterflyend-1 *WR, T *QX+, A *WI+0%, *PX, 16, A, B *QX-, A ST ‖SUB ST B, *PX *PX+ B *QX B, ；初始化蝶形指针；初始化输入数据的索引；以下是每一步的运算过程；PX 指向参加蝶形运算第一个数据；的实部（PR）；QX 指向参加蝶形运算第二个数据；的实部（QR）；AR1 是组个数计数器；以下是每一组的运算过程；将每一组中的蝶形的个数装入 BRC ；重复执行至 butterflyend-1 处；A:=QR*WR‖QX*QI ；A:=QR*WR+QI*WI ；B:=( QR*WR+QI*WI)+PR ；‖QX 指向 QR ；PR’:=( (QR*WR+QI*WI)+PR)/2 ；B:=PR-( QR*WR+QI*WI) ；QR’:=(PR-( QR*WR+QI*WI))/2

‖MPY *QX+, A MAS ADD ST *QX, *WR+0%, *PX, 16, B, A A, B *QX+ ; QI’: =((QR*WI-QI*WR)+PI)/2 ‖SUB LD ST *PX, B *WR, T B, *PX+ ‖MPY butterflyend: *QX+, A PSHM MVDK AR0 d_data_idx, AR0 MAR MAR BANZS *PX+0 *QX+0 group, POPM MAR AR0 *QX- *GROUP_COUNTER- A B=A-1 BUTTERFLY_COUNTER 1, d_data_idx d_grps_cnt A A d_twid_idx *STAGE_COUNTER_ A, -1, d_data_idx, #1, B; B, A, d_grps_cnt, A, d_twid_idx, A, stage, d_twid_idx, bk ar0 st0 -1, LD SUB STLM STL LD STL LD STL BANZD MVDK POPM POPM POPM Fft_end: RET [T=WT] ；A:=QR*WI ；‖QX 指向 QI ；A:= QR*WI+QI*WR ；B:= (QR*WI+QI*WR)+PI ；‖QX 指向 QR ；B:=PI- (QR*WI-QI*WR) ；T:=WR ；PI’:=(PI- (QR*WI-QI*WR))/2 ；‖PX 指向 PR ；A:=QR*WR‖QX 指向 QI ；更新指针以准备下一组蝶形的运算；保存 AR0 ；AR0 中装入在该步运算中每一组所；用的蝶形的数目；增加 PX 准备进行下一组的运算；增加 QX 准备进行下一组的运算；当组计数器减一后不等于零时，延 ;迟跳转至 group 处；恢复 AR0（一字节）；修改 QX 以适应下一组的运算更新指；针和其他索引数据以便进入下一个；步骤的运算；修改蝶形个数计算器；下一步计算的数据索引翻倍；下一步计算的组数目减少一半；下一步计算的旋转因索引减少一半 AR0 ；AR0=旋转因子索引（两字节）；恢复环境变量 3）分离复数 FFT 的输出为奇部分和偶部分分离 FFT 输出为相关的四个序列：RP、RM、IP 和 IM，即偶实数、奇实数、偶虚数和奇虚数四部分，以便第四部形成最终结果。利用信号分析的理论我们把 D(K)通过下面的公式分为偶实数 RP[K]、奇实数 RM[K]、偶

分享到：

赞收藏

资料库

FFT算法的DSP实现.doc

相关推荐

开发技术

热门标签

最新资料