Anyone who started playing with BASIC on personal computers in the 1970s and 1980s will recognize the importance of that language for noodling around. The creators of BASIC intended it to be used to teach computer programming concepts to FORTRAN students, but noodling around made BASIC ideal for teaching other concepts in science, engineering, and math. I am a huge fan of DFT/FFT books (especially these two: Understanding the FFT and Understanding FFT Applications by Anders Zonst of Citrus Press, Titusville, Florida) where authors provided hundreds of demo programs written in PCBASIC (a generic term I am using for this article). I am assuming they did this because BASIC programs are readable by nonprogrammers AND computers with BASIC were so ubiquitous that their readers could play with the example programs. In fact, I am going to assume that anyone with the curiosity to learn DFT/FFT will already own a personal computer with at least one programming language.
BASIC on personal computers of the 1970s was usually implemented in ROM, and every implementation was different (Apple2, TRS80, HeathKitH9). Starting with the IBMPC in 1981, Microsoft, began publishing 16bit software products like GWBASIC (1983), QuickBASIC (1985) and QBASIC (1991) which worked well on 16bit operating systems like MSDOS (1981) up through Windows3.11 (1992). These 16bit language interpreters were also supported on 32bit operating systems starting with Windows95 through to Windows7 via a process known as THUNKING.
The big problem today is that 64bit computers run 64bit operating systems, like Windows10, where 32bit programs are THUNKED but not 16bit programs. Technical workarounds exist including "setting up a virtual machine" on your 64bit OS but why go to all that bother when all you want to do is noodle around? Perhaps it is time to ditch BASIC
There are numerous free Python tutorials on the net along with a plethora of examples so I won't waste your time duplicating them here. But this webpage will show how I converted a few BASIC examples found in Understanding the FFT (Anders Zonst of Citrus Press, Titusville, Florida) into Python3.
BASIC allows floating point variables to be used in FORNEXT loops but Python does not because floats are unpredictable. For example, loss of floating point precision causes the FORNEXT loop in BASIC to terminate earlier than you would first expect (I believe the author knew about this and depended upon it; but he could be forgiven since his intention was to teach a few DFT concepts)
Original PCBASIC source code
1 REM ============================================================= 2 REM title : DFT1_0.bas (fig 1.3 on Page 8) 3 REM book : "Understanding the FFT" by Anders E. Zonst 4 REM : (c) Citrus Press. Titusville, Florida. 5 REM language: GWBASIC/QBASIC/QuickBASIC 6 REM ============================================================= 10 PRINT "*** DFT1.0  GENERATE SQUARE WAVE ***" 12 INPUT "NUMBER OF TERMS";N 20 PI = 3.14159265358# 30 FOR I = 0 TO 2*PI STEP PI/8 32 Y=0 40 FOR J=1 TO N STEP 2: Y=Y+SIN(J*I)/J: NEXT 50 PRINT Y 60 NEXT I 70 END
Output (terms=1 means no harmonics)
*** DFT1.0  GENERATE SQUARE WAVE ***
NUMBER OF TERMS: 1
0.0
.3826835
.7071068
.9238795
1.0
.9238795
.7071068
.3826835
1.509958E07
 .3826832
 .7071066
 .9238795
1.0
 .9238794
 .7071065
 .382683
Program Description (BASIC)  

