The Fourier Series

Outline for today

• Discrete Fourier Transform
• Applications of the DFT

The Fourier series

The Fourier series of a function is an expansion into sines and cosines:

\[ f(x)=a_{0}+a_{1} \cos x+b_{1} \sin x+a_{2} \cos 2 x+b_{2} \sin 2 x+\cdots \]

. . .

How do we find the coefficients \(a_{0}, a_{1}, b_{1}, a_{2}, b_{2}, \ldots\)?

. . .

We use orthogonality. The functions \(\cos x, \cos 2 x, \cos 3 x, \ldots\) are orthogonal on \([0,2 \pi]\)… Same is true for the sines.

. . .

To compute a coefficient like \(b_{1}\), multiply both sides by the corresponding function \(\sin x\) and integrate from 0 to \(2 \pi\). In other words, take the inner product of both sides with \(\sin x\) :

\[ \int_{0}^{2 \pi} f(x) \sin x d x=a_{0} \int_{0}^{2 \pi} \sin x d x+a_{1} \int_{0}^{2 \pi} \cos x \sin x d x+b_{1} \int_{0}^{2 \pi}(\sin x)^{2} d x+\cdots \]

\[ b_{1}=\frac{\int_{0}^{2 \pi} f(x) \sin x d x}{\int_{0}^{2 \pi}(\sin x)^{2} d x}=\frac{(f, \sin x)}{(\sin x, \sin x)} \]

Analogy with projections

The component of the vector \(b\) along the line spanned by \(a\) is \(b^{\mathrm{T}} a / a^{\mathrm{T}} a\).

A Fourier series is projecting \(f(x)\) onto \(\sin x\). Its component \(p\) in this direction is exactly \(b_{1} \sin x\).

Euler’s formula

The complex exponential \(e^{i x}\) is a combination of \(\cos x\) and \(\sin x\):

\[ e^{i x}=\cos x+i \sin x \]

So we can rewrite the Fourier series using complex exponentials:

\[ f(x)=c_{0}+c_{1} e^{i x}+c_{2} e^{2 i x}+\cdots = \sum_k c_k\cdot e^{i k x} \]

pause

The coefficients \(c_{0}, c_{1}, c_{2}, \ldots\) are related to the \(a\)’s and \(b\)’s.

For example, \(c_{1}=a_{1}-i b_{1}\).

The formula for finding the coefficients \(c_{k}\) is the same as before, but now we use the complex exponential functions \(e^{i k x}\) instead of the sines and cosines:

\[ c_{k}=\int_{0}^{2 \pi} f(x) e^{-i k x} d x \]

Discrete Fourier Transform

Discrete Fourier Series

What if we have only a series of N discrete values \(y_{0}, y_{1}, y_{2}, \ldots\), instead of a continuous function \(f(x)\)?

. . .

First, we imagine that the values are periodic, with period N. So it’s as if they are repeating beyond the N values we have.

. . .

Then we can identify a natural period \(N\), and rewrite \(x=2 \pi n / N\), so the discrete Fourier series becomes

\[ \begin{align} y[n] = \sum_k c[k]\cdot e^{i 2\pi k \frac{n}{N}} \end{align} \]

. . .

In matrix form, \(\mathbf{y} = F \mathbf{c}\)

where

\[ \begin{align} \begin{bmatrix} y_0 \\ y_1 \\ y_2 \\ \vdots \\ y_{N-1} \end{bmatrix} = F c = \begin{bmatrix} 1 & 1 & 1 & \ldots & 1 \\ 1 & e^{-i 2\pi/N} & e^{-i 4\pi/N} & \ldots & e^{-i 2\pi (N-1)/N} \\ 1 & e^{-i 4\pi/N} & e^{-i 8\pi/N} & \ldots & e^{-i 4\pi (N-1)/N} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & e^{-i 2\pi(N-1)/N} & e^{-i 4\pi(N-1)/N} & \ldots & e^{-i 2\pi(N-1)(N-1)/N} \end{bmatrix} \begin{bmatrix}c_0 \\ c_1 \\ c_2 \\ \vdots \\ c_{N-1}\end{bmatrix} \end{align} \]

DFT matrix

