8
$\begingroup$

There are a number of nice Wolfram Demonstrations on Fourier Series, but none quite what I seek. The closest I've found is this one, but it has several unnecessary features (e.g., colors) and missing features. Specifically, it addresses just one single target function (a square wave) and forces the user to scroll time forward and back.

What I seek is more closely demonstrated starting at 2:15 in this video. In brief, time keeps running and the graph scrolls to the right. There are multiple circles, each rotating at angular frequency $\omega$, $2 \omega$, ..., and appropriate phases, all added together to approximate the target function.

Specifically I seek a Manipulate with the following features:

  • A pulldown menu with perhaps 10 leading periodic functions (sine, square, triangle, ramp left, ramp right, rectified sine, ...
  • A selector button list for the number $n$ (from $n = 1$ to $n = 15$) of Fourier components to be included in the approximation
  • a scrolling graph of the selected target function (sine, square, ...) and the trace of the Fourier series including just the first (selected) $n$ components
  • a two-dimensional graph at the left showing the trajectory of the tip of the chain of rotating vectors. (All this should be clear from the video clip.)

You will need to know the Fourier series coefficients for each target function in order to create each graph. Ideally, these coefficients will be computed using FourierSeries or other appropriate function.

I started programming this and got bogged down in the graphics and more. Here is a start:

    Manipulate[
        Row[{Graphics[{Darker[Green], Circle[{0, 0}, 2], 
     Black, Line[{{0, 0}, 2 {Cos[2 \[Pi] t], Sin[2 \[Pi] t]}}], 
    Red, Line[{2 {Cos[2 \[Pi] t], Sin[2 \[Pi] t]}, {2, 
        2 Sin[2 \[Pi] t]}}]}, ImageSize -> 300],
    Plot[{-Sin[2 \[Pi] x], -SquareWave[x]}, {x, -t, -t + 1},
    PlotStyle -> {Red, Blue},
    Exclusions -> None,
    PlotRange -> {-1.1, 1.1}, AxesOrigin -> {-t, 0}, 
    ImageSize -> 500]}],
     {t, 0, 20, AnimationRate -> .3, Appearance -> "Open"},
     {{n, 1, "components"}, 1, 15, 1, ControlType -> SetterBar},
     {{f, "sine", "function"}, {"sine", "square", "ramp left", 
   "ramp right", "pulse train"}, ControlType -> PopupMenu}]

I realize this is real work, so I'll post a bounty at the earliest opportunity (in a few hours).

$\endgroup$
4
  • 3
    $\begingroup$ Have you tried contacting the author asking for a code since he explicitly mentions in the video that he used Mathematica? Why "reinventing the wheel" if someone already coded everything :) $\endgroup$
    – Domen
    Jun 24 at 19:54
  • $\begingroup$ @Domen: I have indeed tried to contact the author/programmer, without success. $\endgroup$ Jun 25 at 3:28
  • 2
    $\begingroup$ FourierSeries gives complex coefficients. The demo is, in effect, the Fourier sine series on a odd function. Sines give a nice way to show how the height of the graph corresponds to the sum of the Fourier components. I think a real-value function can be written $\sum a_k \sin(k t + \phi_k)$, in which the components sum to the height of the graph and $a_k$, $\phi_k$ are computed from the $\pm k$ Fourier coefficients. Would that be what seek? $\endgroup$
    – Michael E2
    Jun 25 at 14:52
  • $\begingroup$ @Michael E2: All my target functions are real. And yes... perhaps the simplest way forward is to compute $a_k$ and $\phi_k$ from the $\pm k$ coefficients. Ideally, I'd like to enter an arbitrary periodic function (period = 1) as a potential input, e.g., SquareWave[x] + .5 TriangleWave[2 x]. $\endgroup$ Jun 25 at 21:07

2 Answers 2

11
+50
$\begingroup$
ClearAll[accumulatedRV,linkedWipers, plotGrid]

