Class blog 2/13

ECE 3340 - Han Q. Le (r)

1. Make a Fourier-series coefficient function

We did this in class. To avoid error, we USE variable name Λ for period instead of 2. Why? Even an experienced programmer will MAKE mistakes when having to change numerical values from one problem to the next.

(1/Λ)(∫_(-(Λ/2))^(Λ/2))(x^3) (E^(-((I 2 π m x)/Λ)))\[DifferentialD]x

Substitute

((I (Λ^3) (m π (-6+(m^2) (π^2)) Cos[m π]-3 (-2+(m^2) (π^2)) Sin[m π]))/(8 (m^4) (π^4)))/.Sin[m π]-> 0

Substitute again:

((I (-6+(m^2) (π^2)) (Λ^3) Cos[m π])/(8 (m^3) (π^3))) /. Cos[m π]-> (((-1))^m)

We still have to solve the problem of m=0:

Limit[((I (Λ^3) (m π (-6+(m^2) (π^2)) Cos[m π]-3 (-2+(m^2) (π^2)) Sin[m π]))/(8 (m^4) (π^4))),m-> 0]

We know that it should be zero as expected for an odd function.

Npw, we define this function in class (but here we include period Λ

coefcube[m_,Λ_]:=If[m==0,0,((I (((-1))^m) (-6+(m^2) (π^2)) (Λ^3))/(8 (m^3) (π^3)))];

But to get an array we have to do this:

mymarray={1,5,2,6,9,0,3};
Map[coefcube[#,Λ]&, mymarray]

What we have done here is that we allow Λ to be an unknown variable. The expression:
                            coefcube[#,Λ] &

is known as an anonymous function, a most seminal concept in the history of computer programming development. Read about it in wikipedia here. Anonymous function.

If you are not ready to learn it, it is OK for now, don’t worry. If you alreaydy know about lambda function in Python, then this is the syntax in Mathematica. For those using MATLAB, look it up.

You will see that I use anonymous function every in my codes. Example, see the highlight below:

fouriergraphics2[{funcArtnm_,σ_,styleparm_}, nFouriercomp_,nintv_]:=Block[{coef0
,coef,xpltarr,fexpn,outFS}
,
coef0=Chop[Table[Quiet[NIntegrate[funcArtnm[x,σ]*(E^(-I 2 π n x)),{x,-0.5,0.5}
,AccuracyGoal->3 ]]
,{n,0,nFouriercomp}],(10^(-4)) ];
If[Length[Dimensions[coef0]]≥ 2 && Depth[coef0]≥   3
,coef=Join[{1},ConstantArray[2,nFouriercomp]]* # & /@ Transpose[coef0]
,coef=Join[{1},ConstantArray[2,nFouriercomp]]*coef0;
];
xpltarr=Range[-0.5,0.5,1./nintv];
fexpn=Table[(E^(I 2 π n xpltarr )),{n,0,nFouriercomp}];
outFS=Re[coef . fexpn ] ;
ListLinePlot[outFS,DataRange-> {-0.5,0.5},PlotRange-> All
,PlotStyle->styleparm[[1]]
,Filling->styleparm[[2]]
,AspectRatio->Automatic
]
];

Programmers can’t live without it.

For now, if you don’t understand, just skip and don’t worry about it. Just treat the whole expression
             (coefcube[#,Λ] &)  -> a blackbox function”

One student in the class wanted to know how to put a whole array in coefcube function. We can redesign the function:

coefcube[{m_,Z_}]:=Module[{q,outc},
(* q is a local variable, outc is also a local var  *)
If[!ListQ[m]  (* if m is NOT a list, i. e. just a number,  *)
,q={m}  (* we make q to be the list of ONE element m  *)
,q=m]; (* if mm is a list, then q = m  *)

(*now define output outc to be a list output using Table *)
outc=Table[If[k==0,0,((I (((-1))^k) (-6+(k^2) (π^2)) (Z^3))/(8 (k^3) (π^3)))],{k,q}]; 
(* in m is a list, we return a list, if m is a number, 
we return a number  *)
If[ListQ[m],outc,outc[[1]]]
]

It is a more versatile function, because we can put in just a number or a list of number, it will be OK:

coefcube[{3.,Λ}]  (* just a number, it is 3 *)

Or:

mymarray={1,5,2,6,9,0,3};
coefcube[{mymarray,Λ}]

See the above? input an array, you get an array, input a number, you get a number. Kind for kind.

By now, hope you appreciate the capability of Mathematica to handle symbolic manipulation. We don’t have a value for Λ, but it works just fine. We can’t do that in most other languages.

2. Use Mathematica built-in function

FourierCoefficient[(x^3),x,n]

The problem? you have to be aware that the period is -π to π. Look it up:

The n(^th) coefficient in the Fourier series expansion of f(t) is by default given by (1/(2π))(∫_(-π)^π)f(t) (e^(-i n t))dt.

Too often, people forget to scale their function period to {-π to π}. Just like I forgot to scale our period to 2 in the class. Why, because I had chosen period = 1 (-0.5 to 0.5) for so many examples that when we change to [-1, 1], I forgot to insert the period. This is the reason to always use a variable name, not a numerical value of that variable. I forgot to change from 1 -> 2, it cost some time to find out. If we use Λ from the  beginning, the mistake would not have happened.

3. Try your hand at it

Below are the 4 codes of 4 common waveforms, taken from one of my tutorial files in Fourier “Lect_3340_LTI_FFT_utilities.nb”

Take a shot and create your own function, see if you can modify the code below (just substitute your function Fourier coefficient in function xFcoeff below). and obtain a demo that you enjoy.
In section 4, I show some results with half-moon bridge.

xFcoeff[m_,wfparm_]:=Block[{wf,η}
,{wf,η}=Take[wfparm,2]; 
If[wf=="square"
,(E^(-I m π η))  η Sinc[m π η] (* coef for square *)
,If[wf=="sawtooth"
,(E^(-2 I m π η))If[m==0,0.5 η,(( I (1-(E^(I m π η)) Sinc[m π η]))/(2 m π))] (* coef for sawtooth *)
,If[wf=="triangle"
,(η/2.) (E^(- I m π η))((Sinc[((m π η)/2)])^2) (* coef for triangle *)
,(η/2) (1-(E^(-(2/η)))) (1/(1+ I m π η))  (* coef for exponential *)
]]]
];

5.1.2 For computation - Square wave

Demo only. You SHOULD extract, modify, adapt, do whichever needed for your own work.

Manipulate[
Grid[{
{Item[Style["waveform",18,Black,FontFamily->"Century Gothic"],Background->Hue[0.35,0.4,1]]
,Item[Style["Fourier decomposition",18,Black,FontFamily->"Century Gothic"],Background->Hue[0.15,0.4,1]]}
,
{Plot[(UnitStep[t]-UnitStep[t-η T]) /. t-> Mod[u,T]
,{u,0,2.5},PlotRange-> {All,{-0.1,1.1}}
,  PlotStyle->{Thickness[0.0065],Hue[0.,0.85,1]}
,Filling->{1-> {0,Hue[0.05,0.85,1,0.4]}}
,AspectRatio->0.35, Background->Black,ImageSize->300]
,Column[{Show[fourierSumplot[k,{"square",η,T},tarr]
,PlotRange-> {All,{-0.1,1.1}}
,AspectRatio->0.35, Background->Black,ImageSize->300]
,Show[fourierCompplot1[k,{"square",η,T},tarr]
,AspectRatio->1.5, Background->Black,ImageSize->300]
}]
}},Alignment-> Top,Frame-> True,FrameStyle->White]
,{{η,0.25},0.1,0.9},{{T,1},0.5,2},{{k,16},0,32,1}
,Initialization:>{
tarr=Range[0,500]/200.;
bogusgraphics=Graphics[{Point[{-1,0}]}];
xFcoeff[m_,wfparm_]:=Block[{wf,η}
,{wf,η}=Take[wfparm,2]; 
If[wf=="square"
,(E^(-I m π η))  η Sinc[m π η]
,If[wf=="sawtooth"
,(E^(-2 I m π η))If[m==0,0.5 η,(( I (1-(E^(I m π η)) Sinc[m π η]))/(2 m π))]
,If[wf=="triangle"
,(η/2.) (E^(- I m π η))((Sinc[((m π η)/2)])^2)
,(η/2) (1-(E^(-(2/η)))) (1/(1+ I m π η))
]]]
];

fourierCompplot1[k_,wfparm_,vararr_]:=Block[{wf,η,tper,xLoc}
,wf=wfparm[[1]]; 
If[TrueQ[wf=="square"],{η,tper}=Take[wfparm,{2,3}] ; 
xLoc=vararr/tper;
If[k==0
,Graphics[{Thickness[0.0065],White,Line[{{xLoc[[1]],η},{Last[xLoc],η}}]}]
,Graphics[{Thickness[0.0065]
,Table[{Hue[(((m-1))/(k+1)),0.8,1]
,Line[Transpose[{xLoc
,2η Sinc[m π η]Re[(E^(I π m (2 xLoc- η))) ] +1.5 Log[1.k/m]}]]},{m,k}]
}]]
,bogusgraphics]
];

fourierSumplot[k_,wfparm_,vararr_]:=Block[{wf,η,tper,coefF,xLoc}
,{wf,η,tper}=Take[wfparm,3]; 
xLoc=vararr/tper;
coefF=Table[xFcoeff[m,wfparm],{m,k}];
Graphics[{Thickness[0.0065], White
,Line[Transpose[{xLoc
,2Re[coefF . Table[(E^(I 2 π m xLoc )),{m,k}]]
+ConstantArray[xFcoeff[0,wfparm],Length[xLoc]]
}]]
}]
];

}]

5.1.2 For computation - Sawtooth

Demo only. You SHOULD extract, modify, adapt, do whichever needed for your own work.

Code

Manipulate[
Grid[{
{Item[Style["waveform",18,Black,FontFamily->"Century Gothic"],Background->Hue[0.35,0.4,1]]
,Item[Style["Fourier decomposition",18,Black,FontFamily->"Century Gothic"],Background->Hue[0.15,0.4,1]]}
,
{Plot[If[t≤ η T,(t/(η T)),0.]/. t-> Mod[u,T]
,{u,0,2.5},PlotRange-> {All,{-0.1,1.1}}
, PlotStyle->{Thickness[0.005],Hue[0.15,1,1]}
,Filling->{1-> {0,Hue[0.14,1,1,0.4]}}
,AspectRatio->0.35, Background->Black,ImageSize->300]
,Column[{Show[fourierSumplot[k,{"sawtooth",η,T},tarr]
,PlotRange-> {All,{-0.1,1.1}}
,AspectRatio->0.35, Background->Black,ImageSize->300]
,Show[fourierCompplot2[k,{"sawtooth",η,T},tarr]
,AspectRatio->1.5, Background->Black,ImageSize->300]
}]
}},Alignment-> Top,Frame-> True,FrameStyle->White]
,{{η,0.75},0.1,0.9},{{T,1},0.5,2},{{k,16},0,32,1}
,Initialization:>{
tarr=Range[0,500]/200.;
bogusgraphics=Graphics[{Point[{-1,0}]}];
xFcoeff[m_,wfparm_]:=Block[{wf,η}
,{wf,η}=Take[wfparm,2]; 
If[wf=="square"
,(E^(-I m π η))  η Sinc[m π η]
,If[wf=="sawtooth"
,(E^(-2 I m π η))If[m==0,0.5 η,(( I (1-(E^(I m π η)) Sinc[m π η]))/(2 m π))]
,If[wf=="triangle"
,(η/2.) (E^(- I m π η))((Sinc[((m π η)/2)])^2)
,(η/2) (1-(E^(-(2/η)))) (1/(1+ I m π η))
]]]
];

fourierCompplot2[k_,wfparm_,vararr_]:=Block[{wf,η,tper,xLoc}
,wf=wfparm[[1]]; 
If[TrueQ[wf=="sawtooth"],{η,tper}=Take[wfparm,{2,3}] ; 
xLoc=vararr/tper;
If[k==0
,Graphics[{Thickness[0.0065],White,Line[{{xLoc[[1]],0.5η},{Last[xLoc],0.5 η}}]}]
,Graphics[{Thickness[0.0065]
,Table[{Hue[(((m-1))/(k+1)),0.8,1]
,Line[Transpose[{xLoc
,2Re[(( I (1-(E^(I m π η)) Sinc[m π η]))/(2 m π))(E^(I 2π m (xLoc- η))) ] +1. Log[k/m]}]]},{m,k}]
}]]
,bogusgraphics]
];

fourierSumplot[k_,wfparm_,vararr_]:=Block[{wf,η,tper,coefF,xLoc}
,{wf,η,tper}=Take[wfparm,3]; 
xLoc=vararr/tper;
coefF=Table[xFcoeff[m,wfparm],{m,k}];
Graphics[{Thickness[0.0065], White
,Line[Transpose[{xLoc
,2Re[coefF . Table[(E^(I 2 π m xLoc )),{m,k}]]
+ConstantArray[xFcoeff[0,wfparm],Length[xLoc]]
}]]
}]
];

}]

5.1.2 For computation - Triangle ramp (symmetric ramp)

Demo only. You SHOULD extract, modify, adapt, do whichever needed for your own work.

Code

Manipulate[
Grid[{
{Item[Style["waveform",18,Black,FontFamily->"Century Gothic"],Background->Hue[0.35,0.4,1]]
,Item[Style["Fourier decomposition",18,Black,FontFamily->"Century Gothic"],Background->Hue[0.15,0.4,1]]}
,
{Plot[If[t≤0.5 η T,((2t)/(η T)),If[t≤ η T,((2(η T-t))/(η T)),0]]/. t-> Mod[u,T]
,{u,0,2.5},PlotRange-> {All,{-0.1,1.1}}
, PlotStyle->{Thickness[0.005],Hue[0.35,1,1]}
,Filling->{1-> {0,Hue[0.32,1,1,0.4]}}
,AspectRatio->0.35, Background->Black,ImageSize->300]
,Column[{Show[fourierSumplot[k,{"triangle",η,T},tarr]
,PlotRange-> {All,{-0.1,1.1}}
,AspectRatio->0.35, Background->Black,ImageSize->300]
,Show[fourierCompplot3[k,{"triangle",η,T},tarr]
,AspectRatio->1.5, Background->Black,ImageSize->300]
}]
}},Alignment-> Top,Frame-> True,FrameStyle->White]
,{{η,0.5},0.1,0.9},{{T,1},0.5,2},{{k,16},0,32,1}
,Initialization:>{
tarr=Range[0,500]/200.;
bogusgraphics=Graphics[{Point[{-1,0}]}];
xFcoeff[m_,wfparm_]:=Block[{wf,η}
,{wf,η}=Take[wfparm,2]; 
If[wf=="square"
,(E^(-I m π η))  η Sinc[m π η]
,If[wf=="sawtooth"
,(E^(-2 I m π η))If[m==0,0.5 η,(( I (1-(E^(I m π η)) Sinc[m π η]))/(2 m π))]
,If[wf=="triangle"
,(η/2.) (E^(- I m π η))((Sinc[((m π η)/2)])^2)
,(η/2) (1-(E^(-(2/η)))) (1/(1+ I m π η))
]]]
];
fourierCompplot3[k_,wfparm_,vararr_]:=Block[{wf,η,tper,xLoc}
,wf=wfparm[[1]]; 
If[TrueQ[wf=="triangle"],{η,tper}=Take[wfparm,{2,3}] ; 
xLoc=vararr/tper;
If[k==0
,Graphics[{Thickness[0.0065],White,Line[{{xLoc[[1]],0.5η},{Last[xLoc],0.5 η}}]}]
,Graphics[{Thickness[0.0065]
,Table[{Hue[(((m-1))/(k+1)),0.8,1]
,Line[Transpose[{xLoc
,η ((Sinc[((m π η)/2)])^2) Re[(E^(I π m (2xLoc- η))) ] +1. Log[k/m]}]]},{m,k}]
}]]
,bogusgraphics]
];

fourierSumplot[k_,wfparm_,vararr_]:=Block[{wf,η,tper,coefF,xLoc}
,{wf,η,tper}=Take[wfparm,3]; 
xLoc=vararr/tper;
coefF=Table[xFcoeff[m,wfparm],{m,k}];
Graphics[{Thickness[0.0065], White
,Line[Transpose[{xLoc
,2Re[coefF .Table[(E^(I 2 π m xLoc )),{m,k}]]
+ConstantArray[xFcoeff[0,wfparm],Length[xLoc]]
}]]
}]
];

}]

5.1.2 For computation - Exponential

Demo only. You SHOULD extract, modify, adapt, do whichever needed for your own work.

Code

Manipulate[
Grid[{
{Item[Style["waveform",18,Black,FontFamily->"Century Gothic"],Background->Hue[0.35,0.4,1]]
,Item[Style["Fourier decomposition",18,Black,FontFamily->"Century Gothic"],Background->Hue[0.15,0.4,1]]}
,
{Plot[(E^(-2 t/(η T)))/. t-> Mod[u,T]
,{u,0,2.5},PlotRange-> {All,{-0.1,1.1}}
, PlotStyle->{Thickness[0.005],Hue[0.55,1,1]}
,Filling->{1-> {0,Hue[0.5,1,1,0.4]}}
,AspectRatio->0.35, Background->Black,ImageSize->300]
,Column[{Show[fourierSumplot[k,{"exponential",η,T},tarr]
,PlotRange-> {All,{-0.1,1.1}}
,AspectRatio->0.35, Background->Black,ImageSize->300]
,Show[fourierCompplot4[k,{"exponential",η,T},tarr]
,AspectRatio->1.5, Background->Black,ImageSize->300]
}]
}},Alignment-> Top,Frame-> True,FrameStyle->White]
,{{η,0.25},0.1,0.9},{{T,1},0.5,2},{{k,16},0,32,1}
,Initialization:>{
tarr=Range[0,500]/200.;
bogusgraphics=Graphics[{Point[{-1,0}]}];
xFcoeff[m_,wfparm_]:=Block[{wf,η}
,{wf,η}=Take[wfparm,2]; 
If[wf=="square"
,(E^(-I m π η))  η Sinc[m π η]
,If[wf=="sawtooth"
,(E^(-2 I m π η))If[m==0,0.5 η,(( I (1-(E^(I m π η)) Sinc[m π η]))/(2 m π))]
,If[wf=="triangle"
,(η/2.) (E^(- I m π η))((Sinc[((m π η)/2)])^2)
,(η/2) (1-(E^(-(2/η)))) (1/(1+ I m π η))
]]]
];

fourierCompplot4[k_,wfparm_,vararr_]:=Block[{wf,η,tper,xLoc}
,wf=wfparm[[1]]; 
If[TrueQ[wf=="exponential"],{η,tper}=Take[wfparm,{2,3}] ; 
xLoc=vararr/tper;
If[k==0
,Graphics[{Thickness[0.0065],White,Line[{{xLoc[[1]],0.5η},{Last[xLoc],0.5 η}}]}]
,Graphics[{Thickness[0.0065]
,Table[{Hue[(((m-1))/(k+1)),0.8,1]
,Line[Transpose[{xLoc
,η (1-(E^(-(2/η))))Re[((E^(I 2π m xLoc))/(1+ I m π η)) ] + Log[1.k/m]}]]},{m,k}]
}]]
,bogusgraphics]
];

fourierSumplot[k_,wfparm_,vararr_]:=Block[{wf,η,tper,coefF,xLoc}
,{wf,η,tper}=Take[wfparm,3]; 
xLoc=vararr/tper;
coefF=Table[xFcoeff[m,wfparm],{m,k}];
Graphics[{Thickness[0.0065], White
,Line[Transpose[{xLoc
,2Re[coefF . Table[(E^(I 2 π m xLoc )),{m,k}]]
+ConstantArray[xFcoeff[0,wfparm],Length[xLoc]]
}]]
}]
];

}]

4. Half-moon bridge

Plot[{( (0.5^2) -(x^2)),If[Abs[x]≤ 0.2,(Sqrt[0.2^2-x^2]),0] }
,{x,-0.5,0.5}, AspectRatio->Automatic
,Filling-> {1-> {{2},Hue[0.6,0.3,0.9]}}
,Axes->False]

moonbridge[x_,σ_]:={( ((((1/2)))^2) -(x^2)),If[Abs[x]≤(1/5),(Sqrt[((((1/5)))^2)-x^2]),0]};

(∫_(-1/2)^(1/2))moonbridge[x,1](E^(-I m 2 π x))\[DifferentialD]x

((-m π Cos[m π]+Sin[m π])/(2 (m^3) (π^3))) /. Sin[m π]-> 0  /. Cos[m π]-> (((-1))^m)

It’s nice to learn about BesselJ function. The two coefficients are:

{-((((-1))^m)/(2 (m^2) (π^2))),((BesselJ[1, ((2 m π)/5)])/(10 m))}

Limit[{((-m π Cos[m π]+Sin[m π])/(2 (m^3) (π^3))),((BesselJ[1,((2 m π)/5)])/(10 m))},m-> 0]

coef[m_]:=If[m==0,{(1/6),(π/50)},{-((((-1))^m)/(2 (m^2) (π^2))),((BesselJ[1,((2 m π)/5)])/(10 m))}];

nchan=32; plotstyle={{Thick,Hue[0.5,0.4,0.4]}};
filling= {1-> {{2},Hue[0.6,0.3,0.9]}};
xarr=Range[-0.5,0.5,1/256];
coefnum0=Transpose @ Table[ coef[m],{m,0,nchan}];
coefnum=Join[{1},ConstantArray[2,nchan]] * # & /@ coefnum0 ;
fExpn=Table[(E^(I 2 π m xarr)),{m,0,nchan}];
fouriersynth=Re[coefnum . fExpn]  ;
imgfourier=ListLinePlot[fouriersynth, DataRange->{-0.5,0.5}
,PlotStyle->plotstyle
,Filling-> filling
,Axes->False
,AspectRatio->Automatic
,ImageSize->400]

Code

Manipulate[
Grid[{{imgfourier
,Manipulate[Column[{fourierprogGr[[i]]
,fouriercompGr[[i]]}]
,{i,1,Length[fourierprogGr],1}]
,Column[fouriercompGr]
}}, Alignment->Top]
,Initialization:> {
nchan=32; plotstyle={{Thick,Hue[0.5,0.4,0.4]}};
filling= {1-> {{2},Hue[0.6,0.3,0.9]}};
xarr=Range[-0.5,0.5,1/256];
coefnum0=Transpose @ Table[ coef[m],{m,0,nchan}];
coefnum=Join[{1},ConstantArray[2,nchan]] * # & /@ coefnum0 ;
fExpn=Table[(E^(I 2 π m xarr)),{m,0,nchan}];
fouriersynth=Re[coefnum . fExpn]  ;
imgfourier=ListLinePlot[fouriersynth, DataRange->{-0.5,0.5}
,PlotStyle->plotstyle
,Filling-> filling
,Axes->False
,AspectRatio->Automatic
,ImageSize->300];
fouriercomp=Table[ Re[fExpn[[k]] * # & /@ coefnum[[ ;;   ,k ]]],{k,Length[fExpn]}];
fouriercompGr=Table[ListLinePlot[f, DataRange->{-0.5,0.5}
,PlotStyle->plotstyle
,Filling-> filling
,Axes->False
,AspectRatio->Automatic
,ImageSize->200],{f,fouriercomp}];
fourierprogGr=Table[
fs=Re[coefnum [[;; ,  1 ;;  k ]]. Take[fExpn,k]]  ;
ListLinePlot[fs, DataRange->{-0.5,0.5}
,PlotStyle->plotstyle
,Filling-> filling
,Axes->False
,AspectRatio->Automatic
,ImageSize->200]
,{k,1,Length[fExpn]}];

}]