\[ \begin{align} \begin{bmatrix} y_0 \\ y_1 \\ y_2 \\ \vdots \\ y_{N-1} \end{bmatrix} = F c = \begin{bmatrix} 1 & 1 & 1 & \ldots & 1 \\ 1 & e^{-i 2\pi/N} & e^{-i 4\pi/N} & \ldots & e^{-i 2\pi (N-1)/N} \\ 1 & e^{-i 4\pi/N} & e^{-i 8\pi/N} & \ldots & e^{-i 4\pi (N-1)/N} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & e^{-i 2\pi(N-1)/N} & e^{-i 4\pi(N-1)/N} & \ldots & e^{-i 2\pi(N-1)(N-1)/N} \end{bmatrix} \begin{bmatrix}c_0 \\ c_1 \\ c_2 \\ \vdots \\ c_{N-1}\end{bmatrix} \end{align} \]

Visualization

From Wikipedia, Original upload by en:User:Glogger, CC BY-SA 3.0 http://creativecommons.org/licenses/by-sa/3.0/, via Wikimedia Commons

Finding the Fourier coefficients

How can we find the Fourier coefficients? . . .

By analogy with the continuous case,

\[ c_{k}=\int_{0}^{2 \pi} f(x) e^{-i k x} d x \]

We will replace the integral by a sum over the discrete values \(y_{0}, y_{1}, y_{2}, \ldots\).

\[ c[k]= \frac{1}{N} \sum_{n=0}^{N-1} y[n] e^{-i 2 \pi k n / N}, \quad \text{for } k=0,1,2, \ldots, N-1 \]

This defines a matrix \(F^\prime\) given as

\[ F^\prime = \frac{1}{N}\begin{bmatrix} 1 & 1 & 1 & \ldots & 1 \\ 1 & e^{i 2\pi/N} & e^{i 4\pi/N} & \ldots & e^{i 2\pi (N-1)/N} \\ 1 & e^{i 4\pi/N} & e^{i 8\pi/N} & \ldots & e^{i 4\pi (N-1)/N} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & e^{i 2\pi(N-1)/N} & e^{i 4\pi(N-1)/N} & \ldots & e^{i 2\pi(N-1)(N-1)/N} \end{bmatrix} \]

Note: this is \(\frac{1}{N}\) times the complex conjugate of the matrix \(F\): \(F^\prime = \frac{1}{N} F^*\).

Inverse DFT

We have just shown two things:

1. The DFT can be written as a matrix-vector product \(\mathbf{y} = F \mathbf{c}\).

2. The coefficients \(\mathbf{c}\) can be found by \(\mathbf{c} = F^\prime \mathbf{y}\).

3. But of course it is also true from 1 that if \(F\) is invertible, that \(\mathbf{c} = F^{-1} \mathbf{y}\). Therefore, \(F^{-1} = F^\prime=\frac{1}{N} F^*\).

4. This means that \(F\) is unitary. Its rows and columns are orthogonal.

A concrete example

Let’s take an example with \(N=4\) and \(y=(2, 4, 6, 8)\).

. . .

We want to find coefficients \(c_0, c_1, c_2, c_3\) such that our Fourier series correctly reproduces the input values at each of the four points.

. . .

(We can’t find more than 4 coefficients because we only have 4 points of data.)

. . .

\[ \begin{align} y[n] = \sum_k c[k]\cdot e^{i 2\pi k \frac{n}{N}} \end{align} \]

. . .

\[ \begin{aligned} c_{0}+c_{1} e^{i x}+c_{2} e^{2 i x}+c_{3} e^{3 i x}=\left\{\begin{array}{lll} 2 & \text { at } & x=0 \\ 4 & \text { at } & x=\pi / 2 \\ 6 & \text { at } & x=\pi \\ 8 & \text { at } & x=3 \pi / 2 \end{array}\right) \end{aligned} \]

Make a table. Note that \(i^2 = i^6 = -1\), \(i^3 = i^7 = -i\), and \(i^4 = 1\).

For n = 0:

k	\(i 2\pi k \times 0\)	\(e^{i 2\pi k \times 0}\)
0	0	1
1	0	1
2	0	1
3	0	1

For n = 1:

k	\(i 2\pi k \frac{1}{4}\)	\(e^{i 2\pi k \frac{1}{4}}\)
0	0	1
1	\(i\pi/2\)	\(i\)
2	\(i\pi\)	\(-1=i^2\)
3	\(3i\pi/2\)	\(-i=i^3\)

For n = 2:

k	\(i 2\pi k \frac{2}{4}\)	\(e^{i 2\pi k \frac{2}{4}}\)
0	0	1
1	\(i\pi\)	\(-1=i^2\)
2	\(2i\pi\)	\(1=i^4\)
3	\(3i\pi\)	\(-1=i^6\)

For n = 3:

