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:
- only convert this BASIC stuff if you want to learn low level conceptual details of DFT and FFT (or want to learn why Python is better than every BASIC).
- In production use, a Python programmer would import either numpy or scipy, then would have access to professionally written FFT routines.
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:
- 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 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
DFT-3 (Discrete Fourier Transform - Generate/Analyze)
Second Road Block: array declaration
| PC-BASIC | Python2 | 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
- https://www.python.org/
- https://www.python.org/downloads/ -- get your python interpreter for any platform along with the IDLE integrated development environment
- https://docs.python.org/3/tutorial/
- https://pythonclock.org/ - Python2.7 will not receive any security updates during/after 2020 so make sure you are using Python3
- https://realpython.com/python-scipy-fft/
- https://pythonnumericalmethods.berkeley.edu/notebooks/chapter24.04-FFT-in-Python.html
- OpenVMS Notes: floats (subtitled: fun with floats)
- Why Radians rather than Degrees?
- https://www.zdnet.com/article/python-is-eating-the-world-how-one-developers-side-project-became-the-hottest-programming-language-on-the-planet/
- https://www.zdnet.com/article/microsoft-we-want-you-to-learn-python-programming-language-for-free/
Back to HomeNeil Rieck
Waterloo, Ontario, Canada.