Discrete Signals


By Prof. Seungchul Lee
http://iai.postech.ac.kr/
Industrial AI Lab at POSTECH

Table of Contents

1. Discrete Time Signals¶

A signal $x[n]$ is a function that maps an independent variable to a dependent variable.

In this course, we will focus on discrete-time signals $x[n]$:

  • Independent variable is an integer: $n \in \mathbb{Z}$
  • Dependent variable is a real or complex number: $x[n] \in \mathbb{R}$ or $\mathbb{C}$




1.1. Plot Real Signals¶

  • plot for continuous signals in Python

  • stem for discrete signals in Python

$$x(t) = \sin(2\pi t)$$
In [1]:
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
In [2]:
t = np.linspace(0,2,200)
x = np.sin(2*np.pi*t)

# plot
plt.figure(figsize = (10,4))
plt.plot(t, x)
plt.xlim([np.min(t), np.max(t)])
plt.xlabel('t [sec]')
plt.ylabel('x(t)')
plt.show()
$$x[n] = \sin \left(\frac{2 \pi}{N}n \right)$$
In [3]:
N = 20
n = np.arange(0, N)
x = np.sin(2*np.pi/N*n)

# plot
plt.figure(figsize = (10, 8))
plt.subplot(2,1,1)
plt.plot(n, x, 'o')
plt.xlim([np.min(n), np.max(n)])
plt.title('plot')

plt.subplot(2,1,2)
plt.stem(n, x)
plt.xlim([np.min(n), np.max(n)])
plt.title('stem')
plt.show()

1.2. Signal Sounds¶

In [4]:
from scipy.io import wavfile
from IPython.display import Audio

fs, data = wavfile.read('./data_files/hamlet.wav')
In [5]:
# signal: 'ALas, Poor Yorick!'

plt.figure(figsize = (10, 4))
plt.plot(data)
plt.show()
In [6]:
Audio('./data_files/hamlet.wav')
Out[6]:
$$x[n] = \cos \left(\frac{2 \pi}{N}k n \right)$$
In [7]:
# cosine wave

N = 44100
n = np.arange(0, 2*N)
k = 40 # frequancy
coswave = np.cos(2*np.pi/N*k*n)

# plot
plt.figure(figsize = (10, 4))
plt.plot(n/fs, coswave)
plt.show()
In [8]:
Audio(coswave, rate = fs)
Out[8]:

Chirp

$$x[n] = 0.3 \cos \left(\frac{2 \pi}{3000000}n^2 \right)$$
In [9]:
# chirp (https://en.wikipedia.org/wiki/Chirp)
n = np.arange(0,44100*2)
chirp = 0.3*np.cos(2*np.pi/3000000*n**2)

# plot
plt.figure(figsize = (10, 4))
plt.plot(n/fs, chirp)
plt.xlim([0, 0.4])
plt.show()
In [10]:
Audio(chirp, rate = 44100)
Out[10]:

White noise

In [11]:
# white gaussian noise
N = 44100
noise = 0.1*np.random.randn(2*N,1)

# plot
plt.figure(figsize = (10, 4))
plt.plot(n/fs, noise)
plt.xlim([0, 0.4])
plt.show()
In [12]:
Audio(noise.T, rate = fs)
Out[12]:

2. Complex Sinusoid¶

2.1. Complex Numbers¶

$$z_1 = a_1 + b_1i, \quad \vec{z}_1 = \begin{bmatrix} a_1 \\ b_1 \end{bmatrix} $$


$$z_2 = a_2 + b_2i, \quad \vec{z}_2 = \begin{bmatrix} a_2 \\ b_2 \end{bmatrix} $$

  • Addition
$$ \begin{align*} z &= z_1 + z_2 = (a_1 + a_2) + (b_1 + b_2)i \\ \\ \vec{z} &= \vec{z}_1 + \vec{z}_2 = \begin{bmatrix} a_1 \\ b_1 \end{bmatrix}+ \begin{bmatrix} a_2 \\ b_2 \end{bmatrix} = \begin{bmatrix} a_1 + a_2 \\ b_1 + b_2 \end{bmatrix} \end{align*} $$




  • Multiplication