k	\(i 2\pi k \frac{3}{4}\)	\(e^{i 2\pi k \frac{3}{4}}\)
0	0	1
1	\(3i\pi/2\)	\(-i=i^3\)
2	\(3i\pi\)	\(-1 = i^6\)
3	\(9i\pi/2\)	\(i = i^9\)

Solving for the coefficients

\[ \begin{align*} & c_{0}+c_{1}+c_{2}+c_{3}=2 \\ & c_{0}+i c_{1}+i^{2} c_{2}+i^{3} c_{3}=4 \\ & c_{0}+i^{2} c_{1}+i^{4} c_{2}+i^{6} c_{3}=6 \\ & c_{0}+i^{3} c_{1}+i^{6} c_{2}+i^{9} c_{3}=8 \end{align*} \]

In other words, \(F c = y\) for

\[ F=\left[\begin{array}{cccc} 1 & 1 & 1 & 1 \\ 1 & i & i^{2} & i^{3} \\ 1 & i^{2} & i^{4} & i^{6} \\ 1 & i^{3} & i^{6} & i^{9} \end{array}\right] \]

Check that \(F F^\prime = I\):

\[\begin{aligned} F F^\prime =\frac{1}{4}\left[\begin{array}{cccc} 1 & 1 & 1 & 1 \\ 1 & i & i^{2} & i^{3} \\ 1 & i^{2} & i^{4} & i^{6} \\ 1 & i^{3} & i^{6} & i^{9} \end{array}\right]\left[\begin{array}{cccc} 1 & 1 & 1 & 1 \\ 1 & (-i) & (-i)^{2} & (-i)^{3} \\ 1 & (-i)^{2} & (-i)^{4} & (-i)^{6} \\ 1 & (-i)^{3} & (-i)^{6} & (-i)^{9} \end{array}\right]=I \end{aligned} \]

Check in Sympy:

F:
F prime:
F F prime:

\(\displaystyle \left[\begin{matrix}1 & 1 & 1 & 1\\1 & i & -1 & - i\\1 & -1 & 1 & -1\\1 & - i & -1 & i\end{matrix}\right]\)

\(\displaystyle \left[\begin{matrix}0.25 & 0.25 & 0.25 & 0.25\\0.25 & - 0.25 i & -0.25 & 0.25 i\\0.25 & -0.25 & 0.25 & -0.25\\0.25 & 0.25 i & -0.25 & - 0.25 i\end{matrix}\right]\)

\(\displaystyle \left[\begin{matrix}1.0 & 0 & 0 & 0\\0 & 1.0 & 0 & 0\\0 & 0 & 1.0 & 0\\0 & 0 & 0 & 1.0\end{matrix}\right]\)

A simplification

Notice that our matrix F had many terms which were exponents of \(i\).

In general, we can simplify the matrix \(F\) by defining \(w = e^{i 2\pi/N}\). This is the \(N\)th “primitive root of unity”, so that \(w^N = 1\).

. . .

For N = 4, \(w = e^{i \pi/2} = i\).

. . .

Write the first few terms of this: \[ \left[\begin{array}{ccccc} 1 & 1 & 1 & \cdot & 1 \\ 1 & w & w^{2} & \cdot & w^{n-1} \\ 1 & w^{2} & w^{4} & \cdot & w^{2(n-1)} \\ \cdot & \cdot & \cdot & \cdot & \cdot \\ 1 & w^{n-1} & w^{2(n-1)} & \cdot & w^{(n-1)^{2}} \end{array}\right]\left[\begin{array}{c} c_{0} \\ c_{1} \\ c_{2} \\ \cdot \\ c_{n-1} \end{array}\right]=\left[\begin{array}{c} y_{0} \\ y_{1} \\ y_{2} \\ \cdot \\ y_{n-1} \end{array}\right], \quad w = e^{i 2\pi/n} \]

Simplest form of Fourier matrix

Then we can write the Fourier equation as

\[ \left[\begin{array}{ccccc} 1 & 1 & 1 & \cdot & 1 \\ 1 & w & w^{2} & \cdot & w^{n-1} \\ 1 & w^{2} & w^{4} & \cdot & w^{2(n-1)} \\ \cdot & \cdot & \cdot & \cdot & \cdot \\ 1 & w^{n-1} & w^{2(n-1)} & \cdot & w^{(n-1)^{2}} \end{array}\right]\left[\begin{array}{c} c_{0} \\ c_{1} \\ c_{2} \\ \cdot \\ c_{n-1} \end{array}\right]=\left[\begin{array}{c} y_{0} \\ y_{1} \\ y_{2} \\ \cdot \\ y_{n-1} \end{array}\right], \quad w = e^{i 2\pi/n} \]

