Fourier Series and Transforms

In this Maple worksheet we will explore some of the properties of Fourier series and transforms. In particular, we will plot the Fourier series for a finite number of terms to get some idea of the rate of convergence of the series. As an example, let's start with the Fourier series for the square wave.

> S:=(4/Pi)*Sum(sin((2*n+1)*x)/(2*n+1),n=0..N);

This series will converge to the square wave in the limit that N tends to infinity. Now let's plot the result for various values of N:

> plot({subs(N=4,S),subs(N=8,S),subs(N=16,S)},x=-Pi..Pi,color=[red,green,blue],numpoints=200);

We can see that the approximation improves as we increase N, although rather slowly. Let's try N=50; we'll need to increase the number of points plotted in order to get sufficient resolution.

> plot(subs(N=50,S),x=-Pi..Pi,numpoints=800);

This looks pretty good, but we see that the series still overshoots near the discontinuities. Let's have a closer look near one of these points:

> plot(subs(N=50,S),x=-Pi/10..Pi/10,numpoints=100);

It looks like the series overshoots to about 1.18 near the discontinuity, at x=0.03. Let's try a larger value of N:

> plot(subs(N=100,S),x=-Pi..Pi,numpoints=1600);

> plot(subs(N=100,S),x=-Pi/20..Pi/20,numpoints=800);

Now the maximum is at about x=0.016, with a height of 1.18. In fact, with a little work one can show that the maximum occurs at x=Pi/2*N, and that the maximum is 1.179, so that the series overshoots by about 18%. This happens for any N; the width over which the overshoot occurs decreases as 1/N. This overshoot is known as the ``Gibbs phenomenon,'' and is symptomatic of any attempt to fit a discontinuous function with a finite number of analytic functions.

Next, let's turn to Fourier transforms, which can be thought of as the limit of a Fourier series when the period of the function being transformed tends to infinity. First, the Fourier transform of the function f(x)=1 for -1<x<1, and zero otherwise:

> f(k):=int(exp(-I*k*x),x=-1..1);

> simplify(%);

> plot(f(k),k=-10..10);

We see that the transform oscillates, which is expected since f(x) has sharp discontinuities. We can try something smoother, like the Gaussian:

> f1(x):=exp(-x^2);

> f1(k):=int(f1(x)*exp(-I*k*x),x=-infinity..infinity);

> plot(f1(k),k=-10..10);

The result is also a Gaussian, which is smooth without oscillations. This is a general rule---the smoother f(x) is, the smoother the transform will be. Let's try another-the Lorentzian:

> f2(x):=1/(x^2+a^2);

Now let's use Maple's transform package to compute the transform. We first need to load the package:

> with(inttrans);

Now calculate the transform:

> f2(k):=fourier(f2(x),x,k);

You can see that Maple is careful-it did not assume that a is real, or positive. Plot the result, assuming that a=1:

> plot(subs(a=1,f2(k)),k=-5..5);