Fourier Series

We start with the classical Fourier expansion of a periodic function and then present some important expansions that are widely used in applications. We will assume that all periodic functions have period 2ℓ.

There are two equivalent forms of Fourier series---one is complex form

\[ f(x) \,\sim\, \sum_{k=-\infty}^{\infty} \alpha_k e^{k{\bf j} \pi x/\ell} \]
where
\[ \alpha_k = \frac{1}{2\ell} \int_{-\ell}^{\ell} f(x)\, e^{-k{\bf j} \pi x/\ell} \,{\text d} x , \qquad k=0, \pm 1, \pm 2, \ldots ; \]

and another is a trigonometric form

\[ f(x) \,\sim\, \frac{a_0}{2} + \sum_{k=1}^{\infty} \left[ a_k \cos \frac{k \pi x}{\ell} + b_k \sin \frac{k \pi x}{\ell} \right] , \]
where coefficients are defined according to Euler--Fourier formulas (they were discovered by Leonhard Euler many years before Fourier used them in his research):
\begin{align*} a_0 &= \frac{1}{\ell} \int_{-\ell}^{\ell} f(x)\,{\text d} x , \\ a_k &= \frac{1}{\ell} \int_{-\ell}^{\ell} f(x)\, \cos \frac{k \pi x}{\ell} \,{\text d} x , \qquad k= 1, 2, 3, \ldots ; \\ b_k &= \frac{1}{\ell} \int_{-\ell}^{\ell} f(x)\, \sin \frac{k \pi x}{\ell} \,{\text d} x , \qquad k= 1, 2, 3, \ldots . \end{align*}

The complex Fourier series is more elegant and shorter to write down than the one expressed in term of sines and cosines, but it has the disadvantage that the coefficients might be complex even if the given function is real-valued. The trigonometric form is based on Euler's formulas:

\[ \sin \theta = \frac{1}{2{\bf j}} \,e^{{\bf j}\theta} - \frac{1}{2{\bf j}} \,e^{-{\bf j}\theta} = \Im \,e^{{\bf j}\theta} = \mbox{Im} \,e^{{\bf j}\theta}, \qquad \cos \theta = \frac{1}{2} \,e^{{\bf j}\theta} - \frac{1}{2} \,e^{-{\bf j}\theta} = \Re \,e^{{\bf j}\theta} = \mbox{Re} \,e^{{\bf j}\theta}. \]

Here, j is the unit vector in the positive vertical direction on the complex plane, so \( {\bf j}^2 =-1. \) Mathematica has four default commands to calculate Fourier series:

FourierSeries        (* to calculate complex coefficient expansion *)
FourierTrigSeries  (* to calculate standard Fourier expansion via sine and cosine *)
FourierCosSeries (* to calculate cosine Fourier series *)
FourierSinSeries   (* to calculate sine Fourier series *)

 

Sometimes, it is convenient to use polar form:

\[ f(x) \,\sim\, \frac{a_0}{2} + \sum_{k=1}^{\infty} A_k \sin \left( \frac{k \pi x}{\ell} + \phi_k \right) = \frac{a_0}{2} + \sum_{k=1}^{\infty} A_k \cos \left( \frac{k \pi x}{\ell} - \varphi_k \right) , \]
where \( A_k = \sqrt{a_k^2 + b_k^2} \quad\mbox{and} \quad \varphi_k = \arctan \left( b_k / a_k \right) , \quad \phi_k = \arctan \left( a_k / b_k \right) . \)

Before we get into the topic of convergence, we need to define first a couple of terms that we’ll run into in the rest of the section. First, we say that f(x) has a jump discontinuity at x = a if the limit of the function from the left, denoted f(a-0) or f(a-) and the limit of the function from the right, denoted f(a+0) or f(a+), both exist and \( f(a+0) \ne f(a-0) . \) A function f(x) is called piecewise continuous on an interval [a,b], if it is continuous on [a,b] except for at most finitely many points \( x_1 , x_2 , \ldots , x_n \) at each of them the function has both the left-hand and right-hand limits: \( f(x_k -0) \ \mbox{ and }\ f(x_k +0) , \quad k=1,2,\ldots , n . \) Now we give the formal definition.

