512点FFT算法实现

N=512 points FFT Program ***

******************************************************************

.title "fft.asm"

.mmregs

.copy "coeff.inc"

.def _c_int00

sine: .usect "sine",512

cosine: .usect "cosine",512

fft_data: .usect "fft_data",2048

d_input: .usect "d_input",2048

fft_out: .usect "fft_out",512

STACK .usect "STACK",9

K_DATA_IDX_1 .set 2

K_DATA_IDX_2 .set 4

K_DATA_IDX_3 .set 8

K_TWID_TBL_SIZE .set 128

K_TWID_IDX_3 .set 32

K_FLY_COUNT_3 .set 4

*** N points FFT ***

K_FFT_SIZE .set 256 ;N

K_LOGN .set 9 ;LOG(N) =512 则N=9

PA0 .set 0

PA1 .set 1

.bss d_twid_idx,1

.bss d_data_idx,1

.bss d_grps_cnt,1

.sect "fft_prg" ;

*** 位码倒序Bit Reversal Routine ***

.asg AR2,REORDERED ；AR2中装入第一个位倒序数据指针.asg AR3,ORIGINAL_INPUT ；AR3中放入输入地址

.asg AR7,DA TA_PROC_BUF ；AR7中放入处理后输出的地址

.text

Start:

SSBX FRCT

STM #STACK+9,SP

STM #d_input,AR1

RPT #2*K_FFT_SIZE-1

PORTR PA1,*AR1+

STM #sine,AR1 ;正弦系数表

RPT #255 ;

MVPD sine1,*AR1+

STM #cosine,AR1 ;

RPT #255 ;

MVPD cosine1,*AR1+

STM #d_input,ORIGINAL_INPUT

STM #fft_data,DATA_PROC_BUF

MVMM DATA_PROC_BUF,REORDERED ;RECORDERED point fft_data too STM #K_FFT_SIZE-1,BRC

RPTBD bit_rev_end-1

STM #K_FFT_SIZE,AR0 ;this is double words order MVDD *ORIGINAL_INPUT+,*REORDERED+ ;Bit Reversal

MVDD *ORIGINAL_INPUT-,*REORDERED+

MAR *ORIGINAL_INPUT+0B

bit_rev_end:

* * * * FFT Code * * * * *

.asg AR1,GROUP_COUNTER ；定义FFT计算的组指针

.asg AR2,PX ；AR2指向蝶型运算的第一个数据指针.asg AR3,QX ；AR3指向蝶型运算的第二个数据指针.asg AR4,WR ；AR4指向余弦表的指针

.asg AR5,WI ；AR5指向正弦表的指针

.asg AR6,BUTTERFL Y_COUNTER ；AR6指向蝶形结的指针

.asg AR7,STAGE_COUNTER ；定义数据处理缓冲指针

* * * stage 1级蝶形运算，计算2点的FFT * * *

STM #0, BK ;让BK=0，使*ARn+0%=*ARn+0 why not=16? LD #-1, ASM ;每步输出时右移一位避免溢出

STM #fft_data, PX ；PX指向蝶型运算第一个数的实部

LD *PX,16,A ;AH:=Re[x(0)] AH=PR

STM #fft_data+K_DATA_IDX_1,QX ;QX point 蝶形运算第二个数的实部STM #K_FFT_SIZE/2-1,BRC ;stage 1 is N/2-1设置循环块计数器RPTBD stage1end-1

STM #K_DA TA_IDX_1+1,AR0

SUB *QX,16,A,B ;BH=Re[x(0)]-Re[x(4)] BH=PR-QR

ADD *QX,16,A ;AH=Re[x(0)]+ Re[x(4)] AH=PR+QR

STH A,ASM,*PX+ ;PRV=(PR+QR)/2

ST B,*QX+ ;QRV=(PR-QR)/2

||LD *PX,A ;AH=PI

SUB *QX,16,A,B ;BH=Im[x(0)]-Im[x(4)] BH=PI-QI

ADD *QX,16,A ;AH=Im[x(0)]+Im[x(4)] AH=PI+QI

STH A,ASM,*PX+0 ;PIV=(PI+QI)/2

ST B,*QX+0% ;PIV=(PI-QI)/2;BK=0 why here circle access ?

||LD *PX,A AH为下一个PR

stage1end:

* * * Stage 2 第二级计算4点的FFT* * *

STM #fft_data,PX

STM #fft_data+K_DATA_IDX_2,QX

STM #K_FFT_SIZE/4-1,BRC

LD *PX,16,A ;AH=Re[x(0)]+ Re[x(4)],AH=1

RPTBD stage2end-1

STM #K_DA TA_IDX_2+1,AR0

;1st butterfly

SUB *QX,16,A,B ;BH={Re[x(0)]+ Re[x(4)]}-{Re[x(2)]+ Re[x(6)]},BH=0 ADD *QX,16,A ;AH={Re[x(0)]+ Re[x(4)]}+{Re[x(2)]+ Re[x(6)]},AH=2 STH A,ASM,*PX+

ST B,*QX+

||LD *PX,A

SUB *QX,16,A,B ;BH={Im[x(0)]+ Im[x(4)]}-{Im[x(2)]+ Im[x(6)]},BH=0 ADD *QX,16,A ;AH={Im[x(0)]+ Im[x(4)]}+{Im[x(2)]+ Im[x(6)]},BH=0 STH A,ASM,*PX+ ；PIV=（PI+QI）/2

STH B,ASM,*QX+ ；PIV=（PI-QI）/2

;2nd butterfly

MAR *QX+ ；QX中地址加1

ADD *PX,*QX,A ;AH=PR+QI

SUB *PX,*QX-,B ;BH=PR-QI

STH A,ASM,*PX+ ; PRV=（PI+QI）/2

SUB *PX,*QX,A ;AH=PR-QI

ST B,*QX ;QRV=（PI-QI）/2

||LD *QX+,B ;BH=QR very important, " BH= {Re[x(2)]- Re[x(6)]} ",=1 ST A,*PX ; PIV=(PI-QR)/2

||ADD *PX+0%,A ; AH=PI+QR

ST A,*QX+0% ;QIV=（PI+QR）/2

||LD *PX,A ;AH=PR

Stage2end:

* * * Stage 3 through Stage logN * * *

STM #K_TWID_TBL_SIZE,BK ;为旋转因子表格的大小值