home *** CD-ROM | disk | FTP | other *** search
/ Power-Programmierung / CD1.mdf / basic / library / trubasic / chaos / xfft.tru < prev    next >
Encoding:
Text File  |  1980-01-01  |  8.4 KB  |  275 lines
  1. 1000 !PROGRAM PEND/FFT.TRU
  2. 1001 CLEAR
  3. 1002 PRINT"                ***PENDULUM - FAST FOURIER TRANSFORM OF VARIABLES***"
  4. 1003 PRINT
  5. 1004 PRINT"THIS PROGRAM PROVIDES A PHASE DIAGRAM, TIME SERIES, AND FFT OF THE"
  6. 1005 PRINT"ANGLE OR ANG. VELOCITY OF THE PENDULUM.  A BLINKING CURSOR INDICATES"
  7. 1006 PRINT"THAT THE PROGRAM IS READY FOR THE NEXT STEP.  PRESS ANY KEY TO CONTNUE."
  8. 1007 PRINT"THE HANNING OPTION IS USED TO SMOOTH THE EFFECT OF THE ABRUPT WINDOW"
  9. 1008 PRINT"AND IS RECOMMENDED IN MOST CASES."
  10. 1009 PRINT
  11. 1010 LIBRARY "SGLIB.TRC"
  12. 1020 !
  13. 1030 DIM thetadata(5000),thetadotdata(5000),xreal(0 to 5000),ximag(0 to 10000)
  14. 1040 DIM tpoint(0 to 5000),power(2048),frequency(2048)
  15. 1050 DECLARE DEF accel
  16. 1060 DECLARE DEF bitr
  17. 1070 INPUT prompt"Max frequency (try 0.5): ":maxfreq
  18. 1080 INPUT prompt"Input driving force strength: ":g
  19. 1090 INPUT prompt"Input damping term:":q
  20. 1100 INPUT prompt"Input initial position: ": xint
  21. 1110 INPUT prompt"Input initial velocity: ": vint
  22. 1120 INPUT Prompt"Input min.time:":tmin
  23. 1130 INPUT prompt"No. of FFT points(..256,512,1024,2048,4096) : ":number
  24. 1140 PRINT"Desired power spectrum quantity"
  25. 1150 PRINT"  1)Power spectrum of angular velocity"
  26. 1160 PRINT"  2)Power spectrum of angle"
  27. 1170 INPUT prompt"Choose 1 or 2 :":ps
  28. 1180 
  29. 1190 LET del=.5/maxfreq
  30. 1200 LET tmax=number*del+tmin
  31. 1210 LET w=0.6666667
  32. 1220 LET eps=1.0e-6
  33. 1230 LET tstep=0.5
  34. 1240 
  35. 1250 LET t=0
  36. 1260 LET x=xint
  37. 1270 LET v=vint
  38. 1280 LET sumvel=0
  39. 1290 LET count=0
  40. 1300 LET p=1
  41. 1310 CALL SETTEXT("PENDULUM PHASE DIAGRAM","ANGLE","ANGULAR VELOCITY")
  42. 1320 CALL RESERVELEGEND
  43. 1330 PRINT"computing data"
  44. 1340 !
  45. 1350 FOR i=1 to 10000000
  46. 1360     CALL rk4(x,v,tstep,xnew,vnew,t,w,g,q)   ! Take a step of size tstep
  47. 1370     LET tshalf=tstep/2
  48. 1380     CALL rk4(x,v,tshalf,xnh,vnh,t,w,g,q)    !Take two half steps
  49. 1390     CALL rk4(xnh,vnh,tshalf,xn,vn,t+tshalf,w,g,q)
  50. 1400     LET d1=abs(xn-xnew)
  51. 1410     LET d2=abs(vn-vnew)
  52. 1420     LET delta=max(d1,d2)
  53. 1430     IF delta<eps then
  54. 1440        IF t>tmin then
  55. 1450           LET tnew=t + tstep
  56. 1460           LET w1=mod(-t,del)
  57. 1470           LET w2=mod(tnew,del)
  58. 1480           IF w1<tstep then
  59. 1490              IF w2<tstep then
  60. 1500                 LET ts=w1
  61. 1510                 CALL rk4(x,v,ts,xp,vp,t,w,g,q)
  62. 1520                 IF abs(xp)>pi then LET xp=xp-2*pi*abs(x)/x
  63. 1530                 LET thetadata(p)=xp
  64. 1540                 LET thetadotdata(p)=vp
  65. 1550                 LET p=p+1
  66. 1560                 LET sumvel= sumvel + v
  67. 1570                 LET count = count + 1
  68. 1580              END IF
  69. 1590           END IF
  70. 1600        END IF
  71. 1610        LET x=xnew
  72. 1620        LET v=vnew
  73. 1630        LET t=t+tstep
  74. 1640        LET tstep=tstep*.95*(eps/delta)^.25
  75. 1650        IF abs(x)>pi then
  76. 1660           LET x=x-2*pi*abs(x)/x
  77. 1670        END IF
  78. 1680     ELSE
  79. 1690        LET tstep=tstep*.95*(eps/delta)^.2
  80. 1700     END IF
  81. 1710     IF t>tmax then LET i=10000001
  82. 1720 NEXT i
  83. 1730 LET n=p-1
  84. 1740 LET meanvel=sumvel/count
  85. 1750 CLEAR
  86. 1760 CALL setgraphtype("xy")
  87. 1770 CALL datagraph(thetadata,thetadotdata,1,0,"white")
  88. 1780 CALL ADDLEGEND("G="&STR$(G)&"   Q="&str$(Q),0,1,"WHITE")
  89. 1790 CALL DRAWLEGEND
  90. 1800 GET KEY keyvariable
  91. 1810 !
  92. 1820 !
  93. 1830 !PREPARATION OF THE FFT DATA
  94. 1840 CLEAR
  95. 1850 INPUT prompt" HANNING OPTION Y/N? ": hanning$
  96. 1860 LET tgamma=log2(n)
  97. 1870 IF abs(int(tgamma)-tgamma)=0 then
  98. 1880    LET gamma=tgamma
  99. 1890    GOTO 1920
  100. 1900 END IF
  101. 1910 LET gamma=int(tgamma)+1
  102. 1920 PRINT "gamma= ";gamma
  103. 1930 LET newn=2^gamma
  104. 1940 LET nu=gamma
  105. 1950 FOR i=n+1 to newn
  106. 1960     LET xreal(i)=0
  107. 1970 NEXT i
  108. 1980 LET n=newn
  109. 1990 PRINT"n=";n
  110. 2000 CLEAR
  111. 2010 IF ps=1 then LET title$="ANGULAR VELOCITY"
  112. 2020 IF ps=2 then LET title$="ANGLE"
  113. 2030 CALL settext("TIME SERIES","TIME",title$)
  114. 2040 CALL setxscale(tmin,tmax)
  115. 2050 FOR k=0 to n-1
  116. 2060 
  117. 2070     IF ps = 1 then
  118. 2080        LET xreal(k)=thetadotdata(k+1)
  119. 2090     ELSE IF ps = 2 then
  120. 2100        LET xreal(k)=thetadata(k+1)
  121. 2110     END IF
  122. 2120     IF hanning$="y" then LET xreal(k) = xreal(k)*(.5-.5*cos(2*pi*k/(n-1)))
  123. 2130     LET ximag(k)=0
  124. 2140     LET tpoint(k)=tmin+k*del
  125. 2150 NEXT k
  126. 2160 CALL setgraphtype("")
  127. 2170 CALL datagraph(tpoint,xreal,1,0,"white")
  128. 2180 GET KEY keyvariable
  129. 2190 FOR i= 1 to 100
  130. 2200 NEXT i
  131. 2210 
  132. 2220 !FFT ALGORITHM
  133. 2230 CLEAR
  134. 2240 PRINT "Calculating FFT"
  135. 2250 LET n2=n/2
  136. 2260 LET nu1=nu-1
  137. 2270 LET k=0
  138. 2280 FOR l=1 to nu
  139. 2290     DO while k<(n-1)
  140. 2300        FOR i=1 to n2
  141. 2310            LET argument=k/2^nu1
  142. 2320            LET garbage=int(argument)
  143. 2330            LET p=bitr(garbage,nu)
  144. 2340            LET arg =2*pi*p/n
  145. 2350            LET c=cos(arg)
  146. 2360            LET s=sin(arg)
  147. 2370            LET k1=k+1
  148. 2380            LET k1n2=k1+n2
  149. 2390            LET treal=xreal(k1n2)*c+ximag(k1n2)*s
  150. 2400            LET timag=ximag(k1n2)*c-xreal(k1n2)*s
  151. 2410            LET xreal(k1n2)=xreal(k1)-treal
  152. 2420            LET ximag(k1n2)=ximag(k1)-timag
  153. 2430            LET xreal(k1)=xreal(k1)+treal
  154. 2440            LET ximag(k1)=ximag(k1)+timag
  155. 2450            LET k=k+1
  156. 2460        NEXT i
  157. 2470        LET k=k+n2
  158. 2480     LOOP
  159. 2490     LET k=0
  160. 2500     LET nu1=nu1-1
  161. 2510     LET n2=int(n2/2)
  162. 2520 NEXT l
  163. 2530 
  164. 2540 FOR k=1 to n
  165. 2550     LET i=bitr(k-1,nu)+1
  166. 2560     IF i<=k then GOTO 2630
  167. 2570     LET treal=xreal(k)
  168. 2580     LET timag=ximag(k)
  169. 2590     LET xreal(k)=xreal(i)
  170. 2600     LET ximag(k)=ximag(i)
  171. 2610     LET xreal(i)=treal
  172. 2620     LET ximag(i)=timag
  173. 2630 NEXT k
  174. 2640 
  175. 2650 !GRAPHING THE FFT
  176. 2660 CLEAR
  177. 2670 INPUT prompt"Plot 1)power spectrum, or 2)log power spectrum: ":pps
  178. 2680 INPUT prompt"Frequency variable - 1)linear, or 2)log: ":freqvar
  179. 2690 LET maxfreq=.5/del
  180. 2700 LET minfreq=1/(number*del)
  181. 2710 
  182. 2720 CLEAR
  183. 2730 !Y-AXIS
  184. 2740 IF pps = 1 then
  185. 2750    LET TITLE$="POWER SPECTRUM"
  186. 2760    LET YAXIS$="POWER"
  187. 2770 ELSE
  188. 2780    LET TITLE$="LOG POWER SPECTRUM"
  189. 2790    LET YAXIS$="LOG POWER"
  190. 2800 END IF
  191. 2810 !X-AXIS
  192. 2820 IF freqvar=2 then
  193. 2830    LET XAXIS$="LOG FREQUENCY"
  194. 2840 ELSE
  195. 2850    LET XAXIS$="FREQUENCY"
  196. 2860 END IF
  197. 2870 
  198. 2880 !DRAW AXES
  199. 2890 CLEAR
  200. 2900 CALL setxscale(minfreq,maxfreq)
  201. 2910 CALL setyscale(1e-6,.99)
  202. 2920 CALL SETTEXT(TITLE$,XAXIS$,YAXIS$)
  203. 2930 CALL RESERVELEGEND
  204. 2940 
  205. 2950 !PLOT POINTS
  206. 2960 FOR i=1 to n/2
  207. 2970     LET frequency(i)=i/(n*del)
  208. 2980     LET power(i)=(((xreal(i))^2+(ximag(i))^2))/(n^2)
  209. 2990 
  210. 3000 NEXT i
  211. 3010 !PLOT TEXT
  212. 3020 CALL setaxes(0)
  213. 3030 IF pps=1 then
  214. 3040    IF freqvar=1 then CALL setgraphtype("xy")
  215. 3050    IF freqvar=2 then CALL setgraphtype("logx")
  216. 3060 END IF
  217. 3070 IF pps=2 then
  218. 3080    IF freqvar=1 then CALL setgraphtype("logy")
  219. 3090    IF freqvar=2 then CALL setgraphtype("logxy")
  220. 3100 END IF
  221. 3102 IF NUMBER=4096 THEN
  222. 3110    CALL datagraph(frequency,power,1,1,"white")
  223. 3112 ELSE
  224. 3114    CALL DATAGRAPH(FREQUENCY,POWER,1,0,"WHITE")
  225. 3116 END IF
  226. 3120 CALL ADDLEGEND("N="&STR$(N)&"  MAX FREQ ="&STR$(MAXFREQ)&"  DEL F="&STR$(MINFREQ),0,1,"WHITE")
  227. 3130 CALL ADDLEGEND("  G="&STR$(G)&"  Q="&STR$(Q),0,1,"WHITE")
  228. 3140 CALL drawlegend
  229. 3150 IF hanning$="y" then
  230. 3160    !   CALL ADDLEGEND("  HANNING",0,1,"WHITE")
  231. 3170 END IF
  232. 3180 GET KEY keyvariable
  233. 3190 INPUT PROMPT "Another with Hanning? y/n: ":hann$
  234. 3200 IF hann$= "y" THEN GOTO 1840
  235. 3210 INPUT PROMPT "FFT of different quantity (y/n)? :": diffquant$
  236. 3220 IF diffquant$ = "y" THEN
  237. 3230    IF ps = 2 THEN LET ps = 1
  238. 3240    IF ps = 1 THEN LET ps = 2
  239. 3250    GOTO 2000
  240. 3260 END IF
  241. 3270 INPUT PROMPT "Different presentation of same FFT? (y/n): ":diffplot$
  242. 3280 IF diffplot$ ="y" THEN GOTO 2660
  243. 3290 END
  244. 3300 !
  245. 3310 SUB rk4(x,v,tstep,xnew,vnew,t,w,g,q)
  246. 3320     DECLARE DEF accel
  247. 3330     LET xk1=tstep*v
  248. 3340     LET vk1=tstep*accel(x,v,t,w,g,q)
  249. 3350     LET xk2=tstep*(v+vk1/2)
  250. 3360     LET vk2=tstep*accel(x+xk1/2,v+vk1/2,t+tstep/2,w,g,q)
  251. 3370     LET xk3=tstep*(v+vk2/2)
  252. 3380     LET vk3=tstep*accel(x+xk2/2,v+vk2/2,t+tstep/2,w,g,q)
  253. 3390     LET xk4=tstep*(v+vk3)
  254. 3400     LET vk4=tstep*accel(x+xk3,v+vk3,t+tstep,w,g,q)
  255. 3410     LET vnew=v+(vk1+2*vk2+2*vk3+vk4)/6
  256. 3420     LET xnew=x+(xk1+2*xk2+2*xk3+xk4)/6
  257. 3430 END SUB
  258. 3440 !
  259. 3450 !
  260. 3460 DEF accel(x,v,t,w,g,q)
  261. 3470     LET accel=-sin(x)-(1/q)*v+g*cos(w*t)
  262. 3480 END DEF
  263. 3490 
  264. 3500 !BIT REVERSER FUNCTION
  265. 3510 DEF bitr(j,nu)
  266. 3520     LET j1=j
  267. 3530     LET ibitr=0
  268. 3540     FOR i=1 to nu
  269. 3550         LET j2 = int(j1/2)
  270. 3560         LET ibitr=ibitr*2+(j1-2*j2)
  271. 3570         LET j1=j2
  272. 3580     NEXT i
  273. 3590     LET bitr=ibitr
  274. 3600 END DEF
  275.