Osnovne operacije nad signalima

Digitalna obrada signala, Vladimir Petrović

Ovaj notebook će predstaviti osnovne operacije nad diskretnim signalima: skaliranje i kašnjenje. Na osnovu njih ćemo definisati linearne, vremenski nepromenljive sisteme (LTI - linear time-invariant systems). Na kraju ćemo naučiti šta je impulsni odziv LTI sistema i kako izgleda odziv sistema na proizvoljnu pobudu - uvešćemo pojam konvolucije.

# Set to False only if a HTML or PDF is generated
USE_WIDGETS = False

Skaliranje signala i kašnjenje

Skaliranje

Neka je dat diskretan signal $ x[n] $. On se može skalirati konstantom $ a $ ako se svaki odbirak signala pomnoži tom konstantom.

$$ y[n] = a x[n] $$

Signal praktično predstavlja vektor odbiraka i skaliranje je ništa drugo do množenje vektora skalarom.

$$ \mathbf{y} = a \mathbf{x} $$

Evo jednostavnog primera u Pajtonu.

import numpy as np
if USE_WIDGETS:
    %matplotlib widget
else:
    %matplotlib inline
import matplotlib as mpl
mpl.rc('text', usetex = True)
mpl.rc('font', family = 'serif', size = 18)
import matplotlib.pyplot as plt

from matplotlib import animation, rc

N = 21
n = np.arange(N)
x = np.exp(-0.25*n)

a = 0.3
y = a*x;

fig, ax = plt.subplots()
plt.subplots_adjust(bottom=0.15)

# stem x[n]
plt.stem(n, x, label=r'$x[n]$')

# stem y[n]
markerline, stemlines, baseline = plt.stem(n, y, 'tab:red', label=r'$y[n]$')
plt.setp(markerline, 'markerfacecolor', 'tab:red')

ax.set_xlabel('$n$')
ax.set_ylabel('$x[n], y[n]$')

plt.legend(loc="upper right");

Kašnjenje

Neka je dat diskretan signal $ x[n] $. Kašnjenje predstavlja pomeranje odbiraka signala u desno za određeni broj odbiraka $ n_0 $.

$$ y[n] = x[n - n_0] $$

U slučaju da je $ n_0 $ negativan broj, onda nije prirodno ovu operaciju nazvati kašnjenjem, pa je univerzalan naziv "pomeraj u vremenu". Primetimo da iako je dimenzija vremena nestala prilikom odabiranja i pomeraj je praktično promena indeksa odbiraka, pomeraj nazivamo pomeraj u vremenu radi analogije sa kontinualnim signalima.

Evo jednostavnog primera u Pajtonu. Poigrajmo se i sa različitim vrednostima pomeraja $ n_0 $.

import numpy as np
if USE_WIDGETS:
    %matplotlib widget
else:
    %matplotlib inline
import matplotlib as mpl
mpl.rc('text', usetex = True)
mpl.rc('font', family = 'serif', size = 18)
import matplotlib.pyplot as plt

# definition of shift function
def shiftSignal(x, n0):
    N = len(x)
    y = np.zeros(N)
    if n0 < 0:
        y[:N + n0] = x[-n0:N]
    else:
        y[n0:N] = x[:N-n0]
    return y


N = 31
n = np.arange(N) - 10
x = np.exp(-0.25*n)*(n >= 0)

y1 = shiftSignal(x, 5);
y2 = shiftSignal(x, -5);

fig, ax = plt.subplots(figsize = [8, 4])
plt.subplots_adjust(bottom=0.15)

# stem x[n]
plt.stem(n, x, label=r'$x[n]$')

# stem y1[n]
markerline, stemlines, baseline = plt.stem(n, y1, 'tab:red', label=r'$y_1[n]$')
plt.setp(markerline, 'markerfacecolor', 'tab:red')

# stem y2[n]
markerline, stemlines, baseline = plt.stem(n, y2, 'tab:green', label=r'$y_2[n]$')
plt.setp(markerline, 'markerfacecolor', 'tab:green')

ax.set_xlabel('$n$')
ax.set_ylabel('$x[n], y[n]$')

plt.legend(loc="upper right");

Linearan vremenski nepromenljiv sistem (LTI)

Definicija LTI sistema

