# Python Notes: DFT + FFT

Edit: 2022-05-27

## A few BASIC-to-Python conversion examples

"BASIC hacking noodling" was moved here

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 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 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.

There are numerous free Python tutorials and examples on the web 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 into Python3.

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

```#!/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
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}")
#
#
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
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)

#### 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

```#!/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
# ---------------------------------------------------------------
#
#
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
#
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 forms2) also note the use of three ticks to stop/start a block of remarks====================================================================
'''a9 = 123b9 = 456c9 = 789d9 = "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}")
#``` Back to Home