home *** CD-ROM | disk | FTP | other *** search
/ HAM Radio 1 / HamRadio.cdr / tech / design2 / sabin.asc < prev    next >
Encoding:
Text File  |  1986-10-29  |  11.5 KB  |  289 lines
  1. 10 PRINT:PRINT "WAVEFORM ANALYSIS:"
  2. 20 PRINT "W.E. SABIN - EDN JUNE 1983 - PAGE 243"
  3. 30 INPUT "ENTER EXPONENT M ";M
  4. 40 N=2^M: PI=3.1415928159#:E1=2
  5. 50 DIM X(N+1,5)
  6. 60 PRINT: INPUT "1 TO READ A DATA FILE, 2 TO GO ON";A
  7. 70 ON A GOTO 80,110
  8. 80 INPUT "FILE TO READ";F$
  9. 81 OPEN "I",1,F$
  10. 82 FOR I=1 TO N
  11. 83 INPUT #1,X(I,0)
  12. 84 NEXT I
  13. 85 CLOSE
  14. 100 PRINT "DATA LOADED":PRINT:GOTO 280
  15. 110 REM ****************************************************************
  16. 120 REM * X(I,0) REAL, X(I,1) IMAGINARY                                *
  17. 130 REM * EVALUATE REAL, IMAGINARY IN LINES 190 - 269                  *
  18. 140 REM * FOR AUTOCORRELATION, AUTOSPEC. USE X(I,0), X(I,1)            *
  19. 150 REM * FOR CROSS SPECTRUM, CROSS CORRELATION                        *
  20. 160 REM * AND CONVOLUTION,                                             *
  21. 170 REM * USE X(I,0),X(I,1) AND X(I,2),X(I,3)                          *
  22. 180 REM ****************************************************************
  23. 190 REM
  24. 200 FOR I=0 TO N
  25. 210 IF I<N/2 THEN X(I,0)=1 ELSE X(I,0)=0
  26. 220 X(I,1)=0
  27. 230 NEXT I
  28. 240 REM
  29. 250 REM
  30. 260 REM
  31. 270 PRINT
  32. 280 PRINT: PRINT "ITEMS 1-7 FOR X(I,0), X(I,1) ONLY":PRINT
  33. 290 PRINT "TYPE 1  FOR FORWARD TRANSFORM"
  34. 300 PRINT "TYPE 2  FOR INVERSE TRANSFORM"
  35. 310 PRINT "TYPE 3  TO LIST X REAL, X IMAG"
  36. 320 PRINT "TYPE 4  FOR SINE/COSINE"
  37. 330 PRINT "TYPE 5  FOR MAG AND PHASE"
  38. 340 PRINT "TYPE 6  FOR SMOOTHING"
  39. 350 PRINT "TYPE 7  FOR WINDOWING"
  40. 360 PRINT "TYPE 8  FOR POWER SPECTRUM"
  41. 370 PRINT "TYPE 9  FOR CORRELATION"
  42. 380 PRINT "TYPE 10 FOR CONVOLUTION"
  43. 390 PRINT "TYPE 11 FOR TO SAVE DATA IN FILE"
  44. 400 PRINT "TYPE 12 TO END"
  45. 410 PRINT "TYPE 13 TO EXCHANGE SEQ 1 & 2"
  46. 420 PRINT "TYPE 14 FOR DEQSEQ(1)=SEQ(1)*SEQ(2)"
  47. 430 INPUT X
  48. 440 ON X GOTO 470,1410,500,600,600,1450,1690,1790,2000,2310,1230,2770,460,450
  49. 450 FOR I=1 TO N:X(I,4)=X(I,0)*X(I,2)-X(I,1)*X(I,3)451 X(I,5)=X(I,0)*X(I,3)+X(I,2)*X(I,1):X(I,0)=X(I,4)452 X(I,1)=X(I,5):NEXT:GOTO 270
  50. 460 FOR I=1 TO N:X(I,4)=X(I,0):X(I,5)=X(I,1):X(I,0)=X(I,2):X(I,1)=X(I,3)461 X(I,2)=X(I,4):X(I,3)=X(I,5):NEXT:GOTO 270
  51. 470 REM ******** COMPUTE FORWARD TRANSFORM **********
  52. 480 D=0:GOSUB 2480
  53. 490 GOTO 270
  54. 500 REM ******** DISPLAY EXPONENTIALS **************
  55. 510 PRINT:PRINT "N";TAB(6);"XR";TAB(25);"XI":PRINT
  56. 520 FOR J=O0 TO N STEP 10:PRINT
  57. 530 FOR K=1 TO 10
  58. 540 IF J+K>N GOTO 570
  59. 550 G=(2-E1)*N/2
  60. 560 PRINT J+K-G-1;TAB(6)X(J+K,0);TAB(25)X(J+K,1)
  61. 570 NEXT K
  62. 580 PRINT: PRINT "PRESS 'Q' TO QUIT, ANY OTHER KEY TO CONTINUE"
  63. 581 D$=INPUT$(1): IF D$="Q" THEN GOTO 270
  64. 582 PRINT:NEXT J
  65. 590 GOTO 270
  66. 600 REM ********** CALCULATE SINE, COSINE ****************
  67. 610 D=0:GOSUB 2480
  68. 620 X(1,2)=X(1,0):X(1,4)=X(1,1)
  69. 630 FOR I=2 TO N/2
  70. 640 X(I,2)=X(I,0)+X(N+2-I,0)
  71. 650 X(I,3)=X(N+2-I,1)-X(I,1)
  72. 660 X(I,4)=X(N+2-I,1)+X(I,1)
  73. 670 X(I,5)=X(I,0)-X(N+2-I,0)
  74. 680 NEXT I
  75. 690 REM ******* PRINT SIN/COS OF FIND MAG ***********
  76. 700 IF X=5 GOTO 880
  77. 710 PRINT:PRINT "REAL PART SINE, COSINE":PRINT
  78. 720 FOR J=0 TO N/2 STEP 10
  79. 730 PRINT "N" TAB(6)"COSINE" TAB(25)"SINE":PRINT
  80. 740 FOR K=1 TO 10
  81. 750 IF J+K>N/2 GOTO 270
  82. 760 PRINT J+K-1 TAB(6)X(J+K,2) TAB(25)X(J+K,3) 
  83. 770 NEXT K
  84. 780 GET E$:PRINT:NEXT J
  85. 790 PRINT:PRINT "IMAG PART, SINE, COSINE:P":PRINT
  86. 800 FOR J=0 TO N/2 STEP 10
  87. 810 PRINT "N" TAB(6)"COSINE"TAB(325)"SINE":PRINT
  88. 820 FOR K=1 TO 10
  89. 830 IF J+K>N/2 GOTO 850
  90. 840 PRINT J+K-1 TAB(6)X(J+K),4)TAB(25,)X(J+K,5)
  91. 850 NEXT K
  92. 860 PRINT:INPUT "PRESS <RET> TO CONTINUE:";D$:PRINT:NEXT J:GOTO 270
  93. 870 REM ******** EMD SINE COSINE ROUTINE **********
  94. 880 REM ****** ** START AMPLITUDE, PHASE  **********
  95. 890 FOR I=1 TO N/2
  96. 900 V=SQR(X(I,2)^2+X(I,3)^2)
  97. 910 EW=DQSQR(X(I,4)^2+X(I,5)^2)
  98. 920 IF ABS(X(I,2))<1E-12 THEN X(I,2)=1E-12
  99. 930 IF ABS(X(I,4))<1E-12 THEN X(I,4)=1E-12
  100. 940 Y=ATN(X(I,3)/X(I,2))*57.2958
  101. 950 IF X(I,2)<0 AND X(I,3)>0 THEN Y=Y+180
  102. 960 IF X(I,2)<0 AND X(I,3)<0 THEN Y=Y-180
  103. 970 X(I,2)=V:X(I,3)=Y
  104. 980 Y=ATN(X(I,5)/X(I,4))*57.2958
  105. 990 IF X(I,4)<0 AND X(I,5)>0 THEN Y=Y+180
  106. 1000 IF X(I,4)<0 AND X(I,5)<0 THEN Y=Y-180
  107. 1010 X(I,4)=W:X(I,5)=Y
  108. 1020 NEXT I
  109. 1030 PRINT CHR$(26):PRINT:PRINT "MAGNITUDE AND PHASE":PRINT
  110. 1040 PRINT "REAL PART":PRINT
  111. 1050 FOR J=0 TO N/2 STEP 10
  112. 1060 PRINT "N" TAB(6)"MAG" TAB(25) "PHASE":PRINT
  113. 1070 FOR K=1 TO 10
  114. 1080 IF J+K>N/2 GOTO 1100
  115. 1090 PRINT J+K-1 TAB(6)X(J+K,2) TAB(25)X(J+K,3)
  116. 1100 NEXT K
  117. 1110 PRINT:INPUT "PRESS <RET> TO CONTINUE";D$:PRINT:NEXT J
  118. 1120 PRINT: PRINT "IMAGINARY PART::":PRINT
  119. 1130 FOR J=0 TO N/2 STEP 10
  120. 1140 PRINT "N" TAB(6) "MAG" TAB(25) "PHASE":PRINT
  121. 1150 FOR K=1 TO 10
  122. 1160 IF J+K>N/2 GOTO 1180
  123. 1170 PRINT J+K-1 TAB(6) X(J+K,4) TAB(25) X(J+K,5)
  124. 1180 NEXT K
  125. 1190 PRINT:INPUT "PRESS <RET> TO CONTINUE";D$:PRINT:NEXT J
  126. 1200 GOTO 270
  127. 1210 REM ******* END MAG,PHASE ******
  128. 1220 REM ****** OUTPUT DATA FILE ****
  129. 1230 PRINT CHR$(26):PRINT:PRINT "INSTRUCTIONS TO SAVE DATA IN DATA FILE"
  130. 1240 PRINT "SAVE X OR MAGNITUDE OR PHASE":PRINT
  131. 1250 PRINT  "X IS X(N),X(K) , SPEC.,CONV.,CORR.":PRINT
  132. 1260 PRINT "0=X(I,0) REAL"
  133. 1270 PRINT "1=X(I,1) IMAGINARY"
  134. 1280 PRINT "2=MAGNITUDE, REAL PART"
  135. 1290 PRINT "3=PHASE, REAL PART"
  136. 1300 PRINT "4=MAGNITUDE, IMAGINARY PART"
  137. 1310 PRINT "5=PHASE, IMAGINARY PART"
  138. 1320 INPUT R
  139. 1330 Q=1
  140. 1340 IF R>1 THEN Q=2
  141. 1350 DIM A(N/Q)
  142. 1360 FOR I=1 TO N/Q:A(I)=X(I,R):NEXT
  143. 1370 INPUT "NAME OF FILE TO SAVE DATA";F$
  144. 1371 OPEN "O",1,F$
  145. 1372 FOR I=1 TO N/Q
  146. 1373 PRINT #1,A(I)
  147. 1374 NEXT I
  148. 1375 CLOSE #1
  149. 1380 PRINT CHR$(13):PRINT "DATA FILED"
  150. 1390 GOTO 270
  151. 1400 REM ********* END DATA OUTPUT ************
  152. 1410 REM ********* INVERSE DFS     ************
  153. 1420 D=1:GOSUB 2480
  154. 1430 GOTO 270
  155. 1440 REM ******** SMOOTHING ************
  156. 1450 PRINT CHR$(26):PRINT:PRINT "SEQUENCE SMOOTHING":PRINT
  157. 1460 PRINT "TYPE 1 FOR LOW PASS"
  158. 1470 PRINT "TYPE 2 FOR HIGH PASS"
  159. 1480 INPUT Z
  160. 1490 ON Z GOTO 1500, 1590
  161. 1500 X(1,5)=.25*X(N,0)+.5*X(I1,0)+.25*X(2,0)
  162. 1510 X(N,5)=.25*X(N-1,0)+.5*X(1,0)+.25*X(2,0)
  163. 1520 FOR J=2 TO N-1:X(J,5)=.25*X(J-1,0)+.5*X(J,0)+.25*X(J+1,0):NEXT
  164. 1530 FOR J=1 TO N:X(J,0)=X(J,5):NEXT
  165. 1540 X(1,T5)=.25*X(N,1)+.5*JX(1,1)+.25*X(2,1)
  166. 1550 X(N,5)=.25*X(N-1,1)+.5*X(N,1)+.25*X(1,1)
  167. 1560 FOR J=2 TO N-1:X(J,5)=.25*X(J-1,0)+.5*X(J,1)+.25*X(J+1,1):NEXT J
  168. 1570 FOR J=1 TO N:X(J,1)=X(J,5):NEXT J
  169. 1580 GOTO 270
  170. 1590 X(1,5)=-.25*X(N,0)+.5*X(1,0)-.25*X(2,0)
  171. 1600 X(N,5)=-.25*X(N-1,0)+.5*X(N,0)-.25*X(1,0)
  172. 1610 FOR J=2 TO N-1:X(J,5)=-.25*X(J-1,1)+.5*X(J,0)-.25*X(J+1,0):NEXT
  173. 1620 FOR J=1 TO N:X(J,0)=X(J,5):NEXT J
  174. 1630 X(1,5)=-.25*X(N,1)+.5*X(1,1)-.25*X(2,1)
  175. 1640 X(N,5)=-.25*X(N-1,1)+.5*X(N,1)-.25*X(2,1)
  176. 1650 FOR J=2 TO N-1:X(J,5)=-.25*X(J-1,1)+.5*X(J,1)-.25*X(J+1,1):NEXT J
  177. 1660 FOR J=1 TO N:X(J,1)=X(J,5):NEXT J
  178. 1670 GOTO 270
  179. 1680 REM ******** WINDOWING ************
  180. 1690 PRINT CHR$(26):PRINT:PRINT: PRINT "SEQUENCE WINDOW":PRINT
  181. 1700 PRINT "TYPE 1 FOR HANNING"
  182. 1710 PRINT "TYPE 2 FOR HAMMING"
  183. 1720 INPUT Q
  184. 1730 IF Q=2 THEN Q=.8519
  185. 1740 FOR I=1 TO N
  186. 1750 X(I,0)=X(I,0)*(1-Q*COS(2*PI*(I-1)/N))
  187. 1760 X(I,1)=X(I,1)*(1-Q*COS (2*PI*(I-1)/N))
  188. 1770 NEXT I
  189. 1780 GOTO 270
  190. 1790 REM ******** POWER SPECTRUM ***********
  191. 1800 REM ******* AUTO SPEC. USE X(I,0),X(I,1) ***********
  192. 1810 REM ****** CROSS SPEC. USE X(I,0),X(I,1) AND X(I,2),X(I,3) ********
  193. 1820 PRINT CHR$(26):PRINT:PRINT "TWO-SIDED POWER SPECTRUM":PRINT
  194. 1830 PRINT "TYPE 1 FOR SPECTRUM OF SEQUENCE ONE"
  195. 1840 PRINT "TYPE 2 FOR SPECTRUM OF SEQUENCE TWO"
  196. 1850 PRINT "TYPE 3 FOR CROSS SPECTRUM"
  197. 1860 INPUT F:PRINT
  198. 1870 ON F GOTO 1890,1880,1920
  199. 1880 FOR I=1 TO N: X(I,0)=X(I,2):X(I,1)=X(I,3):NEXT
  200. 1890 D=0:GOSUB 2480
  201. 1900 FOR I=1 TO N:X(I,0)=X(I,0)*X(I,0)+X(I,1)*X(I,1):X(I,1)=0:NEXT
  202. 1910 PRINT:PRINT "AUTOSPECTRUM":PRINT:GOTO 1980
  203. 1920 FOR I=1 TO N:X(I,4)=X(I,0):X(I,5)=X(I,1):X(I,0)=X(I,2):X(I,1)=X(I,3):NEXT
  204. 1930 D=0:GOSUB 2480
  205. 1940 FOR I=1 TO N:X(I,2)=X(I,0):X(I,3)=X(I,1):X(I,0)=X(I,4):X(I,1)=X(I,5):NEXT
  206. 1950 D=0:GOSUB 2480
  207. 1960 FOR I=1 TO N:X(I,4)=X(I,0)*X(I,2)+X(I,1)*X(I,3)1961 X(I,5)=X(I,1)*X(I,2)-X(I,0)*X(I,3):X(I,0)=X(I,4):X(I,1)=X(I,5):NEXT
  208. 1970 PRINT:PRINT "CROSS SPECTRUM":PRINT
  209. 1980 PRINT "TYPE 3 TO LIST POWER SPECTRUM":PRINT:GOTO 290
  210. 1990 REM ******* CORRELATION *************
  211. 2000 PRINT CHR$(26):PRINT:PRINT "CORRELATION"
  212. 2010 PRINT: PRINT "FOR LINEAR CORRELATION, DOUBLE THE VALUE OF N AND FILL IN ";"ZEROS FROM N/2 TO N-1 IN X(N) BEFORE PROCEEDING"
  213. 2020 PRINT:PRINT "TYPE 1 FOR AUTOCORRELATION OF SEQUENCE X(N) IN X(I,0),X(I,1)"
  214. 2030 PRINT "TYPE 2 FOR AUTOCORRELATION OF SEQUENCE X(N) IN X(I,2),X(I,3)."
  215. 2040 PRINT "TYPE 3 FOR CROSS-CORRELATION"
  216. 2050 INPUT C
  217. 2060 PRINT:PRINT "TYPE 1 FOR CORRELATION:"
  218. 2070 PRINT "TYPE 2 FOR COVARIANCE."
  219. 2080 INPUT Q
  220. 2090 PRINT:PRINT "TYPE 1 IF LINEAR"
  221. 2100 PRINT "TYPE 2 IF CIRCULAR"
  222. 2110 INPUT E1
  223. 2120 IF Q=2 THEN GOSUB 2280
  224. 2130 ON C GOTO 2150,2140,2160
  225. 2140 FOR I=1 TO N:X(I,0)=X(I,2):X(I,1)=X(I,3):NEXT
  226. 2150 D=0:GOSUB 2480:GOSUB 2200:D=1:GOSUB 2480:GOSUB 2210:GOTO 2170
  227. 2160 D=0:GOSUB 2480:GOSUB 2260:D=0:GOSUB 2480:GOSUB 2270:D=1:GOSUB 2480: GOSUB 2210
  228. 2170 PRINT: ON Q GOTO 2180,2190
  229. 2180 PRINT "TYPE 3 TO LIST CORRELATION":PRINT:GOTO 290
  230. 2190 PRINT "TYPE 3 TO LIST COVARIANCE":PRINT:GOTO 290
  231. 2200 FOR I=1 TO N:X(I,0)=X(I,0)^2+X(I,1)^2:X(I,1)=0:NEXT
  232. 2210 IF E1=2 GOTO 2250
  233. 2220 FOR I=1 TO N:X(I,0)=X(I,0)*2:X(I,1)=X(I,1)*2:NEXT
  234. 2230 FOR I=1 TO N/2:X(I+N/2,4)=X(I,0):X(I,4)=X(I+N/2,0):X(I+N/2,5)=X(I,1):X(I,5)=X(I+N/2,1):NEXT
  235. 2240 FOR I=1 TO N:X(I,0)=X(I,4):X(I,1)=X(I,5):NEXT
  236. 2250 RETURN
  237. 2260 FOR I=1 TO N:X(I,4)=X(I,0):X(I,5)=X(I,1):X(I,0)=X(I,2):X(I,1)=X(I,3):NEXT:RETURN
  238. 2270 FOR I=1 TO N:X(I,2)=X(I,0)*X(I,4)+X(I,1)*X(I,5):X(I,3)=X(I,0)*X(I,5)-X(I,1)*X(I,4):X(I,0)=X(I,2):X(I,1)=X(I,3):NEXT:RETURN
  239. 2280 U=N/(3-E1):AA=0:BB=0:CC=0:DD=0
  240. 2290 FOR I=1 TO U:AA=AA+X(I,0):BB=BB+X(I,1):CC=CC+X(I,2):DD=DD+X(I,3):NEXT
  241. 2300 FOR I=1 TO U:X(I,0)=X(I,0)-AA/U:X(I,1)=X(I,1)-BB/U:X(I,2)=X(I,2)-CC/U:X(I,3)=X(I,3)=-DD/U:NEXT:RETURN
  242. 2310 REM ********* CONVOLUTION ************
  243. 2320 PRINT CHR$(26):PRINT:PRINT "CONVOLUTION":PRINT
  244. 2330 PRINT "SEQUENCE 1 IN X(I,0),X(I,1)"
  245. 2340 PRINT "SEQUENCE 2 IN X(I,3),X(I,4)":PRINT
  246. 2350 PRINT "FOR LINEAR CONVOLUTION, DOUBLE THE VALUE OF N AND ";"ARGUMENT WITH ZEROS IN BOTH SEQUENCES"
  247. 2360 PRINT:PRINT "TYPE 1 IF LINEAR"
  248. 2370 PRINT "TYPE 2 IF CIRCULAR"
  249. 2380 INPUT QQ
  250. 2390 D=0:GOSUB 2480:GOSUB 2440:GOSUB 2480:GOSUB 2450:D=1:GOSUB 2480:GOSUB 2460
  251. 2400 PRINT:PRINT TYPE 1 TO "TYPE 1 TO MULTIPLY FOR N"
  252. 2410 GET A$:PRINT A$
  253. 2420 IF A$="1" THEN: FOR I=1 TO N:X(I,0)=X(I,0)*N:X(I,1)=X(I,1)*N:NEXT
  254. 2430 PRINT:PRINT "TYPE 3 TO LIST CONVOLUTION OF X(1,N) AND X2(N).":PRINT: GOTO 290
  255. 2440 FOR I=1 TO N:X(I,4)=X(I,0):X(I,5)=X(I,1):X(I,0)=X(I,2):X(I,1)=X(I,3):NEXT:RETURN
  256. 2450 FOR I=1 TO N:X(I,2)=X(I,0)*X(I,4)-X(I,1)*X(I,5):X(I,3)=X(I,0)*X(I,5)+X(I,4)*X(I,1):X(I,0)=X(I,2):X(I,1)=X(I,3):NEXT:RETURN
  257. 2460 IF QQ=1 THEN: FOR I=1 TO N:X(I,0)=2*X(I,0):X(I,1)=2*X(I,1):NEXT:RETURN
  258. 2470 RETURN
  259. 2480 REM ***** FFT ROUTINE, COMPLEX DATA ARRAY ******
  260. 2490 REM ****** X(I,0) REAL, X(I,1) IMAGINARY  ******
  261. 2500 REM ****** D=0, FORWARD. D-=1, REVERSE **********
  262. 2510 N2=N/2:N1=N-1:J=1
  263. 2520 FOR I=1 TO N1
  264. 2530 IF I>=J THEN 2550
  265. 2540 T1=X(J,0):T2=X(J,1):X(J,0)=X(I,0):X(J,1)=X(I,1):X(I,0)=T1:X(I,1)=T2
  266. 2550 K=N2
  267. 2560 IF K>=J THEN 2590
  268. 2570 J=J-K:K=K/2
  269. 2580 GOTO 2560
  270. 2590 J=J+K
  271. 2600 NEXT I
  272. 2610 S1=-1
  273. 2620 IF D=0 THEN 2640
  274. 2630 S1=1
  275. 2640 FOR L=1 TO M
  276. 2650 L1=2^L:L2=L1/2:U1=1:U2=0:W1=COS(PI/L2):W2=S1*SIN(PI/L2)
  277. 2660 FOR J=1 TO L2
  278. 2670 FOR I=J TO N STEP L1
  279. 2680 I1=I+L2
  280. 2690 V1=X(I1,0)*U1-X(I1,1)*U2:V2=X(I1,1)*U1+X(I1,0)*U2
  281. 2700 X(I1,0)=X(I,0)-V1:X(I1,1)=X(I,1)-V2:X(I,0)=X(I,0)+V1:X(I,1)=X(I,1)+V2 
  282. 2710 NEXT I
  283. 2720 U3=U1:U4=U2:U1=U3*W1-U4*W2:U2=U4*W1+U3*W2
  284. 2730 NEXT J,L
  285. 2740 IF D=1 THEN 2760
  286. 2750 FOR I=1 TO N:X(I,0)=X(I,0)/N:X(I,1)=X(I,1)/N:NEXT
  287. 2760 RETURN
  288. 2770 PRINT:PRINT "END":END
  289.