Sistem koji je linearan i vremenski nepromenljiv mora da zadovolji sledeće osobine:

• Linearnost: Ako su $y_1[n]$ i $y_2[n]$ odzivi sistema na pobude $x_1[n]$ i $x_2[n]$, respektivno, i ako je sistem linearan, odziv sistema na pobudu $x[n] = a_1x_1[n]+a_2x_2[n]$ je $y[n] = a_1y_1[n] + a_2y_2[n]$.
• Vremenska invarijantnost: Ako je $y[n]$ odziv sistema na pobudu $x[n]$, i ako je sistem vremenski nepromeljiv (invarijantan), $y[n-n_0]$ je odziv sistema na pobudu $x[n-n_0]$.

U engleskoj literaturi se linearni, vremenski nepromeljivi sistemi skraćeno nazivaju LTI sistemima, od Linear Time-Invariant systems.

LTI sistemi se mogu okrakterisati impulsnim odzivom, tj. odzivom na pobudu u vidu Dirakovog impulsa. Impulsni odziv se najčešće obeležava sa $ h[n] $.

$$ \delta[n] = \begin{cases} 1 &\text{ za } n = 0 \\ 0 &\text{ inače} \end{cases} $$

Drugi način karakterizacije je korišćenjem diferencne jednačine kojom se zadaje veza između odziva i pobude. Na primer:

$$ y[n] = a_1 y[n-1] + b_0 x[n] + b_1 x[n-1], $$

gde su $a_1$, $b_0$ i $b_1$ konstante.

Konačno treći način da se karakteriše LTI sistem je pomoću prenosne funkcije u Z-domenu, o čemu će biti mnogo više reči kod tema koje se budu bavile projektovanjem filtara. Za primer ranije napisane diferencne jednačine, prenosna funckija bi bila:

$$ H(z) = \frac{ b_0 + b_1 z^{-1} }{ 1 - a_1 z^{-1} }. $$

Odziv LTI sistema na proizvoljnu pobudu

Sada kada smo definisali linearan, vremenski nepromenljiv sistem, postavljamo pitanje kakav je odziv tog sistema na proizvoljnu pobudu. Ako je poznata diferencna jednačina koja opisuje LTI sistem, onda se odbirak po odbirak lako može izračunati odziv. Međutim, ako je poznat samo impulsni odziv sistema, postavljamo pitanje da li je moguće na osnovu njega odrediti odziv na pobudu $ x[n] $. Odgovor sledi u narednim redovima.

Impulsni odziv, podsetimo, predstavlja odziv sistema na jedinični Dirakov impuls. LTI sistem je vremenski nepromenljiv, što znači da će na svom izlazu generisati isti impulsni odziv bez obzira na to u kom trenutku se pojavio Dirakov impuls. Dodatno, sistem je linearan, što znači da će odziv na skaliran Dirakov impuls, biti skaliran impulsni odziv. Takođe, odziv na zbir više skaliranih Dirakovih impulsa biće, zbog osobine linearnosti, zbir više skaliranih impulsnih odziva zakašnjenih onoliko koliko su pojedinačni Dirakovi impulsi bili zakašnjeni. Dakle, ako je proizvoljan signal moguće predstaviti kao niz skaliranih Dirakovih impulsa, onda je moguće, na osnovu osobine linearnosti i vremenske nepromenljivosti sistema, odrediti i odziv.

Dekompozicija signala

Zaista, svaki diskretan signal se može predstaviti kao zbir skaliranih Dirakovih impulsa koji su pomereni u vremenu jedan u odnosu na drugi.

$$ x[n] = \sum_{m=-\infty}^{+\infty} x[m]\delta[n-m]. $$

Razmortimo značenje prethodnog izraza:

• S ozirom na to da je $\delta[n-m]=1$ samo za $n=m$, svi članovi sume su jednaki nuli, osim člana uz $n=m$, pa se prema tome lako dobija identitet $x[n]=x[n]$.
• Skup svih zakašnjenih Dirakovih impulsa $\delta[n-m]$ predstavlja bazis prostora svih diskretnih signala. Dakle, svi diskretni signali se mogu predstaviti linearnom kombinacijom elemenata bazisa. Koeficjenti uz bazisne vektore $\delta[n-m]$ se dobijaju skalarnim proizvodom vektora signala i bazisnih vektora $\sum_{n=-\infty}^{+\infty} x[n]\delta[n-m]=x[m]$. Prema tome, izraz predstavlja dekompoziciju vektora signala na njegove projekcije na bazisne vektore (koeficijente $x[m]$).

