Digitalna obrada signala, Vladimir Petrović
U ovom notebook-u ćemo implementirati jednu realizaciju FIR filtra i uvesti osnovne alate za analizu implementacija sa konačnom tačnošću izračunavanja u aritmetici sa fiksnim zarezom.
USE_WIDGETS = False
def importEssentialLibs(USE_WIDGETS):
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
import scipy.signal as signal
return np, mpl, plt, signal
U ovim primerim ćemo uraditi najjednostavniji primer realizacije filtara - direktnu realizaciju FIR filtra. Ona se može realizovati pomoću elemenata za kašnjenje kao na slici ispod.
U narednoj ćeliji ćemo definisati funkciju koja implementira direktnu realizaciju FIR filtra.
np, mpl, plt, signal = importEssentialLibs(USE_WIDGETS)
def firDirect(h, x, state = None):
M = len(h)
if state == None:
delayLine = np.zeros(M)
else:
delayLine = state
y = np.zeros(len(x))
for n in range(len(x)):
delayLine[1:] = delayLine[:M-1]
delayLine[0] = x[n]
for m in range(M):
y[n] += h[m]*delayLine[m]
return y
Iskoristimo ovu funkciju za filtriranje signala:
h = np.array([-0.0136, -0.0139, 0.0254, 0.0523, -0.0124, -0.0880, 0.0252, 0.3169, \
0.4807, 0.3169, 0.0252, -0.0880, -0.0124, 0.0523, 0.0254, -0.0139, -0.0136])
n = np.arange(len(h))
fig = plt.figure(figsize = (10,4))
plt.subplots_adjust(bottom=0.25, wspace = 0.4)
# Stem filter
ax1 = fig.add_subplot(1,2,1)
ax1.stem(n, h)
ax1.set_xlabel(r'$n$')
ax1.set_ylabel(r'$h[n]$')
# Definisanje ulaznog signala, zbir dve sinusoide
F1 = 0.1
F2 = 0.43
N = 256
n = np.arange(N)
x = np.sin(2*np.pi*F1*n) + np.sin(2*np.pi*F2*n)
y = firDirect(h, x)
# Provera ispravnosti funkcije
yScipy = signal.lfilter(h, 1, x)
print(sum(abs(yScipy - y)))
# Stem input
ax2 = fig.add_subplot(1,2,2)
ax2.plot(n, x, label = 'Ulaz')
ax2.plot(n, y, label = 'Izlaz')
ax2.set_xlabel(r'$n$')
ax2.set_ylabel(r'$x[n], y[n]$')
ax2.legend(loc = 'upper right');
1.9816613627821056e-14
Filtri se u hardveru realizuju sa konačnom tačnošću. Ne samo u hardveru, vrlo često se za obradu signala koriste posebni procesori za digitalnu obradu signala (DSP procesori) koji mogu biti optimizovani za malu potrošnju, pa u tom slučaju nemaju floating point jedinicu. Prilikom dizajna bilo kog hardverskog sistema za obradu signala uobičajena je praksa da se algoritam najpre napiše u predstavi sa pokretnim zarezom, a zatim se sve vrednosti prevedu u predstavu sa fiksnim zarezom i algoritam se optimizuje tako da daje zadovoljavajuće performanse sa što manje bita. Tek tada se prelazi na hardversku realizaciju.
Za predstavu brojeva sa fiksnim zarezom u Pajtonu koristićemo paket fxpmath. Pre korišćenja ga, naravno, moramo instalirati iz komandne linije komandom pip install fxpmath, a možemo to uraditi i direkto iz notebook-a. Sažeta dokumentacija za ovaj paket, data je ovde.
pip install fxpmath
Requirement already satisfied: fxpmath in /home/vladimir/oe3dos/venv/lib/python3.9/site-packages (0.4.5)
Requirement already satisfied: numpy in /home/vladimir/oe3dos/venv/lib/python3.9/site-packages (from fxpmath) (1.21.2)
WARNING: You are using pip version 21.2.3; however, version 21.3.1 is available.
You should consider upgrading via the '/home/vladimir/oe3dos/venv/bin/python3.9 -m pip install --upgrade pip' command.
Note: you may need to restart the kernel to use updated packages.
Osnovna klasa iz ovog paketa je Fxp. Kreiranjem objekta ove klase zadajemo vrednost broja ili niza koji konvertujemo u predstavu sa fiksnim zarezom i format izlaznog rezultata. Format podrazumeva definiciju da li su brojevi označeni iili neoznačeni, ukupan broj bita za predstavu broja, zatim broj bita za predstavu razlomljenog dela broja (desno od decimalnog zareza), ponašanje kod prekoračenja ('saturate' ili 'wrap') i način zaokruživanja ('floor', 'trunc', 'around', 'ceil', 'fix'). Podrazumevano, brojevi su označeni. Format se može zadati i stringom, a čitanjem polja dtype možemo dobiti predstavu formata u obliku stringa.
Podaci različitih formata mogu ulaziti u izračunavanja. Pre izračunavanja je neophodno konvertovati podatke u željeni format rezultata. Konverzija se najjednostavnije izvodi tako što se ponovo pozvoe konstruktor klase Fxp sa novim formatom, ali je moguće menjati i pojedinačne opcije.
Množenje će podrazumevano dati rezultat čije su bitske širine jednake zbiru bitskih širina činilaca, sabiranje će uvek biti urađeno tako da se ne dešava prekoračenje. Ako rezultate želimo da zaokružimo na manje bitske širine, to je najlakše uraditi konverzijom.
Treba još napomenuti da su Pajtonove biblioteke za predstavu brojeva sa fiksnim zarezom po pravilu veoma spore. Zbog toga se u profesionalne svrhe za evaluaciju performansi algoritama mnogo češće koristi C, ali to prevazilazi okvire ovog kursa.
Implementirajmo direktnu realizaciju filtra u predstavi sa fiksnim zarezom. U narednoj funkciji se smatra da množači na svom izlazu daju podrazumevane bitske širine (zbir bitskih širina ulaza). Ovo je najčešći slučaj u hardverskim realizacijama, pa je zbog toga tako odabrano. Izlazni format je definisan posebno i prosleđuje se funkciji. Primetiti da se nakon svakog sabiranja bitska širina sume povećava za 1, pa se finalno može dobiti rezultat koji ima mnogo više bita nego što je neophodno. Zbog toga je i uvedeno opciono podešavanje formata izlaza.
from fxpmath import Fxp
np.random.seed(0)
x = 5*np.random.randn(10)
print(x)
xFxp = Fxp(x, True, n_word = 12, n_frac = 8, overflow = 'wrap', rounding = 'trunc')
xFxp
[ 8.82026173 2.00078604 4.89368992 11.204466 9.33778995 -4.8863894 4.75044209 -0.75678604 -0.51609426 2.05299251]
fxp-s12/8([-7.18359375 2. 4.890625 -4.796875 -6.6640625 -4.8828125
4.75 -0.75390625 -0.515625 2.05078125])
np, mpl, plt, signal = importEssentialLibs(USE_WIDGETS)
from fxpmath import Fxp
def firDirectFXP(h, x, state = None, outFxpFormat = None):
M = len(h)
if state == None:
delayLine = Fxp(np.zeros(M), like = x) # like = x, kaze da je format delayLine-a isti kao i format x-a
else:
delayLine = state
if outFxpFormat == None:
y = Fxp(np.zeros(len(x)), signed = True, n_word = x.n_word + h.n_word + int(np.ceil(np.log2(M))), n_frac = x.n_frac + h.n_frac)
else:
y = Fxp(np.zeros(len(x)), dtype = outFxpFormat)
y.overflow = x.overflow
y.rounding = x.rounding
for n in range(len(x)):
delayLine[1:] = delayLine[:M-1]
delayLine[0] = x[n]
sumVal = Fxp(0, like = delayLine[0]*h[0])
for m in range(M):
sumVal += h[m]*delayLine[m]
y[n] = Fxp(sumVal, like = y)
return y
h = np.array([-0.0136, -0.0139, 0.0254, 0.0523, -0.0124, -0.0880, 0.0252, 0.3169, \
0.4807, 0.3169, 0.0252, -0.0880, -0.0124, 0.0523, 0.0254, -0.0139, -0.0136])
nh = np.arange(len(h))
fig = plt.figure(figsize = (12,4))
plt.subplots_adjust(bottom=0.25, wspace = 0.4)
# Prikaz koeficijenata filtra
ax1 = fig.add_subplot(1,3,1)
ax1.stem(nh, h, label = 'float')
ax1.set_xlabel(r'$n$')
ax1.set_ylabel(r'$h[n]$')
# Definisanje ulaznog signala, zbir dve sinusoide
F1 = 0.1
F2 = 0.43
N = 256
n = np.arange(N)
x = np.sin(2*np.pi*F1*n) + np.sin(2*np.pi*F2*n)
# FILTRIRANJE (floating point)
y = firDirect(h, x)
# overflow moze bit 'saturate' ili 'wrap'
overFlowMethod = 'saturate'
# rounding moze biti 'floor', 'trunc', 'around', 'ceil', 'fix'
roundingMethod = 'trunc'
xFxp = Fxp(x, signed = True, n_word = 7, n_frac = 5, overflow = overFlowMethod, rounding = roundingMethod)
hFxp = Fxp(h, signed = True, n_word = 8, n_frac = 6, overflow = overFlowMethod, rounding = roundingMethod)
# Prikaz zaokruzenih koeficijenata
ax1.stem(nh, hFxp, linefmt='C1--',markerfmt = 'C1s', label = 'fixed')
ax1.legend(loc = 'upper right');
# FILTIRARANJE (fixed point)
yFXP = firDirectFXP(hFxp, xFxp)
#yFXP = firDirectFXP(hFxp, xFxp, outFxpFormat = 'fxp-s10/8')
# Prikaz izlaza
ax2 = fig.add_subplot(1,3,2)
ax2.plot(n, y, label = 'float')
ax2.plot(n, yFXP, label = 'fixed')
ax2.set_xlabel(r'$n$')
ax2.set_ylabel(r'$y[n]$')
ax2.legend(loc = 'upper right')
# Prikaz razlike izlaza
ax2 = fig.add_subplot(1,3,3)
ax2.plot(n, yFXP - y)
ax2.set_xlabel(r'$n$')
ax2.set_ylabel(r'$e[n]$');
yFXP.dtype
'fxp-s20/11'