Python Notes: DFT + FFT

  1. The information presented here is intended for educational use.
  2. The information presented here is provided free of charge, as-is, with no warranty of any kind.
Edit: 2019-09-30

A few BASIC-to-Python conversion examples

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 PC-BASIC (a generic term I am using for this article). I am assuming they did this because BASIC programs are readable by non-programmers -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, TRS-80, HeathKit-H9). Starting with the IBM-PC in 1981, Microsoft, began publishing 16-bit software products like GW-BASIC (1983),  QuickBASIC (1985) and QBASIC (1991) which worked well on 16-bit operating systems like MS-DOS (1981) up through Windows-3.11 (1992). These 16-bit language interpreters were also supported on 32-bit operating systems starting with Windows-95 through to Windows-7 via a process known as THUNKING.

The big problem today is that 64-bit computers run 64-bit operating systems, like Windows-10, where 32-bit programs are THUNKED but not 16-bit programs. Technical work-arounds exist including "setting up a virtual machine" on your 64-bit 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.

DFT-1 (Discrete Fourier Transform - Wave Generation )

First Road Block - changes to looping

BASIC allows floating point variables to be used in FOR-NEXT loops but Python does not because floats are unpredictable. For example, loss of floating point precision causes the FOR-NEXT 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 PC-BASIC 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: GW-BASIC/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.509958E-07 
- .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 FOR-NEXT loops to be based upon integers or floats (floating point numbers).
LINES 30 and 40 In BASIC, a FOR-NEXT loop will take you up to and including the limiting value
LINE 30 1) Line 30 of the BASIC program implements a FOR-NEXT loop starting at zero then stepping by "PI/8" radians (0.392699081698724) up to "2*PI" radians (6.28318530717959).
2) A non-technical 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 FOR-NEXT 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 FOR-NEXT 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 16-bit 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.509958E-07" is exponential notation for a very small number (E-07 indicates to move the decimal point seven places to the left)
3) Something like "PRINT USING "-#.#####" would have produced a more uniform output

Two alternative ideas for cyclical looping in Python

Run-Anywhere Python
(this is an incomplete stub)

#!/bin/python3
# Title  : DFT0_0.py
# Author : Neil Rieck
# created: 2019-08-23
# notes  : 1) we want 16 data points (0-15)
#        : 2) we do not want the point associated with 360
#             degrees since that is the start of the next wave
#             (this would mess up trig-based 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: 0-15
    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)

Demo #1: DFT-1.0 (BASIC vs. Python3)

Original PC-BASIC 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: GW-BASIC/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


BASIC Notes:
1) Examples DFT1_0 to DFT1_4  (from the book) are just variations of this program.
2) STEP 2 on Line 40 only produces odd harmonics. Change this to "STEP1" for
    all harmonics
3) "/J" on Line 40 reduces the amplitude for each harmonic
4) Other interesting modifications to Line 40 are described in the author's text

Python Notes:
1) this program appears longer because of my embedded remarks
2) I wasted a few lines at the to to ensure a BLANK input is translated to "1"
    (in BASIC this automatically becomes zero)
3) my "raw output" is easy to code but I've commented that out to show you
   one way to do formatted output

Run-Anywhere Python

#!/bin/python3
# Title  : DFT1_0.py
# Author : Neil Rieck
# created: 2019-08-22
# 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 (0-15)
#   2) we do not want the value associated with 2PI
#
for t in range(0, 16, 1):       # yields: 0-15
    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")

Side-by-side comparison

PC-BASIC: output

 *** DFT1.0 - GENERATE SQUARE WAVE ***
NUMBER OF TERMS: 1000
0.0 
 .7867051 
 .7846908 
 .785939 
 .7848983 
 .7859398 
 .7846898 
 .7867033 
7.510006E-05 
-.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
Comment: 7.510006E-05 is exponential notation where E-05 means move the decimal point 5 places to the left (result: 0.00007510006)

DFT-3 (Discrete Fourier Transform - Generate/Analyze)

Second Road Block - array declaration

PC-BASIC Python2 Python3
  • arrays support is built-in
  • creating arrays is easy
  • arrays are auto initialized
  • array support is not built-in
  • array support is provided by a third-party library called "NumPy"
  • informally known as "numpy.array" on the net
  • "numpy" must be downloaded and installed outside of Python via PIP
  • array initialization can be performed during array creation
  • "numby" is very very popular
  • array support is built-in
  • informally known as "array.array" on the net
  • array initialization left up to the programmer (appears non-intuitive)
  • "numpy" is available here as well (still installed via pip)
  • data processing on very large arrays will be more than twice as fast with "numpy.array" vs "array.array"
  • many scientific demos you see on the net will employ "numpy" with Python3
#!/bin/python3
# Title  : DFT0_1.py
# Author : Neil Rieck
# created: 2019-08-23
# notes  :
# 1) we want to store 16 data points (0-15) 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) right-click 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: 2019-08-29
# notes  :
# 1) we want to store 16 data points (0-15) 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

Demo #2: DFT-3.1

Original PC-BASIC 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

Run-Anywhere Python

#!/bin/python3
# Title  : DFT3_1.py
# Author : Neil Rieck
# created: 2019-08-31
# notes  :
# 1) Converted from DFT3_1.bas (c) Citrus Press
# 2) this program employs array.array which is built-into 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 0-15 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 1-7 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")

External Links


 Back to Home
 Neil Rieck
 Waterloo, Ontario, Canada.