Power spectra of periodic, quasi-periodic, and chaotic signals.

Sinusoidal time-series

Clear[sr,frequency,t,ff,ft,tmax]
tmax=200;
frequency=19;

FourierTransform[Sin[frequency t],t,ω](*fourier transform, analytic*)

sint=Plot[Sin[frequency t],{t,0,10},FrameLabel-> {t,"sin(ωt)"},PlotRange-> All,ImageSize->320,Frame->True,LabelStyle->20];
sr=63;(*Sampling Rate: keep Nyquist limit in mind*)
signal=Table[38 Sin[frequency t],{t,0,tmax,1/sr}];(*sample values*)
ft=Fourier[signal,FourierParameters->{-1,-1}];(*num fourier transform*)
ff=Table[2 π (n-1) sr/Length[signal],{n,Length[signal]}];(*frequency values x-axis*)
sinspecter=ListLinePlot[Transpose[{ff,Log[Abs[ft]^2]}][[1;;Round[0.5 Length[signal]]]],PlotRange->All,Frame->True,FrameLabel->{"Frequency ω","Spectral Power Level"},PlotStyle->{Thick,Red},ImageSize->530,PlotLabel->Style["Power spectrum of sin(ωt)",Large],LabelStyle->20];
Row[{sint,sinspecter},"        "]

Lorenz system, periodic solution (ignore initial transient behaviour)

Clear[d,r,b,sr,t,ff,ft,tmax]
tmax=80;
x0=1;y0=-44;z0=300;
d=10;r=350;b=8/3;(*periodic solution*)

soln=NDSolve[{
x'[t]==d (y[t]-x[t]),
y'[t]==r x[t]-y[t]-x[t] z[t],
z'[t]==x[t] y[t]-b z[t],
x[0]==x0,y[0]==y0,z[0]==z0},
{x,y,z},{t,0,tmax},MaxSteps->Infinity];
diffx=d (y-x);diffy=r x-y-x z;diffz=x y-b z;
fixedpoints=Solve[{diffx==0,diffy==0,diffz==0},{x,y,z}];
Print[Style["Coordinates of fixed points",18,Blue]]
coordinates={x,y,z}/.fixedpoints
Gfp=Graphics3D[{PointSize[0.03],Red,Point[coordinates]}];

pplot=ParametricPlot3D[Evaluate[{x[t],y[t],z[t]}/.soln],{t,0,tmax},PlotPoints->500,AxesLabel->{x[t],y[t],z[t]},BoxRatios->{1, 1, 1},PlotStyle->Thin,ImageSize->400,LabelStyle->20];
px=Plot[Evaluate[x[t]/.soln],{t,0,tmax},FrameLabel-> {t,x[t]},PlotRange-> All,ImageSize->330,Frame->True,LabelStyle->20];
Row[{Show[pplot,Gfp],px},"    "]

sr=63;(*Sampling Rate: keep Nyquist limit in mind*)
signal=Table[Evaluate[x[t]/.soln[[1]]],{t,0,tmax,1/sr}];(*sample values*)
ft=Fourier[signal,FourierParameters->{-1,-1}];(*num fourier transform*)
ff=Table[2 π (n-1) sr/Length[signal],{n,Length[signal]}];(*frequency values x-axis*)
pspecter=ListLinePlot[Transpose[{ff,Log[Abs[ft]^2]}][[1;;Round[0.5 Length[signal]]]],PlotRange->All,Frame->True,FrameLabel->{"Frequency ω","Spectral Power Level"},ImageSize->530,PlotLabel->Style["Power spectrum of x(t), periodic",Large],LabelStyle->20]

Quasi-periodic solution, trajectories on a surface of a torus (see Lecture 8)
Note: No quasi-periodic solution in Lorenz system (see Lecture 9)

Clear[sr,t,ff,ft,tmax]
k1=1.1; k2=2.0;(*coupling constants*)
ω1=5.02; ω2=1.046766352;(*ω natural frequency*)
tmax=300;

Print[Style["Ratio of fundamental frequencies",Red,17]]
ω1/ω2
ww=Abs[ω1-ω2] ;kk=k1+k2;
If[ww>=kk,
Print[Style["Solution is quasi-periodic (coupled cases only)",Red,17]],
Print[Style["Solution is periodic or long term sol'n is (coupled cases only)",Blue,17]]]
integ=NDSolve[{
θ1'[t]==ω1+k1 Sin[θ2[t]-θ1[t]],
θ2'[t]==ω2+k2 Sin[θ1[t]-θ2[t]],
θ1[0]==0,θ2[0]==0},{θ1,θ2},{t,0,tmax}];
sinp1=Plot[Evaluate[Sin[θ1[t]]/.integ],{t,0,tmax},Frame->True,FrameLabel-> {t,Sin[θ1]},PlotRange-> All,ImageSize->290,LabelStyle->20];
sinp2=Plot[Evaluate[Sin[θ2[t]]/.integ],{t,0,tmax},Frame->True,FrameLabel->{t,Sin[θ2]},PlotRange-> All,ImageSize->290,LabelStyle->20];

