This notebook contains code examples from Chapter 7: Discrete Fourier Transform
Copyright 2015 Allen Downey
# Get thinkdsp.py
import os
if not os.path.exists('thinkdsp.py'):
!wget https://github.com/AllenDowney/ThinkDSP/raw/master/code/thinkdsp.py
import numpy as np
import matplotlib.pyplot as plt
from thinkdsp import decorate
PI2 = 2 * np.pi
# suppress scientific notation for small numbers
np.set_printoptions(precision=3, suppress=True)
Here's the definition of ComplexSinusoid, with print statements to display intermediate results.
from thinkdsp import Sinusoid
class ComplexSinusoid(Sinusoid):
"""Represents a complex exponential signal."""
def evaluate(self, ts):
"""Evaluates the signal at the given times.
ts: float array of times
returns: float wave array
"""
print(ts)
phases = PI2 * self.freq * ts + self.offset
print(phases)
ys = self.amp * np.exp(1j * phases)
return ys
Here's an example:
signal = ComplexSinusoid(freq=1, amp=0.6, offset=1)
wave = signal.make_wave(duration=1, framerate=4)
print(wave.ys)
[0. 0.25 0.5 0.75] [1. 2.571 4.142 5.712] [ 0.324+0.505j -0.505+0.324j -0.324-0.505j 0.505-0.324j]
The simplest way to synthesize a mixture of signals is to evaluate the signals and add them up.
from thinkdsp import SumSignal
def synthesize1(amps, freqs, ts):
components = [ComplexSinusoid(freq, amp)
for amp, freq in zip(amps, freqs)]
signal = SumSignal(*components)
ys = signal.evaluate(ts)
return ys
Here's an example that's a mixture of 4 components.
amps = np.array([0.6, 0.25, 0.1, 0.05])
freqs = [100, 200, 300, 400]
framerate = 11025
ts = np.linspace(0, 1, framerate, endpoint=False)
ys = synthesize1(amps, freqs, ts)
print(ys)
[0. 0. 0. ... 1. 1. 1.] [ 0. 0.057 0.114 ... 628.148 628.205 628.262] [0. 0. 0. ... 1. 1. 1.] [ 0. 0.114 0.228 ... 1256.295 1256.409 1256.523] [0. 0. 0. ... 1. 1. 1.] [ 0. 0.171 0.342 ... 1884.443 1884.614 1884.785] [0. 0. 0. ... 1. 1. 1.] [ 0. 0.228 0.456 ... 2512.59 2512.818 2513.046] [1. +0.j 0.995+0.091j 0.979+0.18j ... 0.953-0.267j 0.979-0.18j 0.995-0.091j]
Now we can plot the real and imaginary parts:
n = 500
plt.plot(ts[:n], ys[:n].real)
plt.plot(ts[:n], ys[:n].imag)
decorate(xlabel='Time')
The real part is a mixture of cosines; the imaginary part is a mixture of sines. They contain the same frequency components with the same amplitudes, so they sound the same to us:
from thinkdsp import Wave
wave = Wave(ys.real, framerate)
wave.apodize()
wave.make_audio()
wave = Wave(ys.imag, framerate)
wave.apodize()
wave.make_audio()
We can express the same process using matrix multiplication.
def synthesize2(amps, freqs, ts):
args = np.outer(ts, freqs)
M = np.exp(1j * PI2 * args)
ys = np.dot(M, amps)
return ys
And it should sound the same.
amps = np.array([0.6, 0.25, 0.1, 0.05])
ys = synthesize2(amps, freqs, ts)
print(ys)
[1. +0.j 0.995+0.091j 0.979+0.18j ... 0.953-0.267j 0.979-0.18j 0.995-0.091j]
wave = Wave(ys.real, framerate)
wave.apodize()
wave.make_audio()
To see the effect of a complex amplitude, we can rotate the amplitudes by 1.5 radian:
phi = 1.5
amps2 = amps * np.exp(1j * phi)
ys2 = synthesize2(amps2, freqs, ts)
n = 500
plt.plot(ts[:n], ys.real[:n], label=r'$\phi_0 = 0$')
plt.plot(ts[:n], ys2.real[:n], label=r'$\phi_0 = 1.5$')
decorate(xlabel='Time')
Rotating all components by the same phase offset changes the shape of the waveform because the components have different periods, so the same offset has a different effect on each component.
The simplest way to analyze a signal---that is, find the amplitude for each component---is to create the same matrix we used for synthesis and then solve the system of linear equations.
def analyze1(ys, freqs, ts):
args = np.outer(ts, freqs)
M = np.exp(1j * PI2 * args)
amps = np.linalg.solve(M, ys)
return amps
Using the first 4 values from the wave array, we can recover the amplitudes.
n = len(freqs)
amps2 = analyze1(ys[:n], freqs, ts[:n])
print(amps2)
[0.6 +0.j 0.25-0.j 0.1 +0.j 0.05-0.j]
If we define the freqs
from 0 to N-1 and ts
from 0 to (N-1)/N, we get a unitary matrix.
N = 4
ts = np.arange(N) / N
freqs = np.arange(N)
args = np.outer(ts, freqs)
M = np.exp(1j * PI2 * args)
print(M)
[[ 1.+0.j 1.+0.j 1.+0.j 1.+0.j] [ 1.+0.j 0.+1.j -1.+0.j -0.-1.j] [ 1.+0.j -1.+0.j 1.-0.j -1.+0.j] [ 1.+0.j -0.-1.j -1.+0.j 0.+1.j]]
To check whether a matrix is unitary, we can compute $M^* M$, which should be the identity matrix:
MstarM = M.conj().transpose().dot(M)
print(MstarM.real)
[[ 4. -0. 0. 0.] [-0. 4. -0. 0.] [ 0. -0. 4. -0.] [ 0. 0. -0. 4.]]
The result is actually $4 I$, so in general we have an extra factor of $N$ to deal with, but that's a minor problem.
We can use this result to write a faster version of analyze1
:
def analyze2(ys, freqs, ts):
args = np.outer(ts, freqs)
M = np.exp(1j * PI2 * args)
amps = M.conj().transpose().dot(ys) / N
return amps
N = 4
amps = np.array([0.6, 0.25, 0.1, 0.05])
freqs = np.arange(N)
ts = np.arange(N) / N
ys = synthesize2(amps, freqs, ts)
amps3 = analyze2(ys, freqs, ts)
print(amps3)
[0.6 +0.j 0.25-0.j 0.1 -0.j 0.05-0.j]
Now we can write our own version of DFT:
def synthesis_matrix(N):
ts = np.arange(N) / N
freqs = np.arange(N)
args = np.outer(ts, freqs)
M = np.exp(1j * PI2 * args)
return M
def dft(ys):
N = len(ys)
M = synthesis_matrix(N)
amps = M.conj().transpose().dot(ys)
return amps
And compare it to analyze2:
print(dft(ys))
[2.4+0.j 1. -0.j 0.4-0.j 0.2-0.j]
The result is close to amps * 4
.
We can also compare it to np.fft.fft
. FFT stands for Fast Fourier Transform, which is an even faster implementation of DFT.
print(np.fft.fft(ys))
[2.4+0.j 1. -0.j 0.4-0.j 0.2-0.j]
The inverse DFT is almost the same, except we don't have to transpose $M$ and we have to divide through by $N$.
def idft(amps):
N = len(amps)
M = synthesis_matrix(N)
ys = M.dot(amps) / N
return ys
We can confirm that dft(idft(amps))
yields amps
:
ys = idft(amps)
print(dft(ys))
[0.6 +0.j 0.25-0.j 0.1 -0.j 0.05-0.j]
Let's see what happens when we apply DFT to a real-valued signal.
from thinkdsp import SawtoothSignal
framerate = 10000
signal = SawtoothSignal(freq=500)
wave = signal.make_wave(duration=0.1, framerate=framerate)
wave.make_audio()
wave
is a 500 Hz sawtooth signal sampled at 10 kHz.
hs = dft(wave.ys)
len(wave.ys), len(hs)
(1000, 1000)
hs
is the DFT of this wave, and amps
contains the amplitudes.
amps = np.abs(hs)
plt.plot(amps)
decorate(xlabel='Frequency (unspecified units)', ylabel='DFT')
The DFT assumes that the sampling rate is N per time unit, for an arbitrary time unit. We have to convert to actual units -- seconds -- like this:
N = len(hs)
fs = np.arange(N) * framerate / N
Also, the DFT of a real signal is symmetric, so the right side is redundant. Normally, we only compute and plot the first half:
plt.plot(fs[:N//2+1], amps[:N//2+1])
decorate(xlabel='Frequency (Hz)', ylabel='DFT')
Let's get a better sense for why the DFT of a real signal is symmetric. I'll start by making the inverse DFT matrix for $N=8$.
M = synthesis_matrix(N=8)
And the DFT matrix:
Mstar = M.conj().transpose()
And a triangle wave with 8 elements:
from thinkdsp import TriangleSignal
wave = TriangleSignal(freq=1).make_wave(duration=1, framerate=8)
wave.ys
array([ 1. , 0.5, 0. , -0.5, -1. , -0.5, 0. , 0.5])
Here's what the wave looks like.
wave.plot()
Now let's look at rows 3 and 5 of the DFT matrix:
row3 = Mstar[3, :]
print(row3)
[ 1. -0.j -0.707-0.707j -0. +1.j 0.707-0.707j -1. -0.j 0.707+0.707j 0. -1.j -0.707+0.707j]
row5 = Mstar[5, :]
row5
array([ 1. -0.j , -0.707+0.707j, 0. -1.j , 0.707+0.707j, -1. -0.j , 0.707-0.707j, -0. +1.j , -0.707-0.707j])
They are almost the same, but row5 is the complex conjugate of row3.
def approx_equal(a, b, tol=1e-10):
return np.sum(np.abs(a-b)) < tol
approx_equal(row3, row5.conj())
True
When we multiply the DFT matrix and the wave array, the element with index 3 is:
X3 = row3.dot(wave.ys)
X3
(0.5857864376269053-1.322063213371145e-16j)
And the element with index 5 is:
X5 = row5.dot(wave.ys)
X5
(0.585786437626906-5.534107762827377e-16j)
And they are the same, within floating point error.
abs(X3 - X5)
7.881290833695066e-16
Let's try the same thing with a complex signal:
wave2 = ComplexSinusoid(freq=1).make_wave(duration=1, framerate=8)
plt.plot(wave2.ts, wave2.ys.real)
plt.plot(wave2.ts, wave2.ys.imag)
[0. 0.125 0.25 0.375 0.5 0.625 0.75 0.875] [0. 0.785 1.571 2.356 3.142 3.927 4.712 5.498]
[<matplotlib.lines.Line2D at 0x7f335b48db50>]
Now the elements with indices 3 and 5 are different:
X3 = row3.dot(wave2.ys)
X3
(1.7763568394002505e-15-2.6522243254460964e-16j)
X5 = row5.dot(wave2.ys)
X5
(-2.220446049250313e-16+4.0860826919976536e-16j)
Visually we can confirm that the FFT of the real signal is symmetric:
hs = np.fft.fft(wave.ys)
plt.plot(abs(hs))
[<matplotlib.lines.Line2D at 0x7f335b3f7b50>]
And the FFT of the complex signal is not.
hs = np.fft.fft(wave2.ys)
plt.plot(abs(hs))
[<matplotlib.lines.Line2D at 0x7f335b4b69a0>]
Another way to think about all of this is to evaluate the DFT matrix for different frequencies. Instead of $0$ through $N-1$, let's try $0, 1, 2, 3, 4, -3, -2, -1$.
N = 8
ts = np.arange(N) / N
freqs = np.arange(N)
freqs = [0, 1, 2, 3, 4, -3, -2, -1]
args = np.outer(ts, freqs)
M2 = np.exp(1j * PI2 * args)
approx_equal(M, M2)
True
So you can think of the second half of the DFT as positive frequencies that get aliased (which is how I explained them), or as negative frequencies (which is the more conventional way to explain them). But the DFT doesn't care either way.
The thinkdsp
library provides support for computing the "full" FFT instead of the real FFT.
framerate = 10000
signal = SawtoothSignal(freq=500)
wave = signal.make_wave(duration=0.1, framerate=framerate)
spectrum = wave.make_spectrum(full=True)
spectrum.plot()
decorate(xlabel='Frequency (Hz)', ylabel='DFT')