$$\begin{cases} z_1=r_1e^{i\theta_1}\\ \\ z_2=r_2e^{i\theta_2} \end{cases} \quad\Longrightarrow\quad \begin{cases} z_1\centerdot z_2=r_1r_2e^{i(\theta_1+\theta_2)} \\ \\ {z_1 \over z_2}={r_1 \over r_2}e^{i(\theta_1-\theta_2)} \end{cases} $$



  • Euler's Formula:
$$ e^{i\theta} = \text{cos}\theta + i\text{sin}{\theta}$$
  • Complex number in complex exponential
$$ \begin{align*} \vec{z} &= r \, \text{cos}\theta + i\,r\,\text{sin}\theta \\ &= r\,(\text{cos}\theta + i\text{sin}\theta) \\ &= re^{i\theta}\\ \\ r &: \text{magnitude (length)} \\ \theta &: \text{phase (angle)} \end{align*} $$




  • When $x[n] \in \mathbb{C}$, we can use two signal plots


$$ \begin{align*} x[n] &= \text{Re}\{x[n]\} + j \text{Im} \{x[n]\} & \text{Rectangular form}\\\\ x[n] &= \lvert x[n] \rvert e^{j \angle x[n]} & \text{Polar form} \end{align*} $$


  • For $e^{j \omega t}$




2.2. Geometrical Meaning of $e^{i\theta}$¶

  • $e^{i\theta}: \text{point on the unit circle with angle of } \theta$



  • $\theta = \omega t$
  • $e^{i\omega t}: \text{rotating on an unit circle with angular velocity of } \omega$.



  • Question: what is the physical meaning of $e^{-i\omega t}$?
  • Frequency $\omega$ determines rotation speed and direction of a circular motion
    • $\omega > 0 \implies$ counterclockwise rotation
    • $\omega < 0 \implies$ clockwise rotation

2.3. Sinusoidal Functions from Circular Motions¶



  • A complex sinusoid is a helix
    • Real part ($\cos$ term) is the projection onto the $\text{Re}$ axis.
    • Imaginary part ($\sin$ term) is the projection onto the $\text{Im}$ axis.



$$\text{cos}\;\omega t = \large{e^{i\omega t } + e^{-i\omega t} \over 2}$$

3. Discrete Sinusoids¶

Discrete Sinusoids


$$ \begin{align*}x[n] &= A\cos(\omega_0 n + \phi)\quad \text{or}\\ \\ x[n] &= Ae^{j (\omega_0 n + \phi)} \qquad \text{where} \;\; \omega_0 = \frac{2\pi}{N}k\end{align*}$$


Key question: which one is a higher frequency?

$$\omega_0 = \pi ~~~ \text{or} ~~~ \omega_0 = \frac{3\pi}{2}$$
In [13]:
n = np.arange(0, 8) # 0,1,2,3 ~ 7
x1 = np.cos(np.pi*n)
x2 = np.cos(3/2*np.pi*n)

# plot
plt.figure(figsize = (10, 8))
plt.subplot(2,1,1)
plt.stem(n, x1)
plt.ylim([-1.5, 1.5])
plt.title('$x_1$')

plt.subplot(2,1,2)
plt.stem(n, x2)
plt.ylim([-1.5, 1.5])
plt.title('$x_2$')
plt.xlabel('n')
plt.show()
In [14]:
n = np.arange(0, 8)
x1 = np.cos(np.pi*n)
x2 = np.cos(3/2*np.pi*n)

n_real = np.linspace(0,7,200)
x1_real = np.cos(np.pi*n_real)
x2_real = np.cos(3/2*np.pi*n_real)

# plot
plt.figure(figsize = (10, 8))
plt.subplot(2,1,1)
plt.stem(n, x1, markerfmt='ro')
plt.plot(n_real,x1_real,'--')
plt.ylim([-1.5, 1.5])
plt.title('$x_1$')

plt.subplot(2,1,2)
plt.stem(n, x2, markerfmt='ro')
plt.plot(n_real, x2_real,'--')
plt.ylim([-1.5, 1.5])
plt.title('$x_2$')
plt.xlabel('n')
plt.show()

3.1. Frequency in Discrete Sinusoids¶

$$0 \; \rightarrow \; \frac{1}{8}\pi \; \rightarrow \; \frac{2}{8}\pi \; \rightarrow \; \frac{8}{8}\pi$$
In [15]:
N = 32

