10 ' FFT PROGRAM FROM BYTE MAGAZINE, DECEMBER '78
20 PRINT"		*********************************"
30 PRINT"		*				*"
40 PRINT"		*            F F T 		*
50 PRINT"		*				*"
60 PRINT"		*********************************"
70 PRINT:PRINT"A program to derive and display a 256 point"
80 PRINT"frequency-domain representation of a time-domain"
90 PRINT"function.   The required time-domain function can"
100 PRINT"be either generated within the program by patching"
110 PRINT"in an appropriate routine in lines 200 to 990, or else"
120 PRINT"input and stored from an A/D converter."
130 PRINT"(A sine generator is initially included)."
140 PRINT
150 PRINT"	**** PLEASE PRESS RETURN TO CONTINUE ****"
160 INPUT R
170 DIM X1(256),X2(256)
180 N=256: L=8: PI=3.14159
200 PRINT"	---- GENERATING A TIME FUNCTION ----"
210 PRINT"		(1010Hz sinewave)"
220 T=0
230 FOR Z=0 TO N-1
240 X1(Z)=SIN(2*PI*1000*T)
250 T=T+1.95313E-04
260 NEXT Z
1000 PRINT "Do you want a listing of the generated time function";
1010 INPUT A$
1020 IF A$="NO" THEN 1120
1030 IF A$<>"YES" THEN 1000
1040 B=X1(0)
1050 FOR Z=0 TO N-1
1060 IF ABS(X1(Z))>B THEN B=ABS(X1(Z))
1070 NEXT Z
1080 FOR Z=0 TO N-1
1090 PRINT X1(Z);TAB(41+20*X1(Z)/B);"*"
1100 NEXT Z
1110 '	---- SCALE INPUT FUNCTION ----
1120 FOR Z=0 TO N-1
1130 X1(Z)=X1(Z)/N
1140 NEXT Z
1150 '	---- FFT IN PLACE ALGORITHM ----
1160 PRINT:PRINT"	---- FFT CALCULATION IN PROGRESS ----"
1170 PRINT"	(Please give me a few minutes)""
1180 I1=N/2: I2=1: V=2*PI/N
1190 FOR I=1 TO L
1200 I3=0: I4=I1
1210 FOR K=1 TO I2
1220 X=INT(I3/I1)
1230 GOSUB 1840
1240 I5=Y
1250 Z1=COS(V*I5)
1260 Z2=-SIN(V*I5)
1270 FOR M=I3 TO I4-1
1280 A1=X1(M): A2=X2(M)
1290 B1=Z1*X1(M+I1)-Z2*X2(M+I1)
1300 B2=Z2*X1(M+I1)+Z1*X2(M+I1)
1310 X1(M)=A1+B1: X2(M)=A2+B2
1320 X1(M+I1)=A1-B1: X2(M+I1)=A2-B2
1330 NEXT M
1340 I3=I3+2*I1: I4=I4+2*I1
1350 NEXT K
1360 I1=I1/2: I2=2*I2
1370 NEXT I
1380 '---- OUTPUT RESULTS ----
1390 PRINT"In what form do you want the output?"
1400 PRINT"	Magnitude spectrum plot		(1)"
1410 PRINT"	Table of values			(2)"
1420 INPUT A
1430 IF A=1 THEN 1470
1440 IF A=2 THEN 1660
1450 PRINT"INCORRECT INPUT, 1 OR 2 PLEASE!": GOTO 1390
1460 '	---- OUTPUT MAGNITUDE SPECTRUM PLOT ----
1470 B=0
1480 PRINT:PRINT"	---- MORE CALCULATIONS IN PROGRESS ----"
1490 FOR Z=0 TO N/2
1500 X=Z
1510 GOSUB 1930
1520 IF X3>B THEN B=X3
1530 NEXT Z
1540 FOR Z=0 TO N/2
1550 X=Z
1560 GOSUB 1930
1570 X4=INT(56*X3/B)
1580 C=0
1590 PRINT Z;TAB(5);"!";
1600 C=C+1
1610 IF C<X4 THEN PRINT"=";: GOTO 1600
1620 PRINT ""
1630 NEXT Z
1640 GOTO 1780
1650 '	---- OUTPUT TABLE OF VALUES ----
1660 U=0
1670 Z=0
1680 PRINT"HARMONIC";TAB(14);"REAL";TAB(30);
1690 PRINT"IMAGINARY";TAB(50);"MAGNITUDE"
1700 X=U
1710 GOSUB 1930
1720 PRINT U;TAB(10);X1(Y);TAB(30);X2(Y);TAB(50);X3
1730 U=U+1: Z=Z+1
1740 IF Z>9 THEN 1670
1750 IF U>N/2 THEN 1780
1760 GOTO 1700
1770 '	---- TERMINATE? ----
1780 PRINT"DO you want another output (YES, NO)";
1790 INPUT A$
1800 IF A$="YES"THEN 1390
1810 IF A$<>"NO" THEN 1780
1820 END
1830 '	---- SCRAMBLER ROUTINE ----
1840 Y=0: N1=N
1850 FOR W=1 TO L
1860 N1=N1/2
1870 IF X<N1 THEN 1900
1880 Y=Y+2^(W-1)
1890 X=X-N1
1900 NEXT W
1910 RETURN
1920 '	---- MAGNITUDE (X3) SUBROUTINE ----
1930 GOSUB 1840
1940 X3=SQR(X1(Y)^2+X2(Y)^2)
1950 RETURN
1960 END