Google
 

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