Terminology

We call the matrix \(F\) the DFT (Discrete Fourier Transform) matrix. It is a square unitary matrix of size \(N \times N\).

The matrix \(F^\prime\) is the inverse DFT matrix. It is the complex conjugate of \(F\) divided by \(N\).

Summary

1. The DFT can be written as a matrix-vector product \(\mathbf{y} = F \mathbf{c}\).

2. The coefficients \(\mathbf{c}\) can be found by \(\mathbf{c} = F^\prime \mathbf{y}\).

3. F is easy to compute, with a simple formula for each entry.

4. (There’s actually a really really fast way to compute the DFT using the FFT algorithm.)

Applications of the DFT

Filtering

The DFT is used in many applications, but one of the most common is filtering.

For example, suppose we have a signal that is a sum of two sinusoids:

\[ f(x) = \sin(2\pi x) + 0.5 \sin(20\pi x) \]

What’s the DFT of this signal?

from scipy.fft import fft, fftshift
yf = fft(y)
yfalone = fft(yalone)
# make a non connected plot of yf

plt.plot(np.fft.fftfreq(200),np.abs(yf), color='red')
plt.plot(np.fft.fftfreq(200),np.abs(yfalone), color='blue')
plt.xlim(-0.25, 0.25)
plt.legend(['Full signal', 'low frequency only'])
plt.show()

What’s with the double peaks?

Reconstruct the signal

We can reconstruct the signal by taking the inverse DFT.

Low-pass filtering

What if we just zero out the coefficient corresponding to the second sinusoid:

yf_new = yf.copy()
yf_new[190] = 0
yf_new[10] = 0
plt.plot(np.abs(fftshift(yf_new)))
plt.title('Filtered DFT')
plt.show()

Compute the inverse DFT:

It’s pretty good!

Try one more thing. Can we describe the change in the DFT as a multiplicative vector? We sure can.

So we can describe the new DVD as the old DFT multiplied by a vector of 1s and 0s.

The total process went like this:

1. Take the DFT of the signal.
2. Multiply the DFT by a vector of 1s and 0s to filter out the high frequency.
3. Take the inverse DFT to get the filtered signal.

Filtering as convolution

Suppose we have a signal and we’d like to try to filter out the high frequencies, but we don’t know which ones they are.

We could try with a simple filter like \([1, 1, 1, 1, 1]/5\). This is a simple moving average filter.

That means that we replace each point in the signal with the average of the 5 points before it.

Mathematically, this is:

\[ f_{\text{filtered}}[n] = \sum_{m=0}^{4} f[n-m]/5 \]

OK. We note that this can also be written as a convolution, between the signal and the vector \(h = [1, 1, 1, 1, 1, 0, 0, 0, ...]/5\)., where we have padded the vector with zeros so that it is the same length as the signal vector.

That is, \(h_0 = 1/5\), \(h_1 = 1/5\), \(h_2 = 1/5\), \(h_3 = 1/5\), \(h_4 = 1/5\), and \(h_5 = 0\), all the way to \(h_n=0\).

Then we can write this as:

\[ f_{\text{filtered}}[n] = \sum_{m=0}^{n} f[n-m] h[m] \]

This operation is called convolution.

It’s messy to compute the convolution directly. But we can do it in the Fourier domain…

The convolution theorem

The convolution of two signals \(f\) and \(h\) is the inverse DFT of the product of the DFTs of \(f\) and \(h\).

In other words, if \(f = \text{ifft}(F)\) and \(h = \text{ifft}(H)\), then the convolution of \(f\) and \(h\) is \(\text{ifft}(F \cdot H)\).

That’s much easier to calculate – only a dot product!

Proof

(fill in if I have time)

Example

Let’s take the signal from before, \(f(x) = \sin(2\pi x) + 0.5 \sin(20\pi x)\), and filter it with the moving average filter \(h = [1, 1, 1, 1, 1]/5\).

from scipy.signal import convolve
N = 200
x = np.linspace(0, 1, N)
y = np.sin(2*np.pi*x) + 10*np.sin(3.3*np.pi*x+.2)+ 10*np.random.normal(0, 0.1, N)
h = np.array([1, 1,1,1,1])/5
y_filtered = convolve(y, h, mode='same')
plt.plot(x, y, label='original')
plt.plot(x, y_filtered, label='filtered')
plt.legend()
plt.show()

Now try it again using the convolution theorem.

