Directory of image this file is from
This file as a plain text file
/PAFFT -AN OVERLAY TO DAQUAN MS FOR POWER AVERAGER FFT. / /DEC-8E-APAFA-A-LA / /COPYRIGHT 1972 /DIGITAL EQUIPMENT CORPORATION /MAYNARD MASSACHUSETTS 01754 /
/PAFFT OVERLAY FOR DAQUAN UNDER PS8 /FILE PAFFT.2 /THIS IS A SET OF OVERLAYS FOR DAQUAN WHICH WILL /CREATE PAFFT--DAQUAN+FFT ALA ROTHMAN: /INCLUDED ARE ROUTINES TO CREATE AVERAGED POWER SPECTRUM, /DO HANNING SMOOTHING, AND RETAIN DATA IN A THIRD ARRAY. /PAFFT 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 SWAB= 7431 SWBA= 7447 DAD= 7443 DST= 7445 DCM= 7575 SETM2= 7344 SET3= 7325 PCTR= 125 OFFSET= 77 MQA= 7501 CDF= 6201 CIF= 6202 LSR= 7417 NPTS= 70 NMI= 7411 SCA= 7441 TEM3= 134 TEMP1= 132 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
FIELD 0 /UPDATA COMMAND TABLE. *2103 /MOD SM: EXIT 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 /MOD TO SET UP ARRAYS AND LINKS TO FFT. *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. OVERLAY SQUEEZE. *3716 HDS, TEXT /H,3,11(0,1,2):/ HDPR, TEXT /FACTOR=/ HDSA, TEXT /SAVE?/ /COME HERE FROM FFT TO PRINT SCALE FACTOR /IN AC ON ENTRY. TRET, FLOAT HEDIT1 HDPR SET3 DCA DP /3 PLACE INTEGER. FLOUT DISPLAY /EXIT TO DISPLAY. /NOP OUT PS-8 PAGE 7600 SWAP. *200 7200 *2242 7200 *2451 7000 *2525 7200 *3747 7000 *2223 7200 *2504 7000 *2535 7000 SWBA=7447 *2552 SWBA *602 SWBA *766 SWBA
/FIELD 1 DISPATCHER & NEW CODE WILL OVERLAY 'MODIFY' FUNCTION. *400 /TO 554 IS OPEN ZAPIT, ISZ BLOCK /SET FOR BUFFER 2 CALL ZAP /CLEAR THE ARRAYS DCA BLOCK DFT0, CDF CIF 10 /GO DO FFT JMP I .+1 DFT ZAP, CLRR DIFT0, CDF CIF 10 /GO DO INVERSE FFT JMP I .+1 DIFT STPT, 0 /SUB. TO SET UP NO, OF PTS. AND POEWR 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 /IS NUMBER OF PTS CDF 0 JMP I STPT KP12, 12 PW2, 4 PTS1, 3 /COMPUTE POWER SPECTRUM POWR, HEDIT1 /FACTOR= HDPR GETNO /GET POWER SCALE DOWN FACTOR FIXT TAD M1 DCA PSH /AS SHUFT COUNTER DCA BLOCK CALL SINT /INIT POINTERS,ETC.
/GET REAL AND IMG. VALUES AND SQUARE THEM. /THEN ADD THEM ALL UP FOR ALL VALUES. PLP, TAD I YIND JMS PMUL DCA TEMP TAD I AUTO JMS PMUL TAD TEMP DCA I ZIND ISZ CNTR JMP PLP DISPLAY SINT, SPINIT M1, -1 /SUB TO SQUARE AC AND SCALE IT. PMUL, 0 SPA CIA DCA MPD /ABS. VALUE TAD MPD MQL MUY MPD, 0 LSR PSH, 0 /SET UP BEFORE CALL CLA MQA /12 BIT LOW ORDER PRODUCT IN AC JMP I PMUL /ROUTINE TO SAVE OR RESTORE DATA FROM 3RD ARRAY SAVE, HEDIT1 /SAVE? HDSA DCA BLOCK /SET UP POINTERS. INITAR TAD AR3 DCA AUTO CALL QUERY /GET ANSWER,0=YES,1=NO CDF 10 SNA CLA JMP SVE /SAVE CHANNEL 1. TAD I AUTO /GET STORAGE ARRAY INTO CHANNEL 1 DCA I YIND 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. FILTER TAD CNP /DO 3 PT. FILTER CDF CIF 10 DCA I INS1 /MODIFY INSTRUCTION AND DO IT. JMP I .+1 RLSM SHAV, SMOT11 CNP, NOP SHAN, CDF CIF 10 /EXIT TO HANNING FILTER JMP I .+1 HNSM INS1, INST /ROUTINE TO OFFSET FFT DATA ONLY SETOFF, TAD FTOFST SNA CLA TAD K1000 DCA FTOFST /SET FFT DISPLAY OFFSET. DISPLAY *SCPINT+4 JMP M13+1 *M13+1 TAD FTOFST DCA TEMP1 /ADD OFFSET FOR DISPLAY OF FFT. JMP SCPINT+5 *167 FTOFST, 0 /MODS TO AVERAGER AND ITS DISPLAY *4437 JMP I PSTX *4317 STPWR *4445 CLL CLA IAC RAL *4574 PSTX, PST FSX, FS *5324 PST, TAD OFFSET CLL RAL TAD YIND DCA YIND FS, CDF 10 JMP I .+1 4440 *4475 JMP I FSX *4511 AZN /PREVENT AVERAGING TIME DATA *4252 DCA I YIND 6000 6000 6000 NOP NOP *4264 NOP /PRVENT BOXCAR INTEGRATION *4233 DCA 0 *4251 TAD 0
/PAGE 7600 0F PS-8 IS NOT SWAPPED. /ROUTINE ADDED TO AVERAGER TO GET POWER AVERAGE. *2600 STPWR, ISZ BLOCK /CLEAR AND INIT. FOR FFT ARRAYS CALL ZAPRR DCA BLOCK CDF CIF 10 CALL FPR ISZ PCTR JMP I GOMR AZN, TAD NPTS /SET UP COUNTERS,POINTERS CLL RAR DCA NPTS INITAR TAD OFFSET CLL RAL TAD YIND DCA YIND JMP I .+1 4267 GOMR, 4117 FPR, F1PR ZAPRR, CLRR /SUB. TO SET UP VARIABLES FOR POWER AVERAGE. RSTPTS, 0 CALL STPTX TAD NPTS DCA TEM3 INITAR TAD I AR3X DCA YIND HEDIT1 /FACTOR= HDPR GETNO FIXT TAD MAG /STORE SHIFT COUNTER CDF 10 DCA I SHCX DCA I YIND /CLEAR ARRAYS ISZ CNTR JMP .-2 CDF 0 JMP I RSTPTS SHCX, SHC MAG, -14 STPTX, STPT AR3X, AR3
FIELD 1 /SUB. TO DO HANNING OR 3 PT. FILTER *56 HANM, 0 JMS PSTR /SET UP POINTERS 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 TAD 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 /SUB TO SET POINTERS AND COUNTERS PSTR, 0 TAD RTP /POINTER FOR REAL +AC DCA 13 TAD N CIA DCA CNT /PT COUNTER TAD 13 DCA 14 /SET AIR 14=13 TAD 13 TAD 47 /ADD OFFSET DCA 15 /POINTER TO IMG. DCA TMP JMP I PSTR RTP, XRTAB-1
*200 /SUB TO FIND MEAN AND SUBTRACT IT FROM ARRAY. /TO GET 0 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 HIGH ORDER,TMPX IS LOW DCA TMP ISZ CNT JMP ML1 STA /HAVE SUM. SET A SHIFT COUNTER TO GET AVERAGE TAD NU 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
DFT, JMS MIDSET /ZERO MEAN THE DATA CALL DOFFT /DO FFT SKP DIFT, CALL DOIFFT /DO INVERSE FFT CALL SORT /BIT INVERT STORAGE TAD SCALE CDF CIF 0 JMP I .+1 TRET /GO PRINT SCALE FACTOR HNSM, TAD CIAC /DO HANNING SMOOTH DCA INST /SET INSTRUCTION TAD XLOCDF /SET AC SO HANNING THE IMG. JMS HANM RLSM, JMS HANM /THEN HANNING THE REAL CDF CIF 0 JMP I .+1 BEGDIS /SUB TO FFT A SCAN AND AVERAGE ITS POWER F1PR, 0 JMS MIDSET /ZERO MEAN THE DATA CALL DOFFT /DO FFT CALL SORT /BIT INVERT IT TAD CIAC DCA INST /HANNING THE IMG. TAD XLOCDF JMS HANM JMS HANM /THEN THE REALS JMS PSTR /RESET POINTERS IAC TAD 15 TAD XLOCDF DCA PUTA /POINTER TO POWER AVERAGE ARRAY TAD CNT STL RAR DCA CNT TAD SCALE /NORMALIZE USING FFT SCALE CLL RAL TAD SHC DCA PRSH SWAB /MODE B EAE
LLL, TAD PUTA DCA GETA TAD I 13 /GET REAL AND SQUARE IT JMS PMULX DST DTEM /STORE TEMP CLA TAD I 15 /GET IMG AND SQUARE IT JMS PMULX DAD /ADD REAL**2 DTEM DAD /ADD OLD SUM GETA, 0 DST /STORE IN POWER AVERAGE ARRAY PUTA, 0 CLA ISZ PUTA ISZ PUTA ISZ CNT JMP LLL SWBA /BACK TO MODE A EAE CDF CIF 0 JMP I F1PR SHC, -3 DTEM, 0;0 /SUB TO SQUARE AC AND SCALE IT DOWN PMULX, 0 SPA /MODE B EAE ON ENTRY CIA DCA MPDX /ABS. VALUE TAD MPDX MQL MUY MPDX LSR PRSH, 0 /SET UP BEFORE CALL JMP I PMULX /RETURN RESULTS IN MQ,AC MPDX, 0
/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 P<Q? JMP SWAPED /NO, HAVE ALREADY DONE THIS PAIR TAD P /YES. SWAP ORDER TAD XRLOC /FIRST SET UP SUBSCRIPT POINTERS FOR X(P) AND X(Q). DCA PR TAD Q TAD XRLOC DCA QR TAD PR TAD XLOCDF DCA PI TAD QR TAD XLOCDF DCA QI /EXCHANGE: X(P)<=X(Q) AND X(Q)<=X(P) TAD I PR /EXCHANGE REAL PARTS. GET RE[X(P)] DCA TEMPR /STORE IT. TAD I QR /GET RE[X(Q)] DCA I PR /MAKE IT RE[X(P)] TAD TEMPR /GET RE[X(P)] DCA I QR /MAKE IT RE[X(Q)] TAD I PI /EXCHANGE IMAGINARY PARTS. GET IM[X(P)] DCA TEMPR /STORE IT. TAD I QI /GET IM[X(Q)] DCA I PI /MAKE IT IM[X(P)] TAD TEMPR /GET IM[X(P)] DCA I QI /MAKE IT IM[X(Q)] SWAPED, TAD Q /IS Q=0?, IE:ARE WE AT THE TOP OF THE ARRAY SNA CLA JMP I SORTX /YES. DONE. EXIT CMA /NO, Q<=Q-1.IE: MOVE UP THE ARRAY TAD Q DCA Q JMP REVERS /GO BACK AND CONTINUE
/THIS SUBROUTINE TAKES THE INVERSE FFT (IFFT) OF THE DATA IN THE BUFFER. /IT IS ASSUMED THAT THIS DATA IS STORED IN SEQUENTIAL ORDER. /THE RESULTS ARE STORED IN BIT INVERTED ORDER. /THE ALGORITHM USED IS AS FOLLOWS: / THE NORMAL TRANSFORM IS PERFORMED, EXCEPT: / ON FETCHING THE VALUE FOR IM[W^K], WHICH IS / THE SIN(2*PI*K/N), THIS SIN VALUE IS NEGATED. / /THE REASONING FOR THIS IS AS FOLLOWS: / A WEIGHTING FACTOR OF W^(-K) IS USED IN THE IFFT / AND SINCE W^K AND W^(-K) ARE THE SAME EXCEPT THAT / THEIR IMAGINARY PARTS HAVE OPPOSITE SIGNS, IT FOLLOWS / THAT IM[W^K] SHOULD BE REPLACED BY -IM[W^K]. IFFT, 0 CLA CLL TAD CCIA /NEGATE IM[W^K]. GET CIA INSTRUCTION DCA I SGNADJ /AND PUT AT LOCATION ADJSGN. JMS I DOFFT /DO FFT. TAD CNOP /RE-INSTATE NOP AT ADJSGN FOR FFT. DCA I SGNADJ JMP I IFFT /EXIT. SGNADJ, ADJSGN /POINTER TO SIGN ADJUST INSTRUCTION. CCIA, CIA CNOP, NOP
*1000 /SIGNED SINGLE PRECISION MULTIPLY, USING THE EAE. /ENTRY: AC=MULTIPLIER, C(CALL+1)=ADDRESS OF MULTIPLICAND. EXIT:AC=PRODUCT, /AN 11 BIT SIGNED BINARY FRACTION. MULTIP, 0 CLL /AC=ARG1 (MULTIPLIER) SPA /ARG1>0? 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 $