Demonstration of Fourier Series: Selectable Target Function, Selectable Number of Coefficients, Scrolling Graph

Question

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

Answer 1

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}]

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}]

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}]

The reason behind \[FreakedSmiley]s around sin:

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

Answer 2

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

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

---

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}]

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}]