Fourier Notes

On discrete and continuous Fourier analysis, or how to add up waves to get anything you like.

Jim Mahoney, Marlboro College, Fall 1999

Discrete Fourier Series

Let (Subscript[x, n])be coordinates on a uniform grid with spacing Δx over a length L = N Δx,  e.g.

(Subscript[x, n]) =  n Δx ,  
    where   0 < n ≤ N  and
Δx  =  
      (L/(N.))

Likewise, let (Subscript[y, n])be a discrete set of points defined on the grid,  (Subscript[Y, n]) = f((Subscript[x, n])), like this.

We would like to express this function f as a sum of waves -  sin and cosine functions, or complex exponentials.

The longest wavelength λ we will allow is L, the length of the x interval, since that will always give us solutions which are periodic in L.  Essentially that means we are only considering functions such that f(x + L) = f(x).  (For some problems, such as a rope which is fixed at both ends, one chooses instead (λ/2) = L; however, that is not what we are doing here.)  
	  		Longest wavelength.

The shortest wavelength is 2  Δx, which corresponds to an oscillation up and down as fast as possible on the grid, like this.
	 ...	Shortest wavelength.

Choosing the other possible wavelengths to match the other available grid spacings (3  Δx, 4  Δx, ...), we find (N/2)discrete wavelengths.  Since we also have two independent phases (either sin and cosine or positive and negative complex exponentials), the total number of discrete waves is just  N, the same as the number of independent y values.  More on this in a moment.

It is convenient to label each wave by its wavenumber k = ((2  π)/λ)with a spacing Δk = ((2  π)/L).  Then the possible values of k are

(Subscript[k, j]) =  
      ((2  π)/(Subscript[λ, j]))  = 
        j Δk  = 
          j (( 2  π)/L) = 
            j (( 2  π)/(N Δx)) , 
    where  - (N/2) <  j  ≤  (N/2)  andΔk = ((2  π)/L) .

