/DAFFT -AN OVERLAY TO DAQUAN MS FOR FFT ON AVERAGED DATA. / /DEC-8E-ADAFA-A-LA / /COPYRIGHT 1972 /DIGITAL EQUIPMENT CORPORATION /MAYNARD MASSACHUSETTS 01754 / /DAFFT OVERLAY FOR PS-8,CONTAINS FFT /FILE DAFFT.3 /THIS IS A SET OF OVERLAYS FOR DAQUAN WHICH WILL /CREATE DAFFT--DAQUAN+FFT ALA ROTHMAN: /INCLUDED ARE ROUTINES TO CREATE POWER SPECTRUM, /DO HANNING SMOOTHING, AND RETAIN DATA IN A THIRD ARRAY. /DAFFT REQUIRES ADVANCED LAB 8/E WITH KE8/E (EAE). /LOAD DAQUAN(E) BINARY, FFTS-C INTO FIELD 1, THEN THESE OVERLAYS. /PUNCH FIELD 0, 0-7577 & FIELD 1, 0-1577. FIXMRI CALL= 4400 SETM2= 7344 SET3= 7325 MQA= 7501 CDF= 6201 CIF= 6202 LSR= 7417 NPTS= 70 NMI= 7411 SCA= 7441 TEM3= 134 HEDIT1= 4430 GETNO= 4422 BLOCK= 140 YIND= 11 ZIND= 12 CNTR= 126 DISPLAY=5466 TEMP= 131 MQL= 7421 MUY= 7405 SHL= 7413 CLRR= 3540 SMOT11= 3025 ASR= 7415 BEGDIS= 222 AUTO= 13 INITAR= 4425 QUERY= 32 SPINIT= 3014 TEMFP= 150 FTEM1= 142 FIXT= 4421 FLOUT= 4406 CRLFD= 4427 FLOAT= 4420 DP= 57 SCPINT= 230 M13= 363 K1000= 124 TEMP1= 132 FIELD 0 /UPDATE COMMAND TABLE *2103 SMCHK *2132 -1124 /IT: INVERSE TRANSFORM DIFT0 *2156 -624 /FT: FOURIER TRANSFORM DFT0 -2017 /PO-WER SPECTRUM POWR -3224 /ZT: ZERO IMAG ARRAY AND DO FFT ZAPIT -2310 /SH-IFT ARRAY TO OR FROM STORAGE BUFFER SAVE -1706 /OF-FSET FOR FFT DISPLAY. SETOFF 0 /END TABLE *1321 CALL STPTS /SET UP FFT PTS & POWER OF 2 *3 STPTS, STPT RSTPT, RSTPTS *76 1577 /ARRAY BASE-1 2000 /AND OFFSET *4011 CALL RSTPT /NEW TEXT OVERLAYS SQUEEZE *3716 HDS, TEXT /H,3,11(0,1,2):/ HDPR, TEXT /FACTOR=/ HDSA, TEXT /SAVE?/ /ROUTINE TO PRINT FFT SCALE FACTOR IN AC TRET, FLOAT HEDIT1 HDPR SET3 /3 DIGIT INTEGER DCA DP FLOUT DISPLAY *200 /CANCLE MOVING PAGE 7600 FIELD 1. 7200 *2242 7200 *2451 7000 *2525 7200 *3747 7000 *2223 7200 *2504 7000 *2535 7000 *2552 /INSERT SWBA'S INTO THE PROPER PLACES IN DAQUAN 7447 /SWBA *602 7447 /SWBA *766 7447 /SWBA /FIELD 1 DISPATCHER & NEW CODE WILL OVERLAY 'MODIFY' FUNCTION. *400 /TO 554 IS OPEN ZAPIT, ISZ BLOCK /CLEAR BUFFER 2 CALL ZAP DCA BLOCK DFT0, CDF CIF 10 /DO FFT JMP I .+1 DFT ZAP, CLRR DIFT0, CDF CIF 10 /DO INVERSE FFT JMP I .+1 DIFT STPT, 0 /SUB TO SET UP NO. PTS. AND POWER. DCA NPTS /OF 2 FOR FFT IN FIELD 1 TAD NPTS NMI CLA SCA CIA TAD KP12 CDF 10 DCA I PW2 /IS POWER OF 2 TAD NPTS DCA I PTS1 /AND NO. PT. CDF 0 JMP I STPT KP12, 12 PW2, 4 PTS1, 3 /SUB. SIMILAR TO AVERAGER RSTPTS, 0 JMS STPT TAD NPTS DCA TEM3 JMP I RSTPTS /COMPUTE POWER SPECTRUM POWR, HEDIT1 /FACTOR= HDPR GETNO /GET SCALE DOWN FACTOR FIXT TAD M1 DCA PSH /AS A SHIFT COUNTER DCA BLOCK CALL SINT /INIT POINTERS,ETC. PLP, TAD I YIND /GET REAL,SQUARE IT JMS PMUL DCA TEMP TAD I AUTO /GET IMAG.,SQUARE IT JMS PMUL TAD TEMP DCA I ZIND /STORE SUM ISZ CNTR JMP PLP DISPLAY SINT, SPINIT M1, -1 /SUB. TO SQUARE AC AND SCALE IT DOWN PMUL, 0 SPA CIA DCA MPD /ABS. VALUE TAD MPD MQL MUY MPD, 0 LSR PSH, 0 /SET BEFORE ENTRY CLA MQA JMP I PMUL /ROUTINE TO SAVE OR RESTORE DATA FROM 3RD ARRAY SAVE, HEDIT1 /SAVE? HDSA DCA BLOCK /SET UP POINTERS,ETC INITAR TAD AR3 DCA AUTO CALL QUERY /GET ANSWER,0=YES,1=NO CDF 10 SNA CLA JMP SVE /SAVE CHANNEL TAD I AUTO DCA I YIND /GET STORAGE ARRAY INTO CHANNEL 1 ISZ CNTR JMP .-3 DISPLAY SVE, TAD I YIND /PUT CHANNEL 1 INTO STORAGE ARRAY DCA I AUTO ISZ CNTR JMP .-3 DISPLAY AR3, 5577 /BASE OF STORAGE ARRAY-1 /NEW SMOOTHING ROUTINE. SMCHK, HEDIT1 /H,3,11(0,1,2): HDS GETNO /GET CODE FIXT SNA JMP SHAN /DO HANNING FILTER CLL RAR SZA CLA JMP I SHAV /DO 11 PT. TAD CNP /DO 3 PT CDF CIF 10 /SET INSTRUCTION AND DO IT DCA I INS1 JMP I .+1 RLSM SHAV, SMOT11 CNP, NOP SHAN, CDF CIF 10 /EXIT TO HANNING FILTER JMP I .+1 HNSM INS1, INST SETOFF, TAD FTOFST /SET DISPLAY OFFSET. SNA CLA TAD K1000 DCA FTOFST DISPLAY *SCPINT+4 JMP M13+1 *M13+1 TAD FTOFST DCA TEMP1 JMP SCPINT+5 *167 FTOFST, 0 /PREVENT BOX CAR INTEGRATION. *4233 DCA 0 *4251 TAD 0 FIELD 1 *200 DFT, JMS MIDSET /0 MEAN THE DATA AND DO FFT CALL DOFFT SKP DIFT, CALL DOIFFT /DO INVERSE FFT CALL SORT /DO BIT INVERT STORAGE TAD SCALE CDF CIF 0 JMP I .+1 /PRINT SCALE FACTOR TRET HNSM, TAD CIAC /DO HANNING SMOOTH DCA INST /SET INSTRUCTION TAD XLOCDF /SET AC TO HANN IMAG. JMS HANM RLSM, JMS HANM /THEN REALS CDF CIF 0 JMP I .+1 BEGDIS /EXIT /SUB TO DO HANNING OR 3 PTS FILTER HANM, 0 JMS PSTR /SET UP POINTERS 7447 /SWBA ISZ CNT TAD I 13 /DIVIDE Y(I) BY 2 ASR 0 DCA TMPX TAD I 13 /DIVIDE Y(I+1) BY 4 ASR 0 DCA TMPY TAD TMPY JMP INST+1 SMGO, TAD I 13 ASR 0 DCA TMPY TAD TMPY ASR 0 TAD TMP INST, CIA TAD TMPX DCA I 14 TAD TMPX ASR 0 DCA TMP ISZ CNT SKP JMP .+4 TAD TMPY DCA TMPX JMP SMGO TAD TMPX TAD TMPY DCA I 14 JMP I HANM TMPY, 0 MM1, -1 CNT, 0 TMP, 0 TMPX, 0 CIAC, CIA RTP, XRTAB-1 /POINTER TO REAL ARRAY-1 PSTR, 0 TAD RTP DCA 13 /POINTER TO REALS TAD N CIA DCA CNT /PT COUNT TAD 13 DCA 14 /14=13 DCA TMP JMP I PSTR /SUB TO FIND MEAN AND SUBTRACT IT FROM ARRAY /TO GET ZERO MEAN DATA WITH /NO DC OFFSET IN FFT MIDSET, 0 JMS PSTR DCA TMPX ML1, CLL TAD I 13 SMA JMP ML1A TAD TMPX /ADD WHEN Y IS + DCA TMPX RAL TAD MM1 JMP ML1B ML1A, TAD TMPX /ADD WHEN Y IS - DCA TMPX RAL ML1B, TAD TMP /TMP IS HI ORDER,TMPX IS LOW DCA TMP ISZ CNT JMP ML1 STA /HAVE SUM, SET A SHIFT COUNT TAD NU /TO GET AVERAGE DCA MSH TAD TMPX MQL TAD TMP ASR MSH, 0 CLA MQA CIA DCA TMPX /IS - MEAN VALUE OF DATA JMS PSTR ML2, TAD I 13 /UPDATE ARRAY TAD TMPX DCA I 14 ISZ CNT JMP ML2 JMP I MIDSET /FFTS-COMPLEX : (VERSION D) / /THIS IS A SUBROUTINE FOR CALCULATING THE FAST FOURIER /TRANSFORMATION OF A SEQUENCE OF N COMPLEX TIME SAMPLES /WHICH ARE STORED IN MEMORY. IT IS FOR USE WITH A 4K /PDP-8 OR PDP-8/I COMPUTER EQUIPPED WITH AN ASR33 TELETYPE AND AN /EXTENDED ARITHMETIC ELEMENT OPTION AS MINIMUM HARDWARE. /BY JAMES ROTHMAN -- AUGUST, 1968 /PAGE ZERO *3 /TABLE PARAMETERS N, 0 /NUMBER OF POINTS IN COMPUTATION NU, 0 /POWER OF TWO OF POINTS IN COMPUTATION (N=2^NU) L, 0 /INDEX TO SHOW WHAT ARRAY IS BEING CONSTRUCTED S, 0 /GIVES SPACING BETWEEN NODE PAIRS IN THE LTH ARRAY. F, 0 /USED FOR SCALING NODE POSITION TO GET NUMBERS IN NODES. NOVER4, 0 /STORAGE FOR N/4. MAXNU, BIGSNU /LARGEST TABLE SIZE (POWER OF 2) MNOVR2, 0 /STORAGE FOR -N/2 *20 /INDEXING VARIABLES QR, 0 /POINTER TO REAL PART OF X(Q) QI, 0 /POINTER TO IMAG. PART OF X(Q) PR, 0 /POINTER TO REAL PART OF X(P) PI, 0 /POINTER TO IMAG. PART OF X(P) Q, 0 /NUMERICAL INDEX Q(=0,1,...,N-1) P, 0 /NUMERICAL INDEX P(=0,1,...,N-1) K, 0 /NUMBER IN THE NODE BEING OPERATED ON. /LOOP DELIMITERS C, 0 /INTERRUPTS COMPUTATION OF LTH ARRAY EVERY S PASSES /DATA VARIABLES ADD2, 0 /USED BY SUBROUTINE ADDR AS DATA (ADDEND) TEMPR, 0 /TEMPORARY STORAGE REGISTER FOR REAL PARTS SINE, 0 /TEMP. STORAGE FOR SIN(2*PI*K/N) COSINE, 0 /TEMP. STORAGE FOR COS(2*PI*K/N) GR, 0 /REAL PART OF PRODUCT (W^K)*X(P). TEMP STORAGE GI, 0 /IMAG. PART OF (W^K)*X(P). TEMP STORAGE /SUBROUTINE CALL LIST ADDER, ADDR /ADD C(AC) TO C(ADD2) AND SCALE RIGHT ONE SORT, SORTX /BIT INVERTED BUFFER SORTED. INVERT, INVRT /WORD IN AC OF NU BITS IS BIT INVERTED MULT, MULTIP /SINGLE PRECISION SIGNED MULTIPLY AC=ARG1;C(CALL+1)=ADD OF ARG2 GETRIG, TRIGET /FETCH SIN AND COS OF 2*PI*C(AC)/N. DOFFT, FFT /DO FFT OF THE INPUT BUFFER DOIFFT, IFFT /DO INVERSE OF BUFFER /DATA TABLES SINLOC, SINTAB /TABLE OF SIN(2*PI*I/N) FOR I=0,1,2,...,N-1 XRLOC, XRTAB /INPUT BUFFER AND TABLE OF ARRAYS (REAL PARTS) XLOCDF, XITAB-XRTAB /DIFFERENCE IN ADDRESS OF REAL AND IMAG PART TABLES /PSEUDO FLOATING POINT FORMAT FLAGS. SCALE, 0 /PSEUDO EXPONENT OF FOURIER COEFFICIENTS. SHFLAG, 1 /IF =1,ADD WITH SHIFT; IF=0,ADD WITH OUT SHIFT. SHFCHK, 0 /INDICATES IF ALL X'S IN AN ITERATION ARE <.5 /POINTERS TO SINE TABLE LOOK-UP SHIFTS SHIFT1, SHFT1 /THE NUMBER 10-NU MUST BE PLACED SHIFT2, SHFT2 /IN EACH OF THESE LOCATIONS. SHIFT3, SHFT3 *400 /COMPUTATION OF FIRST COMPLEX ARRAY FROM INPUT DATA /NUMBER OF INPUT POINTS IN "N" .LOG(2)(N)IN"NU". FOR DETAILS OF ALGORITHM, SEE FLOWCHART FFT, 0 CLA IAC CLL DCA L /L<=1 DCA SCALE /INITIALIZE FLOATING POINT FORMAT IAC DCA SHFLAG DCA SHFCHK TAD N CLL RTR /INITIALIZE PROGRAM CONSTANTS DCA NOVER4 TAD NU CIA TAD MAXNU DCA I SHIFT1 TAD I SHIFT1 DCA I SHIFT2 TAD I SHIFT2 DCA I SHIFT3 TAD N CLL RAR DCA S /S<=N/2 IS SPACING OF NODE PAIRS IN FIRST ARRAY TAD S CIA DCA MNOVR2 CMA /AC<=-1 TAD S /AC<=N/2-1 TAD XRLOC /BEGINNING OF TABLE OF REAL PARTS. DCA QR /Q<=N/2-1. QR POINTS TO WORD IN MEMORY, WHILE Q IS ACTUAL INDEX TAD NU CIA IAC DCA F /F<=1-NU (=L-NU SINCE L=1) LOOP1, TAD QR /QR=XRLOC+Q AT ALL TIMES. TAD S DCA PR /P<=Q+N/2 TAD QR /XLOCDF=XILOC-XRLOC (XILOC=BEGIN. OF IMAG. PARTS TABLE) TAD XLOCDF /QR+XLOCDF=(S+XRLOC)+(XILOC-XRLOC)=XILOC+S=QI DCA QI /QI=XILOC+Q AT ALL TIMES. QI POINTS TO IMAG. PART OF X(Q) TAD PR TAD XLOCDF /COMPUTE COMPLEX OPERATIONS X(P)<=X(Q)-X(P) AND X(Q)<=X(Q)+X(P) DCA PI /BY REAL AND IMAGINARY PARTS. TAD I QI /IM(X(Q)) (IM () MEANS IMAGINARY PART) DCA ADD2 /MAKE IT ADDEND. DO IMAG. PARTS FIRST TAD I PI /IM(X(P)) JMS I ADDER /FORM ADDITION IM[X(P)+X(Q)]=IM[X(P)]+IM[X(Q)] AND SCALE RIGHT DCA TEMPR /FOR SCALING, THEN STORE. TAD I QI /FORM DIFFERENCE IM[X(Q)-X(P)]=IM[X(Q)]-IM[X(P)] DCA ADD2 TAD I PI CIA JMS I ADDER DCA I PI /PUT AWAY AT IM[X(P)] TAD TEMPR /GET IM[X(P)+X(Q)] DCA I QI /PUT AT IM[X(Q)]. IMAGINARY PARTS DONE. TAD I QR /ADD REAL PARTS NEXT DCA ADD2 TAD I PR /RE=REAL PART JMS I ADDER /FORM RE[X(P)+X(Q)]=RE[X(P)]+RE[X(Q)] (DIVIDED BY 2) DCA TEMPR /STORE TAD I QR /GET RE[X(Q)] DCA ADD2 TAD I PR /AND RE[X(P)] CIA JMS I ADDER /FORM RE[X(Q)-X(P)] (DIVIDED BY 2) DCA I PR /PUT AT RE[X(P)] TAD TEMPR /GET RE[X(Q)+X(P)] DCA I QR /PUT AT RE[X(Q)]. REAL PARTS DONE TAD XRLOC /Q=QR-XRLOC CIA TAD QR /AC IS Q SPA SNA CLA /IS Q>0? (IE-THE WHOLE ARRAY HAS NOT BEEN COVERED) JMP CHKPT /NO. Q=0. DONE WITH FIRST ARRAY. MOVE ON TO OTHERS. CMA /YES. Q<=Q-1. MOVE UP THIS ARRAY. TAD QR /OR EQUIVALENTLY, QR<=QR-1 DCA QR JMP LOOP1 /DO NEXT NODE PAIR CHKPT, TAD L /L GIVES THE NUMBER OF THE VERTICAL ARRAY JUST BUILT CIA TAD NU /IS L=NU? (IE HAS THE LAST ARRAY BEEN COMPUTED?) SNA CLA JMP I FFT /YES. DONE. RESULTS STORED IN BIT REVERSED ORDER. TAD SHFCHK /GET SCALE FACTOR AND ADJUST FOR PROPER DCA SHFLAG /ADDITION ON NEXT ITERATION. TAD SHFCHK SNA CLA ISZ SCALE DCA SHFCHK ISZ L /L<=L+1. MOVE ON TO NEXT ARRAY TAD S /S GIVES SPACING BETWEEN NODE PAIRS, WHICH IS N/2^L CLL RAR /DIVIDE BY 2 AND PUT BACK, SO THAT ON THE LTH PASS THROUGH DCA S /S WILL=N/2^L, THE SPACING. ISZ F /F<=F+1. ON LTH PASS, F WILL BE F=L-NU, THE SCALE FACTOR FOR K. NOP /NOP FOR WHEN F=-1 TO PREVENT ERROR DUE TO SKIP CMA /AC<=-1 TAD N TAD XRLOC DCA PR /P<=N-1. PR POINTS TO RE[X(P=N-1)] SETC, CLA IAC DCA C /C<=1. C BREAKS BUILD LOOP EVERY S ITERATIONS BUILD, TAD PR /SO AS TO AVOID RE-COMPUTATION. TAD XLOCDF DCA PI TAD XRLOC /PR=XRLOC+P CIA TAD PR DCA P /ACTUAL INDEX IS P:(0,1,...,N-1) TAD F /BUILD ARRAY. F=L-NU. SHIFT "P"-F PLACES RIGHT (=NU-L) SNA /SHIFT ZERO PLACES? JMP NOROT /YES. LEAVE ALONE CMA /F COMPLEMENTED IS -F-1=-(F+1)=PLACES TO BE SHIFTED-1 DCA SHIFCT /CONTAINS -F-1 TAD P /GET NODE INDEX LSR /SHIFT P RIGHT SHIFCT+1=-F-1+1=-F=NU-L PLACES SHIFCT, HLT /STORAGE FOR SHIFT COUNT. SKP /AC<=INTEGER PART [P*2^F] NOROT, TAD P /NO ROTATION. JUST GET P=P*2^0 JMS I INVERT /INVERT BIT ORDER AND PUT IN K (NUMBER IN PTH NODE) TAD MNOVR2 /SUBTRACT N/2 TO GET NUMBER IN Q (=K) (P'S NODE PAIR.) JMS I GETRIG /GET REAL AND IMAGINARY PARTS OF W^K. ADJSGN, NOP /SET TO CIA FOR DOING IFFT, NOP FOR FFT. DCA SINE /SIN(2*PI*K/N)=-IM[W^K]. COS IN REGISTER COSINE. TAD I PR /FORM (W^K)*X(P)-A COMPLEX MULTIPLICATION JMS I MULT /DO REAL PART FIRST=RE[X(P)]*COSINE+IM[X(P)]*SINE COSINE /AC=RE[X(P)]*COSINE=RE[X(P)]*RE[W^K] DCA ADD2 /SAVE FOR ADDITION LATER TAD I PI /GET IM[X(P)] JMS I MULT SINE /AC=IM[X(P)]*SINE=-IM[W^K]*IM[X(P)] TAD ADD2 /AC=RE[W^K]*RE[X(P)]-IM[W^K]*IM[X(P)]=RE[X(P)*W^K] DCA GR /STORE AT GR /DO IMAG. PART NEXT=IM[X(P)]*COSINE-RE[X(P)]*SINE=IM[X(P)]*RE[W^K]+RE[X(P)]*IM[W^K] TAD I PI JMS I MULT /AC=IM[X(P)] COSINE /AC=IM[X(P)]*COSINE=IM[X(P)]*RE[W^K] DCA ADD2 /STORE FOR LATER ADDITION TAD I PR /AC=RE[X(P)] JMS I MULT SINE /AC=RE[X(P)]*SINE=-RE[X(P)]*IM[W^K] CIA /AC=RE[X(P)]*IM[W^K] TAD ADD2 /AC=IM[X(P)]*RE[W^K]+RE[X(P)]*IM[W^K]=IM[X(P)*W^K] DCA GI /STORE AT GI. SO GI=IM[X(P)*W^K] AND GR=RE[X(P)*W^K] G=GR+I*GI. TAD S /LOCATE P'S NODE PAIR Q. LOCATED S=N/(2^L) UP ARRAY. CIA /SO SET Q=P-S=INDEX OF NODE PAIR TAD PR /LOCATE X(Q) IN MEMORY BY FIXING POINTERS QR AND QI DCA QR /TO Q'S REAL AND IMAG. PARTS, RESPECTIVELY TAD QR TAD XLOCDF DCA QI TAD I QR /DO THE COMPLEX OPERATIONS: X(P)<=X(Q)-G;X(Q)<=X(Q)+G DCA ADD2 /FIRST DO REAL PART OF X(P). GET RE[X(Q)] AND STORE TAD GR /GET RE[G] CIA JMS I ADDER /SUBTRACT THEM. DCA I PR /RE[X(P)]<=RE[X(Q)]-RE[G] TAD I QI /COMPUTE IMAG. PART OF X(P). GET IM[X(Q)] DCA ADD2 /AND STORE TAD GI /GET IM[G] CIA JMS I ADDER /AND SUBTRACT THEM. DCA I PI /IM[X(P)]<=IM[X(Q)]-IM[G]. X(P) IS NOW DONE. TAD I QR /NEXT COMPUTE X(Q). FIRST REAL PART DCA ADD2 /GET RE[X(Q)] AND STORE TAD GR /GET RE[G] AND ADD TO FORM JMS I ADDER /RE[X(Q)]+RE[G]. DCA I QR /RE[X(Q)]<=RE[X(Q)]+RE[G]. TAD I QI /NOW COMPUTE IMAG PART OF X(Q). GET IM[X(Q)] DCA ADD2 /AND STORE TAD GI /GET IM[G] AND ADD TO FORM JMS I ADDER /IM[X(Q)]+IM[G] DCA I QI /IM[X(Q)]<=IM[X(Q)]+IM[G]. THE NEW NODE PAIR IS COMPUTED. CMA /MOVE UP ARRAY TO NEXT NODE. SET AC=-1 TAD P /TO FORM P-1 DCA P /P<=P-1 CMA TAD PR /DO THE SAME FOR POINTER PR DCA PR TAD C /CHECK ON SPACING. IS A NODE WHICH HAS ALREADY BEEN COMPUTED CIA /ABOUT TO BE RE-DONE, OR EQUIVALENTLY, TAD S /IS C=S? SZA CLA /YES. JMP CNOTS /NO. DO NEXT NODE PAIR TAD P /YES. BUT ARE WE AT THE TOP OF THE ARRAY? CMA /OR, IS S=P+1? (P COMPLEMENTED=-P-1=-(P+1) TAD S SNA CLA JMP I RECHK /YES. DONE WITH THIS ARRAY. DO NEXT ONE. TAD S /NO. MOVE PAST AREA THAT HAS ALREADY BEEN DONE, OR SET P TO P-S. CIA /BY CHANGING THE POINTER TO RE[X(P)] TAD PR DCA PR JMP I RESETC /REINITIALIZE C TO 1 SINCE AN UNUSED AREA HAS BEEN ENTERED. CNOTS, ISZ C /C<=C+1. ANOTHER NODE PAIR HAS BEEN HANDLED. JMP I RBUILD /DO NEXT NODE PAIR IN THIS AREA. RBUILD, BUILD /POINTERS TO RETURN LOCATIONS. RESETC, SETC /WHICH ARE LOCATED ON RECHK, CHKPT /ANOTHER PAGE. SORTX, 0 /SUBROUTINE THAT CMA /SORTS OUT TRANSFORMS BY TAD N /BIT INVERSION OF ADDRESS. DCA Q /Q<=N-1. START FROM BOTTOM OF BUFFER REVERS, TAD Q /P<=BIT INVERTED Q JMS I INVERT /BIT INVERSION ROUTINE DCA P TAD P /FORM Q-P CIA TAD Q SPA SNA CLA /IS P0? CMA CML IAC /NO. MAKE POSITIVE. SET LINK=1 TO SHOW IT WAS NEGATIVE. MQL /LOAD INTO MQ TAD I MULTIP /GET ADDRESS OF MULTIPLICAND DCA ARG2 /STORE TAD I ARG2 /AND RETRIEVE MULTIPLICAND ITSELF. ISZ MULTIP /(FOR EXIT AT CALL+2) SPA /ARG2>0? CMA CML IAC /NO. MAKE POSITIVE. CHANGE LINK,SINCE-1+-1=1 AND -1+1=-1 DCA ARG2 /PUT AWAY AT ARG2 RAL /SIGN IN LINK. PUT INTO AC11 AND DCA SIGN /PUT AWAY AT SIGN (= 1 IF -; =0 IF +) MUY /DO MULTIPLICATION ARG2, HLT /ARGUMENT 2 (MULTIPLICAND) SHL /NORMALIZE BINARY POINT. 0 DCA ARG2 /SAVE HIGH ORDER. NOW ROUND OFF. SHL /SET AC11=MQ0,AC0-10=0 0 MQL TAD SIGN /RESTORE PROPER SIGN CLL RAR /PUT SIGN IN LINK TAD ARG2 /BRING BACK RESULT MQA /RESULT=(HIGH ORDER) .OR. (BIT 0 OF LOW ORDER) SZL /POSITIVE SIGN? CMA IAC /NO. NEGATE JMP I MULTIP /EXIT. SIGNED RESULT IN AC. SIGN, 0 /BIT INVERSION ROUTINE /ENTRY: AC=WORD TO BE INVERTED; EXIT:AC=RESULT /NU CONTAINS THE NUMBER OF BITS IN THE WORD INVRT, 0 DCA WORD /GET WORD TO BE INVERTED DCA WORDP /ZERO OBJECT REGISTER TAD NU /GET NUMBER OF BITS TO BE CIA /INVERTED AND USE TO LIMIT THE DCA FLIPCT /EXTENT OF LOOP FLIP, TAD WORD /PULL OUT RIGHTMOST BIT OF WORD CLL RAR /(RIGHT MOST BIT NOW IN LINK) DCA WORD /(PUT BACK SO A NEW BIT IS OPERATED ON EACH TIME) TAD WORDP /AND PUSH INTO WORDP FROM LEFT RAL DCA WORDP ISZ FLIPCT /ALL BITS DONE? JMP FLIP /NO. DO NEXT BIT TAD WORDP /YES. PICK UP RESULT JMP I INVRT /AND EXIT WORD, 0 WORDP, 0 FLIPCT, 0 /THIS SUBROUTINE FETCHES THE VALUES OF SIN(2*PI*C(AC)/N) /AND OF COS(2*PI*C(AC)/N) FOR C(AC) < N/2+1 /ENTRY: AC=INDEX OF LOOK UP /EXIT : COS(2*PI*C(AC)/N) STORED AT "COSINE" AND / AC=VALUE OF SIN(2*PI*C(AC)/N). TRIGET, 0 DCA K /STORE C(AC) AT K. MQL /CLEAR MQ TAD K /FORM N/4-K. CLL CIA TAD NOVER4 DCA NO4MIK SZL /IS N/4-K<0? JMP QUAD1 /NO. FIRST QUADRANT ANGLE. QUAD2, TAD NO4MIK /2ND QUADRANT. GET -COS AT K-N/4. CIA LSR /MAKE CORRECTIVE RIGHT SHIFT ON INDEX. 0 SHL /FIND ON SINE TABLE FOR 2^MAXNU BY MULTIPLYING SHFT1, HLT /INDEX BY 2^(MAXNU-NU), WHICH IS STORED HERE. TAD SINLOC /LOCATE IT IN MEMORY. DCA INDEX TAD I INDEX CIA /2ND QUADRANT COS IS NEGATIVE. DCA COSINE TAD NO4MIK /GET SIN AT N/2-K TAD NOVER4 JMP SINRET QUAD1, TAD NO4MIK /GET COS AT N/4-K. LSR 0 SHL SHFT2, HLT TAD SINLOC DCA INDEX TAD I INDEX DCA COSINE TAD K /GET SIN AT K. SINRET, LSR 0 SHL SHFT3, HLT TAD SINLOC DCA INDEX TAD I INDEX /AC= SIN VALUE. JMP I TRIGET NO4MIK, 0 /STORAGE FOR N/4-K INDEX, 0 /POINTER TO SINE TABLE /THIS ROUTINE PERFORMS A SINGLE PRECISION ADD WITH ROUNDING. EACH ARGUMENT IS /SHIFTED RIGHT ONCE TO PREVENT OVERFLOW OF BINARY POINT (IF NECESSARY) /AND THEN CHECKED TO SEE IF IT CAN BE NORMALIZED AFTER ADDITION /ENTRY: AC=ADDEND,C(ADD2)=AUGEND /EXIT : AC=RESULT, DIVIDED BY TWO IF NECESSARY. ADDR, 0 DCA ADD1 TAD SHFLAG /SHOULD ADD BE DONE WITH SHIFT? SNA CLA JMP ADDWOS /NO. DO ADD WITH OUT SHIFT TAD ADD1 /YES. GET ADDEND ASR /DO 1 SIGNED RIGHT SHIFT 0 /MQ0=LOW ORDER (LO) OF ADD1 DCA ADD1 TAD ADD2 ASR /MQ0=LO(ADD2) 0 /MQ1=LO(ADD1) DCA ADD2 MQA /GET MQ RAL /L<=LO(ADD2); AC0<=LO(ADD1) CMA CML /COMPLEMENT BOTH. SMA SNL CLA /IF BOTH WERE=1 (NEITHER=0), INTRODUCE A CARRY. IAC ADDWOS, TAD ADD1 /DO THE ADDITION. TAD ADD2 DCA XSUM /STORE THE RESULT TAD XSUM /CHECK TO SEE IF ALREADY NORMALIZED. SPA /IS IT POSITIVE? CIA /MAKE IT POSITIVE. RAL /GET BIT 1. WAS NORMALIZED IF =1 SMA CLA JMP NOTNOR /NOT NORMALIZED. LEAVE SHFCHK ALONE. IAC DCA SHFCHK /SET SHFCHK=1 NOTNOR, TAD XSUM JMP I ADDR /AND EXIT ADD1, 0 /ADDEND STORAGE. XSUM, 0 /TEMPORARY STORAGE FOR SUM. /TABLE OF VALUES OF SIN (2*3.14159*I/1024) FOR I FROM /0 TO 256 INCLUSIVE. SINTAB, 0000 0015 0031 0046 0062 0077 0113 0130 0144 0161 0176 0212 0227 0243 0260 0274 0311 0325 0342 0356 0373 0407 0424 0440 0455 0471 0505 0522 0536 0553 0567 0603 0620 0634 0650 0664 0701 0715 0731 0745 0762 0776 1012 1026 1042 1056 1072 1106 1123 1137 1153 1166 1202 1216 1232 1246 1262 1276 1312 1325 1341 1355 1370 1404 1420 1433 1447 1462 1476 1511 1525 1540 1554 1567 1602 1616 1631 1644 1657 1672 1705 1720 1734 1747 1761 1774 2007 2022 2035 2050 2062 2075 2110 2122 2135 2147 2162 2174 2207 2221 2233 2246 2260 2272 2304 2316 2330 2342 2354 2366 2400 2411 2423 2435 2447 2460 2472 2503 2515 2526 2537 2551 2562 2573 2604 2615 2626 2637 2650 2661 2672 2703 2713 2724 2734 2745 2755 2766 2776 3007 3017 3027 3037 3047 3057 3067 3077 3107 3117 3126 3136 3145 3155 3164 3174 3203 3212 3222 3231 3240 3247 3256 3265 3274 3302 3311 3320 3326 3335 3343 3351 3360 3366 3374 3402 3410 3416 3424 3432 3440 3445 3453 3460 3466 3473 3501 3506 3513 3520 3525 3532 3537 3544 3551 3556 3562 3567 3573 3600 3604 3610 3614 3621 3625 3631 3635 3640 3644 3650 3653 3657 3662 3666 3671 3674 3700 3703 3706 3711 3713 3716 3721 3724 3726 3731 3733 3735 3740 3742 3744 3746 3750 3752 3754 3755 3757 3761 3762 3764 3765 3766 3767 3770 3771 3772 3773 3774 3775 3776 3776 3777 3777 3777 3777 3777 3777 3777 *1600 XRTAB, 0 /DATA BUFFER FOR REAL PARTS *XRTAB+2000 XITAB, 0 /DATA BUFFER FOR IMAGINARY PARTS DATAHI=XITAB+2000 /FIRST LOCATION AVAILABLE FOR PROGRAMMING /DEFINITIONS FOR EAE DVI=7407 NMI=7411 SHL=7413 ASR=7415 LSR=7417 MQL=7421 MUY=7405 MQA=7501 CAM=7621 SCA=7441 SCL=7403 /ASSEMBLY PARAMETERS BIGSNU=12 /LARGEST TABLE HAS DIMENSION 2^10. /FOLLOWING IS PATCH TO CORRECT BUT IN FFTS-C: *1014 RAR DCA SIGN MUY ARG2, HLT SHL 0 DCA ARG2 TAD SIGN SHL 0 TAD ARG2 SPA CLL STA RAR NOP SZL CIA JMP I MULTIP SIGN, 0 $