LINES 1 to 6  Programmer comments begin with REM (for remark) 
LINE 10  In BASIC, data trailing the PRINT statement is passed to the PRINT statement 
LINE 20  Here the author defines PI (notice the trailing octothorpe which forces double precision floating point) 
LINE 30  BASIC allows FORNEXT loops to be based upon integers or floats (floating point numbers). 
LINES 30 and 40  In BASIC, a FORNEXT loop will take you up to and including the limiting value 
LINE 30  1) Line 30 of the BASIC program implements a FORNEXT loop starting at
zero then stepping by "PI/8" radians (0.392699081698724) up to "2*PI"
radians (6.28318530717959). 2) A nontechnical person might think that this loop will get all the way from start (0) to finish (6.28) producing 17 data points. The author's text indicates he was only shooting for 16 data points and that is what he got because of the precision limitations of floating point. (but I think he was aware of that) 3) No one attending "DATA PROCESSING 101" this side of Y2K would be allowed to hand in an assignment where a FORNEXT loop was based upon a float because some newer floating point implementations, like XFLOAT (IEEE double float) will actually produce 17 data points as I show here 
LINE 40  implements an inner FORNEXT loop all on one line. Here, the
individual code segments are separated with a colon. This will actually
run a tiny bit faster on a 16bit BASIC interpreter. Also note that this loop starts at one and steps by two meaning that this algorithm only generates odd harmonics. 
Test Run  1) In the test run (second panel) I input TERMS=1 to output only the
fundamental wave with no harmonics. Caveat: some texts describe the
fundamental as the first harmonic (???) 2) "1.509958E07" is exponential notation for a very small number (E07 indicates to move the decimal point seven places to the left) 3) Something like "PRINT USING "#.#####" would have produced a more uniform output 
RunAnywhere Python
(this is
an incomplete stub)
#!/bin/python3 # Title : DFT0_0.py # Author : Neil Rieck # created: 20190823 # notes : 1) we want 16 data points (015) # : 2) we do not want the point associated with 360 # degrees since that is the start of the next wave # (this would mess up trigbased calculus) #  import math # PI lives here # # the angle method (ugh!) # print("looping via the angle method") for a_10 in range (0, 3600, 225): a = a_10/10 # reduce range by 10 r = math.radians(a) # v = math.sin(r) # get the trigonometric sine # print("angle:", a, " sin:", v) print("angle: {:>7.2f} sin: {: 2.16f}".format(a,v)) # # the radian method (better) # print("looping via radian method (preferred)") for t in range (0, 16, 1): # 16 items: 015 r = (t/16) * 2 * math.pi # convert to radians v = math.sin(r) # get the trigonometric sine # print("angle:", a, " sin:", v) print("time: {:3} sin: {: 2.16f}".format(t,v)) print("Exit")
Output From
Two Methods
angle method angle: 0.00 sin: 0.0000000000000000 angle: 22.50 sin: 0.3826834323650898 angle: 45.00 sin: 0.7071067811865476 angle: 67.50 sin: 0.9238795325112867 angle: 90.00 sin: 1.0000000000000000 angle: 112.50 sin: 0.9238795325112867 angle: 135.00 sin: 0.7071067811865476 angle: 157.50 sin: 0.3826834323650899 angle: 180.00 sin: 0.0000000000000001 angle: 202.50 sin: 0.3826834323650897 angle: 225.00 sin: 0.7071067811865475 angle: 247.50 sin: 0.9238795325112868 angle: 270.00 sin: 1.0000000000000000 angle: 292.50 sin: 0.9238795325112866 angle: 315.00 sin: 0.7071067811865477 angle: 337.50 sin: 0.3826834323650896 radian method (preferred) time: 0 sin: 0.0000000000000000 time: 1 sin: 0.3826834323650898 time: 2 sin: 0.7071067811865476 time: 3 sin: 0.9238795325112867 time: 4 sin: 1.0000000000000000 time: 5 sin: 0.9238795325112867 time: 6 sin: 0.7071067811865476 time: 7 sin: 0.3826834323650899 time: 8 sin: 0.0000000000000001 time: 9 sin: 0.3826834323650897 time: 10 sin: 0.7071067811865475 time: 11 sin: 0.9238795325112865 time: 12 sin: 1.0000000000000000 time: 13 sin: 0.9238795325112866 time: 14 sin: 0.7071067811865477 time: 15 sin: 0.3826834323650904 Exit
Program Description (Python) 

