four.html

> restart;

Fourier Polynomials and Series

>

Fourier Polynomials

Fourier polynomials provide a way of approximating general periodic functions by sums of very simple periodic functions, namely the familiar sine and cosine functions, shifted and scaled. We will consider functions, f(x), that are 2 -periodic, so that f(x + ) = f(x). For example, any function of the form or , where k is a positive integer while and are constants, is 2 -periodic, although it may have a smaller period. That is, the period of , say, is , so that it repeats after every interval of length , but then of course it repeats after every interval of length 2 . The integer k gives the frequency of the sine or cosine function, so large values of k correspond to very wiggly graphs. The numbers and give the amplitudes . Recall also that the graph of cos(x) is obtained from the graph of sin(x) by a horizontal shift of /2, and the graph of Acos(kx) is likewise obtained from the graph of Asin(kx) by a horizontal shift of length /(2k). What do you get if you shift the graph of Asin(kx) horizontally by some other distance? It turns out that any such horizontal shift of Asin(kx) can be represented in the form

(1)

where the amplitude A is equal to . Functions of the form (1) are convenient to work with. Sums of such functions, with possibly an added constant which we denote by , are called Fourier polynomials. Thus a Fourier polynomial of degree n is a function of the form F(x) =

A Fourier polynomial can be specified by giving the Fourier coefficients , ... . For example, suppose , with all other coefficients being 0. The corresponding Fourier polynomials is then

F(x) = . Let's graph the three functions , , and F(x) on the same axes:

> plot([1+cos(x)/2,cos(8*x)/8+sin(8*x)/8,1+cos(x)/2+cos(8*x)/8+sin(8*x)/8], x=-Pi..Pi, color=[blue, green,black]);

>

Note that we only graphed between and , a practice we will continue throughout the worksheet. Observe that the high frequency term is a shifted sine with amplitude ~ 0.177 and that the graph of F(x) is the superposition of this high frequency graph with that of the vertically shifted cos(x)/2.

It will be convenient to specify the Fourier coefficients by giving the constant term , and then two lists for the coefficents , and . The procedure FPoly defined below, takes the constant term , the two lists and the degree and returns the Fourier polynomial as an expression. For our example, the coefficient list for the a's is [1/2,0,0,0,0,0,0,0,0,1/8] while the coefficient list for the b's is [0,0,0,0,0,0,0,0,0,1/8]. Be sure to execute the next command line.

> FPoly:=proc(a0,a,b,n) local k;
a0+sum('a[k]*cos(k*x)+b[k]*sin(k*x)','k'=1..n); end;

As an example let's plot the degree 9 Fourier polynomial F(x) = .

> plot(
FPoly(.25, [0,0,.5,0,.25,0,0,0,.125],[0,0,0,0,0,0,0,0,0],9),x=-Pi..Pi);

>

Note that this is the graph of an even function of x, and that all the sine coefficients were 0. This is an example of the following general fact you will find useful.

Odd and Even:

A Fourier polynomial is an even function when all its sine coefficients are 0. A Fourier polynomial with its constant term removed is an odd function when all the cosine coefficients are 0.

>

Problem 1. By trial and error, estimate the Fourier coefficients of the degree 4 Fourier polynomials whose graphs are shown below. All plots are from to and (hint) each has no more than 3 nonzero coefficients.

>

The Fourier Polynomials of a Function

If f(x) is a 2 periodic function, we define the Fourier coefficients of f to be

, for k = 1,2,...n.

The degree n Fourier polynomial with these coefficients is called the degree n Fourier polynomial of f . The infinite series whose partial sums are the Fourier polynomials of f is called the Fourier series of f . The degree n Fourier polynomial of f is the Fourier polynomial of degree n that best approximates the function f over the whole interval [- , ] in a certain technical sense we won't go into here.

