A Foray into Fourier analysis

Many coupled oscillators... again.

Setting up the equations for many oscillators, we can solve the system

numerically ...

Remove["Global`*"]

n=7;
A=Table[0,{n},{n}];
A[[1,1]]=-2;
A[[1,2]]=1;
A[[n,n-1]]=1;
A[[n,n]]=-2;

Do[
A[[i,i-1]]=1;
A[[i,i]]=-2;
A[[i,i+1]]=1,
{i,2,n-1}]

A=A k;
A//MatrixForm

vecX=Table[(Subscript[x, i])[t],{i,1,n}]

eqn=Table[m D[vecX[[i]],{t,2}]==A[[i]].vecX,{i,1, n}]

L=1;
f[x_]:=Exp[-((((x-L/2))^2)/0.01)]
Plot[f[x],{x,0,1}]

initpos=Table[f[(i-1)(L/((n-1)))],{i,1,n}]

init=Flatten[Table[{(Subscript[x, i])[0]==f[(i-1)(L/((n-1)))],(Subscript[x, i])'[0]==0},{i,1,n}]]

k=1;
m=1;
T=80;
sol=NDSolve[Flatten[{eqn,init}],vecX,{t,0,T}][[1]];

Xsol=Table[(Subscript[X, i])[t_],{i,1,n}];

Xsol=vecX/.sol;

tabv[t_]=Table[{i-1,Xsol[[i]]},{i,1,n}];

Manipulate[ListPlot[{tabv[t]},PlotRange->{{-1.5,L*n+1.5},{-1.5,1.5}}],{t,0,T}]

or analytically with normal modes

(I tried with DSolve, but it just hangs)

The eigenmode (or normal mode) angular frequency of mode i is (Subscript[ω, i])= ω (Sqrt[-(Subscript[eig, i])]), where (Subscript[eig, i]) is the (i^th)eigenvalue of A.

eigentab=Sqrt[Eigenvalues[-A]]//N

The eigenvectors of A give the relative displacements of all the particles for the different modes.  We see that they are arranged from highest to lowest frequency.

eig=Eigenvectors[A]//N

eig[[1]]

It would be nice to see the eigenmodes in action.  We will use the Mathematica routine Manipulate for this purpose.  If the system starts from zero displacement, the contribution to the displacement of particle j from eigenmode i  is given by (Subscript[x, j])= (Subscript[B, i])sin((Subscript[ω, i])t ) eig[[i,j]], where (Subscript[B, i]) is an overall prefactor giving the strength of eigenmode i.  In the following, we have used the Manipulate command to animate eigenmode chosen by a second slider.

Manipulate[ListPlot[Table[{j,0.35 Sin[eigentab[[i]] t]*eig[[i,j]]},{j,1,n}],PlotStyle->PointSize[0.03],PlotRange->{{0,n+1},{-1,1}}],{{t,0,"Time"},0,50} ,{{i,1,"Eigenmode"},1,n,1}]

To find the analytic solution to our problem, we need to figure out how much of each normal mode is present in our initial condition, ie we need to express our initial positions as a liner combination of normal modes.  The fact that the initial velocities are zero simplifies our analysis.

initpos

eig//MatrixForm

Coeff=Table[(Subscript[c, i]),{i,1,n}]

Coeff.eig//MatrixForm

solc=Solve[Coeff.eig==initpos]//First

c=Coeff/.solc

The above numbers tell us how much of each eigenvalue is present

Particle #1’s motion is then (the sum below is over all eignmodes

nid=1;

Sum[c[[i]]eig[[i,nid]]Cos[eigentab[[i]] t],{i,1,n}]

solAnalytic[nid_]:=Sum[c[[i]]eig[[i,nid]]Cos[eigentab[[i]] t],{i,1,n}]

Let’s visualize the whole chain.

solAnalytic[7]

tabA[t_]=Table[{i-1,solAnalytic[i]},{i,1,n}]

Manipulate[ListPlot[{tabA[t],tabv[t]},PlotStyle->{{Blue,PointSize[0.02]},{Red,PointSize[0.01]}},PlotRange->{{-1.5,L*n+1.5},{-1.5,1.5}}],{t,0,T}]

The smaller red dots are the numerical solution.

Let’s look at the motion of say particle #3, both the numerical and analytic solution.

solAnalytic[3]

Plot[{solAnalytic[3],Xsol[[3]]},{t,0,T},PlotStyle->{{Blue,Thickness[0.015]},{Red,Thickness[0.005]}}]

We see that the time series is composed of waves of different frequencies.  For n coupled oscillators, there are only n possible frequencies.

Representing functions as a sum of sines and/or cosines is what Fourier analysis is all about.

A look with Fourier

The Fourier transform tells us what frequencies are present in a curve.  [The Fourier transform is complex, so we will look at its magnitude.]

dt=0.1;
signal=Table[solAnalytic[3],{t,0,T,dt}];

Length[signal]

ListPlot[signal]

Whereas the original signal is in the time domain, the Fourier transform yields a curve in the freqency domain that is symmetric about the middle.

ft=Abs[Fourier[signal]];

f0=(1/(Length[signal] dt))

famp=Table[{f0 (i-1),ft[[i]]},{i,1,Length[ft]}];

ListPlot[famp,PlotRange->All]

ListPlot[famp[[1;;40]],Joined->True,PlotRange->All]

The Fourier Transform peaks tell us that there are 4 frequencies present.  The peak locations should match our normal mode analysis.

eigentab

c//Chop

Because of the symmetry of our initial conditions, only 4 of the 7 modes are active.

Sort[eigentab[[{1,3,5,7}]]/(2 π)]

Fourier analysis

Remove["Global`*"]

Explanatory material taken from Mathematics for Physicists by Susan M. Lea.
Free resource https://en.wikipedia.org/wiki/Fourier_series

The basic idea here is that periodic functions (or functions over a finite interval that can be extended to be periodic) can be represented, perhaps no always perfectly, by sines and cosines.  Discontinuities in the original function can lead to “overshoots” in the Fourier representation.

Fourier Series

Fourier series for a step function

f[x_]:=Piecewise[{{0,x<1/2},{1,x≥1/2}}]

f1=Plot[f[x],{x,0,1},PlotStyle->{Black,Thickness[0.01]}]

Nterms=11;
L=1;
b0 = (1/L)Integrate[f[x],{x,0,L}]
a=Table[(2/L)Integrate[f[x]Sin[2 π n x/L],{x,0,L}],{n,1,Nterms}]
b=Table[(2/L)Integrate[f[x]Cos[2 π n x/L],{x,0,L}],{n,1,Nterms}]

Note that there are no cosine terms, becasue the function, they way we’ve defined it, is an odd function (recalling that in expressing the function as a Fourier series, it is implicitly periodic).

fseries=b0+Sum[a[[n]] Sin[2 π n x/L],{n,1,Nterms}]+Sum[b[[n]] Cos[2 π n x/L],{n,1,Nterms}]

Show[Plot[fseries,{x,-1,2}],f1,PlotRange->All]

We can also use complex exponentials instead of sines and cosines explicitly.

c=Table[(1/L)Integrate[f[x]Exp[-I 2 π n x/L],{x,0,L}],{n,-Nterms,Nterms}]

cseries=Sum[c[[n]]Exp[I 2 π (n-Nterms-1) x/L],{n,1,2Nterms+1}]

Show[Plot[cseries,{x,-1,2}],f1,PlotRange->All]

ExpToTrig[cseries]

fseries

ExpToTrig[cseries]==fseries

Sine and Cosine Series

g[x_]:=x-1.3 (x^2)+(x^3)

g1=Plot[g[x],{x,0,1},PlotStyle->{Black,Bold}]

Full series

Nterms=13;
L=1;
b0 = (1/L)NIntegrate[g[x],{x,0,L}]
a=Table[(2/L)NIntegrate[g[x]Sin[2 π n x/L],{x,0,L}],{n,1,Nterms}]
b=Table[(2/L)NIntegrate[g[x]Cos[2 π n x/L],{x,0,L}],{n,1,Nterms}]

fseries=b0+Sum[a[[n]] Sin[2 π n x/L],{n,1,Nterms}]+Sum[b[[n]] Cos[2 π n x/L],{n,1,Nterms}]

Show[Plot[fseries,{x,-1,2},PlotPoints->100],g1,PlotRange->All]

Sine series

The sine series assumes that we have constructed an odd function extended the original function down to x=-L.  Thus, the sine series representation of the function will have a period of 2L.

a=Table[(2/L)NIntegrate[g[x]Sin[π n x/L],{x,0,L}],{n,1,Nterms}]

sseries=Sum[a[[n]] Sin[π n x/L],{n,1,Nterms}]

Show[Plot[sseries,{x,-2,2}],g1,PlotRange->All]

Cosine series

This is similar to the sine series, except that we are creating an even function when we extend the function to x=-L.

b0 = (1/L)NIntegrate[g[x],{x,0,L}]
b=Table[(2/L)NIntegrate[g[x]Cos[ π n x/L],{x,0,L}],{n,1,Nterms}]

bseries=b0+Sum[b[[n]] Cos[π n x/L],{n,1,Nterms}]

Show[Plot[bseries,{x,-2,2}],g1,PlotRange->All]

The cosine series gives a much better approximation.

Wave equation solution

See Lea

h=1;
L=2;
y0[x_]:=Piecewise[{{2 (h/L) x,x<(L/2)},{2 (h/L)(L-x),x≥(L/2)}}]

Plot[y0[x],{x,0,L}]

Nterms=3

k=n (π/L)

b=Table[(2/L)Integrate[y0[x] Sin[k x],{x,0,L}],{n,1,Nterms}]

v=1;
y[x_,t_]=Sum[b[[n]] Sin[k x]Cos[k v t],{n,1,Nterms}]

Plot[{y[x,0],y[x,0.5],y[x,1]},{x,0,L}]

Manipulate[Plot[y[x,t],{x,0,L},PlotRange->{All,{-1,1}}],{t,0,20}]

Wave equation solution - different initial conditions

See notes for details

h=1;
L=2;
y0[x_]:=Exp[((-(((x-1.2))^2))/0.02)]

Plot[y0[x],{x,0,L},PlotRange->All]

Nterms=17

k=n (π/L)

b=Table[(2/L)NIntegrate[y0[x] Sin[k x],{x,0,L}],{n,1,Nterms}]

v=1;
y[x_,t_]=Sum[b[[n]] Sin[k x]Cos[k v t],{n,1,Nterms}]

Plot[{y[x,0],y[x,0.5],y[x,1]},{x,0,L},PlotRange->All]

Manipulate[Plot[y[x,t],{x,0,L},PlotRange->{All,{-1,1}}],{t,0,4}]

Fun with Sound

Remove["Global`*"]

SetDirectory["~/Movies"]

FileNames[]

ts=Import["AudioTest.wav"]

Length[ts]

ts[[1]]

ts[[1,1]]

ts[[1,1,1]]

amps=ts[[1,1,1]];

ListPlay[amps,SampleRate->48000]

Length[amps]/48000.

ts[[1,2]]

srate=48000

dt=1./srate

tamp=Table[{dt (i-1),amps[[i]]},{i,1,Length[amps]}];

ListPlot[tamp[[70000;;71000]],Joined->True]

Manipulate[ListPlot[tamp[[s;;s+1000]],Joined->True,PlotRange->{All,{-1,1}}],{{s,1},1,Length[tamp]-1000,10}]

fft=Fourier[amps];

fft[[200]]

fft[[L-200+2]]

fftA=Abs[fft];

L=Length[fftA]

ListPlot[fftA,PlotRange->All]

ListPlot[fftA[[1;;Round[L/30]]],PlotRange->All,Joined->True]

srate

dt

f0=(1/(L dt))

5000*f0

famp=Table[{f0 (i-1),fftA[[i]]},{i,1,Length[fftA]}];

ListPlot[famp[[1;;Round[L/30]]],PlotRange->All,Joined->True]

Play[Sin[440 2Pi t],{t,0,1}]

Play[Sin[880 2Pi t],{t,0,1}]

Play[Sin[550 2Pi t],{t,0,1}]

550/f0

Import["https://en.wikipedia.org/wiki/Discrete_Fourier_transform"];

firstN=2000
(*firstN=3*)

ftFilter=Table[0,{i,1,L}];

ftFilter[[1]]=fft[[1]];
ftFilter[[2;;1+firstN]]=fft[[2;;1+firstN]];
ftFilter[[L-firstN;;L]]=fft[[L-firstN;;L]];

rehash=InverseFourier[ftFilter];

r=Re[rehash];

Length[rehash]

ListPlay[amps,SampleRate->48000]

ListPlay[r,SampleRate->48000]