;WF3.MAC BY EC 1-21-78
;NSWEEP CONTAINS NUMBER OF SWEEPS
;CT CONTAINS SWEEP NUMBER = 1,2,...,NSWEEP
;NPT CONTAINS NUMBER OF POINTS PER SWEEP
;N CONTAINS NUMBER FOR FFT=2*(NPT)
;PTX = STARTING ADDRESS OF REAL PART OF DATA FOR FFT
;PTY = STARTING ADDRESS OF IMAGINARY PART OF DATA FOR FFT
;RATE CONTAINS NUMBER OF CLOCK COUNTS TO OVERFLOW AT 1MICROSEC/CT
;XPT CONTAINS STARTING ADDRESS OF REAL PART OF DATA
;YPT CONTAINS STARTING ADDRESS OF IMAG PART OF DATA
;SUM = STARTING ADDRESS OF ONGOING SAMPLE SUM
;S1 CONTAINS NEG OF NUMBER OF OVERFLOWS IN ONGOING SAMPLE SUM
;MPT CONTAINS NEG OF LOG2(NPT)

	.MCALL ..V2.., .REGDEF,.PRINT,.TTYIN,.TTYOUT
	.MCALL .EXIT,.CSISPC,.WRITW,.CLOSE,.LOCK,.UNLOCK
	.MCALL .ENTER,.SRESET

	..V2..
	.REGDEF

	TKS=177560
	TKB=177562
	TPS=177564
	TPB=177566
	XDAC=170440	;DAC 0
	YDAC=170442	;DAC 1
	DAC3=170446
	ADSR=170400	;A/D STATUS REG
	ADBR=170402	;A/D BUFFER REG
	CSR=170420	;CLOCK STATUS REG
	BPR=170422	;CLK BUFFER PRESET REG
	CR=15
	LF=12

	

START:	MOV @#60,KEYMON
	MOV #4000,DAC3
	MOV #ADINT,@#400	;A/D DONE INTERRUPT VECTOR
	MOV #200,@#402		;PROCESSOR STATUS DURING A/D INTERRUPT=7
	MOV #OVINT,@#440	;CLOCK OVERFLOW INTERRUPT VECTOR
	MOV #200,@#442		;PS=6 
	MOV #ST2INT,@#444	;ST2 INTERRUPT VECTOR
	MOV #200,@#446		;PS=7
	BIS #5000,ADSR		;SAMPLE CHANNEL 10
;ASK NUMBER OF SWEEPS(NSWEEP),
;NUMBER OF POINTS PER SWEEP (NPT),
;TIME BETWEEN POINTS (RATE)
	.PRINT #Q1
	JSR PC,TINO
	MOV R1,NSWEEP
	.PRINT #Q2
	JSR PC,TINO
	MOV R1,NPT
	.PRINT #Q3
	JSR PC,TINO
	MOV R1,RATE
	MOV RATE,R4
	NEG R4
	MOV R4,BPR	;CLK COUNTS TO ZERO
	MOV #MEAN,MNPT		;POINTER FOR SWEEP MEAN
	MOV NPT,N
	ASL N
;COMPUTE MPT =-LOG2(NPT)
	CLR MPT
	MOV #1,R0
LN:	CMP NPT,R0	;(NPT)-(R0)>0--MUL (R0) BY 2
			;(NPT)-(R0)<=0--STOP--R0 CONTAINS THE SMALLEST
			;POWER OF TWO >= NPT
	BLE CH
	ASL R0
	DEC MPT
	BR LN
CH:	CLR H		;H=10000/NPT
	MOV #10000,R0
	ASH MPT,R0
	MOV R0,H
	NEG R0
	MOV R0,H1	;H1=-H
	.PRINT #WFQ	;ASK WHICH FILTER
	CLR R0
AN4:	.TTYIN
	CMPB R0,#15
	BEQ B4
	MOV R0,WFD	;STORE IN WFD:60=NO FILTER,61=W-S,62=D-S,63=W-R,64=D-R
	BR AN4
B4:	.TTYIN		;CLEAR MONITOR'S KEYBOARD BUFFER
	CMPB R0,#12
	BNE B4
INIT:	CLR ADSR
	CLR CSR
	CLR CT
	CLR KEY
	CLR VDS			;VERTICAL DISPLAY SCALE
	CLR SC		;CLEAR ALL SCALE FACTORS
	CLR S1
	CLR KS
	CLR SF
	CLR SFF
	MOV #SUM,R1		;SET TO 0 INITIALLY
	MOV #RTSUM,R2
	MOV #ITSUM,R3
	MOV #PSSUM,R4
	MOV #PTX,R0
	MOV #PTY,R5
	MOV N,TN
ZERO1:	CLR (R1)+
	CLR (R2)+
	CLR (R3)+
	CLR (R4)+
	CLR (R0)+
	CLR (R5)+
	DEC TN
	BNE ZERO1
	MOV N,R0
	MOV #KSUM,R1
	MOV #KER,R2
ZEK:	CLR (R1)+
	CLR (R2)+
	SOB R0,ZEK
;WAIT FOR KEYBOARD TO SET CLK FOR ST2 START--INITIALLY ONLY
	MOV #KEYINT,@#60
	MOV #200,@#62
	BIS #10000,@#44	;CONSOLE MODE 1 FOR COMMANDS
	BIS #100,TKS	;KYBD INT EN
	.PRINT #MSG	;"B" TO BEGIN ADDING
	.PRINT #MMSG	;"M" TO MONITOR
MWAIT:	MOV #PTX,DISPTR
	JSR PC,DISPLY
	JMP MWAIT
CWAIT:	MOV #KER,DISPTR
	JSR PC,DISPLY
	JMP CWAIT
;DISPLAY ADDRESSES PTX,,,,PTX +NPT AND SUM,...,SUM+NPT
;INTERRUPT DISPLAY WHEN A/D CONVERSION IS DONE
;ZERO MEAN--X(I)=X(I)-MX,WHERE MX=(X(0)+...+X(NPT-1))/NPT
;USE DOUBLE PRECISION TO COMPUTE MX
TRANS:	CLR R2		;HIGH WORD OF MEAN
	CLR R3		;LOW WORD OF MEAN
	MOV NPT,R1	;COUNT
	MOV #PTX,R0	;STARTING ADDRESS OF X(I) IN R0
MN1:	MOV (R0)+,R5	;X(I) IN R5
	SXT R4		;EXTEND SIGN OF X(I) INTO R4
	ADD R5,R3	;ADD X(I) INTO LOW WORD OF MEAN
	ADC R4		;ADD CARRY INTO SIGN EXTENSION OF X(I)
			;R4 CONTAINS EITHER 0 OR 1 OR -1
	ADD R4,R2	;ADD INTO HIGH WORD OF MEAN
	SOB R1,MN1	;ALL POINTS ADDED IN? NO,ADD NEXT
	ASHC MPT,R2	;YES--SUM IN R2,R3,HI IN R2
			;DIVIDE BY NPT--RIGHT SHIFT (R2VR3)BY
			;LOG2(NPT)=-MPT-- MEAN IS IN R3
	MOV R3,@MNPT	;SAVE THE MEAN FOR EACH SWEEP
	INC MNPT
	INC MNPT	;INC POINTER BY 2 FOR NEXT SWEEP
	MOV NPT,R1	;COUNT
	MOV #PTX,R0	
MN2:	SUB R3,(R0)+	;X(I)=X(I)-MEAN
	SOB R1,MN2
	TST BEGIN	;MONITORING? 0 IN BEGIN
	BNE ADOVER	;NO,SUMMONG
	BIC #100000,CSR	;YES,CLEAR ST2 FLAG
	BIS #20000,CSR	;SET ST2 GO EN
	RTI
