Discrete Fourier Transform & Epicycles: Scripting

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)