1) Programmer comments begin with an octothorpe (also called sharp;
sometimes called hash) 2) The first line shown is called a shebang and is only relevant on UNIX and Linux platforms 
With Python3, every statement including PRINT is implemented as a function so you must include parentheses (soft brackets) when passing data to the function (same a C/C++) 
Python only supports integers in FOR loops (this restriction is intentional due to the fact that a lack of floating point precision will always make floating loops unpredictable) 
1) In Python, a FOR loop will only take you up to a value LESS than the
limiting value. 2) BTW, Python uses indentation to define loops. The control statement must end with a colon. The statements to be executed must be indented by 4 spaces There are no NEXT statements. The loop ends when indentation decreases. 
1) Almost all science, engineering and math employs RADIANs rather
than DEGREEs 2) I have provided two examples (I prefer method #2) 
1) In the test run (yellow panel) you will notice that both routines
produce similar results to the BASIC test run shown above. 2) Also notice that lines "time: 11" and "time: 15" produce slightly different results than "angle: 247.5" and "angle: 337.5". I suspect (but have not verified) that this is due to a slight loss of precision due to the actions of library function "math.radians()" but this is expected in programs because floats have never been more than an approximation. As an aside, financial applications requiring penny accuracy employ decimal or BCD data types. 3) remember that if we were to multiply this result then we would also multiply that error (so it is best to multiply the angle before running it through that function) 
Original PCBASIC source code
1 REM ============================================================= 2 REM title : DFT1_0.bas (fig 1.3 on Page 8) 3 REM book : "Understanding the FFT" by Anders E. Zonst 4 REM : (c) Citrus Press. Titusville, Florida. 5 REM language: GWBASIC/QBASIC/QuickBASIC 6 REM ============================================================= 10 PRINT "*** DFT1.0  GENERATE SQUARE WAVE ***" 12 INPUT "NUMBER OF TERMS";N 20 PI = 3.14159265358# 30 FOR I = 0 TO 2*PI STEP PI/8 32 Y=0 40 FOR J=1 TO N STEP 2: Y=Y+SIN(J*I)/J: NEXT 50 PRINT Y 60 NEXT I 70 END
RunAnywhere Python
#!/bin/python3 # Title : DFT1_0.py # Author : Neil Rieck # created: 20190822 # Notes : Converted from DFT1_0.bas (c) Citrus Press #  import math print("*** DFT1.0  GENERATE SQUARE WAVE ***") print("NUMBER OF HARMONICS: ",end='') junk = input() if (junk==''): junk=1 N = int(junk) # # loop from 0 to 2PI stepping by PI/8 "with integers" # notes: # 1) we but we do want 16 data points (015) # 2) we do not want the value associated with 2PI # for t in range(0, 16, 1): # yields: 015 R = t/16*2* math.pi # Y = 0 # init accum each pass # # usage: step=2 (odd harmonics) for square wave # step=1 (all harmonics) for a ramp wave # for J in range (1, N+1, 2): Y = Y + math.sin(R*J)/J # # raw output # print("angle:",I," value:",Y) # # formatted output (space after b: leaves room for a sign) print("step: {a:6.2f}, value: {b: 2.9f}".format(a=t,b=Y)) print("Exit")
PCBASIC: output
*** DFT1.0  GENERATE SQUARE WAVE ***
NUMBER OF TERMS: 1000
0.0
.7867051
.7846908
.785939
.7848983
.7859398
.7846898
.7867033
7.510006E05
.7867056
.7846903
.7859396
.7848983
.7859381
.7846909
.7867049
Python: output
*** DFT1.0  GENERATE SQUARE WAVE ***
NUMBER OF HARMONICS: 1000
step: 0.00, value: 0.000000000
step: 1.00, value: 0.786704710
step: 2.00, value: 0.784691059
step: 3.00, value: 0.785939359
step: 4.00, value: 0.784898164
step: 5.00, value: 0.785939359
step: 6.00, value: 0.784691059
step: 7.00, value: 0.786704710
step: 8.00, value: 0.000000000
step: 9.00, value: 0.786704710
step: 10.00, value: 0.784691059
step: 11.00, value: 0.785939359
step: 12.00, value: 0.784898164
step: 13.00, value: 0.785939359
step: 14.00, value: 0.784691059
step: 15.00, value: 0.786704710
Exit
PCBASIC  Python2  Python3 