;FORM KTH PARTIAL SUM OF SAMPLED POINTS
ADOVER:	INC CT		;CT CONTAINS CURRENT SWEEP NUMBER 
	CMP CT,NSWEEP	;CT>NSWEEP,DONE
	BLE SWPCOM	;CT<=NSWEEP,DO SWEEP COMPUTATIONS
	JMP OUT
SWPCOM:	MOV #PTX,R0
	MOV #SUM,R1
	MOV NPT,R4
ADON:	MOV (R0),R2
	MOV (R1),R3
	ASH S1,R2	;SCALE INCOMING VALUES BEFORE ADDING
	ADD R2,R3	;ADD INCOMING SAMPLE VALUE TO ONGOING SAMPLE
			;SUM--A(I)=A(I)+X(I)
	BVC AOK		;CHECK FOR OVERFLOW
	JSR PC,SCALEA	;SCALE ON OVERFLOW
	JMP ADON	;TRY AGAIN
AOK:	MOV R3,(R1)+	;STORE ONGOING SUM IN MEMORY
	TST (R0)+	;INC (R0) FOR NEXT PT
	SOB R4,ADON	;NEXT POINT
	CMPB WFD,#'0	;60 IN WFD?
	BNE DOF		;NO,DO FILTER
	BIC #100000,CSR
	BIS #20000,CSR
	RTI
DOF:	JSR PC,KCAL	;CALCULATE FILTER KERNELS FOR SWEEP JUST 
			;COMPLETED BEFORE NEXT STIMULUS
	MOV N,R4	;SET PTX,PTY ADDRESSES TO 0 AT START OF EACH SWEEP
	MOV #PTX,R0
	MOV #PTY,R1
ZERO:	CLR (R0)+
	CLR (R1)+
	SOB R4,ZERO
;	MOV #100,VCT
;VIEW:	MOV #KER,DISPTR
;	JSR PC,DISPLY
;	DEC VCT
;	BNE VIEW
	CMP (SP)+,(SP)+
	MTPS #0
	BIC #100000,CSR		;CLEAR ST2 FLAG
	BIS #20000,CSR		;ST2 GO EN--ST2 WILL SET (15) ONLY 
			;IF (13) OR (0) IS SET AND (15) IS CLEAR
	JMP CWAIT
OUT:	CLR CSR
	CLR ADSR
	CMP (SP)+,(SP)+
	MTPS #0
INVFT:	.PRINT #MSG2
	.PRINT #M7	;SUM SCALE FACTOR
	MOV S1,R1
	NEG R1
	JSR PC,OCOUT
	.PRINT #M6
	.PRINT #M8	;KER SCALE FACTOR
	MOV SC,R1
	JSR PC,OCOUT
	.PRINT #M6
	.PRINT #M15	;KSUM SCALE FACTOR
	MOV KS,R1
	NEG R1
	JSR PC,OCOUT
	.PRINT #M6
;COMPUTE AVERAGE OF NSWEEPS--A(I)=(S(I))/NSWEEP
	MOV NPT,R4		;COUNT
	MOV #SUM,R0
	MOV #AVG,R1
	MOV NSWEEP,R5
	JSR PC,AVNXT	;DIVIDE BY NSWEEP
	.PRINT #MSG3	;YES,PRINT DONE  MSG
	CMPB WFD,#'0	;FILTERED AVG?
	BNE DOFA	;YES
AVDIS:	MOV #AVG,DISPTR
	JSR PC,DISPLY
	JMP AVDIS
DOFA:	JSR PC,FCAL	;CALCULATE FILTERED SUM
	MOV NPT,R4	
	MOV #PTX,R0
	MOV #FAVG,R1
	MOV NSWEEP,R5
	JSR PC,AVNXT	;DIVIDE BY NSWEEP
	.PRINT #M3	;FAVG COMPUTED
	.PRINT #M5	;LAST FFT SCALE FACTOR
	MOV SFF,R1
	JSR PC,OCOUT
	.PRINT #M6
	.PRINT #M9	;IFT SCALE FACTOR
	MOV SF,R1
	JSR PC,OCOUT
	.PRINT #M6
FAVDIS:	MOV #FAVG,DISPTR
	JSR PC,DISPLY
	JMP FAVDIS

