Periodic Phenomena and Fourier Series

Playing with Familiar Periodic Phenomena: Sound

Sounds will not be available on PDF or HTML versions
Let's begin by "looking" at a familiar periodic phenomena:

We index the notes and write an indexed set of frequency (in Hertz) for each of the notes for one octave above middle-c. We write a function to create a Sound for each note.

c =1; d =2; e=3; f=4; g=5; a=6; b=7;
freq[c] =261.6;
freq[d]  =293.7;
freq[e] =329.6;
freq[f] = 349.2;
freq[g] =392.0;
freq[a] = 440.0;
freq[b] = 493.9;
purenote[note_Integer] := purenote[note]=Play[Sin[2 π freq[note]t],{t,0,1}]

We extend the function to get simultaneous notes from a List. We use Thread which takes f[{a,b,c}] to {f[a],f[b],f[c]}

notes[note_List]:= Sound[Thread[purenote[note]]]

Here are examples of their use.

cnote=purenote[c]
notes[{a,c,e}]

We can play with variable amplitudes for a fixed frequency, here we can hear the increased, but non-singular amplitude through zero.

Plot[Sin[540x]/x,{x,-.1,.1},PlotRange->All,Filling-> Axis]
Play[Sin[540x]/x,{x,-1,1}]

We can vary amplitudes and frequencies. Warning, playing with this function can become addictive...

Play[2 Sin[20x Sin[x Exp[-x/.2]  + Sin[x]/x]] +  Exp[(1+Cos[x])] Sin[x Exp[x/10]] Sin[1500 x],{x,0.01,20}]

Playing with Familiar Periodic Phenomena: Music and Instruments

Let's see if we can play this:    l

twoframes = {e,e,f,g,g,f,e,d,c,c,d,e};

We will play it, but it probably not what was intended...

notes[twoframes]

We create a rest of  a fixed length, and then use Riffle to insert a rest after each note. We use a new function, called SoundNote, with None as the first argument, we get no sound.

purenote[rest] =SoundNote[None,.15];

notes[Riffle[twoframes,rest]]

SoundNote can take general strings for arguments as well, here we enter the musical notes from above (as characters), but one octave lower..

TwoFramesLower = {"E3","E3","F3","G3","G3","F3","E3","D3","C3","C3","D3","E3"};

The default is to use a piano to make the sound; here we ask for a duration of 6/10 of a second.

piano[note_String]:= SoundNote[note,.6]

Sound[Thread[piano[TwoFramesLower]]]

We can use other MIDI instruments as well; here is a bagpipe

bagpipe[note_String]:= SoundNote[note,.6,"Bagpipe"]
Sound[Thread[bagpipe[TwoFramesLower]]]

And, now for the birds.

avian[note_String]:= SoundNote[note,.2,"Bird"]
avian[rest] = SoundNote[None,.4];
Sound[Thread[avian[Riffle[TwoFramesLower,rest]]]]

Now, a little Bach meets Beethoven

