home *** CD-ROM | disk | FTP | other *** search
/ Power-Programmierung / CD1.mdf / basic / library / trubasic / chaos / xexpfft.tru < prev    next >
Encoding:
Text File  |  1990-03-15  |  5.8 KB  |  197 lines
  1. 1000 !PROGRAM XEXPFFT.TRU
  2. 1010 LIBRARY "SGLIB.TRC"
  3. 1011 CLEAR
  4. 1012 PRINT"                ***FFT OF SUPERPOSED SINE WAVES***"
  5. 1013 PRINT
  6. 1014 PRINT"THIS PROGRAM TAKES THE FOURIER TRANSFORM OF A GROUP OF SINE"
  7. 1015 PRINT"WAVES WHOSE AMPLITUDES AND FREQUENCIES ARE INPUTS.  BOTH THE"
  8. 1016 PRINT"TIME SERIES AND THE TRANSFORM ARE GRAPHED.  IF A BLINKING CURSOR"
  9. 1017 PRINT" APPEARS PRESS ANY KEY TO CONTINUE.  THE HANNING OPTION SMOOTHS "
  10. 1018 PRINT" THE ABRUPT EFFECT OF THE WINDOW "
  11. 1019 PRINT" AND SUPPRESSES SPURIOUS COMPONENTS."
  12. 1030 DIM thetadata(5000),thetadotdata(5000),xreal(0 to 5000),ximag(0 to 10000)
  13. 1040 DIM tpoint(0 to 5000),power(2048),frequency(2048),FREQ(10),AMPL(10)
  14. 1050 DECLARE DEF bitr
  15. 1060 INPUT prompt"Max frequency : ":maxfreq
  16. 1070 INPUT Prompt"Input min.time:":tmin
  17. 1080 INPUT prompt"No. of FFT points(..256,512,1024,2048,4096) : ":number
  18. 1090 LET ps=1
  19. 1100 LET del=.5/maxfreq
  20. 1110 LET tmax=number*del+tmin
  21. 1120 LET n=number
  22. 1130 
  23. 1140 LET count=0
  24. 1150 LET p=1
  25. 1160 !DEVELOP PERIODIC TIME SERIES DATA
  26. 1170 INPUT PROMPT "HOW MANY FREQUENCY COMPONENTS?":NUMBFREQ
  27. 1180 FOR NF = 1 TO NUMBFREQ
  28. 1190     INPUT PROMPT" STATE FREQUENCY :":FREQ(NF)
  29. 1200     INPUT PROMPT" COMPONENT AMPLITUDE:":AMPL(NF)
  30. 1210 NEXT NF
  31. 1220 FOR P = 1 TO N
  32. 1230     LET TOTAL = 0
  33. 1240     FOR NF = 1 TO NUMBFREQ
  34. 1250         LET TOTAL = TOTAL + AMPL(NF)*SIN(2*PI*FREQ(NF)*(TMIN+P*DEL))
  35. 1260     NEXT NF
  36. 1270     LET THETADOTDATA(P)=TOTAL
  37. 1280 NEXT P
  38. 1290 !
  39. 1300 !
  40. 1310 !PREPARATION OF THE FFT DATA
  41. 1320 CLEAR
  42. 1330 INPUT prompt" HANNING OPTION Y/N? ": hanning$
  43. 1340 LET tgamma=log2(n)
  44. 1350 IF abs(int(tgamma)-tgamma)=0 then
  45. 1360    LET gamma=tgamma
  46. 1370    GOTO 1400
  47. 1380 END IF
  48. 1390 LET gamma=int(tgamma)+1
  49. 1400 PRINT "gamma= ";gamma
  50. 1410 LET newn=2^gamma
  51. 1420 LET nu=gamma
  52. 1430 FOR i=n+1 to newn
  53. 1440     LET xreal(i)=0
  54. 1450 NEXT i
  55. 1460 LET n=newn
  56. 1470 PRINT"n=";n
  57. 1480 CLEAR
  58. 1490 IF ps=1 then LET title$="WAVE DISPLACEMENT"
  59. 1500 CALL settext("TIME SERIES","TIME",title$)
  60. 1510 CALL setxscale(tmin,tmax)
  61. 1520 FOR k=0 to n-1
  62. 1530 
  63. 1540     IF ps = 1 then
  64. 1550        LET xreal(k)=thetadotdata(k+1)
  65. 1560     ELSE IF ps = 2 then
  66. 1570        LET xreal(k)=thetadata(k+1)
  67. 1580     END IF
  68. 1590     IF hanning$="y" then LET xreal(k) = xreal(k)*(.5-.5*cos(2*pi*k/(n-1)))
  69. 1600     LET ximag(k)=0
  70. 1610     LET tpoint(k)=tmin+k*del
  71. 1620 NEXT k
  72. 1630 CALL SETAXES(0)
  73. 1640 CALL setgraphtype("")
  74. 1650 CALL datagraph(tpoint,xreal,1,0,"white")
  75. 1660 GET KEY keyvariable
  76. 1670 FOR i= 1 to 100
  77. 1680 NEXT i
  78. 1690 
  79. 1700 !FFT ALGORITHM
  80. 1710 CLEAR
  81. 1720 PRINT "Calculating FFT"
  82. 1730 LET n2=n/2
  83. 1740 LET nu1=nu-1
  84. 1750 LET k=0
  85. 1760 FOR l=1 to nu
  86. 1770     DO while k<(n-1)
  87. 1780        FOR i=1 to n2
  88. 1790            LET argument=k/2^nu1
  89. 1800            LET garbage=int(argument)
  90. 1810            LET p=bitr(garbage,nu)
  91. 1820            LET arg =2*pi*p/n
  92. 1830            LET c=cos(arg)
  93. 1840            LET s=sin(arg)
  94. 1850            LET k1=k+1
  95. 1860            LET k1n2=k1+n2
  96. 1870            LET treal=xreal(k1n2)*c+ximag(k1n2)*s
  97. 1880            LET timag=ximag(k1n2)*c-xreal(k1n2)*s
  98. 1890            LET xreal(k1n2)=xreal(k1)-treal
  99. 1900            LET ximag(k1n2)=ximag(k1)-timag
  100. 1910            LET xreal(k1)=xreal(k1)+treal
  101. 1920            LET ximag(k1)=ximag(k1)+timag
  102. 1930            LET k=k+1
  103. 1940        NEXT i
  104. 1950        LET k=k+n2
  105. 1960     LOOP
  106. 1970     LET k=0
  107. 1980     LET nu1=nu1-1
  108. 1990     LET n2=int(n2/2)
  109. 2000 NEXT l
  110. 2010 
  111. 2020 FOR k=1 to n
  112. 2030     LET i=bitr(k-1,nu)+1
  113. 2040     IF i<=k then GOTO 2110
  114. 2050     LET treal=xreal(k)
  115. 2060     LET timag=ximag(k)
  116. 2070     LET xreal(k)=xreal(i)
  117. 2080     LET ximag(k)=ximag(i)
  118. 2090     LET xreal(i)=treal
  119. 2100     LET ximag(i)=timag
  120. 2110 NEXT k
  121. 2120 
  122. 2130 !GRAPHING THE FFT
  123. 2140 CLEAR
  124. 2150 INPUT prompt"Plot 1)power spectrum, or 2)log power spectrum: ":pps
  125. 2160 INPUT prompt"Frequency variable - 1)linear, or 2)log: ":freqvar
  126. 2170 LET maxfreq=.5/del
  127. 2180 LET minfreq=1/(number*del)
  128. 2190 
  129. 2200 CLEAR
  130. 2210 !Y-AXIS
  131. 2220 IF pps = 1 then
  132. 2230    LET TITLE$="POWER SPECTRUM"
  133. 2240    LET YAXIS$="POWER"
  134. 2250 ELSE
  135. 2260    LET TITLE$="LOG POWER SPECTRUM"
  136. 2270    LET YAXIS$="LOG POWER"
  137. 2280 END IF
  138. 2290 !X-AXIS
  139. 2300 IF freqvar=2 then
  140. 2310    LET XAXIS$="LOG FREQUENCY"
  141. 2320 ELSE
  142. 2330    LET XAXIS$="FREQUENCY"
  143. 2340 END IF
  144. 2350 
  145. 2360 !DRAW AXES
  146. 2370 CLEAR
  147. 2380 CALL setxscale(minfreq,maxfreq)
  148. 2390 CALL setyscale(1e-6,.99)
  149. 2400 CALL SETTEXT(TITLE$,XAXIS$,YAXIS$)
  150. 2410 CALL RESERVELEGEND
  151. 2420 
  152. 2430 !PLOT POINTS
  153. 2440 FOR i=1 to n/2
  154. 2450     LET frequency(i)=i/(n*del)
  155. 2460     LET power(i)=((((xreal(i))^2+(ximag(i))^2))^.5)/n
  156. 2480 
  157. 2490 NEXT i
  158. 2500 !PLOT TEXT
  159. 2510 CALL setaxes(0)
  160. 2520 IF pps=1 then
  161. 2530    IF freqvar=1 then CALL setgraphtype("xy")
  162. 2540    IF freqvar=2 then CALL setgraphtype("logx")
  163. 2550 END IF
  164. 2560 IF pps=2 then
  165. 2570    IF freqvar=1 then CALL setgraphtype("logy")
  166. 2580    IF freqvar=2 then CALL setgraphtype("logxy")
  167. 2590 END IF
  168. 2595 IF NUMBER =4096 THEN
  169. 2596   LET SYMBOL=1
  170. 2597 ELSE 
  171. 2598 LET SYMBOL=0
  172. 2599 END IF
  173. 2600 CALL datagraph(frequency,power,1,SYMBOL,"white")
  174. 2610 CALL ADDLEGEND("N="&STR$(N)&"  MAX FREQ="&STR$(MAXFREQ)&"  DEL F="&STR$(MINFREQ),0,1,"WHITE")
  175. 2620 CALL drawlegend
  176. 2630 IF hanning$="y" then
  177. 2640    CALL ADDLEGEND("  HANNING",0,1,"WHITE")
  178. 2650 END IF
  179. 2660 GET KEY keyvariable
  180. 2670 INPUT PROMPT "Another with Hanning? y/n: ":hann$
  181. 2680 IF hann$= "y" THEN GOTO 1320
  182. 2690 INPUT PROMPT "Different presentation of same FFT? (y/n): ":diffplot$
  183. 2700 IF diffplot$ ="y" THEN GOTO 2140
  184. 2710 END
  185. 2720 !
  186. 2730 !BIT REVERSER FUNCTION
  187. 2740 DEF bitr(j,nu)
  188. 2750     LET j1=j
  189. 2760     LET ibitr=0
  190. 2770     FOR i=1 to nu
  191. 2780         LET j2 = int(j1/2)
  192. 2790         LET ibitr=ibitr*2+(j1-2*j2)
  193. 2800         LET j1=j2
  194. 2810     NEXT i
  195. 2820     LET bitr=ibitr
  196. 2830 END DEF
  197.