An adaptation from this lovely paper:[br][br][url=https://aapt.scitation.org/doi/10.1119/1.3254017]https://aapt.scitation.org/doi/10.1119/1.3254017[/url]

[size=150]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.[br][br]It also plots the system of epicycles to trace out a closed loop defined by the data set.[br][br]Finally, it calculates the trigonometric interpolation of the data set by means of the IDFT. [br][br][b]Note:[/b] 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. [br][br]I would like to thank [url=https://www.geogebra.org/u/thijs]Thijs[/url], who help me improve the script. [/size]

[size=150]#==============================[br]# Define discrete data signal and Setup[br]#==============================[br]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)}[br][br]N = Length(Sx)[br]LN = Sequence(N)[br][br]#===========================================[br]# Shift zero-frequency component to center of spectrum.[br]#===========================================[br][br]Lk = LN - 1 - Div(N, 2) [br][br]#================================[br]# Calculate DFT: LX is a list of frequencies [br]#================================[br][br]LX = 1 / N * Zip(Sum(Zip((abs(P); arg(P) - 2 * pi * k * (n-1) / N), P, Sx, n, LN)), k, Lk)[br][br]#===============================[br]# Lj : frequency index sorted by size (abs)[br]#===============================[br][br]LAs = Reverse(Sort(Zip((abs(X), n), X, LX, n, LN)))[br]Lj = Zip(y(As), As, LAs)[br][br]#=============================[br]# Use the first M greatest frequencies[br]#=============================[br][br]M = Slider(1, N, 1, 1, 160, false, true, false, false)[br]SetValue(M, N)[br]Lm = Sequence(1, M)[br][br]Lks = First(Zip(Element(Lk, j), j, Lj), M)[br]LXs = First(Zip(Element(LX, j), j, Lj), M)[br]LRs = First(Zip(x(As), As, LAs), M)[br][br]#==========================[br]# Plot M epicycles[br]#==========================[br][br]speed = 0.5[br]t = Slider(0, 2 * pi, 0.0099, speed, 150, false, true, false, false)[br][br]LC1= Zip(Sum(Zip((abs(X); arg(X) + k * t), X, First(LXs, m), k, Lks)), m, Lm )[br]LC = Join({(0, 0)}, LC1)[br][br]Plast = Last(LC)[br]PolyL = Polyline(LC)[br]Epicycles = Zip(Circle(Element(LC, m), R), m, Lm, R, LRs) [br][br]#===========================================[br]# Parametric curve: Inverse of Discrete Fourier Transform[br]#===========================================[br][br]fx(x) = Sum(Zip( x(X) * cos(k * x) - y(X) * sin(k * x), X, LXs, k, Lks))[br]fy(x) = Sum(Zip( x(X) * sin(k * x) + y(X) * cos(k * x), X, LXs, k, Lks))[br]orbit = Curve(fx(t), fy(t), t, 0, 2 * pi)[br][br]#==========================[br]# Extras: Settings[br]#==========================[br][br]SetValue(t, 0)[br]StartAnimation(t, true)[br][br]SetCaption(M, "n")[br]SetLabelMode(M, 9)[br]SetVisibleInView(M, 1, true)[br]SetVisibleInView(M, 2, false)[br][br]SetVisibleInView(t, 1, true)[br]SetVisibleInView(t, 2, false)[br][br]SetVisibleInView(LX , 1, false)[br]SetVisibleInView(LAs, 1, false)[br]SetVisibleInView(LXs, 1, false)[br]SetVisibleInView(LC1, 1, false)[br]SetVisibleInView(LC , 1, false)[br]SetVisibleInView(fx, 1, false)[br]SetVisibleInView(fy, 1, false)[br][br]ShowLabel(PolyL, false)[br]ShowLabel(orbit, false)[br][/size]