Python Notes: DFT + FFT

Edit: 2023-10-29

Python introduction was moved here

A few BASIC-to-Python conversion examples

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 example programs written in PC-BASIC (a generic term I am using for this article). I am assuming he did this because BASIC programs are readable by non-programmers -AND- computers with BASIC were ubiquitous in the last two decades of the twentieth century meaning almost everyone could play with the demos. In fact, I am going to assume that anyone with the curiosity to learn DFT + FFT already owned a personal computer with at least one programming language.

This page will show how I converted a few examples found in Understanding the FFT from BASIC into Python (python3 to be exact).

caveats:

DFT-1 (Discrete Fourier Transform - Wave Generation )

caveat: DFT = Discreet Fourier Transform. This forms the basis for FFT (Fast Fourier Transform)

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 AND their precision can vary from machine to machine). For example, loss of floating point precision can cause the FOR-NEXT loop in BASIC to terminate earlier, or later, than you might expect.

Notes:

  1. The next example from the book mentions "Square Wave"
  2. To avoid confusion, the author wisely used the prompt "TERMS" rather than "HARMONICS".
  3. 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?).
  4. Notice that at "line 40" we are only employing odd harmonics (stepping from 1 by 2 yields: 1, 3, 5...)
  5. Changing this to step from 1 by 1 will produce all harmonics resulting in a descending triangle wave.

Original PC-BASIC 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 PC-BASIC / GW-BASIC / 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.509958E-07 
- .3826832 
- .7071066 
- .9238795 
-1.0 
- .9238794 
- .7071065 
- .382683 
TERMS? 50
 0 
 .7661327 
 .7851163 
 .7929845 
 .7953941 
 .7929847 
 .7851165 
 .7661327 
 3.742222E-06 
-.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 FOR-NEXT loops to be based upon integers or floats (floating point numbers).
30-60 In BASIC, a FOR-NEXT loop will take you up to and including the limiting value
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
40 1) 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.
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.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 Python ideas for looping with "for"

caveat: unlike BASIC "for-next" loops, Python "for-loops" only support integers.

Run-Anywhere Python
(this is an incomplete stub)

#!/usr/bin/python3
# =============================================
# Title  : DFT0_0.py (generate sine wave)
# Author : Neil Rieck
#          editor: vim8.2 (CentOS-7)
#          plugin: pymode
# Created: 2019-08-23
# notes  :
# 1) Derived from: DFT1_0.bas (page 8, fig 1.3)
#    Understanding the FFT (2nd ed)
#    Anders E. Zonst (c) 1995-2000 Citrus Press
# 2) we require 16 data points (0-15)
# 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: 0-3375
    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: 0-15
    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)
3) Block comments begin-end with three single quotes (obviously cannot be used with the shebang)
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)
 

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

Original PC-BASIC 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 PC-BASIC, 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 "STEP 1" 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 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

#!/usr/bin/python3
# =============================================
# Title  : DFT1_0.py
# Author : Neil Rieck
#          editor: vim8.2 (CentOS-7)
#          plugin: pymode
# created: 2019-08-22
# Notes  :
# 1) Derived from: DFT1_0.bas (page 8, fig 1.3)
#    Understanding the FFT (2nd ed)
#    Anders E. Zonst (c) 1995-2000 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 (0-15)
#   2) we do not want the value associated with 2PI
#
PI2 = 2 * math.pi # for t in range(0, 17, 1): # yields: 0-15 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")
 

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 which is close to zero)
 

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

Second Road Block: array declaration

PC-BASIC Python2 Python3
  • array support is built-in
  • creating arrays is easy
  • arrays are auto initialized
  • array support is not built into python2
  • 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
  • "numpy" is very very popular and can be used with python3
  • array support is available via the non-thirty party library array
  • informally known as array.array on the net
  • array initialization left up to the programmer (appears non-intuitive)
  • numpy is also available here (but installed via pip3)
  • data processing on very large arrays will be more than twice as fast with numpy vs 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) 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) right-click 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: 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) 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
 

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 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 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")
 

Variation in PRINT statements etc.

#!/bin/python3
'''
====================================================================
Title : print.py Author : Neil Rieck created: 2021-10-16 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}") #

External Links


 Back to Home
 Neil Rieck
 Waterloo, Ontario, Canada.