Sound[Riffle[Thread[SoundNote[#,.2]&[TwoFramesLower]],{Thread[SoundNote[#,.1]&[TwoFramesLower]]}]]

Roll over Beethoven: Random Notes and Instruments

Let's hear what random notes sound like: SoundNote[n] will play n semitones above middle C.  Here we make a list of random notes and play them.

RandomNotes = Table[SoundNote[RandomInteger[{-15,20}],.2],{36}];
Sound[RandomNotes,10]

Here, we mix in some rests at random

RanRest[]:= Module[{rdur=.2},If[RandomReal[] >0.5,rdur = .4];SoundNote[None,rdur]]
RandomNotesandRests = Table[If[RandomReal[]≥ .33,SoundNote[RandomInteger[{-15,20}],.2],RanRest[]],{96}];
Sound[RandomNotesandRests,20]

Now, we ask for random instruments as well, with "chords" of up to 5 instruments at each beat.

RandomInstruments = Table[If[RandomReal[]>0.5,
Table[
SoundNote[
Table[RandomInteger[{-15,20}],{RandomInteger[{1,5}]}],Round[RandomReal[{1,2}],.2],RandomInteger[{1,15}]],
{RandomInteger[{1,4}]}],RanRest[]],{48}];
Sound[RandomInstruments,20]

Finally, some random percussion events.

percs= {SoundNote["Clap",1],SoundNote["Sticks",1],SoundNote["Shaker",1],SoundNote["LowWoodblock",1],SoundNote["Castanets",1],SoundNote[None,1]};
perctable=Table[RandomChoice[percs],{50}];
Sound[perctable,20]

Periodic Functions of wavelength λ using the modulus function:

Mod is a very useful function that can be used to force objects to be periodic.  Mod[x,λ] return that part of x that lies within 0 and λ.  Or, in other words if we map the real line x to a circle with circumference λ, then Mod[x,λ] returns were x is mapped onto the circle.

modmatdemo[n_Integer,λ_Integer] := Table[{i,Mod[i,λ]},{i,1,n}]//MatrixForm;
modcircledemo[n_Integer,λ_Integer]:= Module[{xpos,angle,cpos},
Graphics[
Table[xpos= 3 Quotient[i-1,λ]; angle=2 Pi Mod[i,λ]/λ;cpos={Cos[angle],Sin[angle]};
{Circle[{xpos,0}],
Text[i,Flatten[{{xpos,0} +1.2*cpos}]],
Text[Mod[i,λ],Flatten[{{xpos,0} +0.8*cpos}]]},{i,1,n}]]];

GraphicsColumn[{modmatdemo[13,5],modcircledemo[26,5]},ImageSize->Full]

Boomerang uses Mod to force a function, f, with a single argument, x, to be periodic with wavelength λ

Boomerang[f_ , x_ , λ_ ] := f[Mod[x,λ]]

AFunction[x_ ] := ((3-x)^3)/27

The following step uses Boomerang to produce a periodic repetition of AFunction over the range 0 < x < 6:

Plot[Boomerang[AFunction,x,6],{x,-12,12},PlotRange->All]

Demonstration of the orthogonality of the trigonometric functions

Integrating (∫_xo^(xo+λ))cos(2mπx/λ) cos(2nπx/λ)dx

fassume= {Minteger ϵ Integers, Ninteger ∈ Integers,xo ∈ Reals, λ > 0}
coscos = Integrate[Cos[((2 π Minteger x)/λ)] Cos[((2 π Ninteger x)/λ)] ,{x,xo,xo + λ},Assumptions->fassume]

Demonstrating that the integral goes to zero for m≠n and λ/2 if m=n (we demonstrate this by using random choices for m and n)

Table[{{{mrand,nrand}=RandomInteger[{1,50},2]},Simplify[coscos/.{Minteger->mrand , Ninteger->nrand}]},{20}]

If we run the above command enough times, we will see  that results where n=m give indeterminate values, for these we should use a Limit.

Limit[coscos,Minteger->Ninteger,Assumptions->fassume]

Performing the same steps for  (∫_xo^(xo+λ))cos(2mπx/λ) sin(2nπx/λ)dx and  (∫_xo^(xo+λ))sin(2mπx/λ) sin(2nπx/λ)dx

cossin = Integrate[Cos[((2 π Minteger x)/λ)] Sin[((2 π Ninteger x)/λ)] ,{x,xo,xo + λ}, Assumptions->fassume]

Table[{{{mrand,nrand}=RandomInteger[{1,50},2]},Simplify[cossin/.{Minteger->mrand , Ninteger->nrand}]},{20}]

Limit[cossin,Minteger->Ninteger,Assumptions-> fassume]

sinsin = Integrate[Sin[((2 π Minteger x)/λ)] Sin[((2 π Ninteger x)/λ)] ,{x,xo,xo + λ}, Assumptions->fassume]

Table[{{{mrand,nrand}=RandomInteger[{1,50},2]},Simplify[sinsin/.{Minteger->mrand , Ninteger->nrand}]},{20}]

Limit[sinsin,Minteger->Ninteger,Assumptions-> fassume]

Fourier Series: Using The Derived Formulas for Determination of Fourier Foefficients

First we will "do it the hard way" and write short programs that evaluate Fourier coefficients; then we will demonstrate how to make use of built-in functions in Mathematica's FourierTransform package…

Define functions based on the formulas derived for the fourier amplitudes

The constant term:

EvenTerms[0,function_ , λ_] := (1/λ)(∫_0^λ)function[dum]\[DifferentialD]dum

A function that defines each even amplitude individually (this is not very efficient, it would be better to evaluate the integral once and use that result)

EvenTerms[SP_Integer,function_ , λ_] := EvenTerms[SP,function , wavelength]=(2/λ)(∫_0^λ)function[dum] Cos[((2 SP π dum)/λ)]\[DifferentialD]dum

Define the zeroth odd term as zero for symmetry with the even terms:

OddTerms[0,function_ , λ_] := 0

OddTerms[SP_Integer,function_ , λ_] := OddTerms[SP,function , λ]=(2/λ)(∫_0^λ)function[dum] Sin[((2 SP π dum)/λ)]\[DifferentialD]dum

A function to create a vector of amplitudes for the odd terms and one for the even terms

OddAmplitudeVector[NTerms_Integer, function_, λ_] :=
Table[OddTerms[i, function, λ], {i,0,NTerms}]

EvenAmplitudeVector[NTerms_Integer, function_, λ_] :=
Table[EvenTerms[i, function, λ], {i,0,NTerms}]

Fourier Series: Example using Derived Functions

Define a function to use as an example:

myfunction[x_ ] :=(x*(2-x)*(1-x)^2)

OriginalPlot= Plot[myfunction[x],{x,0,2},PlotStyle->{Hue[.66],Thickness[0.015]}]

Now calculate Fourier amplitudes for this function…

EvenAmplitudeVector[6, myfunction,2]

The odd terms should all vanish because the function is even around the center x=(λ/2)

OddAmplitudeVector[6, myfunction,2]

OddBasisVector[NTerms_Integer,var_,  λ_] := Table[Sin[((2 π  i var)/λ)], {i,0,NTerms}]

OddBasisVector[6,x,2]

EvenBasisVector[NTerms_Integer,var_,  λ_] := Table[Cos[((2 π  i var)/λ)], {i,0,NTerms}]

EvenBasisVector[6,x,2]

Make the Fourier series representation by taking the dot product of the AmplitudeVectors with the BasisVectors:

FourierTruncSeries[n_, function_,var_ , λ_]:= EvenAmplitudeVector[n,function,λ].EvenBasisVector[n,var,λ] + 
OddAmplitudeVector[n,function,λ].OddBasisVector[n,var,λ]

For myfunction we write the Fourier series out to six terms as:

FourierTruncSeries[6,myfunction,x,2]

FPlot[n_Integer] :=FPlot[n]=Plot[Evaluate[FourierTruncSeries[n,myfunction,x,2]],{x,-2,4},PlotStyle->{Thick,ColorData[n,"ColorList"]}]

Note the periodicity of the plot above.  The wavelength is 2. Only three terms were calculated in the Fourier series representation of myfunction. Below the Fourier series and the exact function are shown superposed.

Now the calculation is repeated but retaining 6 terms in the Fourier series representation and hopefully obtaining a more accurate representation of myfunction.

Show[OriginalPlot,FPlot[3],FPlot[6],PlotRange->{{-0.5,2.5},{-0.1,0.26}}]

Fourier Series Coefficients using Mathematica's FourierSeries package: Preparing the Input Function

Needs["FourierSeries`"]

AFunction[x_] := (((x-3)^3)/27)

Plot[AFunction[x],{x,0,6}]

Mathematica's Fourier Series functions are defined for function that are periodic in the domain x ∈ (-1/2,1/2).  So we need to map the periodic functions to this domain

ReduceHalfHalf[f_ , x_ , λ_ ] := f[(x+1/2)*λ ]

ReducedFunction = ReduceHalfHalf[AFunction,x,6]//Simplify

ExactPlot=Plot[ReducedFunction,{x,-1/2,1/2},PlotRange->All,PlotStyle->{Red,Opacity[0.5],Thickness[0.01]}]

FourierCosCoefficient[ReducedFunction,x,n]

FourierSinCoefficient[ReducedFunction,x,n]

FourierTrigSeries[ReducedFunction,x,5]

Efficient Calculation of Truncated Fourier Series

ManipulateTruncatedFourierSeries[function_, {truncationstart_, truncationend_, truncjump_}] := Manipulate[Plot[Evaluate[FourierTrigSeries[function,x,itrunc]],{x,-1,1},PlotRange-> {-2,2}],{itrunc,{truncationstart,truncationend, truncjump}}];

The function above will work, but it is horribly inefficient!  Because it asks FourierTrigSeries to calculate one more term each time, it is doing some redundant work.  We can fix this up by having
 it calculate one new term and adding to the sum calculated previously. Here it is:

costerm[function_,x_, n_Integer] := Simplify[FourierCosCoefficient[function,x,n]] Cos[2 π n x]
sinterm[function_,x_, n_Integer] :=  Simplify[FourierSinCoefficient[function,x,n]] Sin[2 π n x]

TruncatedFourierSeries[function_,x_, 0] := TruncatedFourierSeries[function,x, 0]=
costerm[function,x,0] + sinterm[function,x,0]

TruncatedFourierSeries[function_, x_,n_Integer] := TruncatedFourierSeries[function, x,n]=TruncatedFourierSeries[function, x,n-1] + costerm[function,x,n] +sinterm[function,x,n]

Convergence of Fourier Series: Gibbs Phenomenon

The following will demonstrate how convergence is difficult where the function changes rapidly---this is known as Gibbs' Phenomenon

Manipulate[GraphicsRow[{plt=Show[Plot[theapprx[truncation],{x,-0.4999,0.4999},PlotRange->{{-0.55,0.55},{-1.4,1.4}},PlotStyle->{Thick,ColorData[1,truncation]}],ExactPlot],Show[plt,PlotRange->{{0.4,0.5},{0.5,1.2}}]},ImageSize->Full],{{truncation,100},1,100,1},Initialization ->(Table[ theapprx[i]=TruncatedFourierSeries[ReducedFunction,x,i],{i,1,100}];)]