#!/bin/python3 # Title : DFT0_1.py # Author : Neil Rieck # created: 20190823 # notes : # 1) we want to store 16 data points (015) in an array # 2) arrays are native to BASIC and C but not Python # 3) one method of array support comes via the "numpy" library # 4) on the Windows version of Python you must download then # install numpy from the Windows CMD tool like so: # a) rightclick START button # b) click RUN item # c) type: CMD # d) type: pip install numby #  import math # import numpy as py # # N = 16 # desired set size A = py.zeros(N) # create array and init to 0.0 for I in range (0, N): print("a(",I,") = ", A[I]) print("Exit") 
#!/bin/python3 # Title : DFT0_2.py # Author : Neil Rieck # created: 20190829 # notes : # 1) we want to store 16 data points (015) in an array # 2) arrays are native to BASIC and C but not Python # 3) one method of array support comes via the "numpy" library # 4) this method uses the "array" library supplied with Python3 #  import array # N = 16 # desired set size A = array.array('d', range(N)) # create array; init to sub scr B = array.array('d',( 0 for x in range(N))) # create array; init to 0.0 for I in range (0, N): # # raw output # print("a(",I,") = ", A[I],"b(",I,") = ", B[I]) # # formatted output print("a({a:2}) = {b: 6.2f} b({a:2}) = {c: 6.2f}".format(a=I,b=A[I],c=B[I])) # print("Exit") 