Steps:

1. Take the DFT of the signal.
2. Pad the filter with zeros to make it the same length as the signal.
3. Take the DFT of the filter.
4. Multiply the two DFTs.
5. Take the inverse DFT to get the filtered signal.
6. Done!

yf = fft(y)
hf = fft(h, N) # pad h with zeros
y_filtered = ifft(yf*hf)
plt.plot(x, y, label='original')
plt.plot(x, y_filtered, label='filtered')
plt.legend()
plt.show()

Then we can take the inverse DFT to get the filtered signal.

Discrete-time filters in continuous frequency space

Back to your textbook, in Chapter 6: ::: notes Here we will work with discrete signals (using time as the index) and discrete filters, but be interested in their continuous frequency representations.

Imagine we have a signal \(x(t)\) which is defined for \(t\) is the interval \([0, T]\)

We will imagine that \(x(t)\) which is defined over all \(t\) but periodic of period \(T\). So the function has frequency \(f=1 / T\) and angular frequency \(\omega=2 \pi / T\)

A discrete filter is a sequence of complex numbers \(\mathbf{h}=\left\{h_{n}\right\}_{n=-\infty}^{\infty}\).

This will be a finite impulse response (FIR) filter if \(h_n\) is 0 for \(n<0\) or \(n>L\) for some L.

In this case we say \(\mathbf{h}\) has length \(L\) and express it in the form \(\mathbf{h}=\left\{h_{n}\right\}_{n=0}^{L}\).

Then we can take its discrete time Fourier transform (DTFT) as the Fourier series \(X(\mathbf{h})(\zeta)=\sum_{n=-\infty}^{\infty} h_{n} e^{i \zeta}, \zeta \in \mathbb{R}\).

Suppose we have a discrete signal \(x(t)\) in sampling periods \(T_{s}\), i.e., at \(t_{n}=n T_{s}\). Define \(x_{n}=x\left(t_{n}\right)\)

:::

If we use our filter on our signal, we get:

\[ y_{n}=h_{0} x_{n}+h_{1} x_{n-1}+\cdots+h_{L} x_{n-L}, \quad n \in \mathbb{Z} \]

. . .

To understand the effect of the filter, examine its action on a single Fourier mode of the signal, with a possible phase shift \(\phi\):

\(x_{m}(t)=c_{m} e^{i m(\omega t+\phi)}\)

. . .

At individual points in time, this gives

\(x_{m, n}=c_{m} e^{i m \omega\left(n T_{s}+\phi\right)}\)

Thus for fixed \(m\) and arbitrary \(n\), we have

\[ \begin{aligned} y_{m, n} & =h_{0} c_{m} e^{i m \omega\left(n T_{s}+\phi\right)}+h_{1} c_{m} e^{i m \omega\left((n-1) T_{s}+\phi\right)}+\cdots+h_{L} c_{m} e^{i m \omega\left((n-L) T_{s}+\phi\right)} \\ & =\left(h_{0} \cdot 1+h_{1} e^{-i m \omega T_{s}}+\cdots+h_{L} e^{-i L m \omega T_{s}}\right) c_{m} e^{i m \omega\left(n T_{s}+\phi\right)} \\ & =H\left(-m \omega T_{s}\right) x_{m, n} \end{aligned} \]

where

\(H(\zeta)=h_{0}+h_{1} e^{i \zeta}+\cdots+h_{L} e^{i L \zeta}, \quad \zeta \in \mathbb{R}\) is the DTFT of \(\mathbf{h}\).

Gain

\(\left|y_{m, n}\right| \leq\left|H\left(-m \omega T_{s}\right)\right|\left|x_{m, n}\right|\)

So we define gain or attenuation of this transformation as \(G(\zeta)=|H(\zeta)|\)

We define phase rotation as \(\Theta(\zeta)=\theta\), where \(H(\zeta)=|H(\zeta)| e^{i \theta}\)

The filter \(h\) will attenuate the signal at frequencies where \(|H(\zeta)|<1\) and amplify the signal at frequencies where \(|H(\zeta)|>1\). It will phase shift the signal by \(\Theta(\zeta)\).

High-pass and low-pass filters

The FIR filter \(\mathbf{h}=\left\{h_{k}\right\}_{k=0}^{L}\) with discrete time Fourier transform \(H(\zeta)\) is a lowpass filter if \(|H(0)|=1\) and \(|H(\pi)|=0\); \(\mathbf{h}\) is a highpass filter if \(|H(0)|=0\) and \(|H(\pi)|=1\).