Let a and b be real numbers such that a < b. A function \( f\, : \, [a,b] \mapsto \,\mathbb{R} \) is said to be piecewise continuous on \( [a,b] \) if the following conditions are satisfied:
  • there exists a finite set \( \{ x_1 , x_2 , \ldots , x_n \} \subset (a, b) \) such that \( x_1 < x_2 < \cdots < x_n \) and f is continuous on each interval
    \[ \left( a , x_1 \right) , \quad \left( x_k , x_{k+1} \right) , \quad k=1,2,\ldots , n-1, \quad \left( x_n , b \right) ; \]
  • all the following one-sided limits exist
    \[ \lim_{x\downarrow a} f \left( x \right) , \quad \lim_{x\uparrow x_k} f\left( x \right) , \quad \lim_{x\downarrow x_k} f\left( x \right) ,\quad k=1,2,\ldots , n-1, \quad \lim_{x \uparrow b} f\left( x \right) . \]
A function \( f\, : \, \mathbb{R} \mapsto \,\mathbb{R} \) is piecewise continuous on \( \mathbb{R} , \) if it is piecewise continuous on every finite subinterval of \( \mathbb{R} . \)

Next, we say that f(x) is piecewise smooth if the function can be broken into distinct pieces and on each piece both the function and its derivative, f'(x), are continuous. A piecewise smooth function may not be continuous everywhere; however, the only discontinuities that are allowed are at a finite number of jump discontinuities. If derivatives in the above condition are replaced by more weak condition so that function f is piecewise monotonic, then such function is said to satisfy the Dirichlet conditions. A function \( f\, : \, \mathbb{R} \mapsto \,\mathbb{R} \) is piecewise smooth on \( \mathbb{R} , \) if it is piecewise smooth on every finite subinterval of \( \mathbb{R} . \)

Example. The function

\[ f(x) = \begin{cases} 2x , & \ \mbox{ if } 0 < x < 1, \\ 1, & \ \mbox{ if } 1 < x < 2, \end{cases} \]
is piecewise smooth on [0,2], but is not continuous on [0,2]. However, \( f(x) = |x|^{1/2} \) is continuous on [-1,1], but is not piecewise smooth on [-1,1]. ■

The finite sum
\[ F_N (x) = \frac{a_0}{2} + \sum_{k=1}^{N} \left[ a_k \cos \frac{k \pi x}{\ell} + b_k \sin \frac{k \pi x}{\ell} \right] \]
is called a trigonometric polynomial of order N. There are three known (another one will be considered later in Cesàro section) kinds of convergence related to the Fourier series:
  • We say that the series \( \sum_{n\ge 0} f_n (x) \) converges to f(x) pointwise in the interval (a,b), if for each \( x \in (a,b) \)
    \[ \left\vert f(x) - \sum_{n=0}^N f_n (x) \right\vert \,\mapsto \, 0 \quad \mbox{as} \quad N \,\to\,\infty . \]
    That is, for each fixed \( x \in (a,b) , \) the numerical sequence \( \sum_{n\ge 0} f_n (x) \) converges to the number f(x).
  • The infinite sum \( \sum_{n\ge 0} f_n (x) \) converges uniformly in [a,b], if
    \[ \max_{a \le x \le b} \left\vert f(x) - \sum_{n=0}^N f_n (x) \right\vert \,\mapsto \, 0 \quad \mbox{as} \quad N \,\to\,\infty . \]
    That is, the overall “distance” between the function f(x) and the partial sums \( F_N (x) = \sum_{n= 0}^N f_n (x) \) converges to zero. Notice that the endpoints of the interval are included in the definition.
  • We say that the series \( \sum_{n\ge 0} f_n (x) \) converges to f(x) in the mean-square (or L2 sense) in (a,b), if
    \[ \int_a^b \left\vert f(x) - \sum_{n=0}^N f_n (x) \right\vert^2 \,{\text d} x \,\mapsto \, 0 \quad \mbox{as} \quad N \,\to\,\infty . \]
    That is, the the “distance” between f(x) and the partial sums in the mean-square sense converges to zero. ■
