.TITLE FFT ;See FFT.PAS, and FFTEST.MAC ; USE OF THE .FFT MEGAMACRO WITH PASCAL AND FORTRAN ;Subroutine FFT below is for use with Oregon Software PASCAL, FORTRAN, ;and other languages, which use R5 to pass a list of argument addresses. .GLOBL FFT,.FFT .MCALL .COS,.SIN,.FFT FFT: .FFT 2(R5),@4(R5),@6(R5) RETURN ;Example: do forward transform on complex array Z of 1024 = 2**10 points: ;FORTRAN: COMPLEX Z(1024) ; CALL FFT(Z,10,-1) ;PASCAL: TYPE COMPLEX = ARRAY [1..2] OF REAL; ; TYPE ARRAYCOMPLEX = ARRAY [0..1023] OF COMPLEX; ; VAR Z: ARRAYCOMPLEX; ; VAR LAYER, SIGN: INTEGER; ; ... ; PROCEDURE FFT(VAR Z:ARRAYCOMPLEX; VAR LAYER, SIGN:INTEGER); NONPASCAL; ; ... ; LAYER:= 10; {LAYER AND SIGN MUST BE VARIABLES} ; SIGN:= -1; ; FFT (Z, LAYER, SIGN); A0=R0 A1=R1 A2=R2 A3=R3 .MACRO E2,A,B ;MACRO IS EQUIV. TO FUNCTION E2: B:= 2**A, A IS A REG. ASL A ;MULTIPLY OFFSET BY 2 MOV E2TAB(A),B ;DO TABLE LOOKUP, PUT IT IN B .ENDM .MACRO RECPLX,INDEX ;RETRIEVES REAL AND IMAG.PARTS FROM MEMORY. ;A0 = REAL(X[INDEX]); A1 = IMAG(X[INDEX]) MOV INDEX,R3 ;R3 = I; ;CAUTION: IT USES R3 DEC R3 ;DECREMENT I SINCE ARRAY STARTS AT 1,NOT 0 ASH #3,R3 ;MULTIPLY OFFSET BY 8 ADD Z,R3 ;ADD ADDRESS OF Z TO OFFSET LDF (R3)+,A0 ;REAL PART, AUTOINCREMENTED BY 4 LDF (R3),A1 ;IMAGINARY PART .ENDM .MACRO LDCPLX,INDEX ;PERFORMS THE OPPOSITE OF RECPLX: MOV INDEX,R3 ;Z[INDEX] = REAL = A0; Z[INDEX] = IMAG = A1 DEC R3 ASH #3,R3 ADD Z,R3 STF A0,(R3)+ ;LOAD REAL PART STF A1,(R3) ;LOAD IMAGINARY .ENDM SIGN: .BLKW 1 N: .BLKW 1 Z: .BLKW 1 .FFT: MOV R0,-(SP) MOV R1,-(SP) MOV R2,-(SP) MOV R3,-(SP) MOV R4,-(SP) MOV R5,-(SP) STF A0,-(SP) STF A1,-(SP) STF A2,-(SP) STF A3,-(SP) MOV N,R0 ;R0=N E2 R0,LX ;LX := E2(N) LDCIF SIGN,A3 MULF PI2,A3 LDCIF LX,A2 DIVF A2,A3 STF A3,FLXPI2 ;FLXPI2 := 6.28318530*SIGN/LX ;IN ORDER TO EVALUATE INT[I], WE'LL USE E2(N - I) ;FOR LAYER= 1 TO N DO MOV #1,LAYER LOOP4: MOV LAYER,R4 DEC R4 ;R2 = LAYER - 1 E2 R4,NBLOCK ;NBLOCK := E2(LAYER - 1) MOV LX,R1 CLR R0 DIV NBLOCK,R0 ;R0 GETS QUOTIENT, R1 GETS REMAINDER MOV R0,LBLOCK ;R0 = LBLOCK := LX DIV NBLOCK ASR R0 MOV R0,LBHALF ;LBHALF := LBLOCK DIV 2 CLR NW ;NW := 0 ;FOR IBLOCK:= 1 TO NBLOCK BEGIN LOOP3 MOV #1,IBLOCK LOOP3: MOV IBLOCK,R3 DEC R3 MUL LBLOCK,R3 MOV R3,LSTART ;LSTART := LBLOCK*(IBLOCK - 1) LDCIF NW,A0 MULF FLXPI2,A0 ;A0 = A1 = FLXPI2*NW .COS A0,W1 .SIN A0,W2 ;FOR I:= 1 TO LBHALF BEGIN LOOP2 MOV #1,I LOOP2: MOV I,J ADD LSTART,J ;J := I + LSTART MOV J,K ADD LBHALF,K ;K := J + LBHALF RECPLX K ;RETRIEVE COMPLEX NUMBER FROM ARRAY STF A0,A2 ;A2 = X[K,1] STF A1,A3 ;A3 = X[K,2] MULF W1,A0 MULF W2,A1 SUBF A1,A0 STF A0,Q1 ;Q[1] := X[K,1]*W[1] - X[K,2]*W[1] LDF A2,A0 MULF W2,A0 LDF A3,A1 MULF W1,A1 ADDF A0,A1 STF A1,Q2 ;Q[2] := X[K,1]*W[2] + X[K,2]*W[1] RECPLX J ;RETRIEVE COMPLEX NUMBERS FROM ARRAY STF A0,A2 ;NOW LET A2 = X[J,1] STF A1,A3 ; A3 = X[J,2] SUBF Q1,A0 ;X[K,1] := X[J,1] - Q[1] SUBF Q2,A1 ;X[K,2] := X[J,2] - Q[2] LDCPLX K ;PUT THESE IN THE ARRAY LDF A2,A0 ADDF Q1,A0 ;X[J,1] := X[J,1] + Q[1] LDF A3,A1 ADDF Q2,A1 ;X[J,2] := X[J,2] + Q[2] LDCPLX J ;PUT THESE IN THE ARRAY INC I ;I := I + 1 CMP I,LBHALF ;IS I = LBHALF? BGT CONT2 JMP LOOP2 ;TOO FAR TO BRANCH CONT2: CMP N,#1 ;IS N > 1 ? BGT NGRE1 ;IF SO BRANCH TO N > 1 JMP CONT3 ;OTHERWISE CONTINUE NGRE1: MOV #2,I ;I := 2; CALL WHILE CONT3: INC IBLOCK ;IBLOCK := IBLOCK + 1 CMP IBLOCK,NBLOCK ;IS IBLOCK = NBLOCK? BGT CONT5 JMP LOOP3 CONT5: INC LAYER ;LAYER := LAYER + 1 CMP LAYER,N ;UNTIL LAYER > N BGT CONT6 JMP LOOP4 CONT6: CLR NW ;FOR K= 1 TO LX BEGIN LOOP5 MOV #1,K LOOP5: MOV #1,R4 SUB K,R4 ADD NW,R4 ;R4 = NW + 1 - K TST R4 ;IF (NW + 1 - K > 0) THEN BGT SWAP JMP CONT7 ;OTHERWISE CONTINUE SWAP: MOV NW,R2 INC R2 ;[NW + 1] = R2 RECPLX R2 STF A0,A2 STF A1,A3 ;HOLD := X[NW + 1] RECPLX K LDCPLX R2 ;X[NW + 1] := X[K] LDF A2,A0 LDF A3,A1 LDCPLX K ;X[K] := HOLD CONT7: MOV #1,I CALL WHILE ;DO THE BIT MANIPULATION INC K ;K = K + 1 CMP K,LX ;IS K = LX? BGT CONT8 JMP LOOP5 CONT8: CMP SIGN,#-1 ;IS SIGN = -1? BEQ FOWARD BR DONE ;OTHERWISE WE ARE DONE FOWARD: MOV N,R4 ;R4 = N E2 R4,R0 ;R0 = E2(N) LDCIF R0,A2 LDF ONE,A3 DIVF A2,A3 ;A3 = SCALE:= 1.0/E2(N) MOV #1,I ;START LOOP6 LOOP6: RECPLX I ;GET X[I,1] AND X[I,2] MULF A3,A0 ;MULTIPLY BOTH PARTS MULF A3,A1 ;BY THE SCALE LDCPLX I ;LOAD NEW NUMBERS INTO ARRAY INC I CMP I,LX BLE LP6 BR DONE LP6: JMP LOOP6 DONE: LDF (SP)+,A3 LDF (SP)+,A2 LDF (SP)+,A1 LDF (SP)+,A0 MOV (SP)+,R5 MOV (SP)+,R4 MOV (SP)+,R3 MOV (SP)+,R2 MOV (SP)+,R1 MOV (SP)+,R0 RETURN ;SUBROUTINE WHILE REPLACES THESE LINES OF PASCAL. IT USES R0 AND R1. ; FLAG := FALSE; ; WHILE ((I <= N) AND NOT FLAG) DO ; BEGIN ; II := I; ; IF (NW < INT[I]) ; THEN FLAG := TRUE ; ELSE NW := NW - INT[I]; ; I := I + 1 ; END; ; NW := NW + INT[II]; WHILE: CLR FLAG ;FLAG := FALSE; DOWH: CMP I,N BLE NEXTBO ;GO TO NEXT BOOLEAN IF THIS IS TRUE JMP CONT4 ;IF IT'S FALSE, SKIP LOOP NEXTBO: TST FLAG BEQ WHILE1 ;IF NOT FALSE, DO WHILE LOOP JMP CONT4 WHILE1: MOV N,R0 SUB I,R0 ;R0 = N - I E2 R0,INT ;INT := E2(N - I) MOV I,II ;II := I; CMP NW,INT ;IF NW < 2*(N-I) BLT INCFL SUB INT,NW ;MAKE ITH BIT ZERO JMP INCRI INCFL: INC FLAG INCRI: INC I JMP DOWH CONT4: MOV N,R0 ;MAKE ITH BIT ONE SUB II,R0 E2 R0,INT ADD INT,NW RETURN E2TAB: .WORD 1,2,4,8.,16.,32.,64.,128.,256.,512.,1024.,2048.,4096.,8192. .WORD 16384.,32768.;THIS WAS THE TABLE OF POWERS OF 2 PI2: .FLT2 6.28318530 ONE: .FLT2 1.0 FLXPI2: .BLKW 2 LAYER: .BLKW 1 NBLOCK: .BLKW 1 LBLOCK: .BLKW 1 LBHALF: .BLKW 1 NW: .BLKW 1 LX: .BLKW 1 IBLOCK: .BLKW 1 LSTART: .BLKW 1 W1: .BLKW 2 W2: .BLKW 2 I: .BLKW 1 J: .BLKW 1 K: .BLKW 1 Q1: .BLKW 2 Q2: .BLKW 2 FLAG: .BLKW 1 II: .BLKW 1 INT: .BLKW 1 .END