Let's look at the number of (Subscript[x, n])and (Subscript[k, j])values a bit more closely and count the number of independent quantities.  In the first place, since f(x) is periodic, (Subscript[y, 0])is the same as (Subscript[y, N]); therefore, we only need n = 0, 1, 2,  ... , N - 1.  Second, the highest possible frequency is \(TraditionalForm
      \`k\_max = \ \(π\/Δx = Nπ\/L\)\)has only one phase, not two.  That is because it is not possible to shift sideways by a half cycle at that frequency; the grid spacing is not small enough.  Thus j = (N/2)and j = -(N/2)refer to the same wave.  This would appear to imply that there are only N-1 different values of (Subscript[k, j]).  However, there is one more value of k which we have not yet mentioned : k = 0, which allows us to include a "zero frequency" or "DC offset" term.

Furthermore, since x is periodic, we could have instead chosen x symmetric about x = 0 with

\(TextForm\`\(x\_n = \ n\ Δx\ , \ 
    where\ \ \  - N/2\  < n ≤ \ N/2\ \ , 
\)\)

which is the way we'll usually write these limits.  But equation (1) is just as correct, since (Subscript[x, N]) = (Subscript[x, 0])and (Subscript[x, -N/2]) = (Subscript[x, N/2]). 

With those details out of the way, we're now ready to define our discrete Fourier transform :

Please note that different authors put the factors of N and 2  π in different places; check your conventions before plugging into someone else's equations.  I'm choosing conventions which (a) make the transition to the integral forms clear, and (b) are commonly seen among physicists.  (Note that engineers tend to use the symbal "j" for (Sqrt[-1]), rather than "i", and also tend to use f (frequency in units of cycles/sec) rather than the physicist's k or ω (radians/meter or radians/sec).

\(Y\_n = 
      \(∑\+\(j = \(-N\)/2 + 1\)\%\(N/2\)\ 
          \(C\_j\) \(e\^\(i\ \(k\_j\) x\_n\)\) 
            Δk\/\(2  π\ \ \)
C\_j = \ \ \ \ 
        ∑\+\(n = \(-N\)/2 + 1\)\%\(N/2\)\ 
          Y\_n\ \(e\^\(\(-i\)\ \(k\_j\) x\_n\)\) Δx\)\)

The symmetry of these two equations is quite striking.  Essentially what we have is a discrete description of the function f in terms of N coefficients, which we can choose as either the numbers (Subscript[Y, n])or (Subscript[C, j]), which specify f in terms of x coordinate values or sinusoidal waves.  

By inserting our definitions of (Subscript[x, n]), Δx,  (Subscript[k, j]), and Δk, we can also write this as

\(Y\_n = \ 
      \(\(1\/L\) 
          \(∑\+\(j = \(-N\)/2 + 1\)\%\(N/2\)\ 
            \(C\_j\) e\^\(2  Π\ i\ \ \((j\ n)\)/N\)\ 
C\_j\) = \ 
        L\/N\ \(\ \ \ ∑\)\+\(n = \(-N\)/2 + 1\)\%\(N/2\)\ Y\_n\ 
          \(e\^\(\(-2\) Π\ i\ \ \((j\ n)\)/\ N\)\ \ \  . \)\)\)

If we then also re-define the Fourier coefficients to include the length L,

(Subscript[G, j]) =  (Subscript[C, j])/L

we get this even simpler form,

\(\(\(Y\_n = \ \(∑\+\(j = \(-N\)/2 + 1\)\%\(N/2\)\ \(G\_j\) 
          e\^\(2  Π\ i\ \ \((j\ n)\)/N\)\ 

      G\_j = \ \ \ \ 1\/N\ \ \ \(\(\ \ \)\(∑\)\)\+\(n = 1\)\%N\ Y\_n\ e\^\(\(-2\) Π\ i\ \ \((j\ n)\)/N\)\)\)\(,\)\)\)

which is the kind of expression that actually gets coded into computer programs.  It turns out that there's a partiuclarly fast algorithm for calculating this when N is a power of 2, called the "Fast Fourier Transform".

Are these two really inverses of each other?  To check, substitute one inside the other - without forgetting that the sum is over a dummy variable.  Replacing, for example, j with j' in the first expression, and substituting into the second gives

\(G\_j\( = \&?\)\ \ \ \ 1\/N\ \ \ \(\(\ \ \)\(∑\)\)\+\(n = 1\)\%N\ \ \((\ ∑\+\(j' = \(-N\)/2 + 1\)\%\(N/2\)\ \(G\_\(j'\)\) 
              e\^\(2  Π\ i\ \ \((j'\ n)\)/N\))\)\ e\^\(\(-2\) Π\ i\ \ \((j\ n)\)/N\)\)

which should be an identify if all this is really working.  Well, of course it is.  To see this, swap the order of the summation signs, putting everything that depends on n on the inside.  Then

\(G\_j\( = \&?\)\ \ 
      1\/N\ \ \ \ 
        \(∑\+\(j' = \(-N\)/2 + 1\)\%\(N/2\)\ 
          G\_\(j'\)\ \ 
            \((\ \ \(\(\ \ ∑\)\+\(n = 1\)\%N\) 
                e\^\(2  Π\ i\ \ \((j' - j)\)\ n/N\))\)\)\)

The part in the parenthesis is 0 if j ≠ j', since in that case we have the sum of complex numbers spaced equally around the unit circle, which all cancel out.  On the other hand, if j = j', then it's just N, since we're adding N copies of (e^0)which is 1.   So our result checks.

(Subscript[G, j])\( = \&?\)  (1/N)    (∑_(j' = -N/2 + 1)^(N/2)) (Subscript[G, j'])  N (Subscript[δ, j - j'])= (Subscript[G, j])   .

An Example

So now that we have the discrete formulas, let's do an example.  A typical one is the "hat" function,

f(x) = { ({{B   for  - a  ≤  x  ≤  a}, {0   otherwise}})  :

To discretize this, let's set set values like shown in this picture, with B = 1,  Δx = 0.1,  N = 60,  
      a = 1.0.  Applying the formulas above gives the following Fourier Transform:

Here the purple shows the original function (Subscript[Y, n]) while the red and orange show the real and imaginary parts of the Fourier components (Subscript[C, j]), that is, the amplitude of the various sinusoids which when added up give the purple curve.  Cool, eh?  As the purple Y gets wider, the red C gets narrower, and vice-versa.  

You can in fact do this sum analytically - but I'll leave that up to you to try.  (Hint: the imaginary part of the complex exponential cancels out becuase it's odd, while the even parts add up in a series which you can sum with the same formula you use to sum power series sums like  1/2 + 1/4 + 1/8 + ... ).

Linear Algebra

Now all that's well and good, but what's going on and where did those formulas come from anyway?  I'm glad you asked.  We're actually looking at a change of basis here, in just exactly the same way that the co-ordinate description of a vector changes when you rotate your co-ordinate system.  Once we have the equations for that kind of transformation, we'll see that in fact they look just like the ones we just wrote down - and in fact that explains why the Fourier Transform and Inverse Fourier Transform look alike, becuase they're both "rotations" (forward and back) to a new basis.

So let's say we have a vector (Overscript[r, →]) = [x,  y,  z]  in the (OverHat[i]), (OverHat[j]) , (OverHat[k]) basis, which we would like to write in terms of a new basis, say (OverHat[(Subscript[e, 1])]), (OverHat[(Subscript[e, 2])]), (OverHat[(Subscript[e, 3])]).  In other words, we want to find (Subscript[a, 1]), (Subscript[a, 2]), (Subscript[a, 3]) such that (Overscript[r, →])  = (Subscript[a, 1])(OverHat[(Subscript[e, 1])])+  (Subscript[a, 2])(OverHat[(Subscript[e, 2])])+ (Subscript[a, 1])(OverHat[(Subscript[e, 3])]).  How do we get the a's?  Well if the e's are orthonormal - that is, if their dot product is 1 with themself and  0 with each other - then it turns out to be pretty easy: we just dot both sides with each e in turn, like this:

(Overscript[r, →])  = x(OverHat[i]) + y (OverHat[j]) + z (OverHat[k])	(OverHat[(Subscript[e, 1])])	· (Overscript[r, →])  = (OverHat[(Subscript[e, 1])])	· ( x(OverHat[i]) + y (OverHat[j]) + z (OverHat[k]) ) =  (OverHat[(Subscript[e, 1])])	· ( (Subscript[a, 1])(OverHat[(Subscript[e, 1])])+  (Subscript[a, 2])(OverHat[(Subscript[e, 2])])+ (Subscript[a, 1])(OverHat[(Subscript[e, 3])])) = (Subscript[a, 1]) (OverHat[(Subscript[e, 1])])	·(OverHat[(Subscript[e, 1])])+ (Subscript[a, 2]) 0 + (Subscript[a, 3]) 0 = (Subscript[a, 1])    or (Subscript[a, 1])= x  (OverHat[(Subscript[e, 1])])· (OverHat[i]) +  y  (OverHat[(Subscript[e, 1])])· (OverHat[j]) +  z  (OverHat[(Subscript[e, 1])])· (OverHat[k])

How does this apply to the Fourier Series?  Quite simply: the original vector is our set of number [(Subscript[Y, n ])] that describe the function, and the new basis vectors (as written in the x basis) are the sinusoids (e^(i k x)). Thus the formula for the Fourier (Subscript[C, j])'s is just the same as the last line in equation (8) above ; the complex exponentials are the basis dot products, there's a sum over the various components in the original basis, and we get out one coefficient of the vector in the new basis.

(This needs to be fleshed out.)

All right, now on to what this looks like for continuous functions.

Continuous Fourier Integral

Well, the short version is that it's just the same as equation (4) in the limit where \(TraditionalForm
      \`dx → \(0\ and\ L → ∞\)\).  The equations look like this.

\(Y \((x)\) = 
      \(∫\_\(-∞\)\^\(+∞\)\ 
          C \((k)\) \(e\^\(i\ k\ x\)\) dk\/\(2  π\ \ \)
C \((k)\) = 
        \ \ \ ∫\_\(-∞\)\^\(+∞\)\ 
          Y \((x)\)\ \(e\^\(\(-i\)\ k\ x\)\) dx\)\)

Here the Fourier Coefficients C(k) give the (complex) amplitude of a continuous range of possible sinusoidal functions, each with wavenumber k, which when added together give the original function Y(x). 

A common variation on this theme is a mixed continuous/discrete case, in which typically Y(x) is continuous over a finite range of x values with a fixed boundry condition (i.e. a string).  This kind of situation leads to a discrete set of Fourier Coefficients (Subscript[C, j]) which typically give the modes of vibration of the string.  Since this case is well covered in most texts (including ours), I'll leave it out here.

More Examples

One particularly useful function is the Dirac Delta function δ(x) ...

- end -

Scratch Space

<< DiscreteMath`KroneckerDelta`

\(\(j = 0; 
		Sum[\ Exp[2\ Pi\ I\ j\ k\ /N], {k, 1, N}]\ \)\)

\(\(j = 1; \ \ 		Sum[\ Exp[2\ Pi\ I\ j\ k\ /N], {k, 1, N}]\ \)\)

\(\(j =. ; \ 		Sum[\ Exp[2\ Pi\ I\ j\ k\ /N], {k, 1, N}]\ \)\)

a =  
      (1/N) Sum[ G[j] Exp[2  Pi I j n], {j, -N/2 + 1, N/2}]

Sum[ a Exp[ -2 Pi I j n] , {n, 1, N}]

(∑_(n = 1)^N)
        (   (1/N) 
            Sum[ G[j] Exp[2  Pi I j n], {j, -N/2 + 1, N/2}]  )
            Exp[-2 π I j'  n]

\(Doesn' t\ seem\ that\ Mathematica\ can\ convert\ this\ sum\ of\ sums\ 
      into\ a\ single\ sum\ by\ itself\ 
      \((\(! \(?)
		\(Disappointing . \)\)\)

\)\)

dx = L/N

dk  =  2  Pi/L

(N/2) (dk)

dk dx

sigma plot

\(TraditionalForm
    \`sigma[x_]\  := \ \((1 + Sign[a - x])\) \((1 + Sign[a + x])\)/4\)

\(TraditionalForm
    \`\(Plot[sigma[x]\  /. \ a -> 1, \ {x, \(-3\), \ 3}, \ 
      PlotStyle -> Thickness[0.01]]\ \)\)

(Subscript[y, n]) vs (Subscript[x, n])plot

xNums  =  Table[{j, (Subscript[x, j])}, {j, 1, 10}]

yNums  =  Table[{Sin[ .3 + j*0.2], (Subscript[y, j])}, {j, 1, 10}]

yText  =  
      Table[ {Text[yNums[[j, 2]], 
             {xNums[[j, 1]], yNums[[j, 1]]}]}, {j, 1, 10} ]

Show[Graphics[yText],  PlotRange -> {{0, 10}, {0, 1.2}},  
      Axes -> {True, True}, Ticks -> {xNums, None}, AxesLabel -> {x, y}]

Plot[Sin[x], {x, 0, 2 Pi},  Ticks -> False,  ]

ListPlot[Level[Table[{1, -1}, {10}], {2}], Ticks -> False,  
      PlotJoined -> True]

t = Table[{1, 2}, {20}]

t[[1]]

Level[t, {2}]

--- end of scratch -- -