It is obvious from the definition that uniform convergence is the strongest of the three, since uniformly convergent series will clearly converge pointwise, as well as in the L2 sense. The converse is not true, since not every pointwise or L2 convergent series is uniformly convergent. Between the pointwise and L2 convergence, neither is stronger than the other, since there are series that converge pointwise, but not in L2, and vice versa.

Theorem: Let f(x) be a square integrable function on the interval \( [- \ell , \ell ] , \) that is, \( \| f \|^2 = \int_{_\ell}^{\ell} |f(x)|^2 {\text d}x < \infty . \) Then the Fourier series for his function \( f(x) \,\sim \, \frac{a_0}{2} + \sum_{k\ge 0} \left[ a_k \cos \frac{k \pi x}{\ell} + b_k \sin \frac{k \pi x}{\ell} \right] , \) where its coefficients are determined via Euler--Fourier formulas:
\[ \begin{split} a_k &= \frac{1}{\ell} \, \int_{-\ell}^{\ell} f(x)\,\cos \frac{k \pi x}{\ell} \, {\text d} x , \quad k=0,1,2,\ldots , \\ b_k &= \frac{1}{\ell} \, \int_{-\ell}^{\ell} f(x)\,\sin \frac{k \pi x}{\ell} \, {\text d} x , \quad k=1,2,\ldots ; \end{split} \]
converges to f(x) in the mean square sense. ■

Andrey Kolmogorov (1903--1987) from Moscow University (Russia), as a student at the age of 19, in his very first scientific work, constructed an example of an absolutely integrable function whose Fourier series diverges almost everywhere (later improved to diverge everywhere).

The first sufficient conditions for the pointwise convergence of the Fourier series were discovered by a German mathematician Johann Peter Gustav Lejeune Dirichlet (1805--1859). A function that satisfies the Dirichlet conditions is also called piecewise monotone. Such function must have right and left limits at each point of discontinuity.

Theorem (P. Dirichlet, 1829): Suppose that a periodic real-valued function f(x) satisfies the following conditions (known as the Dirichlet conditions): Then the function f(x) is represented by a convergent to \( \frac{1}{2} \left( f(a-0) + f(a+0) \right) \) Fourier series at every point x=a.
Theorem: Assume that
\[ F(x) \,=\, \frac{a_0}{2} + \sum_{k= 1}^{\infty} a_k \cos \left( \frac{k \pi x}{\ell} \right) + b_k \sin \left( \frac{k \pi x}{\ell} \right) \]
is the Fourier series for a piecewise continuous function f(x) over \( [- \ell , \ell ] . \) If its derivative f'(x) is piecewise continuous on the interval \( [- \ell , \ell ] \) and has both a left- and right-hand derivative at each point in this interval, then F(x) is pointwise convergent for all \( x \in [- \ell , \ell ] . \) The relation \( F(x) = f(x) \) holds at all points \( x \in [- \ell , \ell ] , \) where f(x) is continuous. If x = a is a point of discontinuity of f, then
\[ F(x) \,=\, \frac{f(a+0)+ f(a-0)}{2} , \]
where f(a-0) and f(a+0) denote the left- and right-hand limits, respectively.

Theorem: Suppose that \( \left\{ (x_k, y_k ) \right\}_{k=0}^N \) are N+1 points, where \( y_k = f(x_k ) , \) and the abscissas are equally spaced:

\[ x_k = -\ell + \frac{2\ell \,k}{N} \qquad \mbox{for} \quad k=0,1,2,\ldots N. \]
If f(x) is periodic with period \( 2\ell \) and 2M < N, then there exists a trigonometric polynomial
\[ F_M (x) = \frac{a_0}{2} + \sum_{j=1}^M \left[ a_j \cos \left( j\pi x \right) + b_j \sin \left( j\pi x \right) \right] \]
that minimizes the quantity
\[ \sum_{k=1}^N \left\vert f(x_k ) - F_M (x_k ) \right\vert^2 . \]
The coefficients ak and bk of this trigonometric polynomial are computed with the formulas:
\begin{align*} a_j &= \frac{2}{N} \,\sum_{k=1}^N f(x_k ) \,\cos \left( j\,\pi \,x_k \right) \qquad\mbox{for}\quad j=0,1,\ldots , M , \\ b_j &= \frac{2}{N} \,\sum_{k=1}^N f(x_k ) \,\sin \left( j\,\pi\, x_k \right) \qquad\mbox{for}\quad j=1,2,\ldots , M . \end{align*}

