Approved For Release 2009 NEL A'S SPFIETS88R001 800160001-6 (U) and the Stanford mean solar magnetic field are updated on a weekly or monthly basis as the data are available. ELF data from the "old" system are recorded twice daily; ELF data from the "new" system are being recorded at the following intervals: half-hourly averages (48 times a day), and two sets of half-daily averages (twice a day) at 00:00 and 12:00 UT, and at 04:30 and 16:30 UT. Whole-day averages are also being recorded. (U) To summarize the status of the satellite-downlink data- acquisition system: ¢ The downlink-geophysical system is in place and in operation at the Los Altos site. ¢ Partial data acquisition from the downlink is in place. ¢ Software for the final data-acquisition system will be completed and implemented in August. E. (U) Geophysical Data/RV Performance Correlation Analysis (U) Data for RV sessions are to be analyzed for statistical cor- relation with respect to the battery of geophysical data sets listed in Table 1, and those acquired by Time Research Institute via measurement (Figure 3). Recording of RV and geophysical data is now in progress. The overall system for data acquisition and analysis of RV performance/ geophysical data is depicted in Figure 4. When enough data have been collected toward the end of the contract period, analysis will be performed. The tasks described in earlier sections are in preparation for this task, and therefore have received the bulk of the effort. Certain subtasks in the analysis task, which require a longer leadtime, have, however, already been completed in preparation for the analysis. (U) The primary statistical program that will be used to scan the data for possible relationships is called EPOCH ANALYSIS. This program reads two files simultaneously. The first file is an event file, the second a data file. The program first reads an event, then scans the data temporally backward and forward in time around the event. This information is stored, a second event is read, and so forth. When all 14 Approved For Release 200b hk A & SaifPbE-Bh788R001800160001-6
Approved For Releasq PAG COBIRES SIRF PRP 00788R001 800160001-6 LOS ALTOS STATION SRI STATION — TT COATT, WESTAR IV DOWNLINK ELF DATA IBM PC JR. No. 1 IBM PC JR. No. 2 —<——- <a DATA DATA DISKETTE DISKETTE IBM XT COMPUTER BURROUGHS 86900 COMPUTER FIGURE 3. (U) REAL-TIME GEOPHYSICAL DATA-ACQUISITION SYSTEM ELF DATA IBM PC JR. No. 3 DATA DISKETTE UNCLASSIFIED (U) the events and surrounding data have been read, a printout is created that lists appropriate cross-correlation statistics between event and data elements. (U) Preliminary scans of data generated during an approximate 100-site series with one remote viewer have been carried out. The session quality was graded on a scale of 0 to 3, and correlations between solar magnetic field and solar sunspot number were investigated. Some correlation between RV performance and solar sunspot number was found, which, if substantiated by further data, would indicate the possibility that performance might improve immediately after a peak in the sunspot number, and would deteriorate just before the sunspot number peaks in its 27-to-29-day cycle. This result is based on data points that are too small in number to be taken seriously at the point, however, and is mentioned only to give an example of the types of correlations that will be sought out and examined during the course of study. 15 Approved For Releasé2ON Gab AeS SibFHibds-00788R001800160001-6
Approved For Release 200INELA'SSPFIE TS 22RO 1800160001 -6 (GEOPHYSICAL DATA) PERFORM- ANCE DATA GEO- PHYSICAL DATA BURROUGHS B6900 COMPUTER ANALYSIS REQUIRING LARGE MAINFRAME (e.g, MULTIVARIATE ANALYSIS) MANUAL ANALYSIS IBM XT COMPUTER MANAGEMENT SYSTEM STATISTICAL ANALYSIS GRAPHS/ STATISTICS EPOCH ANALYSIS STATISTICS AND GRAPHS GRAPHS/ STATISTICS UNCLASSIFIED FIGURE 4 (U) GEOPHYSICAL/PERFORMANCE DATA-ANALYSIS SYSTEM 16 Approved For Release 200bINMiok Af Sites. Bh788R001800160001-6
Approved For Release ANCLTASSIFIED 800160001-6 IV SUMMARY (VU) (U) Approximately 80 percent of the project's data-acquisition tasks have been completed, and all of the equipment and hardware have been purchased and delivered. (U) System calibration should be completed in the near future, and ELF and downlink-data acquisition will have begun in their final formats. (U) At the above point, the focus of effort will turn to analysis of past and present geophysical and ELF data, soon to be followed by the initiation of correlation studies of these data against RV performance. 17 Approved For Releasd OqQUB/98 § S| AFHIBPES -00788R001800160001-6
Approved For Release 2OPARRSE: RI RPE PE pI788R001800160001 -6 THIS PAGE INTENTIONALLY LEFT BLANK 18 Approved For Release 200)RRO AGS BHEd 0788R001800160001-6
Approved For Releasq FRCP Rs STERE ESO 007 88R001 800160001-6 Appendix A FAST FOURIER ROUTINE FOR ELF DATA (U) , FEE We eee We ee He We Wee We ee Be eee ee ee ee ee ie ee ee ee eee ee ee a 7* * 3* PROGRAM DESCRIPTION » 3* * gE A We He ee He eee We eB We ee ee eee ee eee i ee ea ee ee t ;THIS PROGRAM USES THE FAST FOURIER TRANSFORM (FFT) ALGORITHM TO yCALCULATE THE SPECTRAL MAGNITUDES OF AN ARRAY OF CLOSE-PACKED REAL INPUT ;DATA POINTS. THE PROGRAM IS ORGANIZED AS FOLLOWS: 3 1. AT THE START, DATA POINTS ARE STORED IN THE REAL ARRAY X(14), WHERE «+ RUNS FROM @ TO NX-1. NX MUST BE AN INTEGRAL POWER OF 2, OF THE FORM NX = 24MX (IE, MX = LOG2(NX)). SUCCESSIVE ELEMENTS OF X REPRESENT SUCCESSIVE SAMPLES OF AN INPUT SIGNAL, SAMPLED AT REGULAR INTERVALS OF TIME DT. 2. THE REAL ARRAY X(1) IS TREATED FOR THE FFT AS A COMPLEX ARRAY OF NX/2 PAIRS OF REAL AND IMAGINARY ELEMENTS. THAT IS, THE REAL ELEMENTS OF THE ARRAY ARE THE EVEN-NUMBERED INPUT SAMPLES, AND THE IMAGINARY ELEMENTS OF THE ARRAY ARE THE ODD-NUMBERED INPUT SAMPLES. AFTER THE FFT IS CALCULATED, AN ADDITIONAL STEP IS USED TO EXTRACT THE SPECTRUM OF THE REAL INPUT DATA. TRANSFORMING CLOSE-PACKED DATA IN THIS WAY, EVEN THOUGH IT REQUIRES AN ADDITIONAL STEP FOR THE REAL TRANSFORM EXTRACTION, IS FASTER THAN TRANSFORMING A COMPLEX ARRAY OF NX REAL-IMAGINARY PAIRS (2*NX ELEMENTS), WHERE ALL OF THE INITIAL IMAGINARY VALUES ARE ZERO. t 3. THE FFT 1S CALCULATED AS FOLLOWS: A. THE INPUT PAIRS IN X ARE REARRANGED IN BIT-REVERSED ORDER. B. THE FFT IS CALCULATED USING DECIMATION IN TIME, WITH ANGLE ARGUMENTS IN EACH SUB-DFT APPEARING IN NATURAL (IE, INCREASING) ORDER. C. AFTER THE LAST PASS THROUGH THE FFT ALGORITHM, THE ARRAY X CONTAINS SPECTRAL VALUES IN NORMAL ORDER, WITH EACH EVEN POINT A REAL VALUE AND EACH SUBSEQUENT ODD POINT THE CORRESPONDING IMAGINARY VALUE. 4. FOLLOWING THE FFT THE REAL TRANSFORM IS EXTRACTED. NOTE THAT IF THE SAMPLING INTERVAL (THE TIME BETWEEN SUCCESSIVE SAMPLES) IS DT, THEN THE TOTAL SIGNAL INTERVAL PROCESSED IS NX*DT (THE FFT ASSUMES THAT THE MISSING RIGHT END POINT IS THE SAME AS THE FIRST POINT). THE FFT GENERATES A SPECTRUM CONTAINING VALUES AT INCREMENTS OF DF IN FREQUENCY, WHERE DF = 1/(NX*DT). IF WE WERE TO TRANSFORM AN ARRAY OF NX COMPLEX DATA POINTS (WHOSE IMAGINARY VALUES WERE ZERO, SINCE WE ARE CONCERNED WITH A REAL SIGNAL) WE WOULD GENERATE NX COMPLEX SPECTRAL POINTS RANGING IN FREQUENCY FROM @ TO (NX-1)*DF Hz. HOWEVER, THE POINTS FROM (NX/2+1)*DF Hz TO (NX-1)*DF Hz ARE MERELY THE COMPLEX CONJUGATES OF THE LOWER POINTS, AND CONTAIN NO ADDITIONAL INFORMATION. THIS IS BECAUSE THE NYQUIST RATE, OR THE HIGHEST UNALIASED FREQUENCY THAT CAN BE SAMPLED, IS 1/2*NX*DF. WHEN WE TRANSFORM NX REAL CLOSE-PACKED POINTS AND THEN EXTRACT THE REAL TRANSFORM, WE GET ONLY THE LOWER NX/2 COMPLEX SPECTRAL VALUES, WHICH ARE ALL THAT ARE NEEDED. Oe wh he te ee oe Oh OE OO OO OO we wh OO OO Oe HE OF OE HO CE ee Oe ee en we oe OF Oh ee en Be we we Ot we te te se oe we ee we ee we ee we 2l Approved For ReleaswiAGids AS SAF HED 96-00788R001800160001-6
Approved For Release UNCLASSIEIE 0788R001800160001-6 + ee we we wt te ee we me ee ne es te we we en we we me ee as ee ee we ee we we we we oe 5. NEXT, THE SPECTRUM IS CONVOLVED WITH A SHORT WINDOW FUNCTION. THE REAL AND IMAGINARY (EVEN AND ODD) POINTS ARE CONVOLVED SEPARATELY, SINCE THE WINDOW FUNCTION IS A SEQUENCE OF REAL NUMBERS. WINDOWING IS NECESSARY IF THE SIGNAL FILTERS WE ARE EFFECTIVELY SYNTHESIZING ARE TO HAVE A USEFUL SHAPE. WITHOUT WINDOWING WE WOULD FIND THAT EACH SPECTRAL FILTER WOULD HAVE A NARROW PASSBAND BUT SIGNIFICANT SIDELOBE RESPONSES. THAT IS, THE MAGNITUDE OF A SPECTRAL LINE X(1) WOULD DEPEND NOT ONLY ON SIGNAL COMPONENTS NEAR 1*DF IN FREQUENCY, BUT ALSO SIGNIFICANTLY ON COMPONENTS AT OTHER FREQUENCIES AS WELL. WINDOWING BROADENS THE SHAPE OF THE PASSBAND OF EACH SPECTRAL FILTER BUT DECREASES THE SIDELOBE RESPONSES. THE AMOUNT OF BROADENING AND SUPPRESSION OF SIDELOBES DEPENDS ON THE WINDOW ORDER, OR THE LENGTH OF THE WINDOW FUNCTION WHICH IS CONVOLVED WITH THE RAW SPECTRUM. 6. FINALLY, THE MAGNITUDE OF THE SIGNAL AT EACH SPECTRAL FREQUENCY IS CALCULATED AS THE ROOT SUM OF THE SQUARES OF THE REAL AND IMAGINARY SPECTRAL COMPONENTS (IE, THE VECTOR MAGNITUDE OF EACH COMPLEX SPECTRAL POINT) AND SCALED TO MAKE] THE PROCESSING GAIN INDEPENDENT OF THE SIZE OF THE BLOCK LENGTH NX. WE COULD ALSO CALCULATE THE PHASE OF THE SPECTRAL POINTS BUT THIS INFORMATION ISN'T TOO VALUABLE FOR OUR USE. NOTE THAT THE HIGHEST MEANINGFUL SPECTRAL FREQUENCY DEPENDS ON THE ANTI- ALIASING LOW-PASS FILTER THAT 1S USED WHEN THE SIGNAL IS SAMPLED. THERE IS NOT MUCH MEANING TO SPECTRAL POINTS ABOVE THE FILTER CUTOFF FREQUENCY, WHERE SIGNALS ARE ATTENUATED AND FREQUENCY ALIASING BECOMES A PROBLEM. THUS, IN THIS PROGRAM WE DO NOT CALCULATE SPECTRAL COMPONENTS ABOVE ABOUT 3@ Hz (THE NOMINAL CUTOFF FREQUENCY OF THE LOW-PASS FILTER). sTHIS PROGRAM IS WRITTEN IN BASIC FOR THE IBM PC COMPUTER. HOWEVER, IT 31S EASILY ADAPTED TO OTHER MACHINES. FOUR THINGS TO WATCH OUT FOR WHEN s+ TRANSFERRING THE PROGRAM TO ANOTHER VERSION OF BASIC ARE: Pe ae Oe we we we we we we ee at we we we ee we He et ee ee we 1. ALL ARRAYS START WITH INDEX @. THAT IS, THE ARRAY X(NX) RUNS FROM X(@) TO XCNX-1). ARRAY REFERENCES WILL NEED TO BE CHANGED IF THE PROGRAM IS TO BE USED ON MACHINES WHERE THE: FIRST ELEMENT OF ARRAYS HAS INDEX 1. 2. ALL VARIABLES WHOSE NAMES START WITH THE LETTERS "I" THROUGH "N" ARE IMPLICITLY INTEGERS (DEFINT I-N STATEMENT). THESE VARIABLES MAY NEED TO BE RENAMED. IT IS IMPORTANT) THAT LOOP COUNTERS AND ARRAY INDEX VARIABLES BE INTEGERS BECAUSE INTEGER ARITHMETIC (WHERE AVAILABLE) IS FASTER THAN FLOATING-POINT ARITHMETIC. : 3. THIS PROGRAM USES THE INTEGER DIVISION OPERATOR "\", THIS MAY BE CHANGED TO "/" IN OTHER VERSIONS OF BASIC. 4, THIS PROGRAM USES LOGICAL OPERATORS ( "NOT", "AND", “OR") ON INTEGER VARIABLES (NOT LOGICAL VARIABLES) IN THE BIT-REVERSAL ROUTINE. IF THESE OPERATORS ARE NOT AVAILABLE THE: ALTERNATE VERSION OF STATEMENTS 3928-3198 USING ONLY INTEGER ARITHMETIC MUST BE SUBSTITUTED. 22 Approved For Release 20DPNBOS ASS1PTED °788R001800160001-6
Approved For Release 2000/08/08 : CIA-RDP96-00788R001800160001-6 UNCLASSIFIED F RahehaheReheRehehahahaRahehehehehehehshehshehehehehehehihehedehshehshebabehehahehehehshelshelsdsledelalehehehehaheheheheleheieheheiehehsleiahidel 3* * :* CHOICE OF SAMPLING TIME AND BLOCK SIZE bad * ® ge We IO ee He FW ee ee TC Ce 3 3THE SAMPLING TIME DT AND THE BLOCK SIZE NX TOGETHER DETERMINE THE sNUMBER OF SPECTRAL POINTS CALCULATED AND THEIR SPACING IN FREQUENCY. 3THE CHOICE OF SAMPLING TIME ALSO AFFECTS THE FREQUENCIES AT WHICH POWER- sLINE HARMONIC INTERFERENCE THAT IS PASSED BY THE ANTI-ALIASING LOW-PASS sFILTER WILL APPEAR IN THE SPECTRUM. THE FOLLOWING DISCUSSION ASSUMES :THAT THE 38 Hz LOW-PASS FILTER WILL BE USED TO ATTENUATE HIGH FREQUENCY ;SIGNAL COMPONENTS BEFORE SAMPLING. 3 ;IF THE INPUT SIGNAL THAT IS SAMPLED CONTAINED ONLY FREQUENCIES BELOW 338 Hz IT WOULD BE SUFFICIENT TO SAMPLE AT THE NYQUIST RATE OF 68 SAMPLES ;PER SECOND. HOWEVER, SINCE THE FILTER DOES NOT HAVE INFINITE sATTENUATION ABOVE 3@ Hz IT IS NECESSARY TO SAMPLE AT A SOMEWHAT HIGHER ;RATE, AND THEN DISCARD THOSE SPECTRAL POINTS REPRESENTING SIGNALS ABOVE 338 Hz. 3 1. ONE CONSIDERATION IN CHOOSING THE SAMPLING RATE IS THE ACTUAL RATE AVAILABLE WITH A GIVEN A/D CONVERTER. WITH THE CMC BUSSTER DI6 A/D CONVERTER, SAMPLING TIMES CAN BE AS SMALL AS @.8895 s. HOWEVER, IF MORE THAN ONE SIGNAL IS TO BE DIGITIZED AT ONE TIME, THE MINIMUM SAMPLING TIME IS 8.81 s PER CHANNEL, OR A SAMPLING RATE OF 188 SAMPLES/SECOND. THIS SEEMS LIKE A REASONABLE CHOICE. 2. A SECOND CONSIDERATION IS THE EFFECT OF THE SAMPLING RATE ON THE FREQUENCIES OF ALIASED POWER LINE HARMONICS. ONLY HARMONICS AT 68 AND 129 Hz ARE LIKELY TO BE A PROBLEM. WITH OT = 2.81 S, THE NYQUIST RATE IS 5@ Hz, SO THE SPECTRUM WILL CONTAIN POINTS AT FREQUENCIES FROM 8 TO 58 Hz. ANY 6 Hz SIGNAL THAT IS DIGITIZED WILL APPEAR IN THE SPECTRUM AT 48 Hz, WHICH POINT WILL BE THROWN OUT, SO 6@ Hz INTERFERENCE WON'T BE A PROBLEM. HOWEVER, INTERFERENCE AT 128 Hz WILL APPEAR AT 28 Hz IN THE SPECTRUM, AND THIS FACT MUST BE KEPT IN MIND WHEN ANALYZING THE DATA. (HIGHER HARMONICS APPEAR AS: 188 AT 28 Hz, 248 AT 48 Hz, 388 AT B Hz, 368 AT 48 Hz, AND SO ON.) we es ne we we we we ee ee we ee we we we we we se ee we sGIVEN A SAMPLING TIME DT = 8.81 s, WE CAN CALCULATE THE SPECTRAL SPACING ;FOR DIFFERENT CHOICES OF BLOCK SIZE NX, AND WE FIND THE FOLLOWING: + ; NX DF NUMBER OF POINTS <= 38 Hz 3 64 1.5625 Hz 19 ; 128 8.7813 Hz 38 ; 256 8.3986 Hz 76 ; 512 8.1953 Hz 153 ; 1824 8.8977 Hz 387 ; 2048 8.8488 Hz 614 Approved For Releasy PRERGEIARG & fe725}406-00788R00 800160001-6
Approved For Release "UN CL ASSIE 00788R001800160001-6 iED  Redelelehedehehseheke hehehehe iehek heheh dedehedded th tetetehdsbtehetdeheddehehed ttt tt tt ttt kt ttt tt tt8.t ttle! :* x 7* WINDOWING * 3% * ehelslietieielteh Seid t ttt td tt tt EPEC EC PEPE EP OLE ce recrers U sAFTER THE SPECTRUM OF THE INPUT SIGNAL HAS BEEN FOUND, IT 1S CONVOLVED sWITH A SHORT WINDOW SEQUENCE TO IMPROVE THE SHAPE OF THE SYNTHESIZED FILTERS. THOUGH THE SPECTRAL VALUES ARE COMPLEX, THE WINDOW SEQUENCE IS 3A SEQUENCE OF REAL NUMBERS, SOQ THE CONVOLUTION IES DONE SEPARATELY FOR sTHE REAL AND IMAGINARY SEQUENCES IN THE SPECTRUM. THE CONVOLUTION IS OF sTHE FORM 3 3 XCg) <- W3* EXC 9-3)4+XC 9439] #W2*0K( 9-2 )4X( 39 +2)) 3 #WI* CXC 3-1 I4XC 94193 tWO*XCg) 3 sWHERE X¢3) IS THE REAL OR IMAGINARY SUBSEQUENCE OF THE SPECTRUM, AND THE sCONVOLUTION IS ILLUSTRATED FOR A 4TH ORDER WINDOW SEQUENCE. NOTE THAT sIF THE CONVOLUTION IS TO BE DONE IN-PLACE, A SMALL ARRAY ‘MUST BE USED TO sHOLD THE PREVIOUS VALUES OF X(3-3), X¢g-2), X(g-1), and X()) FOR USE IN SUBSEQUENT CONVOLUTIONS. iTHE VALUES TO USE FOR DATA POINTS BEYOND THE ENDS OF THE SPECTRUM ARE sTHE COMPLEX CONJUGATES OF THE POINTS REFLECTED ABOUT THE ‘ENDS, SO THE sAUXILLIARY ARRAY OF PREVIOUS VALUES MUST BE INITIALIZED AS FOLLOWS: H ; Re(-3) = Re(3) Im(-3) = -Im(3) ; Re(-2) = Re(2) Im(-2) = -Im(2) ; Re(-1) = Re(l) Im(-1) = -Im(1) 3 sPOINTS AFTER THE FAR END OF THE SPECTRUM ARE REFLECTED IN THE SAME WAY, sBUT SINCE WE WON'T CALCULATE SPECTRAL POINTS ABOVE 39 Hz THERE IS NO iNEED TQ WINDOW POINTS THERE. ; THE FOLLOWING WINDOW SEQUENCES ARE THE MOST IMPORTANT FOR OUR USE. sTHESE ARE ALL MINIMUM-SIDELOBE WINDOWS WHICH SUPPRESS SIDELOBE RESPONSES 1AS MUCH AS POSSIBLE FOR A GIVEN WINDOW ORDER, ANO HAVE BEEN NORMALIZED sFOR UNITY GAIN. NOTE THAT THE ORDER-1 WINDOW HAMNENG PUR TO NO WINDOW c ;AT ALL. THE ORDER-2 WINDOW IS SIMILAR TO THE HAMMING FUNCTION. i ORDER 1 2 3 4 : war 1. egeee 1.eoga8 1.2999998 ° 1.9999099 ; wi g -@. 42875 ~8.5895613  -#.6727198 ; wW2 a a @.8922278 | #.1878524 ' W3 a g g -9 8146337 Hy ; HIGHEST -13 dB -43 dB -72 dB -98 dB : SIDELOBE 3 ; 3-dB WIDTH 8.89 DF 1.38 DF 1.61 DF 1.86 DF ; 6-dB WIDTH 1.28 OF 1.81 DF 2.25 DF 2.62 DF sFOR OUR USE, I THINK THE 3RD-ORDER WINDOW IS THE MOST APPROPRIATE, AND ; THE PROGRAM USES THIS WINDOW FUNCTION. SINCE W3 = @ IN THIS CASE, THE sACTUAL CONVOLUTION CONTAINS ONLY 3 TERMS. SINCE W@ = 1, THE LAST TERM 31S JUST X(3), WITHOUT A MULTIPLICATION. THUS THE WINDOWING OF EACH sPOINT INCLUDES 5 ADDITIONS AND 2 MULTIPLICATIONS FOR EACH OF THE REAL iAND IMAGINARY COMPONENTS. IF ONLY VERY SHORT BLOCKS (SAY NX = 64) ARE sTO BE PROCESSED, THE 2ND-ORDER WINDOW MIGHT BE MORE APPROPRIATE SINCE sTHE FILTERS ARE A BIT NARROWER IN THAT CASE. 3 24 Approved For Release 2009 AF CLA SSTFTEE)0788R001800160001-6
Approved For Release,290 /08 ; SPER eo teoots0001-6 UNCIAS gE ee We We ee ee We Wee We We eee We Wee ee a 3* * 3* PROGRAM * 36 * POR OIOIIIIOIIOIOIOR OOOO RI IOIGOIOIISIIOIOUIINOIO II IOUIGIOUIGI OI IOI 3 ;DEFINITIONS: 1388 DEFINT I-N "VARIABLES STARTING “I"“~"N" ARE INTEGERS 1819 NX = 1924 "DATA BLOCK SIZE 1928 NXM1 = NX-1 1938 DIM X¢€NX) "ENPUT DATA ARRAY X(€ 4), 1 © By wee, NX-1 1848 OT = 9.81 "INPUT SAMPLE TIME IN SECONDS 1858 DF = 1./(DT*NX) "SPECTRAL LINE SPACING IN Hz 1968 NMAX = 39./DF "NUMBER OF HIGHEST SPECTRAL LINE <= 38 Hz 1198 NF & NXN\2 "FFT SIZE (NR OF COMPLEX POINTS) 1118 MF = LOG(NF)/LOG(2.) "NF = 24MF 1128 NFD2 = NFN2 1138 NFM1 = NF-1 1148 DIM CACNF ) "COSINE AND SINE ARRAYS, 1159 DIM SACNF) * ANGLES FROM @ TO 18@*(NX-1)/NX 116% PI2 = 6.2831853 1229 Wl = -%.5859613 ‘WINDOW COEFFICIENTS 1218 W2 8 #.8922278 AR ee WW ee WR III III IOI IOI TOR III III Cine iin Oi wr ie ie ie ar 3 sCALCULATE SINE AND COSINE ARRAYS: 3 3THIS NEED BE DONE ONLY ONCE. NOTE THAT THE FFT ROUTINE USED SINES AND sCOSINES FROM & TO <18@ DEGREES IN STEPS OF 36@/NF DEGREES, AND THE REAL 3 TRANSFORM EXTRACTION ROUTINE USES SINES AND COSINES UP TO <9@ DEGREES IN s;STEPS OF 36@/NX OEGREES. THUS THE ODD NUMBERED POINTS IN THESE ARRAYS sFOR ANGLES ABOVE 99 DEGREES ARE NOT USED. IF MEMORY SPACE IS A PROBLEM, +THE ARRAYS CAN BE CONDENSED AND ONLY POINTS FROM @ TO £18@ DEGREES IN 3STEPS OF 368/NF CALCULATED FOR USE BY THE FFT. IN THIS CASE, THE VALUES sUSED BY THE REAL TRANSFORM EXTRACTION WILL HAVE TO BE CALCULATED AS sNEEDED. IT IS ALSO POSSIBLE TO USE A TABLE CONTAINING VALUES FOR ANGLES sONLY FROM @ TO 45 DEGREES, BUT THIS MAKES THE LOOKUP PROCEDURE SOMEWHAT sMORE COMPLICATED, INVOLVING A POSSIBLE CHANGE OF SIGN OR EXCHANGE OF SIN sAND COS VALUES DEPENDING ON THE SEMIQUADRANT OF THE ANGLE. 1522 DA = PI2/NX "ANGLE INCREMENT 1518 FOR J = # TO NFM1 "CALCULATE SINES AND COSINES 1528 A = DA*) 15398 CA(J) = COS(A) 1548 SA(J) = SINCA) 1552 NEXT J JCC IGIGCICIOIIOIOI IO IOIUIIIOIICIUIOIIIOCIOIOIUGIGI IOI OIDIOIDIOIIOISIOI III IOI IIO IO IDIGIOI I OIOIOIOIO IOI 3 sSAMPLE INPUT DATA INTO ARRAY X: ; THIS CODE MUST BE WRITTEN TO SAMPLE THE INPUT SIGNAL EVERY %.@1 SECOND sAND STUFF SUCCESSIVE VALUES INTO X. ONLY NX VALUES NEED BE TAKEN. sENTER HERE TO START PROCESSING A NEW BLOCK OF DATA. 32888 CODE YET TO BE WRITTEN 25 Approved For Releasd 2PGQUB ASS S$ [4 HEI26 -00788R001800160001-6
Approved For Release 200 NEL A'S S PEPE BP 28R001800160001-6 9 EW He eH eR He He Wee We Wee Wee eee Wee ee We We ee ee We ee 3 sBIT REVERSAL: $ 3;THE REAL COMPACT SAMPLES IN X ARE NOW TREATED AS A COMPLEX ARRAY OF NF 3(REAL, IMAGINARY) PAIRS. sPRIOR TO THE FFT ROUTINE. 3 3THE FOLLOWING ROUTINE USES LOGICAL OPERATORS. PAIRS AT BIT-REVERSED ADDRESSES ARE NOW SWAPPED ‘INITIAL BIT-REVERSED ADDRESS ‘SCAN ALL ADDRESSES EXCEPT THE FIRST ‘INITIAL BIT MASK IF NOT(K AND L) GOTO 3978 'THIS BIT SET? ' VES - CLEAR IT, ' SHIFT MASK RIGHT, ' AND CONTINUE ‘BIT WAS CLEAR, SET IT ; ‘SKIP IF ADDRESS SAME OR ALREADY DONE ‘ARRAY ADDRESS OF LOWER REAL COMPONENT sARRAY ADDRESS OF UPPER REAL COMPONENT "SWAP REALS "SWAP IMAGINARIES sTHE FOLLOWING VERSION USES ONLY INTEGER ARITHMETIC BUT IS A BIT LESS 3928 K = 1 3918 FOR J = 1 TO NFM1 3828 L = NFD2 3838 3849 K = K AND NOT L 3852 L = LN\2 3868 GOTO 3438 3878 K = K ORL 3088 IF K <= J GOTO 3198 3998 JZ 2 J+J 3198 K2 = K+K 3118 T= X(J2) 3129 X(J2) = X(K2) 3139 XCK2) = T 3148 T = X(J2+1) 3158 XCJ2+1) = XCK24+1) 3189 XCK2+1) = T 31998 NEXT J sEFFICIENT. 3328 FOR J = 1 TO NFM1 33818 Ki = J 33820 K = 8 33839 FOR I = 1 TO MF 33848 K2 = K1N\2 33958 K = (K-K2)*2+K1 33969 Ki = K2 33878 NEXT I 13889 IF K <= J GOTO 3198 33898 J2 = Jtd 33188 K2 = K+K 33118 T = X(J2) 33128 XCJ2) = X(K2) 33139 X(K2) = T 33149 T= X(J2+1) 33159 MCJ2+1) = X(K24+1) 33189 XCK2+1) = T 33198 NEXT J ‘SCAN ALL ADDRESSES EXCEPT THE FIRST “INITIALIZE DIVISOR INITIALIZE BIT-REVERSED ADDRESS sBIT-REVERSE ALL MF BITS "SKIP IF ADDRESS SAME OR ALREADY DONE "ARRAY ADDRESS OF LOWER REAL COMPONENT sARRAY ADDRESS OF UPPER REAL COMPONENT "SWAP REALS : "SWAP IMAGINARIES 26 Approved For Release 2000 M$G1_ AG SHH £49 0788R001800160001-6
Approved For Release ON CLASSIFIE! 6-00788R001800160001-6 gH MRR RIK RR WR EWR KR HWW 3 sFFT ROUTINE: ; THIS ROUTINE USED DECIMATION IN TIME. THE COMPLEX DATA INPUT IS IN BIT- sREVERSED ORDER. 1. THE OUTER LOOP MAKES MF PASSES THROUGH THE DATA. THE INITIAL NUMBER OF SUB-OFT'S IS 1, AND DOUBLES ON EACH PASS. THE INITIAL INCREMENT IN ANGLE BETWEEN SUB-DFT'S IS 188 DEGREES, AND HALVES ON EACH PASS. 2. THE MIDOLE LOOP PROCESSES THE SUB-DFT'S. AT EACH SUB-OFT ANGLES INCREASE IN NATURAL ORDER AS @ (I31)3 8, 98 (132); 9, 45, 9@, 135 (1=3)3; 8, 22.5, 45, ... €I=4)3; AND SO ON. 3. THE INNER LOOP CALCULATES THE BUTTERFLIES FOR A GIVEN ANGLE IN EACH SUB-DFT. EACH DECIMATION-IN-TIME BUTTERFLY HAS THE FORM JQ) <~ Rely d#0Re Ck *C+#Im(k >*8] 3) <- ImOg +b ImCk *C-Relk *S§] Retk) <- Re jg)-CRe(k )*C+im(k)*SI] k)} <= Im€g)-CImCk*C-ReCk )*S] ae ae ee we ae we ae we we we we we we we we ae we we WHERE C AND S ARE THE COSINE AND SINE OF THE ANGLE. 3288 LM=il ‘INITIAL NUMBER OF SUB-DFT'S = 1 3219 LM2 = LM+LM 3229 IDA = NF ‘INITIAL ANGLE INCREMENT = 189 DEG 3239 FOR I = 1 TO MF ‘QUTER LOOP 3246 IA = @ ‘INITIAL SUB-DFT ANGLE = @ DEGREES 3259 FOR L = 1 TO LM ‘MIDDLE LOOP 3268 C = CACIA) ‘LOOK UP COS AND SIN FOR THIS ANGLE 3278 S = SACTIA) 3288 FOR J=L-1 TO NFM1 STEP LM2 ; INNER LOOP 3298 JR = J+d "ARRAY ADDRESS OF LOWER REAL COMPONENT 3388 JI = JR+1 "ARRAY ADDRESS OF LOWER IMAG COMPONENT 3318 KR = JR+LM2 ‘ARRAY ADDRESS OF UPPER REAL COMPONENT 3322 KI = KR+1 "ARRAY ADDRESS OF UPPER IMAG COMPONENT 3339 TR @ XCKRI*®C+XCKID*S 'BUTTERFLY 3348 TI = XCKID*C-XCKR DNS 3358 XCKR) = XCIRI-TR 3362 XCKI) @ XCO1)-TI 3372 XCIJR) = XCIRD+TR 3388 XCJT) = XCQT+TI 3398 NEXT 9 "END OF INNER LOOP 348 IA = IA+IDA ‘BUMP ANGLE 3419 NEXT L ‘END OF MIDDLE LOOP 3349 LM = LM2 ‘DOUBLE NUMBER OF SUB-DFT'S 3439 LM2 = LM2+LM2 3448 IDA = IDAN2 "HALVE ANGLE INCREMENT 3458 NEXT I "BOTTOM OF OUTER LOOP 27 Approved For Releasd PAQQUB/ASS S [4 REI 5-00788R001800160001 -6
Approved For Release 2000 FACE: AIS SPEPE BP 728R001800160001 -6 ahelehahehehehehehshekehekehshehehsiekebeiebekehehehahehedhehahehslehslabelelkehetehtelehetehtekehelehdeh deletehakhbeletetad teh bohtol $ sEXTRACT CLOSE-PACKED REAL DATA TRANSFORM: 3 sNOW WE MUST EXTRACT THE SPECTRUM OF THE INPUT DATA AS IF IT HAD BEEN sCALCULATED FROM A SEQUENCE OF NX COMPLEX INPUT VALUES ALL OF WHOSE s IMAGINARY COMPENENTS WERE ZERO. OS Be we we ee we we we ee we ue we we te we we we we we we we oe 3599 3518 3529 3539 3549 3558 35698 3578 3589 3598 3698 3619 3629 3639 3648 3658 36698 3678 3689 3699 37988 3718 THE ALGORITHM FOR THIS IS AS FOLLOWS: 1. FOR POINTS 9 = 1, 2, «2+, NF/2-1 CALCULATE WHERE k S = SIN WE ALSO Rely) <- Im(3) <- Re(k) <- = NF ~3, (A). HAVE THE §S Re(®) <- Im(8) <- Re(NF/2) < Im(NF/2) ¢ Ra(jy)+Relk )~fRe(3)- Re(kdI*S+#lIm( yg #Im(k 7 *C Im( 3 }-ImCk )-£Re( jy )-Re(k 1 *C- CimGy d+Im(kd1*S Re(y)+Re(k }+0Re lj )-Relk 1 *S- CIm( gy }+ImekI*C Im€k) <= -ImC3)+ImCk )-CRef) )-Re(k )1*C-CImCg )+im(k 1*S A = jy*360/NX DEGREES, C PECIAL POINTS: CRe(8)+Im(9)) g ~ Re(NF/2)*2 ~ -~Im(NF/2)*2 = COS(A), AND AND (THOUGH THESE VALUES ARE NOT CALCULATED): Re(NF) <- CRe(B)-Im(H)1*2 Im(NF) IS UNKNOWN XCH) = XCHI*2, "SPECIAL CASE FOR X(@) XC1) = 8, FOR J = 1 TO NFD2-1 ‘LOOP FOR GENERAL CASE JR = J+d "ARRAY ADDRESS OF LOWER REAL COMPONENT JI = JR+l "ARRAY ADDRESS OF LOWER IMAG COMPONENT KR = NX-JR ‘ARRAY ADDRESS OF UPPER REAL COMPONENT KI = KR+l "ARRAY ADDRESS OF UPPER IMAG COMPONENT C = CACIJ) "LOOK UP COS AND SIN FOR THIS ANGLE S = SA(J} ; RS = XCIRI+X(KR) "REAL SUM RD = XCJR)-XCKR) "REAL DIFF | QS = XCJT)+X(KI) ‘IMAG SUM QD = X(JT)-X¢(KI> ‘IMAG DIFF TR = QS*C-RD*S ‘TERM FOR REALS TI = RD*C+QS*S ‘TERM FOR IMAGS XCJR) = RS+TR X(J1T) = QS-TI XCKR) = RS-TR X(KI) = -Q5-TI NEXT J "BOTTOM OF LOOP XONF) = XONF )*2, "SPECIAL CASE FOR KONE/2) | XONF+1) = XCNF+19*2. 28 Approved For Release 2009 N@LAG SHH §8)0788R001800160001-6
Approved For Releasg INCLASSIFIED 800160001-6 g RII RRR ; WINDOW SPECTRAL DATA TO IMPROVE FILTER SHAPE: ;SEE THE DISCUSSION OF WINDOWING ABOVE. THIS ROUTINE IS FOR A 3RD-ORDER sWINDOW FUNCTION. 3888 R2 2 X(4} "INITIALIZE PREVIOUS VALUE ARRAY, Re(2) 3818 Rl = X(2) *  Re(l) 3828 R@ = X(H) " Rel @) 3838 Q2 = -X(5) '  Tm(2) 3842 Ql = -X(3) ' Im) 3859 Qg = -XC1) ‘  Im(G) = B. 3369 FOR JM = @ TO NMAX*2 STEP 2 ‘WINDOW POINTS <= 39 Hz 3878 J1 = JH+1 ‘ADDRESS OF IMAG COMPONENT 38389 XC ID = CR24XCIH+4) HW24CRT+XKCIG+2) *W1tRO ‘REAL COMPONENT 3898 MCIL) = CQ24XCI14+4) *W24(QL4xXCIL+2>>*W1+QG ‘IMAG COMPONENT 3998 R2 = Rl "SHIFT PREVIOUS ARRAY 39198 Rl = R@ 3928 RO = X(J9+2) *‘ AND INSERT NEW VALUE 39398 Q2 = Q1 3949 Ql = Qs 3958 Qe = X(Jit2) 3969 NEXT J ‘END OF LOOP gH Wee WR RK RRR RR WK RR KKK RARE ; CALCULATE MAGNITUDES: ; THE MAGNITUDE OF EACH SPECTRAL VECTOR IS CALCULATED AS MAGNITUDE( 3) = SQRT[Re( ))*Re(y)+Im( 3 )*Im(0 3] *SF ; WHERE SF = 1/NX IS A SCALING FACTOR, USED FOR 2 REASONS: 1. THE SCALING FACTOR SETS THE PROCESSING GAIN. WITH SF = 1/NX, THE GAIN FOR A SINGLE~-COMPONENT SIGNAL IS 1. THAT IS, If THE INPUT SIGNAL IS A SINE-WAVE AT FREQUENCY n*DF WITH PEAK AMPLITUDE 1.8, THE SPECTRUM AT FREQUENCY n*DF WILL ALSO BE 1.28. 2. THE SCALING FACTOR MAKES THE PROCESSING GAIN INDEPENDENT OF THE BLOCK SIZE NX, SO THE SAME SIGNAL CAN BE PROCESSED WITH DIFFERENT BLOCK SIZES IN ORDER TO GENERATE DIFFERENT SPECTRA WITH VARYING FREQUENCY RESOLUTION, BUT THE MAGNITUDE OF A NARROW-BAND SIGNAL COMPONENT WILL BE THE SAME IN EACH SPECTRUM. we ee ee ee ee we wn we os se oe 3 sNOTE THAT MAGNITUDES ARE ONLY CALCULATED FOR SPECTRAL POINTS <= 38 Hz. sMAGNITUDES ARE STORED AT THE BOTTOM OF THE ARRAY X, IN THE FIRST NMAX ;POINTS. 3 sTHIS ROUTINE CAN BE INCLUDED IN THE WINDOWING ROUTINE ABOVE WITH A SMALL s INCREASE IN PROGRAM EFFICIENCY, BUT IT IS SHOWN SEPARATELY HERE FOR sCLARITY. 49298 SF = 1./NX "SCALE FACTOR 4919 FOR J = 8 TO NMAX "PROCESS POINTS <= 38 Hz 4228 J2 = J+J "ADDRESS OF REAL COMPONENT 4938 K2 = J2+1 "ADDRESS OF IMAG COMPONENT 4949 XCD) = SQRTCXCI2 "XC 92 +X K29"K(K2))*SF ' CALCULATE MAGNITUDE 4959 NEXT J ‘BOTTOM OF LOOP ;END OF PROGRAM. 3 29 Approved For Releasd 2PQQ@UB/ASS S 4 EBs -00788R001800160001-6
Approved For Release ORIEL AS ETFIED rR” 800160001-6 pI GIORGI III III III IOIIOIIUIIIIOIOUIIIOIOIO—IIIOIOIOIOIIOIOIOIOIOIIOIR OIG IID IOI On ink * * 1* TESTING THE PROGRAM * 7* * We We We We We We ee We We We We We We We Wee ee i eee eee ee eee aT i ieee ;THE BEST WAY TO TEST THE PROGRAM IS TO GENERATE A SYNTHETIC INPUT ARRAY sOF DATA POINTS CONTAINING SIGNALS AT KNOWN FREQUENCIES, AND PROCESS THEM 3AND CHECK THAT THE RESULTING SPECTRUM IS CORRECT. TO GENERATE A SIGNAL sARRAY WHICH CONTAINS A SIGNAL AT A FREQUENCY Q*DF USE TH§ FOLLOWING CODE sBEFORE THE BIT-REVERSAL ROUTINE: 2288 Q=... ‘FREQUENCY IN UNITS OF DF. 2218 OP = PI2*Q/NX "PHASE INCREMENT PER SAMPLE 2829 FOR I = @ TO NXM1 "DATA GENERATION LOOP 2839 XCIT) = SINCOP*I) “CALCULATE A DATA POINT 2848 NEXT I "END OF LOOP 3Q IS A REAL NUMBER, AND SHOULD BE CHOSEN SO THAT @ <= Q <= NX/2. NOTE sTHAT THE ARRAY X WILL CONTAIN @ CYCLES OF A SINEWAVE SIGNAL. sAFTER PROCESSING, THE SPECTRUM SHOULD CONTAIN ONLY COMPONENTS NEAR THE iFREQUENCY Q*DF. THAT IS, SPECTRAL LINES X(3) SHOULD BE SMALL EXCEPT FOR 3] NEAR Q. IF Q@ IS AN INTEGER, THE SPECTRAL LINE X(Q) SHOULD BE 1 AND sALL THE OTHERS SHOULD BE @ (OR OF ORDER 1E-6, DEPENDING ON THE ROUND-OFF sERRORS IN THE PROCESSING). USING NON-INTEGER VALUES OF Q WILL RESULT IN 3A SPECTRUM WITH A CLUSTER OF RESPONSES NEAR XC INT(Q)) BUT LINES FAR AWAY sFROM Q SHOULD BE SMALL. THE FALLOFF IN RESPONSE IN THE SPECTRUM AS ONE sMOVES AWAY FROM THE LINE X¢€INT(Q)) DEPENDS ON THE EFFECTIVENESS OF THE sWINDOW FUNCTION, AND THE SHAPE OF THE SYNTHESIZED FILTER RESPONSE CAN BE sPLOTTED BY USING VARIOUS VALUES OF Q. 5 sNOTE THAT IF NMAX ¢ Q@ <= NX/2, THE SIGNAL GENERATED LIES IN THE UPPER sPART OF THE SPECTRUM, WHICH IS THROWN AWAY (DON'T EXAMINE SPECTRAL LINES 3X(j) WHERE 3 > NMAX, AS THEY HAVEN'T BEEN CALCULATED). ALSO, IF sNX/2 ¢ Q ¢ NX, THE SIGNAL WILL BE ALIASED IN THE SPECTRUM, APPEARING AT 3A FREQUENCY OF NX/2~-Q. ; : OAUIOA III IOIOIOIOIIOIIICI III IOI IOI IOI IOI ICIOIGIO I IGIOIOIOIIOIOR GIO IO IOI IOI I IOI IOI Oia 30 Approved For Release 209 A804 A'S S¥F] 2{)0788R001800160001-6