Discrete Fourier Transform & Epicycles: Scripting
- Author:
- Juan Carlos Ponce Campuzano
What is this script about?
The following script calculates the DFT from a data set of size N. In this example, the points belong to a square of radius 4.
It also plots the system of epicycles to trace out a closed loop defined by the data set.
Finally, it calculates the trigonometric interpolation of the data set by means of the IDFT.
Note: You can define a different data set. Just make sure that your data set forms a closed loop in a cycle and that each point is equally spaced (approximately, at least). If your data set contains more than 300 elements, then I recommend you to use a different mathematical software or programming language.
I would like to thank Thijs, who help me improve the script.
GGB script
#==============================
# Define discrete data signal and Setup
#==============================
Sx = {(2, 2), (2, 3/2), (2, 1), (2, 1/2), (2, 0), (2, -1 / 2), (2, -1), (2,-3/2), (2, -2), (3/2, -2), (1, -2), (1 / 2, -2), (0, -2), (-1/2, -2), (-1,-2), (-3/2, -2), (-2, -2), (-2, -3/2), (-2, -1), (-2, -1 / 2), (-2, 0), (-2,1 / 2), (-2, 1), (-2, 3 / 2), (-2, 2), (-3 / 2, 2), (-1, 2), (-1 / 2, 2), (0,2), (1 / 2, 2), (1, 2), (3 / 2, 2)}
N = Length(Sx)
LN = Sequence(N)
#===========================================
# Shift zero-frequency component to center of spectrum.
#===========================================
Lk = LN - 1 - Div(N, 2)
#================================
# Calculate DFT: LX is a list of frequencies
#================================
LX = 1 / N * Zip(Sum(Zip((abs(P); arg(P) - 2 * pi * k * (n-1) / N), P, Sx, n, LN)), k, Lk)
#===============================
# Lj : frequency index sorted by size (abs)
#===============================
LAs = Reverse(Sort(Zip((abs(X), n), X, LX, n, LN)))
Lj = Zip(y(As), As, LAs)
#=============================
# Use the first M greatest frequencies
#=============================
M = Slider(1, N, 1, 1, 160, false, true, false, false)
SetValue(M, N)
Lm = Sequence(1, M)
Lks = First(Zip(Element(Lk, j), j, Lj), M)
LXs = First(Zip(Element(LX, j), j, Lj), M)
LRs = First(Zip(x(As), As, LAs), M)
#==========================
# Plot M epicycles
#==========================
speed = 0.5
t = Slider(0, 2 * pi, 0.0099, speed, 150, false, true, false, false)
LC1= Zip(Sum(Zip((abs(X); arg(X) + k * t), X, First(LXs, m), k, Lks)), m, Lm )
LC = Join({(0, 0)}, LC1)
Plast = Last(LC)
PolyL = Polyline(LC)
Epicycles = Zip(Circle(Element(LC, m), R), m, Lm, R, LRs)
#===========================================
# Parametric curve: Inverse of Discrete Fourier Transform
#===========================================
fx(x) = Sum(Zip( x(X) * cos(k * x) - y(X) * sin(k * x), X, LXs, k, Lks))
fy(x) = Sum(Zip( x(X) * sin(k * x) + y(X) * cos(k * x), X, LXs, k, Lks))
orbit = Curve(fx(t), fy(t), t, 0, 2 * pi)
#==========================
# Extras: Settings
#==========================
SetValue(t, 0)
StartAnimation(t, true)
SetCaption(M, "n")
SetLabelMode(M, 9)
SetVisibleInView(M, 1, true)
SetVisibleInView(M, 2, false)
SetVisibleInView(t, 1, true)
SetVisibleInView(t, 2, false)
SetVisibleInView(LX , 1, false)
SetVisibleInView(LAs, 1, false)
SetVisibleInView(LXs, 1, false)
SetVisibleInView(LC1, 1, false)
SetVisibleInView(LC , 1, false)
SetVisibleInView(fx, 1, false)
SetVisibleInView(fy, 1, false)
ShowLabel(PolyL, false)
ShowLabel(orbit, false)