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

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.

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

percsRandomRest:= {
SoundNote["Clap",1],
SoundNote["Sticks",1],
SoundNote["Shaker",1],
SoundNote["LowWoodblock",1],
SoundNote["Castanets",1],
SoundNote[None,RandomReal[{0.8,1.2}]]};
perctable=Table[RandomChoice[percsRandomRest],{50}];
Sound[perctable,20]

percsRandomRange[{a_,b_}]:= {
SoundNote["Clap",RandomReal[{a,b}]],
SoundNote["Sticks",RandomReal[{a,b}]],
SoundNote["Shaker",RandomReal[{a,b}]],
SoundNote[SoundNote["Sticks",1],RandomReal[{a,b}]],
SoundNote["Castanets",RandomReal[{a,b}]],
SoundNote[None,RandomReal[{a,b}]]};
perctable=Table[RandomChoice[percsRandomRange[{0.8,1.1}]],{50}];
Sound[perctable,20]

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

Define functions based on the formulas derived for the fourier amplitudes

Although, Mathematica's functions are generally faster, I use these functions because they are more flexible (i.e., man page for FourierTrigSeries is cumbersome to generalize to arbitrary periodic domains)

The constant term:

EvenTerms[0,function_ , λ_, center_] := (1/λ)(∫_(center - λ/2)^(center +λ/2))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_ , λ_, center_] := EvenTerms[SP,function , λ, center]=(2/λ)(∫_(center - λ/2)^(center +λ/2))function[dum] Cos[((2 SP π dum)/λ)]\[DifferentialD]dum

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

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

OddTerms[SP_Integer,function_ , λ_, center_] := OddTerms[SP,function , λ, center]=(2/λ)(∫_(center-λ/2)^(center+λ/2))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_, λ_, center_] :=
Table[OddTerms[i, function, λ, center], {i,0,NTerms}]

EvenAmplitudeVector[NTerms_Integer, function_, λ_, center_] :=
Table[EvenTerms[i, function, λ,center], {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,1,1]

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

OddAmplitudeVector[6, myfunction,1,1]

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_ , λ_, center_]:= EvenAmplitudeVector[n,function,λ, center].EvenBasisVector[n,var,λ] + 
OddAmplitudeVector[n,function,λ, center].OddBasisVector[n,var,λ]

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

temp3 = FourierTruncSeries[3,myfunction,x,1,1]

temp6 = FourierTruncSeries[6,myfunction,x,1,1]

temp12 = FourierTruncSeries[12,myfunction,x,1,1]

Plot[{temp3,temp6,temp12,myfunction[x]},{x,-1,3},PlotRange-> {-.1,.25},PlotStyle-> {Red, Green, Blue,Black},BaseStyle->{ Large,Thick},ImageSize-> Large]

alt2 =FourierTruncSeries[2,myfunction,x,2,1]

alt6 =FourierTruncSeries[6,myfunction,x,2,1]

Plot[{alt2,alt6,myfunction[x]},{x,-1,3},PlotRange-> {-.1,.3},
PlotStyle->{{Thick,Red},{Thick,Blue},{Thick,Darker[Green]}}]

FPlot[n_Integer] :=FPlot[n]=Plot[Evaluate[FourierTruncSeries[n,myfunction,x,2 π,0]],{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.

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

Afunction[x_] := x^3

Clear[afuncappx]

afuncappx[i_] := afuncappx[i] = FourierTruncSeries[i,Afunction,x,4,0]

ExactPlot = Plot[Afunction[x],{x,-2.5,2.5},PlotStyle-> {Thick,Black}]

Simplify[afuncappx[11]]

afuncappx[16]

Plot[Evaluate[afuncappx[17]],{x,-2.5,2.5}]

Manipulate[GraphicsRow[{plt=Show[Plot[Evaluate[afuncappx[truncation]],{x,-2.5,2.5},PlotRange->All,PlotStyle->{Thick,Hue[truncation/50,1,0.25]}],ExactPlot],Show[plt,PlotRange->{{1.5,2.2},{4,8.5}}]},ImageSize->Full,BaseStyle-> Large](*;theappx[[truncation]]*),{{truncation,10},1,40,1},
Initialization-> (theappx = Table[afuncappx[i],{i,1,20}])]