Ilustrujmo ovo jednostavnim primerom u Pajtonu.

import numpy as np
import scipy.signal as signal
if USE_WIDGETS:
    %matplotlib widget
else:
    %matplotlib inline
import matplotlib as mpl
mpl.rc('text', usetex = True)
mpl.rc('font', family = 'serif', size = 18)
import matplotlib.pyplot as plt

N = 11
n = np.arange(N)
x = 0.5 + 2*signal.triang(N)

s = np.zeros((N, N))
for k in range(N):
    s[k][k] = x[k]
    
fig, ax = plt.subplots(N + 2, figsize=[6, 14])
for k in range(N):
    label = '$x[%d]\delta[n - %d]$' % (k, k)
    ax[k].stem(n, s[k][:], label = label)
    
    ax[k].set_ylim([0, 3])
    ax[k].get_yaxis().set_ticks([])
    
    if k != N - 1: 
        ax[k].get_xaxis().set_ticks([])
    
    if k < N/2:
        ax[k].legend(loc="upper right");
    else:
        ax[k].legend(loc="upper left");

ax[N].axis('off')

ax[N + 1].stem(x)

ax[N + 1].set_xlabel('$n$')
ax[N + 1].set_ylabel('$x[n]$')
ax[N + 1].set_ylim([0, 3])
ax[N + 1].get_yaxis().set_ticks([])

ax[N + 1].set_title('Suma elementarnih signala');

Konvolucija

Nakon prethodne analize, ostaje da vidimo kakav je odziv na pojedinačne impulse proizvoljnog signala $ x[n] $, a zatim i odziv na ceo signal $ x[n] $. S obzirom na to da analiziramo linearne, vremenski nepromenljive sisteme sa poznatim impulsnim odzivom $ h[n] $, odziv na pobudu $x[m]\delta[n-m]$ je $x[m]h[n-m]$.

Kao što smo videli, signal $x[n]$ se može zapisati kao suma skaliranih Dirakovih impulsa: $$ x[n] = \sum_{m=-\infty}^{+\infty} x[m]\delta[n-m]. $$

Korišćenjem osobine linearnosti, zaključujemo da je odziv LTI sistema na proizvoljnu pobudu $x[n]$ suma skaliranih impulsnih odziva pomerenih u vremenu:

$$ \boxed{\displaystyle{y[n] = \sum_{m=-\infty}^{+\infty} x[m]h[n-m] = (x*h)[n].}} $$

Ova relacija se zove konvolucija (convolution) signala $x[n]$ i $h[n]$.

Ilustrujmo konvoluciju diskretnih signala primerom u Pajtonu.

import numpy as np
import scipy.signal as signal
if USE_WIDGETS:
    %matplotlib widget
else:
    %matplotlib inline
import matplotlib as mpl
mpl.rc('text', usetex = True)
mpl.rc('font', family = 'serif', size = 18)
import matplotlib.pyplot as plt

N = 21
n = np.arange(N)

# Impulse response
h = 0.5*np.exp(-0.25*n)

fig, ax = plt.subplots(figsize=[6, 3])
plt.subplots_adjust(bottom=0.15)

# stem h[n]
plt.stem(n, h, label=r'$h[n]$')

ax.set_xlabel('$n$')
ax.set_ylabel('$h[n]$')

# Signal x[n]
x = np.linspace(1, 5, 5)/5
fig, ax = plt.subplots(figsize=[3, 3])
plt.subplots_adjust(bottom=0.25)
plt.subplots_adjust(left=0.25)

# stem x[n]
plt.stem(x, label=r'$x[n]$')

ax.set_xlabel('$n$')
ax.set_ylabel('$x[n]$')

# Convolution
y = signal.convolve(x, h, mode = 'full')

fig, ax = plt.subplots(figsize=[7, 3])
plt.subplots_adjust(bottom=0.25)

# stem y[n]
n = np.arange(len(y))
plt.stem(n, y, label=r'$y[n]$')

ax.set_xlabel('$n$')
ax.set_ylabel('$y[n]$');