accumulatedRV[funcs_ : {Cos, Sin}][radii_List, velocities_List] := 
   Accumulate[radii Transpose @ Through[funcs[velocities #]]] &

linkedWipers[funcs_ : {Cos, Sin}][radii_List, 
   velocities_List, styles_ : Automatic] := 
 Module[{acc = Prepend[{0, 0}] @ accumulatedRV[funcs][radii, velocities] @ #, 
    st = PadRight[styles /. Automatic -> {Directive[Thin, Gray]}, 
      Length@radii, "Periodic"]}, 
  {Red, PointSize @ Large, 
    Point @ acc[[-1]], Point @ Scaled[{.999, 0}, acc[[-1]]], 
    Line[{acc[[-1]], Scaled[{1, 0}, acc[[-1]]]}], 
    MapThread[{#, Line @ #2, Circle[#3, Abs @ #4], 
       PointSize @ Medium, Point @ #3} &] @ 
     {st, Partition[acc, 2, 1], Most @ acc, radii}}] &

Manipulate[Graphics[linkedWipers[][Reverse@Range[7], Range[7], 
    ColorData[97] /@ Range[7]]@t, PlotRange -> 30],
 {t, 0, 1}]

enter image description here

fibonacci10 = N@Fibonacci[Range@10]; 

Manipulate[
 Graphics[{linkedWipers[][Reverse @ fibonacci10, fibonacci10, 
   ColorData[97] /@ Range[10]][t] /. {Red, a__, _Line, b__} :> {Red, a, b}, 
   Thick, Blue, 
   Line[Table[Total[Reverse[fibonacci10] Transpose@
      Through[{Cos, Sin}[fibonacci10  k]]], {k, 0, t, 10^-3}]]}, 
  ImageSize -> Large, PlotRange -> 150], {{t, 4}, 0, 2 Pi}]

enter image description here

We start with a list of wave functions:

wnames = {"squarePlusTriangle","squareWave", "triangleWave", "unitBox", 
   "sawtoothWaveDown","sawtoothWaveUp", 
   ToString[Style["\[FreakedSmiley]", 16], StandardForm] <> "sin" <> 
    ToString[Style["\[FreakedSmiley]", 16], StandardForm]};

wlist = {SquareWave[#] + TriangleWave[#]/2 &, SquareWave, 
   TriangleWave, UnitBox[#] /. 0 -> -1 &, 
   SawtoothWave[{-1, 1}, #] &, SawtoothWave[{-1, 1}, 1 - #] &, 
   Sin[Pi/3 + 2 Pi #] &};

We precompute FourierSinSeries for the functions in wlist (so that Manipulate runs smoothly):

fourierSinSeries = AssociationThread[wnames, 
   FourierSinSeries[# @ x, x, 21, FourierParameters -> {1, 2 Pi}] & /@ wlist];

componentLists = MapApply[List] @ fourierSinSeries;

radiusLists = MapApply[ReplaceAll[_Sin -> 1] @* List] @ fourierSinSeries;

velocityLists = Map[Cases[#, Sin[a_ x] :> a, All] &] @ fourierSinSeries;

waves = AssociationThread[wnames, wlist];

Define a function that generates the plots given a wave and the number of components:

plotGrid[wave_, n_][t_] := 
 Module[{sseries = Take[fourierSinSeries @ wave, n], 
   radiuslist = Take[radiusLists @ wave, n], 
   velocitylist = Take[velocityLists @ wave, n], , fun = waves @ wave},
  Grid[{{Show[Graphics[linkedWipers[][radiuslist, velocitylist, 
         ColorData[97] /@ Range[Length@radiuslist]] @ t], 
      ParametricPlot[Evaluate@{sseries /. {Sin -> Cos, Cos -> Sin}, sseries}, 
        {x, 10^-5, t}, 
      PerformanceGoal -> "Quality"],
      ImageSize -> 1 -> {80, 80},
      PlotRangeClipping -> False, ImageMargins -> 0,
      PlotRange -> {{-1.1, 1.1} Total[radiuslist], {-2, 2}}],
     ParametricPlot[{{t - x, fun[x]}, {t - x, sseries}}, {x, 0, t + 10^-5}, 
      Exclusions -> "Discontinuities",
      ExclusionsStyle -> Directive[Gray, Dashed],
      ImageSize -> 1 -> {400, 80},
      Epilog -> {Red, PointSize[Large], Point[{0, sseries /. x -> t}]},
       PlotRange -> {{0, 1.05}, {-2, 2}}, Axes -> True,
      ImageMargins -> 0, PlotRangePadding -> 0, 
      PlotRangeClipping -> False]}}, Spacings -> {0, 0}]]

... and Manipulate:

Manipulate[Labeled[plotGrid[w, n][t], Style[w, 16], {Top, Left}], 
 {{w, "squareWave", "wave"}, wnames, 
    PopupMenu[Dynamic[w, (w = #; nl = Length[componentLists @ #]; &)], 
    wnames] &}, 
 {{n, nl, "components"}, Range[1, nl], SetterBar}, 
 {{t, .5}, 0, 1}, 
 {{nl, Length[componentLists @ "squareWave"]}, None}]

enter image description here

enter image description here

enter image description here

The reason behind \[FreakedSmiley]s around sin:

enter image description here

edit by Nasser (added movies to make it come alive)

enter image description here

enter image description here

$\endgroup$
6
  • $\begingroup$ Awww gee... you're just too good. It doesn't seem fair. My deepest thanks. As promised, I'll post a bounty in a few hours (when permitted) and gratefully accept this superb answer. $\endgroup$ Jun 25 at 22:27
  • $\begingroup$ @DavidG.Stork, my pleasure. I am hoping someone will figure out (and post) the tricks behind the cool examples in Doga Kurkcuoglu's blog for a more general solution. $\endgroup$
    – kglr
    Jun 25 at 22:44
  • $\begingroup$ @Nasser, thank you so much for your edit. $\endgroup$
    – kglr
    Jun 25 at 22:45
  • $\begingroup$ kglr: When I try to add functions (with names), I often get errors. Try adding (SquareWave + TriangleWave)/2 (or (SquareWave[#] + TriangleWave[#])/2 &) and notice the mismatch between the rotating "wipers" at the left and their trajectory. How do I add new functions while avoiding such problems? $\endgroup$ Jun 25 at 23:28
  • $\begingroup$ David, updated with (SquareWave[#] + TriangleWave[#])/2 added to the list of waves. $\endgroup$
    – kglr
    Jun 25 at 23:50
7
$\begingroup$

Something for one of my classes. It may have "several unnecessary features." And be missing some. Can delete if not a propos.

Mathematica graphics

Some features:

  • Parses variables so students may enter functions of t or of x or of whatever.
  • Syntax errors don't clobber the demo. In fact, an edit window preserves the mis-written code for the user to emend. An alternative would be to allow the error-ridden code, but prevent the Fourier analysis and visualization.
  • Allows for a parameter, so one can see how the Fourier series changes as the parameter moves. (I use it to show phase shift in $f(t-a)$ for instance.)
  • An Animator[] allows the time to increase indefinitely while the phasor linkage is animated and the graph scrolls to the left.
  • An OpenerView[] toggles the animated Fourier phasor components.
  • One can manually enter a high order for a slowly convergent Fourier series.
  • Numerical Fourier[] (FFT) for speed and robustness (can handle functions whose Fourier coefficients cannot be symbolically determined), high base order for accuracy.

Code:

The DEBUG lines and comments could have been eliminated to make the code look shorter.

(*
 * Fourier[]
 *   with componentGrid
 *   and autovar
 *   and parameters
 *   and compiled function
 *)
debug = <||>;
DEBUG // ClearAll;
DEBUG // Attributes = {HoldAllComplete};
$debug = False;
DEBUG[x_] /; TrueQ@$debug := x;
$demoSize = 550;
Manipulate[
 With[{t = var, a = param},
  With[{
    line = ReIm[ (* Fourier phasor sum *)
      I*Accumulate[
         Prepend[ConstantArray[2, order], 1]*
          fc*
          Exp[2 Pi I Range[0, order] (-tt)]
         ] +
       tt + $circleoffset/2], pt = {tt, series /. t -> tt}
    ,
    plot = 
     Plot[ (* function and partial series, 
      with room for rotating phasor sum *)
      {functionC[Mod[t, 1], aa], series},
      {t, tt - 1.25, tt},
      PlotStyle -> {AbsoluteThickness[3], AbsoluteThickness[1.5]},
      Exclusions -> None,
      PlotRange -> {{tt - 1.25, 
           tt + Max[$circleoffset, 0.5/Pi]}, $yrange}, 
      PlotRangePadding -> {{0, Scaled[.05]}, {Scaled[0.1], 
         Scaled[0.1]}},
      ImagePadding -> 25, Frame -> True, AspectRatio -> Automatic,
      GridLines -> {Range[Ceiling[tt - 1.25], 
         Floor[tt + $circleoffset]], None},
      PlotLabel -> Row@{"var: ", var, ", param: ", param}
      ]},
   OpenerView[{ (* plot of function, series, phasor sum; 
     individual phasors k=1..10 *)
     Show[(* plot of function, series, phasor sum *)
      (* line, circle, points for phasor sum *)
      Graphics[ 
       {Lighter@Darker@Yellow,
        (* rolling circle has diam. 1/(2Pi) => unit circumference == 
        one period *)
        Circle[First[line], 1/(2 Pi)], 
        Point@CirclePoints[
          First[line], {1/(2 Pi), -2 Pi*(0 tt + line[[1, 1]])}, 8],
        (* circle rolls on this line *)
        InfiniteLine[First@line - {0, 1/(2 Pi)}, {1, 0}],
        Point[ (* ticks on line *)
         Thread[{
           Range[Floor[tt - 1.25, 0.125], 
            Ceiling[tt + $circleoffset, 0.125], 0.125],
           line[[1, 2]] - 1/(2 Pi)}]]
        }
       ]
      ,
      (* plot of function, series *)
      plot 
      ,
      (* phasor sum *)
      Graphics[ 
       {Line@line, {Thin, 
         MapThread[Circle, {Most@line, Norm /@ Differences@line}]},
        Red, {Opacity[0.5], Line[{Last@line, pt}]},
        PointSize@Medium, Point@pt, PointSize@Small, Point@Last@line}
       ]
      , 
      (* Show[] options *)
      ImageSize -> $demoSize,
      Options@plot
      ]
     ,
     (* individual phasors k=1..10 *)
     componentGrid /. t -> tt 
     },
    Dynamic@open, (* OpenerView[] state *)
    Alignment -> Center]
   ]]
 , 
 (* controls & options *)
 {{tt, 0, Dynamic@var}, 0, 4 (* independent time variable *)
  , Row[{ (* slider: *)
     Manipulator[Dynamic[Clip[tt, #2], (tt = #) &], ##2],
     "  ", (* animator for continuous running: *)
     Animator[#, {0, Infinity}, AnimationRunning -> False, 
      AnimationRate -> 1/5]}] &},
 {{aa, 0., Dynamic@param}, -1, 1,(* parameter value *)
  TrackingFunction -> (updateSeries[Null, order, #] &)},
 {{order, 3}, 1, 10, 1, (* truncation order of series *)
  TrackingFunction -> (updateSeries[Null, #, Null] &)},
 {{function, SquareWave[t - a]/4 + t^2},
  Row[{ (* input expression: *)
     InputField[Dynamic[#, updateSeries[#, order, Null] &]],
     " ",(* setter for independent variable: *)
     SetterBar[
      Dynamic[var, (var = #; {param} = DeleteCases[syms, var]; 
         updateSeries[function, order, aa]) &], syms],
     " edit:",
     InputField[
      Dynamic[functionText, updateSeries[#, order, Null] &]]}] &},
 {{error, None}, Pane@error &}, (* show parse error *)
 (* local variables and functions *)
 {{functionText, function}, 
  None}, (* last edit of function, even if syntax error *)
 {functionC, None}, (* compiled function *)
 {fc, None}, (* Fourier coefficient up to order *)
 {fcStore, None}, (* Fourier coefficients up to high order *)
 {series, None}, (* truncated Fourier series *)
 {$yrange, None}, (* y plot range *)
 {$circleoffset, 
  None}, (* offset from graph to give room for Fourier phasor sum *)
 {componentGrid, None}, (* grid of Fourier components/phasors *)
 {{open, False}, 
  None}, (* OpenerView state (True => component grid is showing *)
 {{var, t}, None}, (* function arg symbol/indep. variable *)
 {{param, a}, None}, (* parameter symbol *)
 {{syms, {t, a}}, None}, (* symbols in the input expression *)
 {{updateSeries, updateSeries}, None}, (* update function *)
 {maxFC, None}, (* max abs val Fourier coefficient *)
 {componentGridFN, 
  None}, (* function that generates the Fourier component grid *)
 (*" initializes update function, phasor grid function, 
     initial values for Manipulate-user variables "*)
 Initialization :> (
   updateSeries[newFunction_, newOrder_, newParamValue_] := 
    Module[{newSyms, fvals},
     DEBUG[debug["time"] = {AbsoluteTime[]}];
     (* Parse function *)
     error = None;
     DEBUG[debug["newFN"] == newFunction];
     If[newFunction =!= Null,
      functionText = newFunction;
      If[Head@newFunction =!= RawBoxes,
       (* parse variable/parameter *)
       newSyms = Variables@Level[newFunction, {-1}];
       If[Length@newSyms == 0, {param, var} = newSyms = {a, t}]; (* 
       defaults *)
       If[
        Length@newSyms == 1, {param, var} = 
         newSyms = Prepend[newSyms, a]];
       If[Length@newSyms == 2,
        If[! MemberQ[newSyms, var],
         {param, var} = newSyms,
         {param} = DeleteCases[newSyms, var] 
         ];
        If[NumericQ[newFunction /. {var -> tt, param -> aa}],
         syms = newSyms;
         function = newFunction;
         With[{fargs = {{var, param}, PiecewiseExpand@function}},
          functionC = 
           Compile @@ 
            Join[fargs, {RuntimeAttributes -> {Listable}, 
              Parallelization -> True}];
          If[! FreeQ[functionC, {46 | 47, __}], (* not compilable *)
           functionC = Function @@ Join[fargs, {Listable}]
           ]
          ],
         error = Pane[
           Row@{newFunction, " did not evaluate to a number at ", var,
              "=", tt, ", ", param, "=", aa},
           $demoSize]
         ],
        error = Pane[
          Row@{newFunction, 
            " should have only one parameter, one variable (", 
            Sequence @@ Riffle[newSyms, ", "], ")"},
          $demoSize]
        ],
       error = Pane[
         Row@{newFunction, " incomplete syntax or syntax error"},
         $demoSize]
       ]
      ];
     DEBUG[debug["time"] = {debug["time"], AbsoluteTime[]}];
     (* update param value *)
     If[newParamValue =!= Null,
      aa = newParamValue
      ];
     (* update order value *)
     If[newOrder =!= Null,
      order = newOrder
      ];
     (* If function parsed ok... *)
     If[error === None,
      If[newFunction =!= Null || newParamValue =!= Null,
       (* calculate Fourier coefficients via FFT *)
       (* currently only called when control is active *)
       fvals = 
        functionC[Subdivide[0., 1., ControlActive[2^12, 2^16]], aa];
       fcStore = fvals // Fourier // #/Sqrt[Length@# - 1] &;
       DEBUG[debug["fcStore"] = fcStore];
       maxFC = Max@Abs@fcStore[[;; 11]]; (* 
       for relative phasor size in componentGrid *)
       (* calculate y range for plotting *)
       (* overestimates Gibbs; 
       does not always contain Fourier component line *)
       $yrange = MinMax@fvals;
       DEBUG[debug["yrange1"] = $yrange];
       $yrange = $yrange - {-1, 1} 0.2 Subtract @@ $yrange;
       If[$yrange[[1]] > Re@First@fcStore - 0.6/Pi, (* 
        rolling circle has diam. 1/(2Pi) *)
        $yrange[[1]] = Re@First@fcStore - 0.6/Pi];
       If[$yrange[[2]] < Re@First@fcStore + 0.6/Pi,
        $yrange[[2]] = Re@First@fcStore + 0.6/Pi];
       DEBUG[debug["yrange2"] = $yrange];
       ];
      DEBUG[debug["time"] = {debug["time"], AbsoluteTime[]}];
      fc = Take[fcStore, order + 1]; (* Fourier coefficients *)
      DEBUG[debug["fc"] = fc];
      series = Total[ (* Fourier series *)
        Prepend[ConstantArray[2, order], 1]*
         Abs[fc]*
         Sin[2 Pi Range[0, order] (-var) + Arg[fc] + Pi/2]
        ];
      DEBUG[debug["series"] = series];
      DEBUG[debug["time"] = {debug["time"], AbsoluteTime[]}];
      $circleoffset = 4.4 Total@Abs@Rest@Take[fcStore, 11];
      componentGrid = componentGridFN @@ fcStore[[2 ;; 11]];
      DEBUG[debug["time"] = {debug["time"], AbsoluteTime[]}];
      DEBUG[debug["parsed"] = {var, param, syms, $yrange}]
      ];
     DEBUG[debug["time"] = AbsoluteTime[] - debug["time"]];
     ];
   componentGridFN = (* grid of rotating Fourier component phasors *)
    Block[{order, maxFC, var},
     (*" Block initial values.
         Circle radii relative to largest Fourier coefficient.
         Grid is precomputed to save time. "*)
     Function@
      Evaluate@
       With[{$fc = 
          I Array[Slot, 10]*Exp[2 Pi I Range[1, 10] (-var)]/maxFC},
        Labeled[
         GraphicsGrid[
          Partition[
           MapIndexed[
            Graphics[{ (* Fourier phasor rep. by circle/angle line *)
               If[First@#2 > order, Opacity[0.3], Nothing],
               Circle[], Blue,
               Circle[{0, 0}, Abs[#1]],
               Line[{{0, 0}, ReIm@#1}]},
              PlotRange -> {{-1.5, 1.5}, {-1.1, 1.1}}
              ] &,
            $fc (* scaled Fourier components *)
            ],
           5],
          ImageSize -> $demoSize
          ],
         Style["Fourier Components", "Label"],
         Top
         ]
        ]
     ];
   updateSeries[SquareWave[t - a]/4, 3, 0.]
   ),
 TrackedSymbols :> {tt, functionC, order, aa, var}
 ]

edit by Nasser (we must have some movies!)

enter image description here

enter image description here


2D example, kinda like in the video

But because it's simple silhouette, it doesn't quite have the loopiness of the video's 2D example.

img = Import["https://i.stack.imgur.com/G4eFq.jpg"];

SeedRandom[1];
(paths = FindCurvePath[
    pts = 
     Nest[DeleteCases[#, 
        Alternatives @@ 
         Replace[Rest@Nearest[#, RandomChoice[#], {5, 2}], {} -> 
           Pi]] &, PixelValuePositions[EdgeDetect[img], 1], 1300]
    ]) // Map@Length

Graphics[
 GraphicsComplex[pts, {FaceForm[Red, Blue], paths // Polygon}], 
 ImageSize -> Small]
{xInt, yInt} = 
  Interpolation[#, InterpolationOrder -> 1] & /@ 
   Transpose@pts[[First@paths]];

nn = 2^12;
xFC = Fourier[xInt[Subdivide[1., 380., nn]]]/Sqrt[nn];
yFC = Fourier[yInt[Subdivide[1., 380., nn]]]/Sqrt[nn];

linkage[fc_, tt_, offset_, pt_] := With[{
    line = ReIm[
      I*Accumulate[
         Prepend[ConstantArray[2, Length@fc - 1], 1]*
          fc*
          Exp[2 Pi I Range[0, Length@fc - 1] (-(tt - 1)/379)]
         ] +
       offset]
    },
   Graphics[
    {Line@line, {Thin, 
      MapThread[Circle, {Most@line, Norm /@ Differences@line}]},
     Red, {Opacity[0.5], Line[{Last@line, pt}]},
     PointSize@Medium, Point@pt, PointSize@Small, Point@Last@line}
    ]
   ];

Manipulate[
 With[{xFCtrunc = Take[xFC, order + 1], 
   yFCtrunc = Take[yFC, order + 1]},
  With[{xT = Total[
      Prepend[ConstantArray[2, order], 1]*
       Abs[xFCtrunc]*
       Sin[
        2 Pi Range[0, order] (-(t - 1)/379) + Arg[xFCtrunc] + Pi/2]
      ],
    yT = Total[
      Prepend[ConstantArray[2, order], 1]*
       Abs[yFCtrunc]*
       Sin[
        2 Pi Range[0, order] (-(t - 1)/379) + Arg[yFCtrunc] + Pi/2]
      ]},
   Show[
    ParametricPlot[{{xInt[t], yInt[t]}, {xT, yT}}, {t, 1, 380}, 
     PlotStyle -> {AbsoluteThickness[3], AbsoluteThickness[1.6]}],
    (* some abracadabra to rotate things around to
     * the desired positions next to the curve *)
    linkage[-I * Conjugate@xFCtrunc, -tt, -100 I - 
      I * Conjugate@xFCtrunc[[1]], {xT, yT} /. t -> tt],
    linkage[yFCtrunc, tt, -550 + yFCtrunc[[1]], {xT, yT} /. t -> tt],
    PlotRange -> {{-600, 350}, {-400, 500}}]
   ]],
 {{tt, 1., HoldForm[t]}, 1., 380.},
 {{order, 2}, 1, 100, 1}]

enter image description here

Discussion: Note the first Manipulate (and the first examples in the video and in @kglr's answer) are 1D Fourier series. The plane (i.e. 2D) curve comes from the inclusion of the independent variable t as a coordinate-dimension. In the 2D example above, t is omitted as a coordinate and serves only as parameter. Moving to 3D poses no problems for computing Fourier series, but as @DavidG.Stork points out in the comments below, the difficulty is deciding how to visualize the interactions of the linkages effectively. Likewise, moving to 4D and higher is mathematically easy, but the problem of visualization grows in complexity. The Fourier series for each coordinate $u$ in this approach in effect determines a (hyper)plane $u=u(t_1)$, in which the $u$-component of the end of the linkage must lie. The linkage itself lies in a plane. One dimension corresponds to $u$ and the plane must be parallel to or contain the $u$ axis. The other dimension seems to have no mathematical significance and can be arranged to correspond to any direction perpendicular to the direction of the $u$ axis.

3D, since David asked about it

For an example, we approximate a random polygonal path. It takes until the order 6 series to get the basic shape, and order 7 unexpectedly starts to approximate the straight lines.

Order 6 Order 7
SeedRandom[1];
path = Append[#, First@#] &@RandomReal[{-0.25, 1.8}, {10, 3}];
pathIFN = 
  Interpolation[MapIndexed[{(#2 - 1.)/(Length@path - 1), #1} &, path],
    PeriodicInterpolation -> True, InterpolationOrder -> 1];

nn = 2^12;
samplePts = pathIFN[Subdivide[0., 1., nn]];
xFC = Fourier[samplePts[[All, 1]]]/Sqrt[nn];
yFC = Fourier[samplePts[[All, 2]]]/Sqrt[nn];
zFC = Fourier[samplePts[[All, 3]]]/Sqrt[nn];

ClearAll[circle3D, disk3D, circlePoints3D];
circlePoints3D[center_List, radii_List, angles_List : {0, 2 \[Pi]}] :=
  Table[center + {Cos[theta], Sin[theta]} . radii, {theta, 
    N@First@angles, N@Last@angles, (Last@angles - First@angles)/
    Round[60 (Last@angles - First@angles)/(2 Pi)]}];
circle3D[center_List, radii_List, angles_List : {0, 2 \[Pi]}] := 
  Line@circlePoints3D[center, radii, angles];
disk3D[center_List, radii_List] := 
  Polygon@Most@circlePoints3D[center, radii];

linkage3D // ClearAll;
linkage3D[fc_, tt_, offset_, plane_, pt_] := With[{
    line = ReIm[
       I*Accumulate[
          Prepend[ConstantArray[2, Length@fc - 1], 1]*
           fc*
           Exp[2 Pi I Range[0, Length@fc - 1] (-tt)]
          ] +
        offset] . plane (* plane should be orthonormal *)
    (*,radii=Orthogonalize[N@plane]*)
    },
   {Line@
     line, {MapThread[{ColorData[97][2 + #1], 
        circle3D[##2]} &, {Range[Length@line - 1], Most@line, 
       plane # & /@ Norm /@ Differences@line}]},
    Magenta, {Opacity[0.5], 
     Line[{Last@line, Last@line + Projection[pt, Cross @@ plane], 
       pt}]},
    PointSize@Medium, Point@pt, PointSize@Small, Point@Last@line}
   ];

plot = ParametricPlot3D[pathIFN[t], {t, 0, 1}, 
   PlotStyle -> Directive[ColorData[97, 2], AbsoluteThickness[1.]], 
   AxesLabel -> {x, y, z}];
Manipulate[
 With[{
   xFCtrunc = Take[xFC, order + 1],
   yFCtrunc = Take[yFC, order + 1],
   zFCtrunc = Take[zFC, order + 1]},
  With[{xT = Total[
      Prepend[ConstantArray[2, order], 1]*
       Abs[xFCtrunc]*
       Sin[2 Pi Range[0, order] (-t) + Arg[xFCtrunc] + Pi/2]
      ],
    yT = Total[
      Prepend[ConstantArray[2, order], 1]*
       Abs[yFCtrunc]*
       Sin[2 Pi Range[0, order] (-t) + Arg[yFCtrunc] + Pi/2]
      ],
    zT = Total[
      Prepend[ConstantArray[2, order], 1]*
       Abs[zFCtrunc]*
       Sin[2 Pi Range[0, order] (-t) + Arg[zFCtrunc] + Pi/2]
      ]},
   Show[
    plot,
    ParametricPlot3D[{xT, yT, zT}, {t, 0, 1}, 
     PlotStyle -> AbsoluteThickness[2.5], AxesLabel -> {x, y, z}]
    ,
    Graphics3D[Dynamic@{
       AbsoluteThickness[1.6],
       linkage3D[-I*xFCtrunc, tt, -1.5 I + I*xFCtrunc[[1]], 
        N@{{1, 0, 0}, {0, 1, 0}}, {xT, yT, zT} /. t -> tt],
       linkage3D[-I*yFCtrunc, tt, -1.2 I + I*yFCtrunc[[1]], 
        N@{{0, 1, 0}, {0, 0, 1}}, {xT, yT, zT} /. t -> tt],
       linkage3D[-I*zFCtrunc, tt, -1.4 I + I*zFCtrunc[[1]], 
        N@{{0, 0, 1}, {1, 0, 0}}, {xT, yT, zT} /. t -> tt]
       }]
    ,
    PlotRange -> 1.6, ViewPoint -> {2.2, 2.6, 1.8}]
   ]],
 {{tt, 0., HoldForm[t]}, 0., 1.},
 {{order, 6}, 1, 20, 1}]
$\endgroup$
8
  • 1
    $\begingroup$ Oh very very nice... especially the static graphical representation of the Fourier components. ($+1$) $\endgroup$ Jun 30 at 14:15
  • $\begingroup$ A more general question: Is there code for synthesizing a curve in three-dimensions, presumably by chaining together rotating spheres? $\endgroup$ Jun 30 at 14:54
  • $\begingroup$ @DavidG.Stork If you want a Fourier approximation to a given closed path in 3D, then one could find the Fourier series of each of the three (periodic) coordinates. Then one could have three circles instead of two circles as in the video you linked. The plane of each circle would have to contain the point on the path, so each plane would have to rotate as the point moves. So I guess each circle could be represented by a sphere. I don't know of any code for that. Nor have I ever seen such an animation. Is that what you had in mind? $\endgroup$
    – Michael E2
    Jun 30 at 22:02
  • $\begingroup$ Thanks for the animations, @Nasser. :) $\endgroup$
    – Michael E2
    Jun 30 at 22:03
  • $\begingroup$ Michael E2: Yes... that is basically what I had in mind. I indeed know how to take three Fourier transforms of a three-dimensional curve—one in each dimension. The question is how to best visualize the total, integrated result. The single (complex) Fourier transform and attendant circles work in two spatial dimensions because the horizontal and vertical dimensions are naturally represented by the real and the imaginary parts of a complex function. Generalizing this idea one dimension further is a bit tricky. I've never seen it done, but I think it would be very cool indeed. $\endgroup$ Jun 30 at 22:12

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.