a( 0 ) = 0.0 a( 1 ) = 0.0 a( 2 ) = 0.0 a( 3 ) = 0.0 a( 4 ) = 0.0 a( 5 ) = 0.0 a( 6 ) = 0.0 a( 7 ) = 0.0 a( 8 ) = 0.0 a( 9 ) = 0.0 a( 10 ) = 0.0 a( 11 ) = 0.0 a( 12 ) = 0.0 a( 13 ) = 0.0 a( 14 ) = 0.0 a( 15 ) = 0.0 Exit 
a( 0) = 0.00 b( 0) = 0.00 a( 1) = 1.00 b( 1) = 0.00 a( 2) = 2.00 b( 2) = 0.00 a( 3) = 3.00 b( 3) = 0.00 a( 4) = 4.00 b( 4) = 0.00 a( 5) = 5.00 b( 5) = 0.00 a( 6) = 6.00 b( 6) = 0.00 a( 7) = 7.00 b( 7) = 0.00 a( 8) = 8.00 b( 8) = 0.00 a( 9) = 9.00 b( 9) = 0.00 a(10) = 10.00 b(10) = 0.00 a(11) = 11.00 b(11) = 0.00 a(12) = 12.00 b(12) = 0.00 a(13) = 13.00 b(13) = 0.00 a(14) = 14.00 b(14) = 0.00 a(15) = 15.00 b(15) = 0.00 Exit 
Original PCBASIC source code
6 REM ****************************************** 8 REM *** (DFT3.1) GENERATE/ANALYZE WAVEFORM *** 10 REM ****************************************** 12 PI=3.141592653589793#:P2=2*PI:K1=PI/8:K2=1/PI 14 DIM Y(16),FC(16),FS(16),KC(16),KS(16) 16 CLS:FOR J=0 TO 16:FC(J)=0:FS(J)=0:NEXT 20 GOSUB 108: REM  PRINT COLUMN HEADINGS 30 GOSUB 120: REM  GENERATE FUNCTION Y(X) 40 GOSUB 200: REM  PERFORM DFT 60 GOSUB 140: REM  PRINT OUT FINAL VALUES 70 PRINT:PRINT "MORE (Y/N)? "; 72 A$ = INKEY$:IF A$="" THEN 72 74 PRINT A$:IF A$ = "Y" THEN 16 80 END 100 REM ****************************************** 102 REM * PROGRAM SUBROUTINES * 104 REM ****************************************** 106 REM * PRINT COLUMN HEADINGS * 107 REM ****************************************** 108 PRINT:PRINT 110 PRINT "FREQ F(COS) F(SIN) Y(COS) Y(SIN)" 112 PRINT 114 RETURN 118 REM ****************************** 120 REM *** GENERATE FUNCTION F(X) *** 121 REM ****************************** 122 FOR I=0 TO 15:K3=I*K1 124 Y(I)=COS(K3)+COS(3*K3)/9+COS(5*K3)/25+COS(7*K3)/49 126 NEXT 128 FOR I = 1 TO 7 STEP 2: KC(I)=1/I^2:NEXT 130 RETURN 132 REM ****************************** 138 REM * PRINT OUTPUT * 139 REM ****************************** 140 FOR Z=0 TO 15 142 PRINT Z;" "; 144 PRINT USING "##.#####_ ";FC(Z),FS(Z),KC(Z),KS(Z) 146 NEXT Z 148 RETURN 200 REM ************************** 202 REM * SOLVE FOR COMPONENTS * 204 REM ************************** 206 FOR J=0 TO 16:REM SOLVE EQNS FOR EACH FREQUENCY 208 FOR I = 0 TO 15:REM MULTIPLY AND SUM EACH DATA POINT 210 FC(J)=FC(J)+Y(I)*COS(J*I*K1):FS(J)=FS(J)+Y(I)*SIN(J*I*K1) 212 NEXT I 214 FC(J)=FC(J)/16: FS(J)=FS(J)/16:REM FIND MEAN VALUE 216 NEXT J 218 RETURN
Python Output
DFT3.1 Generate/Analyze Waveform FREQ F(COS) F(SIN) Y(COS) Y(SIN) 0 0.00000 0.00000 0.00000 0.00000 1 0.50000 0.00000 1.00000 0.00000 2 0.00000 0.00000 0.00000 0.00000 3 0.05556 0.00000 0.11111 0.00000 4 0.00000 0.00000 0.00000 0.00000 5 0.02000 0.00000 0.04000 0.00000 6 0.00000 0.00000 0.00000 0.00000 7 0.01020 0.00000 0.02041 0.00000 8 0.00000 0.00000 0.00000 0.00000 9 0.01020 0.00000 0.00000 0.00000 10 0.00000 0.00000 0.00000 0.00000 11 0.02000 0.00000 0.00000 0.00000 12 0.00000 0.00000 0.00000 0.00000 13 0.05556 0.00000 0.00000 0.00000 14 0.00000 0.00000 0.00000 0.00000 15 0.50000 0.00000 0.00000 0.00000 Exit
RunAnywhere Python
#!/bin/python3 # Title : DFT3_1.py # Author : Neil Rieck # created: 20190831 # notes : # 1) Converted from DFT3_1.bas (c) Citrus Press # 2) this program employs array.array which is builtinto Python3 # 3) the numpy.array alternative must be installed via pip #  # # display column headings # def display_column_headings(): print("\nFREQ F(COS) F(SIN) Y(COS) Y(SIN)") print # # generate data # def generate_data(): for I in range(N): # step 015 by 1 K3=I*K1 Y[I]=(math.cos(K3) + math.cos(3*K3)/9.0 + math.cos(5*K3)/25.0 + math.cos(7*K3)/49) for I in range(1,7+1,2): # step 17 by 2 KC[I]=1/(I**2) # # display data # def display_data(): for Z in range(N): print("{a:2} {b: 7.5f} {c: 7.5f} {d: 7.5f} {e: 7.5f}" .format(a=Z,b=FC[Z],c=FS[Z],d=KC[Z],e=KS[Z])) # # solve for components # def solve_via_dft(): for J in range(N): # solve for each frequency for I in range(N): # solve for each data point FC[J]=FC[J]+Y[I]*math.cos(J*I*K1) FS[J]=FS[J]+Y[I]*math.sin(J*I*K1) FC[J]=FC[J]/16 # find mean value FS[J]=FS[J]/16 # find mean value # ======================================= # main program # ======================================= print("DFT3.1 Generate/Analyze Waveform") import math import array # P2=2*math.pi K1=math.pi/8.0 K2=1/math.pi N = 16 # desired set size # # the following lines were coded for clarity, not speed # Y = array.array('d', range(N)) FC = array.array('d', range(N)) FS = array.array('d', range(N)) KC = array.array('d', range(N)) KS = array.array('d', range(N)) for I in range(N): Y[I] = 0.0 FC[I] = 0.0 FS[I] = 0.0 KC[I] = 0.0 KS[I] = 0.0 # # let's get going # display_column_headings() generate_data() solve_via_dft() display_data() print("Exit")