U prethodnom primeru smo iskoristili ugrađenu funkciju iz scipy biblioteke. Pokažimo sada da rezultat konvolucije odgovara sumi zakašnjenih i skaliranih impulsnih odziva, tako što ćemo uraditi dekompoziciju signala $x[n]$ na skalirane Dirakove impulse.

import numpy as np
import scipy.signal as signal
if USE_WIDGETS:
    %matplotlib widget
else:
    %matplotlib inline
import matplotlib as mpl
mpl.rc('text', usetex = True)
mpl.rc('font', family = 'serif', size = 18)
import matplotlib.pyplot as plt

N = 21
n = np.arange(N)

# Impulse response
h = 0.5*np.exp(-0.25*n)

# Signal x[n]
L = 5
x = np.append(np.linspace(1, L, L), np.zeros(N - L))/5

y = np.zeros((N + L - 1,))

# Decomposition
s = np.zeros((L, L))
for k in range(L):
    s[k][k] = x[k]

# Scaled impulse responses
hs = np.zeros((L, len(y)))
for k in range(L):
    hs[k][k:k+N] = x[k]*h    

n = np.arange(N + L - 1)

fig, ax = plt.subplots(L + 2, figsize=[6, 10])
for k in range(L):
    label = '$x[%d]\delta[n - %d]$' % (k, k)
    ax[k].stem(s[k][:], label = label)
    
    hk = hs[k][:]
    label = '$x[%d]h[n - %d]$' % (k, k)
    markerline, stemlines, baseline = ax[k].stem(n, hk, 'tab:red', label=label)
    plt.setp(markerline, 'markerfacecolor', 'tab:red')
    
    ax[k].set_ylim([0, 1.2])
    ax[k].get_yaxis().set_ticks([0, 1])
    
    if k != L - 1: 
        ax[k].get_xaxis().set_ticks([])
    
    ax[k].legend(loc="upper right", fontsize = 14);
    
    y += hk

ax[L].axis('off')

ax[L + 1].stem(y)

ax[L + 1].set_xlabel('$n$')
ax[L + 1].set_ylabel('$y[n]$')
ax[L + 1].set_ylim([0, 1.2])
ax[L + 1].get_yaxis().set_ticks([0, 1])

ax[L + 1].set_title('Suma svih odziva');

diff = y - signal.convolve(x, h, mode = 'full')[:N + L - 1]
print(f'Razlika dvaju konvolucija je: {diff}')
print(f'Srednja kvadratna razlika je: {np.mean(np.square(diff))}')

Razlika dvaju konvolucija je: [0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00
 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00
 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00
 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00
 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00
 0.00000000e+00 1.73472348e-18 0.00000000e+00 0.00000000e+00
 0.00000000e+00]
Srednja kvadratna razlika je: 1.2037062152420224e-37

Srednja kvadratna razlika reda 10^-37 zanemarljiva je i nastaje zbog različitog redosleda po kome se izvršavaju komutativne računske operacije.

Konvolucija kontinualnih signala

Slični zaključci vezani za odziv diskretnih sistema se mogu izvesti i za kontinualne signale i sisteme, a sume će zameniti integrali: $$ y(t) = \int_{\tau=-\infty}^{+\infty} x(\tau)h(t-\tau) = (x*h)(t). $$

Pogledajmo animaciju konvolucije kontinualnih signala.

from IPython.display import Video
display(Video("videos/conv_rects.mp4"))

display(Video("videos/conv.mp4"))

Izvor: DSP Illustations

Vežba

Napišite program koji po definiciji računa konvoluciju dva diskretna signala. Najpre napisati funkciju koja računa konvoluciju dva niza

def conv(x, h):
        ...
        return y

Na kraju uporediti rezultate sa konvolucijama izračunatim korišćenjem numpy i scipy paketa.

import numpy as np
import scipy.signal as signal
if USE_WIDGETS:
    %matplotlib widget
else:
    %matplotlib inline
import matplotlib as mpl
mpl.rc('text', usetex = True)
mpl.rc('font', family = 'serif', size = 18)
import matplotlib.pyplot as plt