n = np.arange(0, 32)
x1 = np.cos(0*np.pi*n)
x2 = np.cos(1/8*np.pi*n)
x3 = np.cos(2/8*np.pi*n)
x4 = np.cos(1*np.pi*n)

# plot
plt.figure(figsize = (10, 8))
plt.subplot(4,1,1)
plt.stem(n, x1)
plt.xlim([0, 31])
plt.ylim([-1.5, 1.5])
plt.ylabel('0')

plt.subplot(4,1,2)
plt.stem(n, x2)
plt.xlim([0, 31])
plt.ylim([-1.5, 1.5])
plt.ylabel('1/8 $\pi$')

plt.subplot(4,1,3)
plt.stem(n, x3)
plt.xlim([0, 31])
plt.ylim([-1.5, 1.5])
plt.ylabel('2/8 $\pi$')

plt.subplot(4,1,4)
plt.stem(n, x4)
plt.xlim([0, 31])
plt.ylim([-1.5, 1.5])
plt.ylabel('8/8 $\pi$')
plt.show()
$$2\pi \; \rightarrow \; \frac{15}{8}\pi \; \rightarrow \; \frac{14}{8}\pi \; \rightarrow \; \frac{8}{8}\pi$$
In [16]:
n = np.arange(0, 32)
x5 = np.cos(2*np.pi*n)
x6 = np.cos(15/8*np.pi*n)
x7 = np.cos(14/8*np.pi*n)
x8 = np.cos(1*np.pi*n)

# plot
plt.figure(figsize = (10, 8))
plt.subplot(4,1,1)
plt.stem(n, x5)
plt.xlim([0, 31])
plt.ylim([-1.5, 1.5])
plt.ylabel('2')

plt.subplot(4,1,2)
plt.stem(n, x6)
plt.xlim([0, 31])
plt.ylim([-1.5, 1.5])
plt.ylabel('15/8 $\pi$')

plt.subplot(4,1,3)
plt.stem(n, x7)
plt.xlim([0, 31])
plt.ylim([-1.5, 1.5])
plt.ylabel('14/8 $\pi$')

plt.subplot(4,1,4)
plt.stem(n, x8)
plt.xlim([0, 31])
plt.ylim([-1.5, 1.5])
plt.ylabel('1 $\pi$')
plt.show()

3.2. Aliasing¶

Consider two sinusoids with two different frequencies


$$ \begin{align*} \omega \quad &\implies \quad x_1[n] = e^{j(\omega n + \phi)}\\ \omega + 2 \pi \quad &\implies \quad x_2[n]= e^{j((\omega+2\pi) n + \phi)}=e^{j(\omega n + \phi)}e^{j2\pi n} \end{align*} $$
  • But note that
$$x_2[n] = x_1[n]$$
  • The signal $x_1$ and $x_2$ have different frequencies but are identical
  • We say that $x_1$ and $x_2$ are aliases; this phenomenon is called aliasing

Aliasing occurs when a system is measured at an insufficient sampling rate. It is perhaps best explained through an example.

In [17]:
t = np.linspace(0, 10*2*np.pi, 300)
x = np.sin(t)
y = np.cos(t)

# plot
plt.figure(figsize = (10, 8))
plt.subplot(2,1,1)
plt.plot(t, x, 'o')
plt.xlim([np.min(t), np.max(t)])
plt.ylabel('sin(t)')

plt.subplot(2,1,2)
plt.plot(t, y ,'o')
plt.ylabel('cos(t)')
plt.xlim([np.min(t), np.max(t)])
plt.show()
In [18]:
ts = np.linspace(0, 10*2*np.pi, 12)
xs = np.sin(ts)
ys = np.cos(ts)

# plot
plt.figure(figsize = (10, 8))
plt.subplot(2,1,1)
plt.plot(t, x)
plt.plot(ts, xs, 'o--')
plt.xlim([np.min(t), np.max(t)])
plt.ylabel('sin(t)')

plt.subplot(2,1,2)
plt.plot(t, y)
plt.plot(ts, ys, 'o--')
plt.ylabel('cos(t)')
plt.xlim([np.min(t), np.max(t)])
plt.show()
In [19]:
ts = np.linspace(0, 10*2*np.pi, 11)
xs = np.sin(ts)
ys = np.cos(ts)