The calculation of the integrals that give the Fourier coefficients is usually very tedious. We provide Maple procedures to do these calculations. The procedures acoefs and bcoefs take the function f and the degree n and return lists of the coefficients , resp. for k = 1..n. A separate procedure called azero returns to constant coefficient . We also provide similar procedures called Nacoefs, Nbcoefs, Nazero that do the same things with all integrals computed numerically. You may find it useful to try both kinds on the problems. Be sure to execute the following commands.

> acoefs:=proc(f,n) local k;
[seq(int(f(x)*cos(k*x),x=-Pi..Pi)/Pi,k =1..n)];
end;
bcoefs:=proc(f,n) local k; [seq(int(f(x)*sin(k*x),x=-Pi..Pi)/Pi,k=1..n)];
end;
azero:=proc(f) int(f(x),x=-Pi..Pi)/(2*Pi); end;
Here are the purely numerical versions of these procedures
Nacoefs:=proc(f,n) local k;
[seq(evalf(Int(f(x)*cos(k*x),x=-Pi..Pi))/evalf(Pi),k=1..n)];
end;
Nbcoefs:=proc(f,n) local k;
[seq(evalf(Int(f(x)*sin(k*x),x=-Pi..Pi))/evalf(Pi),k=1..n)];
end;
Nazero:=proc(f) evalf(Int(f(x),x=-Pi..Pi))/evalf(2*Pi); end;

Let's try this out on the simple function below, whose graph is shown.

> f:=x->piecewise(x>0,1,0);

> plot(f(x),x=-3..3);

>

The part of the graph of f between -3 and 0 does not show as the function is zero there. We compute the Fourier coefficients up to degree 20. We will call the coefficient lists a and b. We print out the Fourier polynomial of degree 9 and we plot the graph of f together with the graphs of the Fourier polynomials of degrees 5,10,15 and 20. Note that for k 's except k = 0 . This illustrates an extension of the odd/even principle we discussed above.

If a function is even, all its sine coefficients will be 0. If a function can be made into an odd function by some vertical shift, then all its cosine coefficients will be 0.

For our present function f , shifting vertically down by 1/2 produces an odd function. It is often very useful to observe even or odd symmetry of the function you are dealing with and to simply set the appropriate coefficient list equal to a list of zeros without doing any calculation.

> a:=acoefs(f,20): b:=bcoefs(f,20): a0:=azero(f):

> FPoly(a0,a,b,9);

> plot([f(x),seq(FPoly(a0,a,b,5*n),n=1..4)],x=-Pi..Pi);

>

Problem 2 Consider the function g(x) = for x between and repeated periodically.

a) Compute and have Maple print the Fourier polynomials of degrees 5,10, and 15.

b) Plot them together with the graph of g(x) on the interval [- , ] .

c) Write a paragraph discussing the Fourier coefficients and your results are to be expected in view of the symmetry of the function g.

>

The Energy Spectrum

We have seen above that the Fourier series of the function f, is a sum of functions of the form

, where the coefficients are computed as integrals of cosines and sines times f(x). We call this function, , the k-th harmonic of f . As noted in the introductory discussion the k-th harmonic is a shifted sine with frequency k and amplitude , for k >0. We will call the square of the amplitude, , the energy of the k-th harmonic of f, if k > 0, and we define . By analogy with some ideas from physics, we say that a -periodic function has energy given by . For functions f(x) with just a finite number of

points of discontinuity, we have the Energy Theorem, whose proof is beyond this course:

Energy Theorem. .

If we interpret the function f as a wave or a signal, this says that the total energy of the wave or signal is the sum of the energies of its harmonics. A graph of as a function of k, is called the energy spectrum of f.

Let's apply this to a new function, h(x), defined as follows. (Execute the next line.)

> h:=x->x-floor(x);

Here is the graph of h , as usual from .

> plot(h(x), x=-Pi..Pi);

>

Let's compute the coefficients of the degree 20 Fourier polynomial for this function. We will use a0, a, and b as the names of the constant term, the cosine coefficient list, and the sine coefficient list.

