Trailing-Edge
-
PDP-10 Archives
-
decuslib10-04
-
43,50336/fft.mac
There are no other files named fft.mac in the archive.
TITLE FFT - FAST FOURIER TRANSFORM SUBROUTINE
SUBTTL DEFINITIONS
;A.R. BALDOCK. -ENGINEERING SCHOOL.
;UNIVERSITY OF WESTERN AUSTRALIA.
REPEAT 0,< FAST FOURIER TRANSFORM SUBROUTINE WAS WRITTEN
TO OPTIMIZE THE PROCEDURES AND SELECT FAST ALGORITHMS
FOR SUBSECTIONS.
THE SUBROUTINE IS DESIGNED TO BE CALLED FROM A FORTRAN
MAINLINE, WITH THE FOLLOWING CALLING SEQUENCE:
CALL FFT(SIGN,IGAMMA,XR,XI)
WHERE SIGN IS THE SIGN OF THE EXPONENTIAL FUNCTION
IGAMMA IS THE POWER OF TWO NUMBER OF SAMPLES
(NSAMPLES=2^IGAMMA)
XR IS THE REAL VECTOR TO BE TRANSFORMED
XI IS THE IMAGINARY VECTOR TO BE TRANSFORMED.
*********NOTE*************
THE ORIGINAL VECTORS ARE DESTROYED BY THE TRANSFORM
BECAUSE IT IS DONE "IN PLACE" TO CONSERVE
CORE REQUIREMENTS.
>
;ACCUMULATOR ASSIGNMENTS
FLAG=0 ;ALWAYS
T1=1
T2=2
T3=3
T4=4
I=5
J=6
K=7
L=10
;FLOATING POINT ACCUMULATORS
FA=11
FB=12
FC=13
SP=17 ;STACK POINTER
OPDEF PJRST [JRST] ;RETURNS WITH A POPJ
ENTRY FFT
DEFINE LOTRAN(A)<
HRRZ T2,K
ADD T2,A
HRRZ T3,J
ADD T3,A
MOVE FA,(T2)
FADR FA,(T3)
MOVE FB,(T3)
FSBRM FB,(T2)
MOVEM FA,(T3)
>
DEFINE HITRAN(A,B,C)<
HRRZ T1,J
ADD T1,A
MOVE FA,(T1)
FADR FA,B
MOVEM FA,(C)
FSB FA,B
FSBR FA,B
MOVEM FA,(T1)
>
DEFINE SHUF(A)<
HRRZ T2,I
ADD T2,A
HRRZ T3,J
ADD T3,A
MOVE FA,(T3)
FDVR FA,FN
MOVE FB,(T2)
FDVR FB,FN
MOVEM FB,(T3)
MOVEM FA,(T2)
>
SUBTTL - SET UP PROCEDURES
FFT: CAIA ;STORE JSA ARGS HERE (SKIP ALWAYS FOR F-10)
PUSH SP,16 ;SAVE P.C. ON STACK FOR F-40
MOVE SAVWD ;PRESERVE
BLT SAV16 ;IMPORTANT AC'S
MOVE T2,@1(16) ;IGAMMA
CAIGE T2,1 ;DINKUM ?
JRST FFEXIT ;NO JUST FOOLING
MOVEI T3,2 ;MY WORD- SO QUIT STALLING
ASH T3,-1(T2) ;HOW MANY SAMPLES
MOVEM T3,N ;SAVE IT
HRRZ T2,2(16) ;ADDRESS OF XR
MOVEM T2,XR ;SAVE MAINLINE ADDRESS
HRRZ T2,3(16) ;ADDRESS OF "XI"
MOVEM T2,XI ;SAVE IT ALSO
PUSH SP,T3 ;ARGS ON STACK
PUSHJ SP,FLOAT ;FLOAT IT
MOVE FA,(SP) ;GRAB FLOATED VERSION
MOVEM FA,FN ;SAVE
MOVSI FC,(1.0) ;BUT NOT IF
LDB T2,[POINT 7,@(16),6] ;GET SIGN
SUBI T2,54 ;REMOVE OFFSET
MOVNS T2 ;MAKE SIGNED UNIT
MOVM T3,T2 ;CHECK VALIDITY
CAIE T3,1 ;OF SIGN
JRST NOTSGN ;GRUMBLE
SKIPL T2 ;INVERSE
MOVEM FC,FN ;MAKE UNITY IN THIS CASE
PUSH SP,T2
PUSHJ SP,FLOAT ;FLOAT IT
POP SP,FB ;RESTORE SIGNED FLOATING UNITY
POP SP,FA ;RECOVER ORIGINAL FN
FMPR FB,TWOPI ;*(2.*PI)
FDVR FB,FA ;/N
MOVEM FB,TPN ;SAVE
MOVEI L,1 ;CURRENT ORDER
MOVE T1,N
ASH T1,-1 ;COMPUTE KS=N/2
MOVEM T1,KS
SUBTTL - ROUTINE TO PERFORM ORDER 2 TRANSFORM
MOVN I,KS ;SET UP THE INDEX
FWLOOP: MOVE J,T1 ;COMPUTE THE
ADD J,I ;THE J INDEX
MOVE K,T1 ;AND THE
ASH K,1 ;K INDEX
ADD K,I
LOTRAN(XR) ;PERFORM ORDER 2 REAL TRANSFORM
LOTRAN(XI) ;AND ORDER 2 IMAGINARY TRANSFORM
AOJL I,FWLOOP ;AND LOOP TILL ALL DONE
SUBTTL - MAIN LOOP FOR HIGHER ORDER TRANSFORMS.
FINTST: CAMN L,@1(16) ;DOES L.EQ.IGAMMA ?
JRST REORDR ;YES, SO JUST REORDER AND EXIT
AOS L ;NO, UPDATE TO HIGHER ORDER
SETZM M ;ENSURE IT IS CLEAR FIRST
MOVE T1,KS ;ADJUST
ASH T1,-1 ;KS
MOVEM T1,KS
MOVE T2,@1(16) ;IGAMMA
SUB T2,L
MOVE I,N ;FOR AN INDEX
SOS I
MOVEI T4,1 ;COMPUTE
ASH T4,(T2) ;A DIVISOR
MAINLP: MOVN K,KS ;INDEX FOR OUTER LOOP
MNLOP1: MOVE T1,M ;OLD VALUE
MOVE T2,I ;COMPUTE
IDIV T2,T4 ;A NEW VALUE OF M
CAMN T2,M ;SAME AS OLD ?
JRST SAMARG ;YES, SKIP OVER
MOVEM T2,M ;NO, SAVE NEW VALUE
PUSH SP,M
PUSHJ SP,BITSWP ;SWAP BITS
PUSHJ SP,FLOAT ;AND FLOAT IT
POP SP,T1 ;RETREIVE
FMPR T1,TPN ;COMPUTE PHI=TPN*K
PUSH SP,T1
PUSHJ SP,GCOS ;COMPUTE COS PHI
PUSH SP,T1
PUSHJ SP,GSIN ;AND SIN PHI
SUB SP,[XWD 2,2] ;PREVENT STACK CREEPAGE
;RETURN WITH COS(PHI) AND SIN(PHI) ON THE STACK
SUBTTL - MAIN INNER LOOP FOR HIGHER ORDER TRANSFORMS
SAMARG: MOVE T2,I
ADD T2,XR ;XR(I)
MOVE T3,I
ADD T3,XI ;XI(I)
MOVN FA,2(SP) ;SINE
FMPR FA,(T3) ;SINE*XI(I)
MOVE FB,1(SP) ;COS
FMPR FB,(T2) ;COS*XR(I)
FADR FB,FA ;GR
MOVE FA,2(SP) ;SINE
FMPR FA,(T2) ;SINE*(XR(I)
MOVE FC,1(SP)
FMPR FC,(T3)
FADR FC,FA
MOVN J,KS
ADD J,I
HITRAN(XR,FB,T2) ;PERFORM HI ORDER REAL TRANSFORM
HITRAN(XI,FC,T3) ;AND HI ORDER IMAGINARY TRANSFORM
SOS I ;DECREMENT "I"
AOJL K,MNLOP1 ;IF MORE TO DO LOOP ROUND
MOVEI T1,1 ;SEE IF
ADD T1,I ;"I" WAS
CAMG T1,KS ;EQUAL TO "KS"
JRST FINTST ;YES, SEE IF WE'RE THRU YET
SUB I,KS ;NOT YET, SO RESET "I"
JRST MAINLP ;AND KEEP GOING
SUBTTL - REORDERING PROCEDURE
REPEAT 0,< THIS SECTION REORDERS THE SHUFFLED DATA
AND IF NECESSARY NORMALISES IT.
>
REORDR: MOVN I,N ;SET
HRLZS I ;UP INDEX AND COUNTER
MOVSI T4,(1.0) ;COMPUTE THE RECIPROCAL
FDVR T4,FN ;ONCE ONLY
REORD1: PUSH SP,I
PUSHJ SP,BITSWP ;MAKE BIT SWAPPED WORD
POP SP,J ;GRAB IT
MOVNI T1,(I) ;TEST
ADD T1,J ;FOR COMPARISON
JUMPE T1,REORD2 ;AND HANDLE
JUMPL T1,REORD3 ;ACCORDINGLY
SHUF(XR) ;SHUFFLE REAL DATA
SHUF(XI) ;AND IMAGINARY DATA
JRST REORD3 ;CKECK IF DONE
REORD2: HRRZ T2,I
ADD T2,XR
HRRZ T3,J
ADD T3,XI
MOVE FA,T4 ;1.0/RN
FMPRM FA,(T2) ;ONLY NORMALIZE
MOVE FB,T4
FMPRM FB,(T3) ;IMAGINARIES ALSO
REORD3: AOBJN I,REORD1 ;LOOP ROUND TILL DONE
FFEXIT: MOVS SAVWD ;RESTORE
BLT 16 ;IMPORTANT AC'S
POPJ SP, ;AND RETURN TO MAINLINE
SUBTTL - MISCELLANEOUS ROUTINES
;ROUTINE TO REVERSE ORDER BITS OF A NUMBER
BITSWP: MOVE T3,@1(16) ;IGAMMA
MOVE T1,-1(SP) ;DATA FROM STACK
SETZ T2,
LSHC T1,-1
ROT T2,2
SOJG T3,.-2 ;LOOP TILL DONE
ROT T2,-1 ;DON'T OVERDO IT
MOVEM T2,-1(SP) ;SAVE RESULT ON STACK
POPJ SP, ;AND RETURN
;ROUTINE TO FLOAT A SINGLE PRECISION NUMBER
FLOAT: MOVE FA,-1(SP) ;GET ARG FROM STACK
HLRE FB,FA ;COPY HI PART OF FA INTO LOW PART OF FB
HLL FA,FB ;FILL HI PART OF FA WITH SIGN
FSC FA,233 ;FLOAT LOW PART
SKIPGE FA ;FOR NEGATIVES
AOJE FB,TPOPJ ;TWOS COMPLEMENT HI PART
FSC FB,255 ;FLOAT HI PART
FADR FA,FB ;COMBINE
PJRST TPOPJ ;STORE FA ON STACK AND EXIT
;ERROR MESSAGE ROUTINE FOR CARELESS MAINLINES
NOTSGN: OUTSTR [ASCIZ !ILLEGAL SIGN FOR EXPONTIAL FUNCTION
!]
CALLI 12 ;EXIT
SUBTTL - TRIGONOMETRIC ROUTINES
;ROUTINES TO CALCULATE SINE AND COSINE FUNCTIONS OF RADIAN ARGUMENTS
GCOS: MOVE FB,-1(SP) ;GET ARG FROM STACK
FADR FB,PIOV2 ;ADD PI/2.
SKIPA ;AND FALL INTO SINE ROUTINE
GSIN: MOVE FB,-1(SP) ;ARG FROM STACK
MOVEM FB,SAVARG ;SAVE IT
MOVMS FB ;GET MAGNITUDE OF ARG
CAMG FB,[XWD 170000,0] ;SIN(X)=X IF X<2^-9
POPJ SP, ;SO EXIT IMMEDIATELY
FDV FB,PIOV2 ;DIVIDE X BY PI/2.
CAMG FB,ONE ;IS X/(PI/2)LESS THAN 1.0 ?
JRST SIN2 ;YES, IN FIRST QUADRANT ALREADY
MULI FB,400 ;NO, SEPARATE EXPONENT AND FRACTION
LSH FC,-202(FB) ;GET X MODULO 2*PI
TLZ FC,(1B0) ;SUPPRESS ERROR MESSAGES
MOVEI FB,200 ;PREPARE FLOATING FRACTION
ROT FC,3 ;SAVE THREE BITS FOR QUADRANT
LSHC FB,33 ;AG NOW IN THE RANGE (-1,1)
FAD FB,[0.0] ;NORMALIZE THE ARGUMENT
JUMPE FC,SIN2 ;REDUCED TO FIRST QUADRANT IF FC=000
TLCE FC,1000 ;SUBTRACT 1.0 FROM ARG IF
FSB FB,ONE ;BITS ARE 001 OR 011
TLCE FC,3000 ;CHECK FOR SECOND QUADRANT - BITS 001
TLNN FC,3000 ;OR THIRD QUADRANT - BITS 010
MOVNS FB ;THIRD OR FOURTH QUADRANT
SIN2: SKIPGE SAVARG ;CHECK FOR ORIGINAL SIGN
MOVNS FB ;SINE IS AN ODD FUNCTION
MOVEM FB,SAVARG ;STORE QUADRANT ADJUSTED ARGUMENT
FMPR FB,FB ;SQUARE ARG FIRST
MOVE FA,SICST9 ;THE ORDER 9 CONSTANT
FMP FA,FB ;MULTIPLY BY X SQUARED
FAD FA,SICST7 ;ADD SEVENTH ORDER CONSTANT
FMP FA,FB ;BUILD UP
FAD FA,SICST5 ;NESTED SINE
FMP FA,FB ;ROUTINE
FAD FA,SICST3 ;USING A
FMP FA,FB ;TAYLORS SERIES
FAD FA,PIOV2 ;REPRESENTATION
FMPR FA,SAVARG ;MULTIPLY BY THE ARG(SIN(X)HAS ODD POWERS)
TPOPJ: MOVEM FA,-1(SP) ;RETURN RESULT ON THE STACK
POPJ SP, ;AND RETURN
SUBTTL - IMPURE AND DATA AREA
SAVWD: XWD 1,SAV1
TWOPI: 6.28318530
PIOV2: 201622077325
ONE: 1.0
SICST3: 577265210372
SICST5: 175506321276
SICST7: 606315546346
SICST9: 164475536722
SAVARG: 0
N: 0
KS: 0
M: 0
XR: 0
XI: 0
FN: 0
TPN: 0
SAV1: BLOCK 20
SAV16=.-1
END