def conv(x, h):
    '''Implementation of 1D convolution of two signals by defintion.

    ARGS:
        - x [1d numpy array] : First vector of input samples.
        - h [1d numpy array] : Second vector of input samples.

    RETURNS:
        - y [1d numpy array] : Vector of output samples. 
    '''
    
    Nx = len(x)
    Nh = len(h)
    Ny = Nx + Nh - 1
    y = np.zeros(Ny)
       
    ## OVDE VAS KOD!
    
    return y

Nx = 10
Nh = 15
# test it:
# ramp signals
x = np.exp(-0.25*np.linspace(1, Nx, Nx))
h = signal.triang(Nh)

y = conv(x, h)

fig, ax = plt.subplots()
plt.stem(y)

# calculate using np.convolve
yn = np.convolve(x, h)
fig, ax = plt.subplots()
plt.stem(yn)
# calculate using signal.convolve
ys = signal.convolve(x, h)
fig, ax = plt.subplots()
plt.stem(ys)
# compare results
print(y - yn)
print(y - ys)

[-0.0973501  -0.27051653 -0.50272878 -0.78092596 -1.09493624 -1.43683779
 -1.80046108 -2.18100128 -2.38001619 -2.43765906 -2.37721029 -2.22479166
 -2.00074683 -1.72091947 -1.39764861 -1.04054394 -0.75443956 -0.52363026
 -0.35186677 -0.22608822 -0.13612277 -0.07404859 -0.03369615 -0.01026062]
[-0.0973501  -0.27051653 -0.50272878 -0.78092596 -1.09493624 -1.43683779
 -1.80046108 -2.18100128 -2.38001619 -2.43765906 -2.37721029 -2.22479166
 -2.00074683 -1.72091947 -1.39764861 -1.04054394 -0.75443956 -0.52363026
 -0.35186677 -0.22608822 -0.13612277 -0.07404859 -0.03369615 -0.01026062]

Analizirajmo izraz za konvoluciju malo detaljnije. Sa tim ciljem, ponovićemo izraz:

$$ y[n] = \sum_{m=-\infty}^{+\infty} x[m]h[n-m] = (x*h)[n] $$

Primetimo da n-ti odbirak konvolucije predstavlja sumu dve pomnožene sekvence: jedne originalne sekvence ($x[m]$) i druge ($h[m]$) koja je invertovana u vremenu ($h[-m]$) i pomerena za n ($h[n-m]$). Ilustrujmo ove pomeraje na primeru:

import numpy as np
import scipy.signal as signal
if USE_WIDGETS:
    %matplotlib widget
else:
    %matplotlib inline
import matplotlib as mpl
mpl.rc('text', usetex = True)
mpl.rc('font', family = 'serif', size = 18)
import matplotlib.pyplot as plt

Nx = 15
Nh = 10

# ramp signals
h = np.exp(-0.25*np.linspace(1, Nh, Nh))
x = signal.triang(Nx)

fig, ax = plt.subplots(2, figsize = [6, 4])
plt.subplots_adjust(bottom=0.15)
ax[0].stem(x)
ax[0].get_xaxis().set_ticks([])
ax[0].set_ylabel('$x[n]$')

ax[1].stem(h)
ax[1].set_ylabel('$h[n]$')
ax[1].set_xlabel('$n$')

nx = np.arange(len(x))

n = np.arange(-Nh, Nx+7, 1)
x = np.append(np.zeros(Nh), x)
x = np.append(x, np.zeros(7))

fig, ax = plt.subplots(6, figsize = [6, 10])
plt.subplots_adjust(hspace = 0.5)

ax[0].stem(n, x)
ax[0].set_ylabel('$x[n]$')
ax[0].axvline(x=-0.5, linestyle=':', color='tab:green')
ax[0].axvline(x=14.5, linestyle=':', color='tab:green')

for i in range(5):
    hi = np.append(np.zeros(i*5 + 1), h[-1::-1])
    hi = np.append(hi, np.zeros(len(n)-i*5-1-Nh))
    ax[i+1].stem(n, hi)
    ylabel = '$h[%d - n]$' % (i*5)
    ax[i+1].set_ylabel(ylabel)
    ax[i+1].axvline(x=-0.5, linestyle=':', color='tab:green')
    ax[i+1].axvline(x=14.5, linestyle=':', color='tab:green')

ax[5].set_xlabel('$n$');

