Fourier Analysis

Eric R. Bittner

•Fourier Series

You should have already met Fourier series (see e.g., http://mathworld.wolfram.com/FourierSeries.html) on a number of occasions. Fourier series are expansions of any periodic function, f(x), in terms of an infinite of the form

f(x) = a _ 0/2 + Underoverscript[∑, n = 1, arg3] (a _ n cos(n x) + b _ n sin(n x)) .

• Orthogonality

In this section the the basic orthogonality integrals will be examined.

•∫ _ 0^(2 π) sin(n x) sin(m x) d x

In this section you will prove that

∫ _ 0^(2 π) sin(n x) sin(m x) d x = π δ _ (n, m)

Compute ∫ _ 0^(2 π) sin(n x) sin(m x) d x. Note that Mathematica does not assume that n or m are integers.   

n != m: Use Simplify with appropriate conditions to show that this integral vanishes.

n = m: In this case take the Underscript[lim, n -> m] (you may have to look up Limit) and Simplify the result.

You should now have the desired result.

•∫ _ 0^(2 π) sin(n x) cos(m x) d x

Prove that ∫ _ 0^(2 π) sin(n x) cos(m x) d x = 0.

•∫ _ 0^(2 π) cos(n x) cos(m x) d x

Compute ∫ _ 0^(2 π) cos(n x) cos(m x) d x for general n and m.

•1/L ∫ _ (-L)^L cos((m π x)/L) cos((n π x)/L) d x

Compute 1/L ∫ _ (-L)^L cos((m π x)/L) cos((n π x)/L) d x for general n and m.

•4.1.2 Square wave

The Fourier series of a function assumes that the function is periodic outside its specified domain. Consider the function defined by

f(x) = { 1, 0 <= x <= π , ó© -1, π <= x <= 2 π

with the additional (periodicity) condition

f(x) = f(x mod (2 π)),

which causes the function to repeat itself outside the domain [0, 2 π].

Complete the following

f(x_ /; 0 <= x <= π) := ...

f(x_ /; ...) := ...

f(x_) := f(Mod[..., ...])

to define the square wave of equations (4.1.0) and (4.1.0) above.

To verify your definition visualize the function.  

•4.1.3 Fourier coefficients

From the orthogonality integrals computed above it can be shown that the Fourier coefficients are:

a _ 0 = 1/π ∫ _ 0^(2 π) f(x) d x, a _ n = 1/π ∫ _ 0^(2 π) f(x) cos(n x) d x,    n >= 1, <br /> b _ n = 1/π ∫ _ 0^(2 π) f(x) sin(n x) d x,    n >= 1.

For the square wave f(x) = { 1, 0 <= x <= π ó© -1, π <= x <= 2 π, compute a _ 0 by splitting the integration into two parts and compute the exact integral for each domain. (Hint: the first integral to be computed is ∫ _ 0^π 1 d x). Save the result as a _ 0 = ···. Is your answer reasonable? (What value do you expect for the DC-term or average of a square wave)?

Similarly compute a _ n and simplify your result using the methods of the orthogonality section above and save the result as a _ n_ = ···. Why is your answer expected?

Compute and simplify b _ n. Save the result as b _ n_ = ···.

•4.1.4 Partial sum of the Fourier series to m terms

Using a _ 0, a _ n, and b _ n computed in §4.1.3 the sum to m terms found in the previous section can be written as

FourierSeries[m_] := a _ 0/2 + Underoverscript[∑, n = 1, arg3] (a _ n cos(n x) + b _ n sin(n x))

Produce a Table of Plots of FourierSeries[m] to see the partial sums for m = 1, 2, ..., 10. Set PlotRange -> {-1.5, 1.5} in each of the plots. Select all the graphical output of the previous section and use Animate Selected Graphics (under the Cell menu) on your result to visualize the partial sums.

It is a good idea to wrap FourierSeries[m] with Evaluate so that the sum is not re-computed for each value of x.

•4.2 DFT and FFT

In this section you will be introduced to the basics of the discrete Fourier transform (DFT) and fast Fourier transform (FFT).

Use of the FFT is important in practice because experimentally you can only measure a quantity finitely often at a discrete set of times. Hence you always have a discrete dataset which you can analyse using the FFT. It is important to understand the relationship between the results returned by the FFT and the amplitudes and frequencies present in the data you are analysing.

As an aside, there are a number of other transform functions relevant including the discrete cosine transform (DCT), used in JPEG image compression, the Hartley transform, and the discrete wavelet transform (DWT). Wavelets are finitely restricted waves (wave-packets) which permit orthogonal expansion of functions which are both spatially and temporally restricted. Recently, wavelets have been adopted by the FBI as their image compression standard for fingerprints.

•4.2.1 DFT and FFT

Read the discussion on fast transform algorithms to see why the FFT is faster than the DFT.

•4.2.2 One-dimensional Example

•Analysis of a Periodic Signal

Evaluate the linear superposition of 3 harmomic waves, sin((π t)/128), sin((3 π t)/128), cos((5 π t)/128), for t = 1, 2, ..., 256 and call the result öÄ.

Visualise öÄ using ListPlot. Join up the data values using the PlotJoined option and use PlotStyle and RGBColor to show signal in red. Call your plot ÷‚.

Compute the DFT (using Fourier) of öÄ and call the result ö·. Display ListPlots of the real (Re) and imaginary (Im) values of ö· over a suitably chosen PlotRange. (Hint: it might be useful to divide ö· by Length[ö·]^(1/2)).  Interpret this (graphical) result in light of the data. What do the largest values in each graph correspond to?  Where is the “DC-term” (i.e., the average) located and what is its value?

Consider the following superposition of 3 signals:
  s(t_) := a e^(- α t) sin((π t)/128) + b e^(-β t ) sin(3 π t/128) + c cos(5 π/128 t)
  where {a, b, c} = {2, 5, 4}and {α, β, γ} = {π/128, 4 π/128, 7 π/128}.
Plot the signal as a function of time.
Now, discretize the signal and add a small random number to each component to simulate a noisy signal
  noisy(t_) := a e^(- α t)(sin((π t)/128) + Random[]) + b e^(-β t ) ( sin(3 π t/128) + Random[]) + c    c e^(-γ t) ( cos(5 π/128 t) + Random[])

s[t_] := {a e^(- α t) sin((π t)/128) , b e^(-β t ) sin(3 π t/128) , c e^(-γ t) cos(5 π/128 t)} /. {a -> 2, b -> 5, c -> 4, α -> 0.1 π/128, β -> .4 π/128, γ -> .7 π/128}

Plus @@ s[t]

4 e^(-0.01718058482431918` t) Cos[(5 π t)/128] + 2 e^(-0.002454369260617026` t) Sin[(π t)/128] + 5 e^(-0.009817477042468103` t) Sin[(3 π t)/128]

Plot[Plus @@ s[t], {t, 0, 2^11}, PlotRange -> All] ;

[Graphics:HTMLFiles/index_60.gif]

data = Table[4 e^(- t/200) (Cos[π t/50] + δ Random[]) + 2 e^(- t/200) (Sin[   π t/30] + δ Random[]) + 5 e^(- t/1000) (Sin[ π    t/100] + δ Random[]) /. {δ -> 0}, {t, 0, 10000, 10}] ; ListPlot[data, PlotRange -> All, PlotJoined -> True] ; ft = Fourier[data] ;

[Graphics:HTMLFiles/index_62.gif]

ListPlot[Take[Abs[ft], 200], PlotRange -> All]

[Graphics:HTMLFiles/index_64.gif]

-Graphics -

2^11 70/2^8 // N

560.`

•Now for somefun!!

Mathematica can be used to do image processing of a complex data set.  For example, here is a picture of the sun.

g = Import["/Library/WebServer/Documents/Chem3111/greensun.gif"]

Show[g]

[Graphics:HTMLFiles/index_70.gif]

Here's the structure of the picture.

Shallow[InputForm[g]]

Graphics[Raster[<< 4 >>], Rule[<< 2 >>], Rule[<< 2 >>], Rule[<< 2 >>]]

Basically, the first element is a raster array, So, we can pull that out

d = g[[1, 1]] ;

Dimensions[d]

{512, 512}

Now we can do a false color plot!

ListDensityPlot[d, ColorFunction -> ( Hue[1 - #/2] &)] ;

[Graphics:HTMLFiles/index_77.gif]

ListDensityPlot[d, ColorFunction -> (CMYKColor[.1, 0.1, 1 - #, 1 - #] &)] ;

[Graphics:HTMLFiles/index_79.gif]

ListDensityPlot[d, ColorFunction -> (RGBColor[0, #, 1] &)] ;

[Graphics:HTMLFiles/index_81.gif]

ListContourPlot[d, ContourShading -> Hue]

[Graphics:HTMLFiles/index_83.gif]

-ContourGraphics -

Here's a plot of the intensity variation:

ListPlot[Sort[Flatten[d]]]

[Graphics:HTMLFiles/index_86.gif]

-Graphics -

Show[Graphics[Raster[d^2 / Max[d^2]]],
    AspectRatio->Automatic]

[Graphics:HTMLFiles/index_88.gif]

-Graphics -

•4.4 Sound

And now for some fun ...  Mathematica has the ability not only to Plot functions but to Play them as well.

•4.4.1 Simple Sounds

For those that are curious, the lowest note on the piano is A sounding at 27.50 Hz and the highest is the top C sounding at 4186.009 Hz. A440 concert pitch refers to tuning the A above middle C to 440.000 Hz. See http://www.pianoworld.com/pitch.htm.

Look up Play. Play a note of frequency 261.626 Hz (i.e., middle-C) for 1 second. Use Play to generate a (linear) superposition of two tones, middle-C and A440, and play this for 1 second.

Consider

Play [ sin ( 2 π 256 t^2 ) , {t, 0, 2} ] ;

[Graphics:HTMLFiles/index_91.gif]

Cell[]At what rate is the frequency of this sound increasing (Hint: think about derivatives)?

To make question 0 clearer, consider the following plot:

f(t_) = 2 t^2 - t + 1 ;

Plot[sin(f(t)), {t, 0, 5}, PlotStyle -> RGBColor[1, 0, 0]] ;

[Graphics:HTMLFiles/index_95.gif]

What is the instantaneous frequency of sin(f(t))? One way of answering this is to consider the periodic function, sin(f^'(t _ 0) t - (f^'(t _ 0) t _ 0 - f(t _ 0))) ≡ sin(ω _ t _ 0 t - φ _ t _ 0), that has the same frequency as sin(f(t)) at t = t _ 0. For example, with t _ 0 = 3:

Plot[Evaluate[sin(f^'(3) t - (f^'(3) 3 - f(3)))], {t, 0, 5}, PlotStyle -> RGBColor[0, 0, 1]] ;

[Graphics:HTMLFiles/index_102.gif]

Show[%, %%] ;

[Graphics:HTMLFiles/index_104.gif]

From this it can be seen that the (angular) frequency of sin(f(t)) at t = t _ 0 is ω = 2 πν = f^'(t _ 0) and the phase is φ = f(t _ 0) - f^'(t _ 0) t _ 0. Hence the instantaneous frequency of sin(2 π f t^2) is

∂ (2 π f t^2)/∂ t/(2 π)

2 f t

and the rate of increase of frequency is 2 f t   Hz per second.  See how many times you can play this before the lab TA gets irrate.  

•4.4.2 Touch-tone Phone

It turns out that all the 12 tones on a touch-tone phone (i.e., 1 - 9, 0, * and #) are generated by forming the superposition of two tones of differing frequencies, played for 0.1 seconds.

•Forming pairs

Use Outer to generate a list of all the possible pairs of frequencies from the lists {697, 770, 852, 941} and {1209, 1336, 1477}.  You should get a total of 4 × 3 = 12 pairs. You will need to appropriately Flatten this result to get rid of one level of brackets. Call the result tones.

Write a simple procedure called PlayTouchTone[{f1_, f2_}] which takes two frequencies and play a superposition of them for 0.1 seconds. (Hint:You have already done this above for middle C and A440 in question 0.

•The numbers 0-9

Play each pair of frequencies generated in question 0. Go to a touch-tone phone and check which sound corresponds to each member of your array. Order your array so that the “numerical” sounds are in ascending order, i.e., from 0 to 9.

Define TouchTone by entering

TouchTone(n_Integer /; 0 <= n <= 9) := PlayTouchTone(tones [[ n + 1 ]])

Note that this code can only handle single-digits.

•Multiple digits

Generalize to multiple digits by making TouchTone have the Attribute Listable and use IntegerDigits to turn a single integer into a list of digits.  Check your definition of TouchTone using ??.

Dial the Chemistry Office using TouchTone[32701] ;  You can animate this sequence of graphics/sounds using Animate Selected Graphics under the Cell menu.

Converted by Mathematica  (August 29, 2002)