Warning: a0, a, and b will now have values throughout the worksheet associated with our specific function h, i.e. a0, a and b are GLOBAL variables. When finished with this section it's a good idea to restore them to being just letters again, by the commands a0:='a0', a:='a', b:='b'.

This time, we will use the numerical integration versions of the procedures to calculate the Fourier coefficients.

>

> a0:=Nazero(h); a:=Nacoefs(h, 20); b:=Nbcoefs(h,20);

>

Look at the results. The list a has entries that are 0 except for numerical integration error. Do you see why this makes sense? We now plot the energy spectrum using the listplot command as follows. (Note:We are plotting against k for k = 1..20. We don't plot .)

> with(plots): listplot([seq(a[k]^2+b[k]^2,k=1..20)],style=POINT,symbol=BOX);

>

The spectrum tells us that there is a strong harmonic at k=6, with strong, but less energetic harmonics at k=7,13 and 19.

Let's plot the function h(x) and the degree 20 Fourier polynomial.

> plot([h(x), FPoly(a0,a,b,20)],x=-Pi..Pi, color = [red,black]);

>

Now let's plot h(x) and just the constant term with the four dominant harmonics we identified from the energy spectrum.

> plot([h(x), a0+ a[6]*cos(6*x)+b[6]*sin(6*x) + a[7]*cos(7*x)+b[7]*sin(7*x)+ a[13]*cos(13*x) + b[13]*sin(13*x)+a[19]*cos(19*x)+b[19]*sin(19*x)],x=-Pi..Pi, color=[red,blue]);

>

These harmonics captured some of the shape, but of course didn't do as well as the full Fourier polynomial, which, on the other hand, requires more computation.

Let us see how much of the total energy of the function h is contained in these harmonics. We calculate the total energy, E, in the function and then the total energy, EH, in the harmonics from k=0 to 20 (remember that k=0 is a special case), and finally, the total energy, DHE contained in the dominant k-th harmonics, namely those with k=0, k=6,7,13,19.

> E:= evalf(Int(h(x)^2,x=-Pi..Pi)/Pi);
EH:= 2.0*a0^2+sum(a[k]^2+b[k]^2, k=1..20);
DHE:= 2.0*a0^2 +a[6]^2+b[6]^2+a[7]^2+b[7]^2+a[13]^2+b[13]^2+a[19]^2+b[19]^2;
EH/E; DHE/E;

>

Thus we see that in this case, about 91.9% of the total energy is contained in the constant term and the 4 dominant harmonics we identified from the energy spectrum.

Problem 3. Consider the function Sq defined below, with graph as shown. (Sq is short for "square wave.")

> Sq:=x-> signum(sin(Pi*(x-1/4)));

> plot(Sq(x),x=-Pi..Pi);

>

a) Using the numerical versions of the functions to find coefficients e.g. Nacoefs, plot the 10th and the 20th degree Fourier polynomials for Sq together with the graph of Sq itself.

b) Use the energy spectrum to identify the 3 most energetic harmonics and plot the sum of these and the constant term along with the graph of Sq.

c) Find the percentage of the total energy of Sq contained in the constant term and the 3 most energetic harmonics you found in part b).

Problem 4. Consider the function q(x) = for x in the interval .

a) Compute the 10th degree Fourier polynomial and display it. (Use the acoefs, bcoefs procedures for this.)

b) Use a) to guess a formula for the energy of the k-th harmonic.

c) * Verify your guess in part b) either by hand calculation or by using Maple. If you use Maple you should execute the following command first

> assume(k, integer);

>

d) Apply the Energy Theorem to this function and your guess for the energy in the k-th harmonic to find a formula, which will involve , for the sum of the infinite series

+ .......

e) Use your result in d) to find an estimate for using only simple arithmetic and a single square root. You should, of course, have Maple do all the arithmetic! How could you get a better estimate, using only the same kind of simple arithmetic?

Homework Problems

Your homework for this worksheet consists of Problems 1-4 above.

Last modified August 1999.

>