Na osnovu prethodne analize, možemo napisati kod lepše. Funkcija napisana bez ugnježdenih for petlji je data u narednoj ćeliji.

import numpy as np
import scipy.signal as signal
if USE_WIDGETS:
    %matplotlib widget
else:
    %matplotlib inline
import matplotlib as mpl
mpl.rc('text', usetex = True)
mpl.rc('font', family = 'serif', size = 18)
import matplotlib.pyplot as plt

def conv(x, h):
    '''Implementation of 1D convolution of two signals by defintion.

    ARGS:
        - x [1d numpy array] : First vector of input samples.
        - h [1d numpy array] : Second vector of input samples.

    RETURNS:
        - y [1d numpy array] : Vector of output samples. 
    '''
    
    Nx = len(x)
    Nh = len(h)
    Ny = Nx + Nh - 1
    
    # Replace so that h has smaller length
    if Nx < Nh:
        p = h
        h = x
        x = p
        Nx = len(x)
        Nh = len(h)
    
    # init
    y = np.zeros(Ny)
    
    # n < Nh
    for n in np.arange(0, Nh):
        y[n] = np.sum(x[0:n + 1] * h[n::-1])
    # n >= Nh and n <= Nx
    for n in np.arange(Nh, Nx):
        y[n] = np.sum(x[n + 1 - Nh:n + 1] * h[-1::-1])
    # n >= Nx
    for n in np.arange(Nx, Ny):
        y[n] = np.sum(x[n + 1 - Nh:Nx] * h[-1:n-Nx:-1])
    
    return y

Nx = 10
Nh = 15
# test it:
# ramp signals
x = np.exp(-0.25*np.linspace(1, Nx, Nx))
h = signal.triang(Nh)

y = conv(x, h)

fig, ax = plt.subplots()
plt.stem(y)

# calculate using np.convolve
yn = np.convolve(x, h)
fig, ax = plt.subplots()
plt.stem(yn)
# calculate using signal.convolve
ys = signal.convolve(x, h)
fig, ax = plt.subplots()
plt.stem(ys)
# compare results
print(y - yn)
print(y - ys)

[ 0.00000000e+00  0.00000000e+00  0.00000000e+00  0.00000000e+00
  0.00000000e+00  0.00000000e+00  0.00000000e+00  0.00000000e+00
 -4.44089210e-16  0.00000000e+00  0.00000000e+00  0.00000000e+00
  0.00000000e+00  0.00000000e+00  0.00000000e+00 -2.22044605e-16
  1.11022302e-16  0.00000000e+00  0.00000000e+00  0.00000000e+00
  0.00000000e+00  0.00000000e+00  0.00000000e+00  0.00000000e+00]
[ 0.00000000e+00  0.00000000e+00  0.00000000e+00  0.00000000e+00
  0.00000000e+00  0.00000000e+00  0.00000000e+00  0.00000000e+00
 -4.44089210e-16  0.00000000e+00  0.00000000e+00  0.00000000e+00
  0.00000000e+00  0.00000000e+00  0.00000000e+00 -2.22044605e-16
  1.11022302e-16  0.00000000e+00  0.00000000e+00  0.00000000e+00
  0.00000000e+00  0.00000000e+00  0.00000000e+00  0.00000000e+00]

Primena konvolucije u konvolucionom reverbu

Reverberacija je zvučni efekat koji nastaje kao rezultat interakcije zvuka sa okruženjem slušaoca. Najlakše ga je objasniti kao superpoziciju zvučnih signala koji nastaju odbijanjem od raznih površina unutar prostorije. Često se ovaj efekat postiže i veštački, obradom signala, čime se postiže da signali snimljeni u anehoičnoj sobi studija zvuče kao da su snimljeni u sportskoj hali, amfiteatru pozorišta i sl. Zanimljivo je da je za postizanje ovog efekta dovoljno uraditi jednu konvoluciju. Ako su na raspolaganju impulsni odzivi prostorija, konvolucijom originalnog signala i impulsnih odziva dobijaju se signali koji zvuče kao da su snimljeni u tim prostorijama. Impulsni odzivi se najčešće snimaju tako što se proizvede impulsni zvuk (pljesak, pucanj iz pištolja, udar filmske klapne) koji predstavljaju impulsnu pobudu.