;DISPLAY ADDRESSES PTX,...,PTX+NPT AND SUM,...SUM+NPT
;WRITE ON DISK
;REINITIATE FOR NEXT CURVE
;SUBROUTINE TO DISPLAY (I,(R2(I)),I=0,1,...,(NPT-1)
DISPLY:	MOV NPT,DISCT
	MOV H1,XDAC
NDIS:	MOV @DISPTR,R1
	INC DISPTR
	INC DISPTR
	ASH VDS,R1	;UP OR DOWN THE DISPLAY--NO CHANGE IN ARRAY
			;VALUES IN MEMORY
	ADD #4000,R1	;4000=0 VOLT IN D/A CONVERSION
POINT:	MOV R1,YDAC
	ADD H,XDAC	;INCREMENT X-COORD BY H
	BIS #1,DAC3
	DEC DISCT
	BNE NDIS
	RTS PC

;ST2 INTERRUPT SERVICE ROUTINE
ST2INT:	BIS #112,CSR	;CLK OVFW INT EN,1MHZ,MODE 1
	MOV #PTX,R0	;DO NXT SWEEP--REINITIALIZE FOR SAMPLING
	MOV NPT,R3	;RESTART SAMPLING PT CT
	BIS #5100,ADSR	;AD DONE INT EN,CHANNEL 10
	RTI

;CLOCK OVERFLOW INTERRUPT SERVICE ROUTINE
OVINT:	BIC #200,CSR	;CLEAR CLK OVFLW FLAG
	INC ADSR	;START A/D CONVERSION
	RTI

;A/D DONE INTERRUPT SERVICE ROUTINE
ADINT:	SUB #4000,ADBR	;0 VOLT=4000 IN A/D CONVERSION
	MOV ADBR,(R0)+	;PUT VALUE IN PTX ARRAY
	DEC R3		;SAMPLED ALL PTS?
	BNE RET		;NO
	BIC #100,ADSR	;YES,DISABLE AD INT
	BIC #113,CSR	;STOP CLK,DISABLE CLK OVFW INT AND MODE 1
CADB:	TST ADBR	;CLEAR AD BUFFER
	TSTB ADSR	;(7) RESETS IF ANY WORDS STILL IN AD BUFFER
	BMI CADB	;READ ADBR UNTIL ALL CLEAR
	JMP TRANS	;BETWEEN SWEEP CALCULATIONS
RET:	RTI

;KEYBOARD INTERRUPT SERVICE ROUTINE
KEYINT:	MOVB TKB,KEY
ECHO:	TSTB TPS	;READY TO PRINT?
	BPL ECHO	;NO
	MOVB KEY,TPB	;YES,ECHO
	BIS #100,TKS
	BIC #200,KEY
	CMPB KEY,#'B	;BEGIN SWEEPS?
	BNE MR		;NO
	CLR ADSR
ADB:	TST ADBR	;READ ADBR UNTIL ALL CLEAR
	TSTB ADSR
	BMI ADB
	INC BEGIN
	CMP(SP)+,(SP)+
	MTPS #0
	MOV #60000,CSR	;ST2 INT EN,ST2 GO EN
	JMP CWAIT
MR:	CMPB KEY,#'M	;MONITOR?
	BNE C		;NO
	CLR BEGIN
	CLR ADSR
	MOV #60000,CSR	;YES
	RTI
C:	CMPB KEY,#'C	;CONTRACT SUM?
	BNE XP		;NO
	ASL NPT
	ASR H
	ASR H1
	RTI
XP:	CMPB KEY,#'X	;EXPAND SUM?
	BNE U		;NO
	ASR NPT
	ASL H
	ASL H1
	RTI
U:	CMPB KEY,#'U	;UP DISPLAY?
	BNE D		;NO
	INC VDS
	RTI
D:	CMPB KEY,#'D	;DOWN DISPLAY
	BNE S		;NO
	DEC VDS
	RTI
S:	CMPB KEY,#'S	;SUM?
	BNE A		;NO
	JMP SDIS	;YES
A:	CMPB KEY,#'A	;AVG?
	BNE W
	JMP ADIS
W:	CMPB KEY,#'W	;FILTERED AVERAGE?
	BNE V		;NO
	JMP WDIS	;YES
V:	CMPB KEY,#'V	;RTSUM
	BNE J	;NO
	JMP RTDIS
J:	CMPB KEY,#'J	;ITSUM
	BNE HP
	JMP ITDIS
HP:	CMPB KEY,#'H	;KER
	BNE Q
	JMP KDIS
Q:	CMPB KEY,#'Q	;KSUM?
	BNE O
	JMP KSDIS
O:	CMPB KEY,#'O	;PSSUM
	BNE MN
	JMP PSDIS
MN:	CMPB KEY,#'M	;MEAN
	BNE L
	JMP MDIS
L:	CMPB KEY,#'L	;PLOT?
	BNE Z		;NO
	JMP LDIS
Z:	CMPB KEY,#'Z	;START PLOTTING?
	BNE F		;NO
	CMP (SP)+,(SP)+	;YES
	MTPS #0
	JMP PLOT
F:	CMPB KEY,#'F	;FFT OF AVG? RESULT IN PTX,PTY ADDRESSES
	BNE NV		;NO
	JMP DOFT	;YES
NV:	CMPB KEY,#'N	;INVERSE FFT?
	BNE R		;NO
	JMP DOIFT	;YES
R:	CMPB KEY,#'R	;REAL PART OF FT?	
	BNE IM		;NO
	JMP RDIS	;YES
IM:	CMPB KEY,#'I	;IMAGINARY PART OF FT?
	BNE P		;NO
	JMP IDIS	;YES
P:	CMPB KEY,#'P	;POWER SPECTRUM?
	BNE K		;NO
	JMP PDIS	;YES
K:	CMPB KEY,#'K	;KEEP?
	BNE T		;NO
	CMP (SP)+,(SP)+
	MTPS #0		;ALLOW INTERRUPTS
	JSR PC,RITE
BACK:	WAIT
	BR BACK
T:	CMPB KEY,#'T	;RESTART WITH SAME PARAMETERS?
	BNE G
	JMP GO
G:	CMPB KEY,#'G	;RESTART?
	BNE E		;NO
	JMP GO
E:	CMPB KEY,#'E	;EXIT TO MONITOR?
	BNE NOM		;NO
	BIC #10000,@#44		;CONSOLE MODE 0
	MOV KEYMON,@#60
	MOV #340,@#62
	BIS #100,TKS
	MOV #1000,SP
	.EXIT		;YES
NOM:	RTI		;NO MATCH
SDIS:	CMP (SP)+,(SP)+	;RESTORE STACK POINTER
	MTPS #0		;PS=0
	MOV #SUM,PLTPTR
SDISL:	MOV #SUM,DISPTR
	JSR PC,DISPLY
	JMP SDISL
ADIS:	CMP (SP)+,(SP)+
	MTPS #0
	MOV #AVG,PLTPTR
ADISL:	MOV #AVG,DISPTR
	JSR PC,DISPLY
	JMP ADISL
WDIS:	CMP (SP)+,(SP)+
	MTPS #0
	MOV #FAVG,PLTPTR
WDISL:	MOV #FAVG,DISPTR
	JSR PC,DISPLY
	JMP WDISL
RDIS:	CMP (SP)+,(SP)+
	MTPS #0
	MOV #PTX,PLTPTR
RDISL:	MOV #PTX,DISPTR
	JSR PC,DISPLY
	JMP RDISL
IDIS:	CMP (SP)+,(SP)+
	MTPS #0
	MOV #PTY,PLTPTR
IDISL:	MOV #PTY,DISPTR
	JSR PC,DISPLY
	JMP IDISL
RTDIS:	CMP (SP)+,(SP)+
	MTPS #0
	MOV #RTSUM,PLTPTR
RTDISL:	MOV #RTSUM,DISPTR
	JSR PC,DISPLY
	JMP RTDISL
ITDIS:	CMP (SP)+,(SP)+
	MTPS #0
	MOV #ITSUM,PLTPTR
ITDISL:	MOV #ITSUM,DISPTR
	JSR PC,DISPLY
	JMP ITDISL
KDIS:	CMP (SP)+,(SP)+
	MTPS #0
	MOV #KER,PLTPTR
KDISL:	MOV #KER,DISPTR
	JSR PC,DISPLY
	JMP KDISL
KSDIS:	CMP (SP)+,(SP)+
	MTPS #0
	MOV #KSUM,PLTPTR
KSDISL:	MOV #KSUM,DISPTR
	JSR PC,DISPLY
	JMP KSDISL
PSDIS:	CMP (SP)+,(SP)+
	MTPS #0
	MOV #PSSUM,PLTPTR
PSDISL:	MOV #PSSUM,DISPTR
	JSR PC,DISPLY
	JMP PSDISL
MDIS:	CMP (SP)+,(SP)+
	MTPS #0
	MOV #MEAN,PLTPTR
MDISL:	MOV #MEAN,DISPTR
	JSR PC,DISPLY
	JMP MDISL
LDIS:	CMP (SP)+,(SP)+
	MTPS #0
	MOV PLTPTR,R2
	MOV NPT,R3
	DEC R3
	MOV #0,XDAC
	MOV #4000,YDAC
FIRST:	BIS #1,DAC3	;WAIT FOR "Z" TO START PLOTTING
	BR FIRST
PLOT:	MOV (R2)+,R1
	ASH VDS,R1
	ADD #4000,R1
	MOV #2000,R0	;HOLD
	MOV R1,YDAC
	ADD H,XDAC
PLOTPT:	BIS #1,DAC3
	SOB R0,PLOTPT
	SOB R3,PLOT
	JMP FIRST	;HOLD AT LAST POINT
PDIS:	CMP(SP)+,(SP)+
	MTPS #0
;COMPUTE X(I)^2+Y(I)^2 AND PUT IN PS ARRAY
	MOV N,R5
	MOV #PTX,R0
	MOV #PTY,R1
	MOV #PS,R4
	CLR PSF		;POWER SCALE FACTOR
PSLP:	MOV (R0),R2	;X(I) IN R2
	MUL R2,R2	;X(I)^2 IN R2VR3--HI IN R2
	ASHC #4,R2	;LS4 FOR MUL WITH PT AFTER (12)
	BVC P1		;OVERFLOW?
	JSR PC,PSCALE	;YES,SCALE
	JMP PSLP
P1:	MOV R2,TX2	;X(I)^2 IN TX2
	MOV (R1),R2	;Y(I) IN R2
	MUL R2,R2	;Y(I)^2 IN R2VR3,HI IN R2
	ASHC #4,R2
	BVC P2
	JSR PC,PSCALE
	JMP PSLP
P2:	ADD TX2,R2	;X(I)^2+Y(I)^2 IN R2
	BVC P3
	JSR PC,PSCALE
	JMP PSLP
P3:	MOV R2,(R4)+	;STORE IN PS ARRAY
	CMP (R0)+,(R1)+	;INC FOR NXT PT
	SOB R5,PSLP	;THRU ALL PTS?
	MOV #PS,PLTPTR
PDISL:	MOV #PS,DISPTR
	JSR PC,DISPLY
	JMP PDISL
GO:	CMP (SP)+,(SP)+
	MTPS #0
	BIC #10000,@#44		;CONSOLE MODE 0
	MOV KEYMON,@#60
	MOV #340,@#62
	BIS #100,TKS
	CLR BEGIN
	MOV #1000,SP
	CMPB KEY,#'T
	BNE G1
	JMP INIT
G1:	JMP START
DOFT:	MOV #AVG,R0
	MOV #PTX,R1
	MOV #PTY,R2
	MOV NPT,R3
LP:	MOV (R0)+,(R1)+	;TRANSFER AVG,0 TO PTX,PTY ADDRESSES
	CLR (R2)+
	SOB R3,LP
	CLR DIR		;FORWARD FT
	CMP (SP)+,(SP)+
	MTPS #0
	JSR PC,FFT
	MOV #PTX,PLTPTR
	JMP RDISL
DOIFT:	MOV #1,DIR	;INVERSE FT
	CMP (SP)+,(SP)+
	MTPS #0
	JSR PC,FFT
	MOV #PTX,PLTPTR
	JMP RDISL

;SUBROUTINE TO SCALE ON OVERFLOW IN ONGOING KERNEL SUM
KSCALE:	MOV R1,-(SP)
	MOV R4,-(SP)
	MOV N,R4
	MOV #KSUM,R1
KBAC:	ASR (R1)+
	SOB R4,KBAC
	DEC KS
	MOV (SP)+,R4
	MOV (SP)+,R1
	RTS PC

;SUBROUTINE TO SCALE ON OVERFLOW IN ONGOING SAMPLE SUM
SCALEA:	MOV R1,-(SP)	;SAVE (R1)
	MOV R4,-(SP)	;SAVE (R4)
	MOV NPT,R4
	MOV #SUM,R1
BAC:	ASR (R1)+	;DIVIDE BY 2
	SOB R4,BAC
	DEC S1		;DECREASE S1 BY 1 EACH TIME SCALING IS DONE
	MOV (SP)+,R4	;RESTORE (R4)
	MOV (SP)+,R1	;RESTORE (R1)
	RTS PC

;SUBROUTINE TO DIVIDE ARRAY POINTED TO BY R0 BY NUMBER IN R5
;QUOTIENT STORED IN ARRAY POINTED TO BY R1
;R4 CONTAINS NUMBER OF POINTS IN ARRAY
AVNXT:	MOV (R0)+,R3
	SXT R2
	DIV R5,R2	;(R2VR3)/R5--QUOTIENT IN R2,REMAINDER IN R3
	BVC OK
	MOV R0,TEMP
	.PRINT #E1
	MOV TEMP,R0
OK:	MOV R2,(R1)+
	SOB R4,AVNXT
	RTS PC


;SUBROUTINE TO SCALE DURING POWER SPECTRUM CALCULATION
PSCALE:	MOV R0,-(SP)
	MOV R1,-(SP)
	MOV R2,-(SP)
	MOV R3,-(SP)
	MOV #PTX,R0
	MOV #PTY,R1
	MOV #PS,R2
	MOV N,R3
SL1:	ASR (R0)+	;X(I)/2
	ASR (R1)+	;Y(I)/2
	ASR (R2)	;(X(I)^2+Y(I)^2)/2
	ASR (R2)+	;(X(I)^2+Y(I)^2)/4
	SOB R3,SL1
	INC PSF
	MOV (SP)+,R3
	MOV (SP)+,R2
	MOV (SP)+,R1
	MOV (SP)+,R0
	RTS PC
;SUBROUTINE TO CALCULATE FILTER KERNELS
;TAKE FOURIER TRANSFORM OF INCOMING DATA
KCAL:	CLR DIR		;FORWARD FFT
	JSR PC,FFT	;TRANSFORM N=2*NPT ARRAY IN PTX,PTY ADDRESSES
;Z(I)=X(I)+J*Y(I),X(0),...,X(NPT-1),0,...0
		;0,...0,0...,0
;RETURNS WITH Z^(I)=X^(I)+J*Y^(I) IN PTX,PTY ADDRESSES


;FORM KTH PARTIAL SUM OF TRANSFORMS
;S^(I)=S^(I)+Z^(I)/2^S,WHERE S=CURRENT # OVERFLOWS,I=0,1,2,...,N-1
;U(I)=U(I)+X(I)/2^S
;V(I)=V(I)+Y(I)/2^S

;FORM KTH PARTIAL SUM OF POWER SPECTRA
;W(I)=W(I)+ABS(Z^(I))^2/2^2*S=W(I)+(X^(I)^2+Y^(I)^2)/2^2*S

;SCALE Z,S,AND W ARRAYS ON OVERFLOW

;FIRST DIVIDE X(I) AND Y(I) BY 2^S,I=0,1,...,N-1
;S=0 INITIALLY
;S INCREASES BY 1 EACH TIME SCALING IS DONE
	TST SC		;SC=0,NO DIVISION--START ADDING INTO ON-
			;GOING SUM
	BEQ ALOOP
	MOV N,R5
	MOV #PTX,R0
	MOV #PTY,R1
SCA:	MOV SC,R4
SCL:	ASR (R0)
	ASR (R1)
	SOB R4,SCL	;SHIFTED RT S TIMES?--NO,SHIFT AGAIN
	CMP (R0)+,(R1)+	;SCALE NEXT POINT
	SOB R5,SCA	;SCALED ALL POINTS?--NO,SCALE NXT PT
			;YES--START ADDING INTO ONGOING SUMS
ALOOP:	MOV N,TN	;CT WITH TN DURING CALC
	MOV #PTX,XPT	;REINITIALIZE POINTERS
	MOV #PTY,YPT
	MOV #RTSUM,PTU
	MOV #ITSUM,PTV
	MOV #PSSUM,PTP
	MOV #KER,PTK
NXLP:	MOV @XPT,TX1	;X(I) IN TX1--MADE IT THRU THE SCALING TO 
			;NEXT I
	MOV @YPT,TY1	;Y(I) IN TY1
	MOV @PTU,TU	;U(I) IN TU
	MOV @PTV,TV	;V(I) IN TV
	MOV @PTP,TW	;W(I) IN TW
;W(I)=W(I)+X(I)^2+Y(I)^2
	MOV TX1,R2	;X(I) IN R2
	MUL R2,R2	;X(I)^2 IN R2,R3--HI IN R2
	ASHC #4,R2	;LEFT SHIFT 4 FOR FRACTIONAL MUL WITH POINT
			;AFTER BIT 12
	BVC V1		;CHECK FOR OVERFLOW
	JSR PC,SCALE1	;SHIFTS ALL ARRAYS RIGHT 1 BIT
	JMP NXLP	;TRY AGAIN
V1:	MOV TY1,R4	;Y(I) IN R4
	MUL R4,R4	;Y(I)^2 IN R4,R5--HI IN R4
	ASHC #4,R4
	BVC V2
	JSR PC,SCALE1
	JMP NXLP
V2:	ADD R2,R4	;X(I)^2+Y(I)^2 IN R4
	BVC V3
	JSR PC,SCALE1
	JMP NXLP
V3:	ADD R4,TW	;W(I)+(R4) IN TW
	BVC V4
	JSR PC,SCALE1
	JMP NXLP
V4:	ADD TX1,TU	;U(I)+X(I) IN TU
	BVC V5
	JSR PC,SCALE1
	JMP NXLP
V5:	ADD TY1,TV	;V(I)+Y(I) IN TV
	BVC V6
	JSR PC,SCALE1
	JMP NXLP
;FOR CT=2,...NSWEEP,COMPUTE KERNEL IF WFD CONTAINS 3 OR 4. IF WFD CONTAINS 1 OR 2,NOT UNTIL CT=NSWEEP
;FOR CT=1,KERNEL=1--MOVE TW,TU,TV TO ADDRESSES AND DO NEXT PT,
;THEN NEXT SWEEP
V6:	CMP CT,#1	;CT=1?
	BNE DECIDE	;NO
NOKER:	MOV TU,R2	;SET KER=10000 (1 WITH PT AFTER BIT 12)
	MOV TV,R0
	MOV #10000,R4	;KER IN R4
	JMP V13
DECIDE:	CMPB WFD,#'3	;WFD>=63,KER
	BGE V61
	CMP CT,NSWEEP	;CT=NSWEEP	;CT=NSWEEP,KER
	BEQ V61
	JMP NOKER
V61:	MOV TU,R2	;U(I) IN R2
	MUL R2,R2	;U(I)^2 IN R2,R3--HI IN R2
	ASHC #4,R2
	BVC V7
	JSR PC,SCALE1
	JMP NXLP
V7:	MOV TV,R4	;V(I) IN R4
	MUL R4,R4	;V(I)^2 IN R4,R5--HI IN R4
	ASHC #4,R4
	BVC V8
	JSR PC,SCALE1
	JMP NXLP
V8:	ADD R2,R4	;U(I)^2+V(I)^2 IN R4
	BVC V9
	JSR PC,SCALE1
	JMP NXLP
V9:	CMPB WFD,#'1	;W-S?
	BNE WF2
	JMP WSR
WF2:	CMPB WFD,#'2	;D-S?
	BNE WF3		;NO
	JMP DSR	;YES
WF3:	CMPB WFD,#'3	;W-R?
	BNE WF4
	JMP WSR	;YES
WF4:	CMPB WFD,#'4	;D-R?
	BNE WF5		;NO
	JMP DSR
WF5:	.PRINT #E2	;NO MATCH
	JMP GO
DSR:	TST R4		;ABS(S(I))^2=0?
	BNE V91		;NO,CALC KER
	JMP NOKER
V91:	MOV R4,R1	;ABS(S(I))^2 IN R1
	SUB TW,R4	;(ABS(S(I))^2-W(I) IN R4
	BVC V10
	JSR PC,SCALE1
	JMP NXLP
V10:	CLR R5
	ASHC #-4,R4	;RT SHIFT FOR DIVISION WITH PT AFTER BIT 12
			;IN NUMERATOR AND DENOMINATOR
	DIV R1,R4	;(R4VR5)/(R1) = 
			;ABS(S(I))^2-W(I)/ABS(S(I))^2 IN R4 WITH 
			;PT AFTER BIT 12
			;REMAINDER IN R5
	BVC GOON	;OVERFLOW CLEAR--OK
	JMP NOKER	;SET KER=1
;WIENER KERNEL = (R4)*CT/(CT-1),CT=2,...,NSWEEP
;COMPUTE CT/(CT-1)
GOON:	MOV CT,R2	;2<=K<512
	MOV R2,R0
	DEC R0		;(K-1) IN R0
	CLR R3
	ASHC #-4,R2	;PT AFTER BIT 6--6 L S +10 R S = 4 R S 
	ASH #6,R0
	DIV R0,R2	;K/(K-1) IN R2 WITH PT AFTER BIT 6
	ASH #6,R2	;LEFT SHIFT 3 TO PUT PT AFTER BIT 12
	MUL R2,R4	;KERNEL IN R4,R5--HI IN R4
	ASHC #4,R4	;LS 4 FOR FRACT MUL WITH PT AFTER BIT 12
	BVC V11
	JSR PC,SCALE1
	JMP NXLP
	JMP V11		;D-SR COMPUTED
WSR:	TST TW		;TW=0?
	BNE WV91	;NO,CALC KER
	JMP NOKER
WV91:	SUB TW,R4	;(ABS(S(I))^2-TW IN R4
	BVC WV10
	JSR PC,SCALE1
	JMP NXLP
WV10:	CLR R5
	ASHC#-4,R4	;RT SHIFT 4 FOR DIVISION WITH PT AFTER BIT 12 IN NUMERATOR AND DENOMINATOR
	DIV TW,R4	;(R4VR5)/TW=ABS(S(I))^2-TW/TW IN R4 WITH PT AFTER BIT 12,REMAINDER IN R5
	BVC WGOON
	JMP NOKER	;SET KER=1
;WIENER KERNEL2 = (R4)/(CT-1),CT=2,...,NSWEEP
WGOON:	MOV CT,R1	;CT>=2
	DEC R1		;CT-1 IN R1
	MOV R4,R5	;(ABS(S(I))^2-TW)/TW IN R5--PT AFTER BIT 12
	SXT R4		;SIGN EXTEND INTO R4
	DIV R1,R4	;(R4VR5)/(R1)--KER2=QUOTIENT IN R4 WITH PT AFTER BIT 12
			;(HAVE DIVIDED BY AN INTEGER IN R1 WITH PT AFTER BIT 0),REMAINDER IN R5
	BVC  V11
	.PRINT #E1
;S(I) = (KERNEL)*S(I)
;KERNEL IS REAL
;U(I)=KER*U(I), V(I)=KER*V(I)
V11:	MOV TU,R2	;U(I) IN R2
	MUL R4,R2	;KER*(R2) IN R2,R3--HI IN R2
	ASHC #4,R2
	BVC V12
	JSR PC,SCALE1
	JMP NXLP
V12:	MOV TV,R0	;V(I) IN R0
	MUL R4,R0	;KER*(R0) IN R0,R1--HI IN R0
	ASHC #4,R0
	BVC V13
	JSR PC, SCALE1
	JMP NXLP
;PUT NUMBERS IN RIGHT ADDRESSES AND DO NEXT POINT
V13:	MOV R2,@PTU	;KER*U(I) IN RTSUM +2*I
	MOV R0,@PTV	;KER*V(I) IN ITSUM +2*I
	MOV TW,@PTP	;W(I) IN PSSUM +2*I
	MOV R4,@PTK	;KER(I) IN KER+2*I
	DEC TN	;ADDED IN ALL N PTS?
	BEQ UPK		;YES,DO NEXT SWEEP
;INCREMENT POINTERS
	INC XPT
	INC XPT
	INC YPT
	INC YPT
	INC PTU
	INC PTU
	INC PTV
	INC PTV
	INC PTP
	INC PTP
	INC PTK
	INC PTK
	JMP NXLP	;NO,DO NEXT POINT
;ADD KER(I) TO KSUM(I),I=0,1,...,(N-1)
UPK:	MOV #KER,R0
	MOV #KSUM,R1
	MOV N,R4
KAD:	MOV (R0),R2
	MOV (R1),R3
	ASH KS,R2	;PRESCALE
	ADD R2,R3
	BVC KOK
	JSR PC,KSCALE
	JMP KAD
KOK:	MOV R3,(R1)+
	TST (R0)+
	SOB R4,KAD
	RTS PC
;SUBROUTINE TO CALCULATE FILTERED AVERAGE--TAKES INVERSE FT OF 
	;FILTERED FT SUM
FCAL:	MOV N,R5	;MOVE ARRAY TO BE TRANSFORMED TO PTX
			;PTY ADDRESSES
	MOV #PTX,R0
	MOV #PTY,R1
	MOV #RTSUM,R2
	MOV #ITSUM,R3
RPL:	MOV (R2)+,(R0)+
	MOV (R3)+,(R1)+
	SOB R5,RPL
	MOV #1,DIR	;DIR CONTAINS 1 FOR INVERSE TRANSFORM
	JSR PC,FFT
;RETURNS WITH FILTERED SUM S=IFT(KER*S^) IN PTX,PTY ADDRESSES,I=0,...,(NPT-1)
	RTS PC
;DISPLAY ADDRESSES PTX,...,PTX+NPT AND SUM,...SUM+NPT
;WRITE ON DISK
;REINITIATE FOR NEXT CURVE

;SUBROUTINE TO SCALE ON OVERFLOW IN FILTER CALCULATIONS
SCALE1:	MOV #PTX,R0	;INCOMING REAL
	MOV #PTY,R1	;INCOMING IMAGINARY
	MOV #RTSUM,R2	;ONGOING REAL
	MOV #ITSUM,R3	;ONGOING IMAG
	MOV #PSSUM,R4	;ONGOING POWER SUM
	MOV N,R5	;COUNT
DAGN:	ASR (R0)+
	ASR (R1)+
	ASR (R2)+
	ASR (R3)+
	ASR (R4)	;SHIFT W RT 2 ,NOT 1, SINCE W IS SUM OF SQUARES
	ASR (R4)+
	SOB R5,DAGN	;THRU ALL PTS?--NO, GO BACK
	INC SC		;YES,INCREMENT S BY 1
	RTS PC
;SUBROUTINE TO READ TYPED DECIMAL NUMBER INTO R1
TINO:	CLR R1
	CLR R0
AN1:	.TTYIN
	CMPB R0,#15
	BEQ B1
	SUB #60,R0
	MUL #12,R1		;MULTIPLY BY 10
	ADD R0,R1
	BR AN1
B1:	.TTYIN
	CMPB R0,#12
	BNE B1
	RTS PC

;SUBROUTINE TO WRITE ON DISK
RITE:	BIC #10000,@#44		;CM0
	MOV #153520,@#60	;ADDRESS OF MONITOR'S KEYBOARD INTERRUPT SERVICE ROUTINE
	MOV #340,@#62
	BIS #100,TKS
	.PRINT #FQ		;FILE NAME?
	.LOCK
	.CSISPC #DBLK,#DEXT
	.ENTER #AREA,#0,#DBLK
	.UNLOCK
	.PRINT #FD		;FILE ENTERED
	.WRITW #AREA,#0,PLTPTR,NPT,#0	;CHANNEL 0,ADDRESS OF MEMORY TO BE WRITTEN IS IN PLTPTR,NUMBER OF WORDS 
					;TO BE WRITTEN IS IN NPT,START EACH FILE AT BLK 0
	BCS BADWRT		;WRITE ERROR
	.LOCK
	.CLOSE #0
	.UNLOCK
	.PRINT #GMSG		;WRITE OK MSG
BER:	BIS #10000,@#44		;CM1
	MOV #KEYINT,@#60
	MOV #200,@#62
	BIS #100,TKS
	RTS PC
BADWRT:	.PRINT #WMSG
	JMP BER
GMSG:	.ASCIZ /WRITE OK/
WMSG:	.ASCIZ /WRITE ERROR/
FQ:	.ASCIZ <CR><LF> /FILE NAME?/
FD:	.ASCIZ <CR><LF> /FILE ENTERED/
WFQ:	.ASCIZ <CR><LF> /WHICH FILTER?/
MSG:	.ASCIZ <CR><LF>	/PRESS "B" TO BEGIN ADDING/ 
MMSG:	.ASCIZ	<CR><LF>	/PRESS "M" TO MONITOR SWEEPS/
MSG2:	.ASCIZ <CR><LF>	/NSWEEPS DONE/
Q1:	.ASCIZ	/NUMBER OF SWEEPS?/
Q2:	.ASCIZ	/NUMBER OF POINTS PER SWEEP?/
Q3:	.ASCIZ	/TIME BETWEEN POINTS?/
MSG3:	.ASCIZ	<CR><LF>	/AVERAGE COMPUTED/
M3:	.ASCIZ	<CR><LF> 	/FILTERED AVERAGE COMPUTED/
E1:	.ASCIZ <CR><LF> /DIVISION ERROR/
E2:	.ASCIZ <CR><LF> /NO FILTER MATCH?/
M1:	.ASCIZ	/TRANSFORM DONE/
M2:	.ASCIZ	/WHAT NEXT?/
M5:	.ASCIZ <CR><LF> /FFT SCALE FACTOR/
M6:	.ASCIZ	<CR><LF>
M7:	.ASCIZ <CR><LF> /SUM SCALE FACTOR/
M8:	.ASCIZ <CR><LF> /KERNEL SCALE FACTOR/
M9:	.ASCIZ <CR><LF> /IFT SCALE FACTOR/
M15:	.ASCIZ <CR><LF> /KSUM SCALE FACTOR/
	.EVEN

;FFT SUBROUTINE
;N CONTAINS NUMBER OF POINTS = 2^M
;PTX = STARTING ADDRESS OF REAL PART OF DATA
;PTY = STARTING ADDRESS OF IMAGINARY PART OF DATA
;DELT CONTAINS DISTANCE BETWEEN POINTS IN A PAIR =N/2,N/4,...2,1
;PR CONTAINS NUMBER OF PAIRS PER CELL = N/2,N/4,...,2,1 = DELT
;CN CONTAINS # CELLS = 1,2, 4,...,N/2 = 2^(L-1),WHERE L = PASS #
;2*PR*CN =N
;DISTANCE BETWEEN CELLS = N/CN = 2*PR
;SF CONTAINS SCALE FACTOR--AT END OF TRANSFORM,EACH POINT HAS BEEN MUL-
	;TIPLIED BY 2^(-SF)
;DIR CONTAINS 0 FOR FORWARD TRANSFORM,1 FOR INVERSE TRANSFORM


FFT:	MOV N,DELT
	CLR SF
;COMPUTE LOG2N
	CLR R0
	MOV #1,R1
TW1:	ASL R1
	INC R0
	CMP R1,N
	BLT TW1
	MOV R0,M	;M=LOG2(N) --USE M TO COUNT PASSES
	MOV R0,LOGN	;SAVE LOG2N

	MOV #31104,PI	;PI=3.11037 (OCTAL)--PT AFTER BIT 12
PASS:	MOV DELT,TPR	;TPR =2*PR=#POINTS PER CELL =N,N/2,...,2
	ASR DELT	;DELT =N/2,N/4,...2,1=PR
	MOV DELT,HPR	;HPR=PR/2=N/4,...1,0
	ASR HPR
	DEC M		;M=# PASSES--M=M-L,WHERE L=PASS #
	MOV M,R0	;COUNT = M-L,WHERE L=PASS #
	TST R0
	BEQ FI		;M=0=>ON LAST PASS--DON'T USE ARG
	MOV PI,R2
DIV:	ASR R2		;PI/2^(M-L) IN R2
	SOB R0,DIV
	MOV R2,TARG	;PI/2^(M-L) IN TARG
	CLR ARG		;0 IN ARG INITIALLY ON EACH PASS
FI:	CLR I		;I IS THE INDEX
NXTI:	MOV #2,ROT	;FOR PI/2 SHIFT
	MOV I,TEMPI	;STORE FIRST VALUE OF I
;Z(I)=Z(I)+Z(I+DELT)
;Z(I+DELT)=W^K(Z(I)-Z(I+DELT)),WHERE K=I*2^(L-1),W=EXP(+-2*PI*J/N)
	;+FOR FORWARD TRANSFORM,-FOR INVERSE TRANSFORM
;COMPUTE SINE AND COSINE--ARG=2*PI*K/N=PI*I/2^(M-L)
	TST I		;SIN(0)=0,COS(0)=1
	BNE XI
	MOV #0,SIN
	MOV #10000,COS		;PT AFTER BIT 12
	BR NEXT
;I*PI/2^(M-L) IN TARG
XI:	ADD TARG,ARG	;(PI/2^(M-L)IN TARG)+((PI/2^(M-L))*(I-1)
			;IN ARG)=((PI/2^(M-L))*I IN ARG)
	MOV ARG,R1	;ARG IN R1 FOR SINE SUBROUTINE
;SINE
	JSR PC,SIN1	;X IN R1--RETURNS SIN(X),0<=X<=PI/2
	MOV R1,SIN	;SIN(X) IN R1
;COSINE
COMCOS:	MOV PI,R5	;PI/2-X
	ASR R5		;PI/2 IN R5
	SUB ARG,R5	;	PI/2-(I*PI/2^(M-L)) IN R5
	MOV R5,R1
	JSR PC,SIN1	;RETURNS SIN(PI/2-X)=COS(X),0<=X<=PI/2
	MOV R1,COS
;COMPUTE ADDRESSES OF POINTS IN PAIR
NEXT:	MOV I,R0	;Z(I)
	ASL R0		;ADDRESS IS PT+I*2--2*I IN R0
	MOV I,R1	;Z(I+DELT)
	ADD DELT,R1	;I+DELT IN R1
	ASL R1		;2*(I+DELT) IN R1
	
;SCALE ON OVERFLOW
;PUT X(I),X(I+DELT),Y(I),Y(I+DELT) IN TEMPORARY ADDRESSES TO RESTORE
;IF SCALING IS NECESSARY

A1:	MOV PTX(R0),TX1	;STORE X(I) IN TX1 FOR REPLACEMENT IN 
			;PTX+2*I IF SCALING IS NECESSARY
	MOV PTX(R1),TX2
	MOV PTY(R0),TY1
	MOV PTY(R1),TY2

;DO ARITHMETIC
;FIRST POINT OF PAIR
	MOV PTX(R0),R2		;PUT X(I) IN R2
	ADD PTX(R1),PTX(R0)	;X(I+DELT)+X(I) TO PTX+2*I
	BVC OK1
	JSR PC,SCALE
	BR A1		;TRY AGAIN WITH SCALED DATA
OK1:	MOV PTY(R0),R4	;PUT Y(I) IN R4
	ADD PTY(R1),PTY(R0)	;Y(I+DELT)+Y(I) TO PTY+2*I
	BVC OK2
	JSR PC,SCALE
	BR A1
;SECOND POINT OF PAIR
OK2:	SUB PTX(R1),R2		;X(I)-X(I+DELT) IN R2
	BVC OK3
	JSR PC,SCALE
	BR A1
OK3:	MOV R2,TEMPX
	MUL COS,R2	;COS*(R2) IN R2,R3--HI IN R2
	ASHC #4,R2	;LEFT SHIFT 4 FOR FRACTIONAL MUL WITH PT AFTER
			;BIT 12
	SUB PTY(R1),R4	;Y(I)-Y(I+DELT) IN R4
	BVC OK4
	JSR PC,SCALE
	BR A1
OK4:	MOV R4,TEMPY
	MUL SIN,R4	;SIN*(R4) IN R4,R5--HI IN R4
	ASHC #4,R4	;L S FOR FRACTIONAL MUL
	TST DIR
	BEQ FWX
	NEG R4		;INVERSE
FWX:	SUB R4,R2	;COS*(X(I)-X(I+DELT))-SIN*(Y(I)-Y(I+DELT)) IN R2
			;COS(X(I)-X(I+DELT))+SIN(Y(I)-Y(I+DELT)) FOR INVERSE
	BVC OK5
	JSR PC,SCALE
	BR A1
OK5:	MOV R2,PTX(R1)	;REAL PART TO PTX+2*(I+DELT)
	MOV TEMPX,R2	;X(I)-X(I+DELT) IN R2
	MUL SIN,R2	;SIN*(R2) IN R2,R3--HI IN R2
	ASHC #4,R2
	MOV TEMPY,R4	;Y(I)-Y(I+DELT) IN R4
	MUL COS,R4	;COS*(R4) IN R4,R5--HI IN R4
	ASHC #4,R4
	TST DIR
	BEQ FWY
	NEG R2		;INVERSE
FWY:	ADD R4,R2	;SIN*(X(I)-X(I+DELT))+COS(Y(I)-Y(I+DELT))IN R2
			;-SIN(X(I)-X(I+DELT))+COS(Y(I)-Y(I+DELT)) INVERSE
	BVC OK6
	JSR PC,SCALE
	BR A1
OK6:	MOV R2,PTY(R1)	;IMAGINARY PART TO PTY+2*(I+DELT)


;PAIR INDEXED BY I+TPR USES SAME SIN AND COS AS I
	ADD TPR,I	;I=I+TPR
	CMP I,N		;I-N<0,DO CORRESPONDING PAIR IN NEXT CELL--
			;SAME ARG
			;I-N>=0,DO CORRESPONDING PAIR IN SAME CELL--
			;ARGS DIFFER BY PI/2
	BLT NEXT
	TST HPR
	BEQ SHUF	;HPR=0=>HAVE FINISHED LAST PASS HERE
	DEC ROT		;2 CORRESPONDING PAIRS IN 1 CELL
	BEQ NXTPR
	MOV TEMPI,I	;RESTORE FIRST VALUE OF I
;I+HPR USES PI/2 ROTATION
	ADD HPR,I	;I=I+HPR
	MOV SIN,R0	;ROTATE BY PI/2
	MOV COS,SIN	;SIN(X+PI/2)=COS(X)
	NEG R0		;-SIN(X) IN R0
	MOV R0,COS	;COS(X+PI/2)=-SIN(X)
	BR NEXT		;DO PI/2-SHIFTED PAIR IN SAME CELL AND 
			;CORRESPONDING PAIRS (SAME ARG) IN FOLLOWING CELLS
NXTPR:	MOV TEMPI,I	;RESTORE I
	INC I		;DO NEXT PAIR--I=0,1,2,...,(HPR-1)
	CMP I,HPR	;I-HPR<0,DO NEXT PAIR
			;I-HPR>=0,NEXT PASS
	BGE PASTST
	JMP NXTI
PASTST:	TST M		;M>0,DO NEXT PASS
			;M=0,TRANSFORM DONE
	BEQ SHUF
	JMP PASS

;SHUFFLE
SHUF:	CLR I		;I=0,1,2,...,N-1
INVT:	CLR R1		;R1 WILL CONTAIN BIT-INVERTED INDEX
	MOV I,R0	;BIT INVERT R0 LOG2(N) PLACES TO GET INDEX 
			;OF THE POINT TO PUT IN THE ITH PLACE
			;AND WHERE TO PUT THE PT NOW IN THE ITH PLACE
	MOV LOGN,R2
	CLC		;CLEAR CARRY--ROTATE BIT FROM RIGHT END OF
			;R0 INTO CARRY,OUT OF CARRY INTO RIGHT END 
			;OF R1
SHIFT:	ROR R0
	ROL R1
	SOB R2,SHIFT		;ROTATE LOG2(N) TIMES
	MOV I,R0
	CMP R0,R1	;R1 CONTAINS BIT-INVERTED VALUE OF I,REL LOG2N
			;(R0)-(R1)<0,INTERCHANGE
			;(R0)-(R1)>=0, NO INTERCHANGE
	BGE NXTPT
;INTERCHANGE PTX+2*(R0) AND PTX+2*(R1)
	ASL R0		;2*I IN R0
	ASL R1		;2*I' IN R1,WHERE I'=BIT INVERTED I
	MOV PTX(R0),R3	;(PTX+2*I) IN R3
	MOV PTX(R1),PTX(R0)	;(PTX+2*I') TO PTX+2*I
	MOV R3,PTX(R1)	;(PTX+2*I) TO PTX+2*I'
	MOV PTY(R0),R3
	MOV PTY(R1),PTY(R0)
	MOV R3,PTY(R1)

NXTPT:	INC I
	CMP I,N		;I-N<0,NEXT POINT
			;I-N>=0,DONE
	BMI INVT
;SAVE SF IN SFF IF FORWARD FT
;DIVIDE BY 2^(LOGN-SF-SFF IF INVERSE FT (THIS ASSUMES SF IS SAME 
;FOR FORWARD FT FOR EACH SWEEP)
	TST DIR		;IFT?
	BNE FSCL	;YES
	MOV SF,SFF	;NO,FORWARD FT
	JMP DONE
FSCL:	MOV LOGN,R3	;SCALE IFT
	SUB SF,R3	
	SUB SFF,R3
	SUB SC,R3
	NEG R3
	MOV #PTX,R0
	MOV #PTY,R1
	MOV N,R2
OLOOP:	MOV (R0),R4
	MOV (R1),R5
	ASH R3,R4	;X(I)/2^(LOGN-SF-SFF) IN R4
	ASH R3,R5	;Y(I)/2^(LOGN-SF-SFF) IN R5
	MOV R4,(R0)+
	MOV R5,(R1)+
	SOB R2,OLOOP
DONE:	RTS PC
;SUBROUTINE TO SCALE ON OVERFLOW
;RESTORE DATA POINTS TO ORIGINAL ADDRESSES
SCALE:	MOV TX1,PTX(R0)
	MOV TX2,PTX(R1)
	MOV TY1,PTY(R0)
	MOV TY2,PTY(R1)
;SHIFT ALL POINTS RIGHT 1 BIT
	MOV N,R5	;COUNT
	MOV #PTX,R2	;REAL
	MOV #PTY,R3	;IMAGINARY
LOOP:	ASR (R2)+
	ASR (R3)+
	SOB R5,LOOP
	INC SF 		;INCREASE SF BY 1 WHENEVER SCALING IS DONE
	RTS PC
	
;SUBROUTINE TO COMPUTE SIN(X),0<-X<=PI/2
;USES FIRST 4 TERMS OF THE TAYLOR SERIES,WITH BINARY POINT AFTER BIT 12
SIN1:	MOV R1,X	;X IN R1
	MOV R1,R2
	MUL R1,R2	;X^2 IN R2,R3--HI IN R2
	ASHC #4,R2	;L S 4 FOR PT AFTER BIT 12;L S 1 FOR PT AFTER BIT 15
	MUL R1,R2	;X^3 IN R2,R3
	ASHC #4,R2
	MUL #1253,R2	;X^3/3! IN R2,R3--12525 FOR PT AFTER BIT 15
	ASHC #4,R2
	SUB R2,R1	;(R1)-(R2)=X-X^3/3! IN R1
	MUL @#X,R2	;X^4/3! IN R2,R3
	ASHC #4,R2
	MUL @#X,R2	;X^5/3! IN R2,R3
	ASHC #4,R2
	MUL #315,R2	;X^5/5! IN R2,R3--3146 FOR PT AFTER BIT 15
	ASHC #4,R2
	ADD R2,R1	;(X-X^3/3!)+X^5/5! IN R1
	MUL @#X,R2	;X^6/5! IN R2,R3
	ASHC #4,R2
	MUL @#X,R2	;X^7/5! IN R2,R3
	ASHC #4,R2
	MUL #141,R2	;X^7/7! IN R2,R3--1414 FOR PT AFTER BIT 15
	ASHC #4,R2
	SUB R2,R1	;(X-X^3/3!+X^5/5!)-X^7/7! IN R1
	RTS PC

;PRINT 6 DIGIT OCTAL NUMBER IN R1
OCOUT:	MOV R1,R0
	ASH #-17,R0	;(15)->(0) IN R0
	BIC #177776,R0	;CLEAR ALL BUT (0) IN R0
	BIS #60,R0	;ASCII BIAS
	.TTYOUT		;PRINT (15)
	MOV #5,DIGCT		;CT DIGITS
	MOV #-14,RTSHF
PD:	MOV R1,R0
	ASH RTSHF,R0	;(14,13,12)->(2,1,0),...
	JSR PC,TY	;TYPE DIGIT
	ADD #3,RTSHF
	DEC DIGCT	;THRU ALL DIGITS?
	BNE PD		;NO,DO NEXT DIGIT
	RTS PC		;YES,RETURN

;SUBROUTINE TO TYPE DIGIT IN (2,1,0) OF R0
TY:	BIC #177770,R0	;CLEAR ALL BUT (2,1,0) IN R0
	BIS #60,R0	;ASCII BIAS
	.TTYOUT
	RTS PC

AREA:	.BLKW 10
DEXT:	.WORD 0,0,0,0
DBLK:	.BLKW 60
PTK:	.WORD 0
KS:	.WORD 0
BEGIN:	.WORD 0
DISCT:	.WORD 0
DISPTR:	.WORD 0
PLTPTR:	.WORD 0
KEYMON:	.WORD 0
WFD:	.WORD 0
VCT:	.WORD 0
CT:	.WORD 0
H:	.WORD 0
H1:	.WORD 0
VDS:	.WORD 0
MNPT:	.WORD 0
TEMP:	.WORD 0
KEY:	.WORD 0
PSF:	.WORD 0
SFF:	.WORD 0
DIGCT:	.WORD 0
RTSHF:	.WORD 0
NSWEEP:	.WORD 0
NPT:	.WORD 0
RATE:	.WORD 0
MPT:	.WORD 0
S1:	.WORD 0
SC:	.WORD 0
N:	.WORD 0
DELT:	.WORD 0
PR:	.WORD 0
SF:	.WORD 0
DIR:	.WORD 0
M:	.WORD 0
LOGN:	.WORD 0
PI:	.WORD 0
TPR:	.WORD 0
HPR:	.WORD 0
I:	.WORD 0
ROT:	.WORD 0
TEMPI:	.WORD 0
TARG:	.WORD 0
ARG:	.WORD 0
SIN:	.WORD 0
COS:	.WORD 0
TX1:	.WORD 0
TX2:	.WORD 0
TY1:	.WORD 0
TY2:	.WORD 0
TEMPX:	.WORD 0
TEMPY:	.WORD 0
X:	.WORD 0
TN:	.WORD 0
TU:	.WORD 0
TV:	.WORD 0
TW:	.WORD 0
XPT:	.WORD 0
YPT:	.WORD 0
PTU:	.WORD 0
PTV:	.WORD 0
PTP:	.WORD 0
PTX:	.BLKW 4000	;RESERVE SPACE FOR 2048 REAL POINTS
PTY:	.BLKW 4000	;2048 IMAGINARY POINTS
SUM:	.BLKW 2000	;N0 OF PTS IN SAMPLE SUM IS 1/2
AVG:	.BLKW 2000	;NO OF PTS FOR FFT
RTSUM:	.BLKW 4000
ITSUM:	.BLKW 4000
PSSUM:	.BLKW 4000
FAVG:	.BLKW 2000
PS:	.BLKW 4000
KER:	.BLKW 4000
KSUM:	.BLKW 4000
MEAN:	.BLKW 2000
	.END START