sr=11;(*Sampling Rate: keep Nyquist limit in mind*)
signal=Table[Evaluate[Sin[θ1[t]]/.integ[[1]]],{t,0,tmax,1/sr}];(*sample values*)
signal2=Table[Evaluate[Sin[θ2[t]]/.integ[[1]]],{t,0,tmax,1/sr}];(*sample values*)
ft=Fourier[signal,FourierParameters->{-1,-1}];(*num fourier transform*)
ft2=Fourier[signal2,FourierParameters->{-1,-1}];(*num fourier transform*)
ff=Table[2 π (n-1) sr/Length[signal],{n,Length[signal]}];(*frequency values x-axis*)
qpspecter1=ListLinePlot[Transpose[{ff,Log[Abs[ft]^2]}][[1;;Round[0.5 Length[signal]]]],PlotRange->All,Frame->True,FrameLabel->{"Frequency ω","Spectral Power Level"},ImageSize->530,LabelStyle->20];
qpspecter2=ListLinePlot[Transpose[{ff,Log[Abs[ft2]^2]}][[1;;Round[0.5 Length[signal]]]],PlotRange->All,Frame->True,FrameLabel->{"Frequency ω","Spectral Power Level"},PlotLabel->Style["Power spectrum of sin(ω(t)), quasi-periodic",Large],ImageSize->530,LabelStyle->20];
Row[{sinp1,qpspecter1},"   "]
Row[{sinp2,qpspecter2},"   "]

Lorenz system, chaotic solution

Clear[d,r,b,sr,tmax,ff,ft]
tmax=100;
x0=6.1;y0=6;z0=10;
d=10;r=28;b=8/3;

soln=NDSolve[{
x'[t]==d (y[t]-x[t]),
y'[t]==r x[t]-y[t]-x[t] z[t],
z'[t]==x[t] y[t]-b z[t],
x[0]==x0,y[0]==y0,z[0]==z0},
{x,y,z},{t,0,tmax},MaxSteps->Infinity];
diffx=d (y-x);diffy=r x-y-x z;diffz=x y-b z;
fixedpoints=Solve[{diffx==0,diffy==0,diffz==0},{x,y,z}];
Print[Style["Coordinates of fixed points",18,Blue]]
coordinates={x,y,z}/.fixedpoints
Gfp=Graphics3D[{PointSize[0.03],Red,Point[coordinates]}];

pplot=ParametricPlot3D[Evaluate[{x[t],y[t],z[t]}/.soln],{t,0,tmax},PlotPoints->1000,AxesLabel->{x[t],y[t],z[t]},PlotStyle->Thin,ImageSize->300,LabelStyle->20];
px=Plot[Evaluate[x[t]/.soln],{t,0,tmax},FrameLabel-> {t,x[t]},PlotRange-> All,ImageSize->330,Frame->True,LabelStyle->20];
(*py=Plot[Evaluate[y[t]/.soln],{t,0,tmax},AxesLabel-> {t,y[t]},PlotRange-> All,ImageSize->230];
pz=Plot[Evaluate[z[t]/.soln],{t,0,tmax},AxesLabel-> {t,z[t]},PlotRange-> All,ImageSize->230];*)
Row[{Show[pplot,Gfp],px},"     "]

sr=20;(*Sampling Rate: keep Nyquist limit in mind*)
signal=Table[Evaluate[x[t]/.soln[[1]]],{t,0,tmax,1/sr}];(*sample values*)
ft=Fourier[signal,FourierParameters->{-1,-1}];(*num fourier transform*)
ff=Table[2 π (n-1) sr/Length[signal],{n,Length[signal]}];(*frequency values x-axis*)
cspecter=ListLinePlot[Transpose[{ff,Log[Abs[ft]^2]}][[1;;Round[0.5 Length[signal]]]],PlotRange->All,Frame->True,FrameLabel->{"Frequency ω","Spectral Power Level"},PlotStyle->Thickness[0.002],ImageSize->530,PlotLabel->Style["Power spectrum of x(t), chaotic",Large],LabelStyle->20]

Comparison of above spectra

Column[{Show[pspecter,sinspecter],qpspecter2,cspecter},Spacings->4]

Creation of this notebook was supported by IT Academy program of Information Technology Foundation for Education (HITSA).