Sledeći Pajton kod ilustruje konvolucioni reverb. Za više primera posetiti sajt.

import numpy as np
import scipy.signal as signal
if USE_WIDGETS:
    %matplotlib widget
else:
    %matplotlib inline
import matplotlib as mpl
mpl.rc('text', usetex = True)
mpl.rc('font', family = 'serif', size = 18)
import matplotlib.pyplot as plt

import IPython
from IPython.display import Markdown, Image
from scipy.io import wavfile

def plotTimeSequence(fs, x, xlabel, ylabel):
    fig, ax = plt.subplots(figsize=[10, 3])
    plt.subplots_adjust(bottom=0.25)

    # stem x[n]
    n = np.arange(len(x))
    t = n/fs
    plt.plot(t, x, label=ylabel)

    ax.set_xlabel(xlabel)
    ax.set_ylabel(ylabel);
    if not USE_WIDGETS:
        plt.show()

             
fs, x = wavfile.read('sounds/singing.wav')
x = x.astype('float64') 

display(Markdown("Originalni signal snimljen u anehoičnoj sobi"))
IPython.display.display(IPython.display.Audio(x, rate = fs))

plotTimeSequence(fs, x, '$t (\mathrm{s})$', '$x(t)$')

#### IR 1 ############################################################################
fs, h = wavfile.read('sounds/ir1_1st_baptist_nashville_balcony.wav')
h = h.astype('float64') 

display(Markdown("__1.__ Impulsni odziv unutrašnjosti katedrale sa slike"))
IPython.display.display(Image(filename='images/ir1_church_inside.jpg'))
IPython.display.display(IPython.display.Audio(h, rate = fs))

plotTimeSequence(fs, h, '$t (\mathrm{s})$', '$h(t)$')

y = signal.convolve(x, h)

display(Markdown("Konvolucija originalnog signala i impulsnog odziva"))
IPython.display.display(IPython.display.Audio(y, rate = fs))

plotTimeSequence(fs, y, '$t (\mathrm{s})$', '$y(t)$')

#### IR 2 ############################################################################
fs, h = wavfile.read('sounds/ir2_central_hall_university_york.wav')
h = h.astype('float64') 

display(Markdown("__2.__ Impulsni odziv unutrašnjosti amfiteatra sa slike"))
IPython.display.display(Image(filename='images/ir2_central-hall-university-york.jpg'))
IPython.display.display(IPython.display.Audio(h, rate = fs))

plotTimeSequence(fs, h, '$t (\mathrm{s})$', '$h(t)$')

y = signal.convolve(x, h)

display(Markdown("Konvolucija originalnog signala i impulsnog odziva"))
IPython.display.display(IPython.display.Audio(y, rate = fs))

plotTimeSequence(fs, y, '$t (\mathrm{s})$', '$y(t)$')

#### IR 3 ############################################################################
fs, h = wavfile.read('sounds/ir3_hoffmann-lime-kiln.wav')
h = h.astype('float64') 

display(Markdown("__3.__ Impulsni odziv unutrašnjosti tunela sa slike"))
IPython.display.display(Image(filename='images/ir3_hoffmann-lime-kiln.jpg'))
IPython.display.display(IPython.display.Audio(h, rate = fs))

plotTimeSequence(fs, h, '$t (\mathrm{s})$', '$h(t)$')

y = signal.convolve(x, h)

display(Markdown("Konvolucija originalnog signala i impulsnog odziva"))
IPython.display.display(IPython.display.Audio(y, rate = fs))

plotTimeSequence(fs, y, '$t (\mathrm{s})$', '$y(t)$')

Originalni signal snimljen u anehoičnoj sobi

1. Impulsni odziv unutrašnjosti katedrale sa slike

Konvolucija originalnog signala i impulsnog odziva

2. Impulsni odziv unutrašnjosti amfiteatra sa slike

Konvolucija originalnog signala i impulsnog odziva

3. Impulsni odziv unutrašnjosti tunela sa slike

Konvolucija originalnog signala i impulsnog odziva

Literatura

• J.-F. Bercher, A Journey in Signal Processing with Jupyter, dostupno na linku
• M. Popović, Signali i sistemi, Akademska misao, 2006.
• Interaktivni web sajt DSP Illustrations
• OpenAir baza impulsnih odziva za reverb