Edit: 20220414
Anyone who played 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 programming concepts to FORTRAN students, but noodling around made BASIC ideal for teaching other concepts in science, technology, engineering, and math (STEM). 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 the author 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 almost everyone 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 8bit 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 for Python.
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 page will show how I converted a few BASIC examples found in Understanding the FFT (Anders E. Zonst, Citrus Press, Titusville, Florida) into Python3.
caveat: DFT = Discreet Fourier Transform. This forms the basis for FFT (Fast Fourier Transform)
BASIC allows floating point variables to be used in FORNEXT loops but Python does (not because floats are unpredictable AND their precision can vary from machine to machine). For example, loss of floating point precision can cause the FORNEXT loop in BASIC to terminate earlier (or later) than you would first expect.
Notes: the next example from the book mentions "Square Wave". To avoid confusion, the author wisely used the prompt
"TERMS" rather than "HARMONICS". Remember that "TERMS=1" always produces a sine wave (the first harmonic really means the
fundamental frequency) while a perfect square wave is only possible with "TERMS=infinity" (so 50 might be good enough here?).
Notice that at "line 40" we are only employing odd harmonics (stepping from 1 by 2 yields: 1, 3, 5...)
Changing this to step from 1 by 1 will produce all harmonics resulting in a descending triangle wave.
Original PCBASIC source code
1 REM ========================================== 2 REM DFT1_0.bas (fig 1.3 on Page 8) 3 REM "Understanding the FFT" by Anders E. Zonst 4 REM (c) Citrus Press. Titusville, Florida. 5 REM PCBASIC / GWBASIC / QBASIC / QuickBASIC 6 REM ========================================== 10 PRINT "*** DFT1.0  GENERATE SQUARE WAVE ***" 12 INPUT "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
TERMS? 1
0.0
.3826835
.7071068
.9238795
1.0
.9238795
.7071068
.3826835
1.509958E07
 .3826832
 .7071066
 .9238795
1.0
 .9238794
 .7071065
 .382683
TERMS? 50
0
.7661327
.7851163
.7929845
.7953941
.7929847
.7851165
.7661327
3.742222E06
.7661327
.7851161
.7929845
.7953941
.7929843
.785116
.7661321
LINES  Program Description (BASIC) 

1 to 6  Programmer comments begin with REM (for remark) 
10  In BASIC, data trailing the PRINT statement is passed to the PRINT statement 
20  Here the author defines PI (notice the trailing octothorpe which forces double precision floating point) 
30  BASIC allows FORNEXT loops to be based upon integers or floats (floating point numbers). 
3060  In BASIC, a FORNEXT loop will take you up to and including the limiting value 
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 
40  1) 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. 2) TERM=1 always produces a sine wave (think: Transverse Flute) 3) Stepping from one by two only employs odd harmonics (square wave; sounds like most reed instruments) 4) Stepping from one by one employs all harmonics (triangle wave; sounds like most stringed instruments) 
Test Run  1) In the test run (second panel) TERMS=1 only produces a fundamental wave with no harmonics. Caveat: many references
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 
caveat: unlike BASIC "fornext" loops, Python "forloops" only support integers.
RunAnywhere Python
(this is an incomplete stub)
#!/bin/python3 # ============================================= # Title : DFT0_0.py (generate sine wave) # Author : Neil Rieck # editor: vim8.2 (CentOS7) # plugin: pymode # Created: 20190823 # notes : # 1) Derived from: DFT1_0.bas (page 8, fig 1.3) # Understanding the FFT (2nd ed) # Anders E. Zonst (c) 19952000 Citrus Press # 2) we require 16 data points (015) # 3) do not include 360 degrees which is really # zero degrees of the next wave # ============================================= import math # PI + sin() # # the angle method (ugh!) # print("looping via the angle method") for A in range(0, 3600, 225): # 16 items: 03375 a = A/10 # reduce range by 10 r = math.radians(a) # angle to radian v = math.sin(r) # trigonometric sine # print("angle:", a, " sin:", v) # print("angle: {:>7.2f} sin: {: 2.16f}".format(a, v)) print(f"angle: {a:>7.2f} sin: {v: 2.16f}") # # 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) # trigonometric sine # print("angle:", a, " sin:", v) # print("time: {:3} sin: {: 2.16f}".format(t, v)) print(f"time: {t:3} sin: {v: 2.16f}") 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) Comments begin with an octothorpe (also called sharp; sometimes called hash) 2) The first line shown is called a shebang (programmer slang for sharp followed by bang) and is only relevant on UNIX and Linux platforms (this system has two python languages installed so a shebang is required to start python3) 
With Python3, 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 DFT1_0.bas (fig 1.3 on Page 8) 3 REM "Understanding the FFT" by Anders E. Zonst 4 REM (c) Citrus Press. Titusville, Florida. 5 REM PCBASIC, 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 # editor: vim8.2 (CentOS7) # plugin: pymode # created: 20190822 # Notes : # 1) Derived from: DFT1_0.bas (page 8, fig 1.3) # Understanding the FFT (2nd ed) # Anders E. Zonst (c) 19952000 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 want 16 data points (015) # 2) we do not want the value associated with 2PI #
PI2 = 2 * math.pi # for t in range(0, 17, 1): # yields: 015 R = t / 16 * PI2 # 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("point:",t," value:",Y) # formatted output (space after : leaves room for a sign) # print("point: {a:6.2f}, value: {b: 2.9f}".format(a=t, b=Y))
print(f"point: {t:6.2f}, value: {Y: 2.9f}") 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) array support can come 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 np #
#
N = 16 # desired set size
A = np.zeros(N) # create array and init to 0.0
for I in range (0, N):
print("a(",I,") = ", A[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
#!/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) array support can come 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 N
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.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 pip3 #  # # 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")
#!/bin/python3
'''
====================================================================
Title : print.py Author : Neil Rieck created: 20211016 notes :
1) PRINT comes in many forms
2) also note the use of three ticks to stop/start a block of remarks
==================================================================== '''
a9 = 123
b9 = 456
c9 = 789
d9 = "hello world"
# older positional formatted print example print("optional {0} stuff {1} goes {2} here {3}".format(a9, b9, c9, d9)) # older named formatted print example print("optional {w} stuff {x} goes {y} here {z}".format(w=a9, x=b9, y=c9, z=d9)) # newer literal formatted print example (be sure to prefix with an 'f') print(f"optional {a9} stuff {b9} goes {c9} here {d9}") #