Example: Consider the function \( f (x) = x^2 \) on the interval [-1,1]. Let us take 16 equally spaced points \( x_k = -1 + \frac{2k}{15} , \quad k=0,1,2,\ldots , 15. \) Suppose we want to find the trigonometric polynomial approximation for M = 6 to the 15 data points \( \left\{ (x_k, y_k ) \right\}_{k=1}^{15} . \)

Since the periodic extension is assumed, at a point of discontinuity x = 1, the function value f(1) must be computed using the formula

\[ f(1) = \frac{f(1-0) + f(1+0)}{2} = \frac{1+1}{2} =1. \]
The function \( f (x) = x^2 \) is an even continuous function that is extended periodically; hence the coefficients for the sine terms are all zero (i.e. \( b_j =0 \) for all j). The trigonometric polynomial of degree M = 6 involves only the cosine terms, and we get
\[ a_j = \frac{2}{15} \,\sum_{k=1}^{15} \left( \frac{2k}{15} -1 \right)^2 \cos \left( j\,\pi\, x_k \right) \qquad\mbox{for}\quad j=0,1,2,\ldots , 6 . \]
Now we ask Mathematica for help:
Do[a[j_] := (2/15)* Sum[(2*k/15 - 1)^2 *Cos[j*Pi*(2*k/15 - 1)], {k, 1, 15}], {j, 0, 6}]
Then we calculate the trigonometric polynomial and plot it.
C6[x_] := a[0]/2 + Sum[a[j]*Cos[j*Pi*x], {j, 1, 6}]
Plot[{x^2 , FC6[x]}, {x, -1, 1}, PlotRange -> {-0.1, 1}, PlotStyle -> {{Thick, Blue}, {Thick, Orange}}]
As we see from the above graph, six term approximation gives relatively good representation of the given function. ■

Computing an integral in Mathematica is fairly painless, and it's tempting to simply use a partial Fourier sum depending to the number of terms n.

Let's consider the following half-wave rectifier:

f[t_] = Piecewise[{{Cos[t], -Pi/2 < t < Pi/2}, {0, -Pi < t < -Pi/2 && Pi/2 < t < Pi}}]
Out[1]= { Cos[t] - Pi/2 < t < Pi/2
             { 0 True
Plot[f[t], {t, -Pi, Pi}, PlotStyle -> Thick]

Mathematica will treat the function f[t] outside the interval [-π , π ] as zero

Plot[f[t], {t, -3*Pi, 3*Pi}, PlotStyle -> Thick]

As we see, Mathematica rescales vertical and horizontal lines so that the picture becomes almost symmetric. You can return to the proper picture by using AspectRatio option:

Plot[f[t], {t, -3*Pi, 3*Pi}, PlotStyle -> Thick, AspectRatio -> 1/3/Pi, Ticks -> {True, None}]

So we get exactly the same graph. Now we are ready to expand this function into Fourier series by finding its Fourier coefficients:

\begin{align*} a_0 &= \frac{1}{\pi} \int_{-\pi /2}^{\pi /2} \cos x\,{\text d} x = \frac{2}{\pi} , \\ a_k &= \frac{1}{\pi} \int_{-\pi /2}^{\pi /2} \cos x\, \cos \left( k\, x \right) {\text d} x = -\frac{2}{\pi}\, \frac{1}{k^2 -1} \, \cos \frac{k\pi}{2} , \qquad k= 2, 3, \ldots ; \\ b_k &= \frac{1}{\pi} \int_{-\pi /2}^{\pi /2} \cos x\, \sin \left( k\, x \right) {\text d} x =0, \qquad k= 1, 2, 3, \ldots . \end{align*}
Then the corresponding Fourier series becomes
\[ f(t) \,=\,\frac{1}{\pi} + \frac{1}{2}\, \cos t - \frac{2}{\pi} \, \sum_{k\ge 1} \frac{1}{k^2 -1} \,\cos \left( k\, x \right) . \]
A partial sum with N = 20 terms gives a very good approximation:

We show how these commands work with other examples.

Example: We consider a piecewise continuous function that cannot be extended into Taylor series because it is identically zero on the interval (1,2). However, we can find its Fourier series.

f[x_] = Piecewise[{{1 - x, 0 < x < 1}, {0, 1 < x < 2}}]
Out[2]= { 1-x, if 0<x<1 and 0 for 1<x<2.

The standard Mathematica command FourierTrigSeries provides you the Fourier series of the function that is extended periodically from the standard interval (-π,π):

FourierTrigSeries[f[x], x, 3]
Out[2]= 1/(4 \[Pi]) + (2 Cos[x] Sin[1/2]^2)/\[Pi] + (Cos[2 x] Sin[1]^2)/(2 \[Pi]) + (2 Cos[3 x] Sin[3/2]^2)/(
9 \[Pi]) + ((2 - 2 Sin[1]) Sin[x])/(2 \[Pi]) - ((-2 + Sin[2]) Sin[2 x])/(4 \[Pi]) - ((-3 + Sin[3]) Sin[3 x])/(9 \[Pi])

With option FourierParameters, it is equivalent to

FourierTrigSeries[f[x], x, 3, FourierParameters -> {1, 1}]
We can plot the partial sum with say 50 terms. Since Fourier partial sums oscillate near points of discontinuity exhibiting so called Gibbs phenomenon, we indicate maximum wiggle (overshoot) and minimum one:

curve = FourierTrigSeries[f[x], x, 50];
Plot[curve, {x, -1.5, 13.5}, PlotRange -> {-0.3, 1.3}, Ticks -> {{Pi, 2*Pi, 3*Pi}, {1.09, -0.09}}, PlotStyle -> Thick]

However, it is not what we want: we need the Fourier series of the function f[x] extended periodically from the interval of length 2:

FourierTrigSeries[f[x], x, 3, FourierParameters -> {1, Pi}]
Out[3]= 1/4 + (2 Cos[\[Pi] x])/\[Pi]^2 + (2 Cos[3 \[Pi] x])/(9 \[Pi]^2) + Sin[\[Pi] x]/\[Pi] + Sin[2 \[Pi] x]/(2 \[Pi]) + Sin[3 \[Pi] x]/( 3 \[Pi])

To check, we find the Fourier coefficients manually:

Integrate[f[x]*Cos[n*Pi*x], {x, 0, 2}, Assumptions -> Element[n, Integers]]
Out[4]= (1 - Cos[n \[Pi]])/(n^2 \[Pi]^2)
So
\[ a_0 = \frac{1}{2} , \qquad a_n = \frac{1- (-1)^n}{n^2 \pi^2} , \quad n=1,2,\ldots . \]
Integrate[f[x]*Sin[n*Pi*x], {x, 0, 2}, Assumptions -> Element[n, Integers]]
Out[5]= (n \[Pi] - Sin[n \[Pi]])/(n^2 \[Pi]^2)

which is 1/(n*π). Hence, its Fourier series becomes

\[ f(x) = \frac{1}{4} + \sum_{n\ge 1} \left[ \frac{1- (-1)^n}{n^2 \pi^2} \, \cos \left( n \pi x \right) + \frac{1}{n\pi}\, \sin \left( n \pi x \right) \right] \]
We can also represent f(x) in other forms:
\[ f(x) = \frac{1}{4} + \sum_{n\ge 1} \sqrt{ \left( \frac{1- (-1)^n}{n^2 \pi^2} \right)^2 + \frac{1}{n^2 \pi^2}} \, \sin \left( n \pi x + \arctan \frac{1- (-1)^n}{n \pi} \right) = \frac{1}{4} + \sum_{n\ge 1} \sqrt{ \left( \frac{1- (-1)^n}{n^2 \pi^2} \right)^2 + \frac{1}{n^2 \pi^2}} \, \cos \left( n \pi x - \mbox{arccot} \frac{1- (-1)^n}{n \pi} \right) . \]
Now we plot some partial sums, starting with the given function itself:

Plot[f[x], {x, -1, 3}, PlotRange -> {0, 1}, AspectRatio -> Automatic, Axes -> {True, False}, PlotStyle -> Thick]

Then we plot some finite sums:

curve5Pi = FourierTrigSeries[f[x], x, 5];
Plot[%, {x, -1, 3}]

Then we repeat calculations with 50 terms:


Next we plot the correct Fourier series partial sums:

curve5=FourierTrigSeries[f[x], x, 50, FourierParameters -> {1, Pi}];
Plot[%, {x, -1, 3}]
■

 

The option FourierParameters has two parameters and when applied, it looks as FourierParameters->{a,b}
In trigonometric form, with setting \( {\bf FourierParameters} \,-> \,\{ a, b \} \) the following series is returned:

\[ \left\vert \frac{b}{\pi} \right\vert^{(1+a)/2} \left[ \frac{a_0}{2} + \sum_{k=1}^n \left( a_k \,\cos (bkt) + b_k \,\sin (bkt) \right) \right] , \]
with
\[ a_k = \left\vert \frac{b}{\pi} \right\vert^{(1+a)/2} \,\int_{-\pi/|b|}^{\pi/|b|} f(t)\,\cos (bkt) \,{\text d}t , \qquad k=0,1,2,\ldots , \]
\[ b_k = \left\vert \frac{b}{\pi} \right\vert^{(1+a)/2} \,\int_{-\pi/|b|}^{\pi/|b|} f(t)\,\sin (bkt) \,{\text d}t , \qquad k=1,2,\ldots . \]

Syntax for the FourierSeries command is: FourierSeries [ function , variable , number of terms ]

Changing the FourierParameters setting allows control over the limits of integration on the coefficients, and therefore, the base frequency of the series. Some work will go into calculating what value for “b” will give the proper limits of integration. There is ordinarily no reason to change “a” from its default setting, which is 1.

Here are two important applications of Fourier series:
\[ \cos \left( a\,\sin x \right) = J_0 (a) + \sum_{k\ge 1} J_{2k} (a) \,\cos \left( 2kx \right) . \]
\[ \delta (x-a) = \frac{1}{2\pi} \, \int_{-\infty}^{\infty} e^{{\bf j} (x-a)t} \, {\text d} t = \lim_{n\to \infty} \, \frac{1}{2\pi} \, \sum_{k=-n}^n \, e^{{\bf j} k(x-a)} \quad \left( = \frac{\sin \left[ \left( n+ \frac{1}{2} \right) (x-a) \right]}{2\pi\, \sin \frac{x-a}{2}} \right) , \]
where j is a unit vector in positive vertical direction on the complex plane, so \( {\bf j}^2 =-1 . \) 

assume(n,Type::NonNegInt);
f1 := piecewise([-2<x<0,-2],[0<<x<2,x]);
L := 2:;
a0 :=int(f1(x),x=-L..L)/(2*L);
an := simplify(int(f1(x)*cos(n*PI*x/L),x=-L..L)/L);
bn := simplify(int(f1(x)*sin(n*PI*x/L),x=-L..L)/L);
\left\{\begin{array}{cl} x & \text{\ if\ \ }x\in \left(0,2\right)\\ -2 & \text{\ if\ \ }x\in \left(-2,0\right) \end{array} -\frac{1}{2} -\frac{{\left(-1\right)}^{n+1}\,{\left({\left(-1\right)}^n-1\right)}^2}{n^2\,\pi ^2} -\frac{2\,\left(2\,{\left(-1\right)}^n-1\right)}{\pi \,n} 

a:=m->subs(an,n=m);
a(0):=a0;
b:=m->subs(bn,n=m);
S:=N->a(0) + sum(a(n)*cos(n*PI*x/L)+b(n)*sin(n*PI*x/L),n=1..N);
plot(f1,S(5),x=-5..5)
picture 

a(0):=a0;
b:=m->subs(bn,n=m);
S:=N->a(0) + sum(a(n)*cos(n*PI*x/L)+b(n)*sin(n*PI*x/L),n=1..N);
plot(f1,S(15),x=-5..5)
plot::Rectangle(-3..0, -1..1) picture
a:=m->subs(an,n=m); a(0):=a0; b:=m->subs(bn,n=m); S:=N->a(0) + sum(a(n)*cos(n*PI*x/L)+b(n)*sin(n*PI*x/L),n=1..N); plot(f1,S(30),x=-5..5)
plot::Rectangle(-3..0, -1..1) picture