10 COM X[103,22],M[19,19],U[19,19],Q[20],V[20],P[78] 30 COM M$[60],N$[72] 40 COM N,K,N8,K8,N9,K9,Q9,Q7,Q5,Q4,Q3,Q2,Q1 45 COM I3,I4,U9,X$[20] 46 REM: UPDATED JUNE 1, 1974 48 REM: NUMBER OF FREQUENCY POINTS RESTRICTED TO 50 50 REM:4JUN74 52 REM FEBRUARY 19,1974 VERSION BY WM.WECKER 54 REM:>SPEC 60 IF Q[1]#0 THEN 70 62 GOTO 9929 70 D2=65 72 P1=3.1416 74 P2=LOG(10) 76 T1=1 78 S1=0 80 DIM L$[72],G$[5],C[65],W[65] 82 DIM Y[505,1] 84 G$="*+0#$" 90 PRINT "DATA IS IN WHICH COLUMN"; 100 INPUT I0 180 PRINT '10,"IS THE INPUT RAW DATA (D) OR AUTOCOVARIANCES (C)"; 190 INPUT T$ 200 IF T$="D" THEN 245 210 IF T$#"C" THEN 180 220 PRINT '10,"HOW MANY OBSERVATIONS IN ORIGINAL DATA"; 230 INPUT N0 245 I1=N-Q4+1 340 N1=I1-1 350 IF T$="C" THEN 367 360 N0=I1 365 PRINT '10,N0;"ACTIVE OBSERVATIONS" 366 GOTO 370 367 PRINT '10,I1;"ACTIVE COVARIANCES" 370 PRINT '10,"HOW MANY DIFFERENT WINDOWS [1 - 5 POSSIBLE]"; 380 INPUT N2 390 IF N2<1 OR N2>5 THEN 370 400 PRINT '10,"TYPE";N2;"LAGS FOR TUKEY LAG WINDOWS"; 405 M1=0 410 MAT INPUT A[N2] 420 FOR J=1 TO N2 430 M1=M1 MAX A[J] 440 NEXT J 450 IF M1 <= N1 MIN D2 THEN 480 460 PRINT '10,"NOT ENOUGH DATA OR SPACE FOR LAG";M1 470 GOTO 370 480 IF 3*M1 <= N0 THEN 520 490 PRINT '10,"SOME LAGS ARE GREATER THAN N/3 --- DO YOU WANT" 491 PRINT "TO ANSWER LAST 2 QUESTIONS AGAIN (Y OR N)"; 500 INPUT A$ 510 IF A$="Y" THEN 370 520 REM: 530 REM: FREQUENCY POINTS RESTRICTED TO 50 540 N3=50 550 REM: 600 IF T$="D" THEN 650 610 FOR K0=1 TO M1 620 C[K0]=X[K0,I0] 630 NEXT K0 640 GOTO 800 650 X1=0 660 FOR I=Q4 TO N 670 X1=X1+X[I,I0] 680 NEXT I 690 X1=X1/N0 695 X9=X1 700 FOR I=Q4 TO N 710 X[I,I0]=X[I,I0]-X1 720 NEXT I 730 FOR K0=1 TO M1 740 L=K0-1 745 X1=0 750 FOR I=Q4 TO N-L 760 X1=X1+X[I,I0]*X[I+L,I0] 770 NEXT I 780 C[K0]=X1/N0 790 NEXT K0 791 FOR I=Q4 TO N 792 X[I,I0]=X[I,I0]+X9 793 NEXT I 800 PRINT LIN(5);"LAG";TAB(10);"AUTOCOVARIANCE";TAB(25);"AUTOCORRELATION";LIN(1) 810 FOR K0=1 TO M1 820 PRINT K0-1;TAB(10);C[K0];TAB(25);C[K0]/C[1] 830 NEXT K0 840 MAT Y=ZER[N3+1,N2] 890 PRINT LIN(5);"SPECTRAL ESTIMATES FOR TUKEY WINDOW" 900 FOR I1=1 TO N2 910 FOR K0=2 TO A[I1] 920 W[K0]=.5*(1+COS(P1*(K0-1)/A[I1]))*C[K0] 930 NEXT K0 940 FOR I=0 TO N3 950 X1=COS(P1*I/N3) 960 V0=V1=0 970 FOR K0=A[I1] TO 2 STEP -1 980 V2=2*X1*V1-V0+W[K0] 990 V0=V1 1000 V1=V2 1010 NEXT K0 1020 Y[I+1,I1]=V0=2*T1*(C[1]+2*(V1*X1-V0)) 1030 S1=S1 MAX V0 1035 NEXT I 1040 NEXT I1 1105 IMAGE"WINDOW LAG VALUE",4X,5(5D.DD,3X) 1110 MAT PRINT USING 1105;A 1115 MAT R=ZER[N2] 1120 FOR I=1 TO N2 1130 R[I]=4/(3*A[I]) 1140 NEXT I 1150 MAT PRINT USING 1151;R 1151 IMAGE"BANDWIDTH",11X,5(5D.DD,3X) 1160 IMAGE"DEGREES OF FREEDOM",2X,5(5D,6X) 1170 MAT R=(2*N0)*R 1180 MAT PRINT USING 1160;R 1183 PRINT "DO YOU WANT DUMP OF SPECTRAL NUMBERS (Y OR N)"; 1184 INPUT A$ 1185 IF A$#"Y" THEN 1200 1190 MAT PRINT Y 1200 PRINT LIN(5);"AUTOSPECTRA PLOTTED ON LOGARITHMIC SCALE" 1210 PRINT LIN(2);"WINDOW LAG VALUE PLOT CHARACTER" 1220 FOR I=1 TO N2 1230 PRINT USING 1235;A[I],G$[I,I] 1235 IMAGE 9D,16X,A 1240 NEXT I 1250 D=INT(LOG(S1)/P2+1) 1260 PRINT LIN(2);TAB(31);"POWER SCALE" 1270 PRINT "FREQ.1.0 1.8 3.2 5.6 1.0 1.8 3.2 5.6 1.0 1.8 3.2 5.6 1.0 1.8 3.2 5.6 1.0" 1280 PRINT USING 1305;D-4,D-4,D-4,D-4,D-3,D-3,D-3,D-3,D-2,D-2,D-2,D-2,D-1,D-1,D-1,D-1,D 1305 IMAGE "SCALE",17("E",SD,1X) 1310 PRINT USING 1315 1315 IMAGE 6X,"1",16("...1") 1316 D=16*D 1317 J1=5 1318 F3=1/(2*N3) 1319 F2=J=0 1320 FOR I=0 TO 50 1330 F1=I/100 1340 READ L$ 1350 RESTORE 1360 IF F2>F1+.005 THEN 1460 1370 F2=F2+F3 1380 J=J+1 1390 FOR I1=N2 TO 1 STEP -1 1400 IF Y[J,I1] <= 0 THEN 1430 1410 K0=INT(16*LOG(Y[J,I1])/P2+66.5-D) 1420 IF K0>0 THEN 1440 1430 K0=1 1440 L$[K0,K0]=G$[I1,I1] 1450 NEXT I1 1460 IF J1=5 THEN 1500 1470 PRINT USING 1475;L$ 1475 IMAGE 5X,67A 1480 J1=J1+1 1490 GOTO 1520 1500 PRINT USING 1505;F1+.001,L$ 1505 IMAGE .DD,2X,67A 1510 J1=1 1520 NEXT I 1540 DATA ". " 1550 GOTO 9998 9929 CHAIN "$IDA29",210 9998 CHAIN "$IDA",150 9999 END