# plot
plt.figure(figsize = (10, 8))
plt.subplot(2,1,1)
plt.plot(t, x)
plt.plot(ts, xs, 'o--')
plt.xlim([np.min(t), np.max(t)])
plt.ylabel('sin(t)')

plt.subplot(2,1,2)
plt.plot(t, y)
plt.plot(ts, ys, 'o--')
plt.ylabel('cos(t)')
plt.xlim([np.min(t), np.max(t)])
plt.show()
In [20]:
ts = np.linspace(0, 10*2*np.pi, 21)
xs = np.sin(ts)
ys = np.cos(ts)

# plot
plt.figure(figsize = (10, 8))
plt.subplot(2,1,1)
# plt.plot(t, x)
plt.plot(ts, xs, 'o--')
plt.xlim([np.min(t), np.max(t)])
plt.ylim([-1,1])
plt.ylabel('sin(t)')

plt.subplot(2,1,2)
# plt.plot(t, y)
plt.plot(ts, ys, 'o--')
plt.ylabel('cos(t)')
plt.xlim([np.min(t), np.max(t)])
plt.show()
In [21]:
%%html
<center><iframe 
width="560" height="315" src="https://www.youtube.com/embed/jHS9JGkEOmA?rel=0" frameborder="0" allowfullscreen>
</iframe></center>

4. Fourier Series¶

4.1. Delta Dirac Function¶

$$ \sum_{n = 1,2,3,\cdots}\cos(n \theta) = \cos \theta + \cos 2\theta + \cos 3\theta + \cdots $$
In [22]:
from scipy.fftpack import fft, ifft, fftshift

theta = np.linspace(0, 6*np.pi, 2**8)
x1 = np.cos(theta)

# plot
plt.figure(figsize = (10, 4))
plt.plot(theta, x1)
plt.ylim([-1.1, 1.1])
plt.xlim([np.min(theta), np.max(theta)])
plt.xlabel(r'$\theta$', fontsize = 15)
plt.show()
In [23]:
x2 = x1 + np.cos(2*theta)
x3 = x2 + np.cos(3*theta)
x4 = x3 + np.cos(4*theta)

# plot
plt.figure(figsize = (10,12))
plt.subplot(3,1,1)
plt.plot(theta, x2/2)
plt.xlim([np.min(theta), np.max(theta)])
plt.ylim([-1.1, 1.1])

plt.subplot(3,1,2)
plt.plot(theta, x3/3)
plt.xlim([np.min(theta), np.max(theta)])
plt.ylim([-1.1, 1.1])

plt.subplot(3,1,3)
plt.plot(theta, x4/4)
plt.xlim([np.min(theta), np.max(theta)])
plt.ylim([-1.1, 1.1])
plt.xlabel(r'$\theta$', fontsize = 15)
plt.show()
In [24]:
x = np.zeros(theta.shape)

N = 11000
for n in range(N):
    x = x + np.cos(n*theta)

# plot
plt.figure(figsize = (10, 4))
plt.plot(theta, x/N)
plt.xlabel(r'$\theta$', fontsize = 15)
plt.ylim([-1.1, 1.1])
plt.show()

4.2. Square Wave¶


$$ \sum_{n = 1,3,5,\cdots}\frac{4}{\pi n} \sin (n\theta) = \frac{4}{\pi} \sin \theta + \frac{4}{3 \pi} \sin 3\theta + \frac{4}{5 \pi} \sin 5\theta + \cdots $$


In [25]:
N = 2**8
theta = np.linspace(0, 2*np.pi, N)
x = np.zeros(theta.shape)

for n in range(1, 5000, 2):
    x = x + 4/(np.pi*n)*np.sin(n*theta)

# plot    
plt.figure(figsize = (10, 4))
plt.plot(theta/np.pi, x)
plt.xlabel(r'$\theta/\pi$', fontsize = 15)
plt.show()

4.3. How to Decompose a Signal into Fundamental Harmonic Sinusoids¶

  • Signal representation by harmonic sinusoids

  • Questions

    • how to decompose it to harmonic sinusoids
    • how to find the coefficients of harmonic sinusoids
    • how to find the frequency components
  • DFT or FFT can do such a job !!!