Filtri sa konačnim impulsnim odzivom (FIR) filtri

Digitalna obrada signala, Vladimir Petrović

U ovom notebook-u ćemo se baviti sintezom sistema sa konačnim impulsnim odzivom. Najpre će biti predstavljen metod ograničavanja impulsnog odziva idealnog filtra, a zatim i optimizacioni postupak projektovanja filtara.

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

Uvod

Filtri sa konačnim impulsnim odzivom (__F__inite __I__mpulse __R__esponse - FIR) u opštem obliku imaju funkciju prenosa izraženu sa: $$ H(z) = \sum_{n=0}^{N} b_n z^{-n}, $$ gde su $b_n$ koeficijenti filtra, tj. koeficijenti impulsnog odziva filtra. $$ H(z) = \mathcal{Z} \left\lbrace h[n] \right\rbrace = \sum_{n=0}^{N} h[n] z^{-n} \implies b_n = h[n] $$

Prenosna funkcija se može zapisati i na sledeći način

$$ H(z) = \sum_{n=0}^{N} b_n z^{-n} = \frac{1}{z^N}\sum_{n=0}^{N} b_n z^{N-n} = k\frac{\prod_{n=0}^{N} (z - z_{on})}{z^N}, $$

što znači da su svi polovi funkcije prenosa u tački $z = 0$ i da oblik funkcije prenosa određuje samo raspored nula. Zbog ove činjenice, FIR filtri moraju biti većeg reda nego IIR filtri koji zadovoljavaju iste specifikacije.

Projektovanje filtara ograničavanjem idealnog impulsnog odziva

Najjednostavniji, ali ne i optimalan, način da se dobiju koeficijenti funkcije prenosa je odsecanjem impulsnog odziva idealnog filtra. Impulsni odziv idealnog filtra je, naravno, beskonačan, ali se može analitički izračunati korišćenjem inverzne Furijeove transformacije.

$$ h_{ID}[n] = \frac{1}{2 \pi} \int_{-\pi}^{\pi} H_{ID} \left( e^{j\Omega} \right) e^{j\Omega n} d\Omega $$

Na primer, za NF filtar važi:

$$ H_{ID} \left( e^{j\Omega} \right) = \begin{cases}0, & -\pi < \Omega < -\Omega_c \\ 1, & -\Omega_c < \Omega < \Omega_c\\ 0, & \Omega_c < \Omega < \pi \end{cases}, $$

pa je impulsni odziv idealnog NF filtra

$$ h_{ID}[n] = \frac{1}{2 \pi} \int_{-\pi}^{\pi} H_{ID} \left( e^{j\Omega} \right) e^{j\Omega n} d\Omega = \frac{1}{2 \pi} \int_{-\Omega_c}^{\Omega_c} 1 \cdot e^{j\Omega n} d\Omega = \frac{\sin(n\Omega_c)}{n\pi}. $$

Ograničavanje impulsnog odziva podrazumeva množenje prozorskom funkcijom koja mora biti simetrična jer se u tom slučaju očuvava simetrija koeficijenata, što je neophodan uslov za linearnost faze FIR filtra. Ako je prozorska funkcija dužine $M$ odbiraka, onda se pre odsecanja impulsni odziv idealnog filtra mora zakasniti za $(M-1)/2$ odbiraka da bi dobijeni sistem bio kauzalan i da bi se očuvala simetrija koeficijenata pri množenju sa prozorskom funkcijom kojoj odbirci počinju od indeksa $n = 0$.

U sledećem primeru ćemo projektovati NF FIR filtar korišćenjem više različitih prozorskih funkcija. Treba obratiti pažnju na deljenje nulom, tj. na neodređeni izraz $$ \frac{\sin(0\cdot\Omega_c)}{0\cdot\pi}. $$

Ovo se dešava kada je red filtra paran, tj. kada je broj odbiraka impulsnog odziva neparan. Tada se središnji odbirak mora izračunati kao: $$ h_{ID}[0] = \lim_{n \rightarrow 0} \frac{\sin(n\Omega_c)}{n\pi} = \frac{\Omega_c}{\pi}. $$

np, mpl, plt, signal = importEssentialLibs(USE_WIDGETS)

import ipywidgets as widgets

Wc = np.pi/4

def lpFIRfilterWindow(windowname, N):
    fig, axs = plt.subplots(1, 3, figsize = [12, 4])
    plt.subplots_adjust(bottom=0.15, wspace = 0.4)
    
    if windowname in "Pravougaona":
        w = signal.boxcar(N+1)
    if windowname in "Trougaona":
        w = signal.triang(N+1)
    if windowname in "Hanova":
        w = signal.hann(N+1)
    if windowname in "Hemingova":
        w = signal.hamming(N+1)
    if windowname in "Blekmenova":
        w = signal.blackman(N+1)
    if windowname in "Bartletova": #vrlo slična trougaonoj samo sa nulama na krajevima
        w = signal.bartlett(N+1)
    
    na = np.arange(-N/2, N/2+1, 1)
    n = np.arange(N+1)
    
    hID = np.zeros(N+1)
    if N % 2 == 0:
        na = np.delete(na, N//2) # remove index 0
        hID[n != N//2] = np.sin(Wc*na)/(na*np.pi)
        hID[N//2] = Wc/np.pi # lim in na->0
    else:
        hID = np.sin(Wc*na)/(na*np.pi)
    
    h = hID * w
    
    # Stem filter
    ax = axs[0]
    ax.stem(n, h)
    ax.set_xlabel(r'$n$')
    ax.set_ylabel(r'$h[n]$')
    
    # Plot spectra
    W, H = signal.freqz(h, 1, worN = 1024)
    
    ax = axs[1]
    ax.plot(W, abs(H))
    ax.set_xlabel(r'$\Omega$')
    ax.set_ylabel(r'$|H(e^{j\Omega}))|$');
    ax.set_xlim([0, np.pi])
    
    ax = axs[2]
    ax.plot(W, np.unwrap(np.angle(H)))
    ax.set_xlabel(r'$\Omega$')
    ax.set_ylabel(r'$H_p(e^{j\Omega}))$');
    ax.set_xlim([0, np.pi])
    
widgets.interact(lpFIRfilterWindow, 
                 N = (5, 30, 1), 
                 windowname=["Pravougaona",
                             "Trougaona",  
                             "Hanova", 
                             "Hemingova", 
                             "Blekmenova",
                             "Bartletova"]);

Primetimo da množenjem različitim prozorskim funkcijama dobijamo različitu frekvencijsku karakteristiku. Naravno, ovo je očekivano, s obzirom na to da će proizvod u vremenskom domenu dati konvoluciju u frekvencijskom domenu. Dakle, frekvencijska karakteristika filtra će biti konvolucija frekvencijske karakteristike idealnog filtra i frekvencijske karakteristike prozorske funkcije. S obzirom na ovu činjenicu, lako je zaključiti da će slabljenje u nepropusnoj zoni i oscilacije u propusnoj zoni biti određeni slabljenjem bočnih lukova prozorske funkcije - što je slabljenje bočnih lukova veće, to su oscilacije u propusnom i nepropusnom opsegu manje. Takođe, širina prelazne zone će biti određena širinom glavnog luka prozorske funkcije - što je glavni luk uži, to je prelazna zona uža.

Primena Kajzerove prozorske funkcije u projektovanju FIR filtara

Od svih prozorskih funkcija, za projektovanje filtara se najčešće koristi Kajzerova prozorska funkcija. Glavni razlog je što se parametrom $\beta$ može podešavati slabljenje bočnih lukova, a brojem odbiraka $N$ se može podešavati širina glavnog luka. Postoje empirijske formule koje daju procenu ovih parametara na osnovu specifikacija filtra. Od njih treba krenuti prilikom projektovanja filtara. Međutim, treba naglasiti da one ne garantuju da će specifikacije biti zadovoljene. Neophodno je proveriti da li su specifikacije zadovoljene i iterativno podešavati parametre prozorske funkcije sve dok one ne budu zadovoljene. Najlakše je podešavati red filtra, tj. povećavati ga sve dok specifikacije ne budu zadovoljene.

Empirijsko određivanje parametara filtra se postiže na sledeći način:

1. Najpre se iz specifikacija graničnih učestanosti odredi širina prelazne zone, na primer, za NF filtar je: $$ B_t = \Omega_a - \Omega_p. $$ Ukoliko ima više prelaznih zona, onda se za parametar $B_t$ uzima ona najuža, na primer za PO filtar je: $$ B_t = \min(\Omega_{p1} - \Omega_{a1}, \Omega_{a2} - \Omega_{p2}). $$

2. Zatim se odrede vrednosti dozvoljenih amplituda oscilacija u propusnom i nepropusnom opsegu: $$ \delta_a = 10^{-0,05\alpha_a}, \; \; \delta_p = \frac{10^{0,05\alpha_p} - 1}{10^{0,05\alpha_p} + 1}. $$ Ukoliko ima više različitih propusnih/nepropusnih opsega, izračunaju se za sve opsege. Zatim se usvoji globalni parametar $\delta$ tako da je: $$ \delta = \min(\delta_{a1}, \delta_{p1}, \delta_{a2}, \delta_{p2}, ...). $$ U slučaju da je $\delta < \delta_a$, onda se izračuna novo slabljenje u nepropusnom opsegu prema izrazu: $$ \alpha_a = -20\log \delta. $$

3. Parametar $\beta$ se odredi na osnovu izraza: $$ \beta = \begin{cases} 0, & \alpha_a < 21\, \mathrm{dB}\\ 0{,}5842(\alpha_a - 21)^{0,4} + 0{,}07886(\alpha_a - 21), & 21\, \mathrm{dB} \leq \alpha_a \leq 50\, \mathrm{dB}\\ 0{,}1102(\alpha_a - 8{,}7), & \alpha_a > 50\, \mathrm{dB}\\ \end{cases}. $$

4. Dužina impulsnog odziva filtra, tj. broj odbiraka prozorske funkcije, se odredi iz sledećih izraza: $$ D = \begin{cases} 0{,}9222, & \alpha_a \leq 21\, \mathrm{dB}\\ \frac{\alpha_a - 7{,}95}{14{,}36}, & \alpha_a > 21\, \mathrm{dB}\\ \end{cases}, $$ $$ M \geq \frac{2\pi D}{B_t} + 1. $$

Nakon što je određena dužina impulsnog odziva i parametar $\beta$, potrebno je odrediti impulsni odziv idealnog filtra. Granične učestanosti idealnog filtra treba odrediti tako da se one nalaze na sredini prelaznih zona, tj: $$ \Omega_{c,i} = \frac{\Omega_{a,i} + \Omega_{p,i}}{2}. $$

NF filtar

np, mpl, plt, signal = importEssentialLibs(USE_WIDGETS)

import ipywidgets as widgets

fig = plt.figure(figsize = (13,3.5))
plt.subplots_adjust(bottom=0.15, wspace = 0.4)

def lpFIRfilterKaiserWindow(Wp, Wa, Ap, Aa):
    
    plt.clf()
    
    Bt = Wa - Wp
    
    # Delta
    delta_p = (10**(0.05*Ap)-1)/(10**(0.05*Ap)+1)
    delta_a = 10**(-0.05*Aa)
    
    delta = min([delta_p, delta_a])
    
    # Alpha_a
    if delta < delta_a:
        AaSpec = -20*np.log10(delta)
    else:
        AaSpec = Aa
        
    # Beta
    beta = 0;
    if AaSpec < 21:
        beta = 0
    elif AaSpec <= 50:
        beta = 0.5842*(AaSpec-21)**0.4 + 0.07886*(AaSpec-21)
    else:
        beta = 0.1102*(AaSpec-8.7)
    
    # Filter order
    if AaSpec <= 21:
        D = 0.9222
    else:
        D = (AaSpec-7.95)/14.36
    M = int(np.ceil(2*np.pi*D/Bt + 1))
    N = M - 1
    
    w = signal.kaiser(M, beta)
    
    Wc = (Wa + Wp)/2
    
    na = np.arange(-N/2, N/2+1, 1)
    n = np.arange(N+1)
    
    hID = np.zeros(N+1)
    if N % 2 == 0:
        na = np.delete(na, N//2) # remove index 0
        hID[n != N//2] = np.sin(Wc*na)/(na*np.pi)
        hID[N//2] = Wc/np.pi # lim in na->0
    else:
        hID = np.sin(Wc*na)/(na*np.pi)
    
    h = hID * w
    
    # Stem filter
    ax1 = fig.add_subplot(1,3,1)
    ax1.stem(n, h)
    ax1.set_xlabel(r'$n$')
    ax1.set_ylabel(r'$h[n]$')
    
    # Plot spectra
    W, H = signal.freqz(h, 1, worN = 1024)
    
    ax2 = fig.add_subplot(1,3,2)
    ax2.plot(W, 20*np.log10(abs(H)))
    ax2.set_xlabel(r'$\Omega$')
    ax2.set_ylabel(r'$|H(e^{j\Omega}))| [dB]$');
    ax2.set_xlim([0, np.pi])
    ax2.set_ylim([-70, 5])
    
    ax2.add_patch(mpl.patches.Polygon([[0, Ap], [Wp, Ap], [Wp, Ap+3], [0, Ap+3]], closed=True,
                                      fill=False, hatch='////', color = 'red'))
    ax2.add_patch(mpl.patches.Polygon([[0, -Ap-3 ], [Wp, -Ap-3], [Wp, -Ap], [0, -Ap]], closed=True,
                                      fill=False, hatch='////', color = 'red'))
    ax2.add_patch(mpl.patches.Polygon([[Wa, -Aa ], [np.pi, -Aa], [np.pi, -Aa+3], [Wa, -Aa+3]], closed=True,
                                      fill=False, hatch='////', color = 'red'))
    
    ax3 = fig.add_subplot(1,3,3)
    ax3.plot(W, np.unwrap(np.angle(H)))
    ax3.set_xlabel(r'$\Omega$')
    ax3.set_ylabel(r'$H_p(e^{j\Omega}))$');
    ax3.set_xlim([0, np.pi])
    
widgets.interact(lpFIRfilterKaiserWindow, 
                 Wp = (0, np.pi/2, 0.1),
                 Wa = (np.pi/2, np.pi, 0.1),
                 Ap = (0.2, 2, 0.2),
                 Aa = (20, 60, 10));

<Figure size 936x252 with 0 Axes>

Za neke kombinacije parametara, možemo primetiti da filtar ne zadovoljava specifikacije. Zbog toga je neophodno uvećavati red filtra sve dok specifikacije ne budu zadovoljene. Proveru specifikacija moramo raditi u programskog kodu i to ćemo uraditi u narednom primeru.

PO filtar

Za idealni PO filtar važi: $$ H_{ID} \left( e^{j\Omega} \right) = \begin{cases} 0, & -\pi < \Omega < -\Omega_{c2} \\ 1, & -\Omega_{c2} < \Omega < -\Omega_{c1}\\ 0, & -\Omega_{c1} < \Omega < \Omega_{c1}\\ 1, & \Omega_{c1} < \Omega < \Omega_{c2}\\ 0, & \Omega_{c2} < \Omega < \pi \end{cases}, $$

pa je impulsni odziv idealnog VF filtra

$$ h_{ID}[n] = \frac{1}{2 \pi} \int_{-\pi}^{\pi} H_{ID} \left( e^{j\Omega} \right) e^{j\Omega n} d\Omega = \frac{1}{2 \pi} \int_{-\Omega_{c2}}^{-\Omega_{c1}} 1 \cdot e^{j\Omega n} d\Omega + \frac{1}{2 \pi} \int_{\Omega_{c1}}^{\Omega_{c2}} 1 \cdot e^{j\Omega n} d\Omega = \frac{\sin(n\Omega_{c2}) - \sin(n\Omega_{c1})}{n\pi}. $$

np, mpl, plt, signal = importEssentialLibs(USE_WIDGETS)

import ipywidgets as widgets

Wa1 = 0.2*np.pi
Wp1 = 0.35*np.pi
Wp2 = 0.65*np.pi
Wa2 = 0.75*np.pi

Ap = 1
Aa = 30

Bt = min([Wp1 - Wa1, Wa2 - Wp2])

Wc1 = (Wa1 + Wp1)/2
Wc2 = (Wa2 + Wp2)/2

# Delta
delta_p = (10**(0.05*Ap)-1)/(10**(0.05*Ap)+1)
delta_a = 10**(-0.05*Aa)

delta = min([delta_p, delta_a])

# Alpha_a
AaSpec = -20*np.log10(delta)

# Beta
beta = 0;
if AaSpec < 21:
    beta = 0
elif AaSpec <= 50:
    beta = 0.5842*(AaSpec-21)**0.4 + 0.07886*(AaSpec-21)
else:
    beta = 0.1102*(AaSpec-8.7)

# Filter order
if AaSpec <= 21:
    D = 0.9222
else:
    D = (AaSpec-7.95)/14.36
M = int(np.ceil(2*np.pi*D/Bt + 1))
N = M - 1
Norig = N 

filterFitsInSpecs = False
first = True
while not filterFitsInSpecs:
    w = signal.kaiser(N+1, beta)

    na = np.arange(-N/2, N/2+1, 1)
    n = np.arange(N+1)

    hID = np.zeros(N+1)
    if N % 2 == 0:
        na = np.delete(na, N//2) # remove index 0
        hID[n != N//2] = np.sin(Wc2*na)/(na*np.pi) - np.sin(Wc1*na)/(na*np.pi)
        hID[N//2] = Wc2/np.pi - Wc1/np.pi  # lim in na->0
    else:
        hID = np.sin(Wc2*na)/(na*np.pi) - np.sin(Wc1*na)/(na*np.pi)

    h = hID * w
    
    numPoints = 1024
    W, H = signal.freqz(h, 1, worN = numPoints)
    Hdb = 20*np.log10(abs(H))
    
    if first:
        first = False
        HdbOriginal = Hdb
        
    DeltaW = np.pi/numPoints
    kSB1 = int(np.ceil(Wa1/DeltaW))
    kPB1 = int(np.floor(Wp1/DeltaW))
    kPB2 = int(np.ceil(Wp2/DeltaW))
    kSB2 = int(np.floor(Wa2/DeltaW))
    
    if np.all(Hdb[:kSB1] < -Aa) and np.all(Hdb[kPB1:kPB2] > -Ap) and \
        np.all(Hdb[kPB1:kPB2] < Ap) and np.all(Hdb[kSB2:] < -Aa):
        filterFitsInSpecs = True 
    else:
        N += 1

fig = plt.figure(figsize = (12,5))
plt.subplots_adjust(bottom=0.15, 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]$')

# Plot spectra
ax2 = fig.add_subplot(1,2,2)
ax2.plot(W, HdbOriginal, label = f'Original: N = {Norig}')
ax2.plot(W, 20*np.log10(abs(H)), label = f'Modified: N = {N}')
ax2.set_xlabel(r'$\Omega$')
ax2.set_ylabel(r'$|H(e^{j\Omega}))| [dB]$');
ax2.set_xlim([0, np.pi])
ax2.set_ylim([-60, 5])
ax2.legend(loc = 'lower center')

# First stop band lim
ax2.add_patch(mpl.patches.Polygon([[0, -Aa ], [Wa1, -Aa], [Wa1, -Aa+2], [0, -Aa+2]], closed=True,
                                  fill=False, hatch='////', color = 'red'))
# Pass band lims
ax2.add_patch(mpl.patches.Polygon([[Wp1, -Ap-2 ], [Wp2, -Ap-2], [Wp2, -Ap], [Wp1, -Ap]], closed=True,
                                  fill=False, hatch='////', color = 'red'))
ax2.add_patch(mpl.patches.Polygon([[Wp1, Ap ], [Wp2, Ap], [Wp2, Ap+2], [Wp1, Ap+2]], closed=True,
                                  fill=False, hatch='////', color = 'red'))
# Second stop band lim
ax2.add_patch(mpl.patches.Polygon([[Wa2, -Aa ], [np.pi, -Aa], [np.pi, -Aa+2], [Wa2, -Aa+2]], closed=True,
                                  fill=False, hatch='////', color = 'red'));

Stepeničasti filtar

FIR filtre možemo projektovati i tako da imaju složeniju frekvencijsku karakteristiku. One ne moraju nužno biti NF, VF, PO i NO oblika. Stoga projektujmo filtar kome su zadate sledeće specifikacije:

• Granične učestanosti nepropusnog opsega su: $\Omega_{a1} = 0$ i $\Omega_{a2} = 0{,}23\pi$, a minimalno slabljenje u ovom opsegu iznosi $\alpha_a = 35\, \mathrm{dB}$.
• Granične učestanosti prvog propusnog opsega su: $\Omega_{p1} = 0{,}3\pi$ i $\Omega_{p2} = 0{,}55\pi$, a maksmimalno slabljenje u ovom opsegu iznosi $\alpha_{p1} = 1\, \mathrm{dB}$. Amplitudska karakteristika u ovom opsegu treba da ima vrednost $H_1 = 0{,}45$.
• Granične učestanosti drugog propusnog opsega su: $\Omega_{p3} = 0{,}65\pi$ i $\Omega_{p4} = \pi$, a maksmimalno slabljenje u ovom opsegu iznosi $\alpha_{p2} = 0{,}65\, \mathrm{dB}$. Amplitudska karakteristika u ovom opsegu treba da ima vrednost $H_2 = 1$.

Idealni filtar ima frekvencijsku karakteristiku kao na slici ispod.

Inverznom Furijevom transformacijom, dobijamo da je njegov impulsni odziv: $$ h_{ID}[n] = \frac{1}{2 \pi} \int_{-\pi}^{\pi} H_{ID} \left( e^{j\Omega} \right) e^{j\Omega n} d\Omega = H_1\frac{\sin(n\Omega_{c2}) - \sin(n\Omega_{c1})}{n\pi} + \frac{\sin(n\pi) - \sin(n\Omega_{c2})}{n\pi}. $$

Dodatna stvar o kojoj treba povesti računa u ovom primeru je određivanje vrednosti slabljenja i varijacija amplitude u propusnom opsegu kada vrednost idealne amplitudske karakteristike nije jednaka 1, već nekoj vrednosti $H$. Izrazi su vrlo slični, ali postoji oređena razlika. $$ \alpha_p = 20\log \left( \frac{H+\delta_p}{H-\delta_p} \right) $$ Odatle je: $$ 10^{0{,}05\alpha_p} = \frac{H+\delta_p}{H-\delta_p} \implies \delta_p = H \frac{10^{0{,}05\alpha_p} - 1}{10^{0{,}05\alpha_p} + 1} $$

np, mpl, plt, signal = importEssentialLibs(USE_WIDGETS)

import ipywidgets as widgets

Wa1 = 0
Wa2 = 0.23*np.pi
Wp1 = 0.3*np.pi
Wp2 = 0.55*np.pi
Wp3 = 0.65*np.pi
Wp4 = np.pi

Ap1 = 1
Ap2 = 0.65
Aa = 35
AaOrig = Aa

H1 = 0.45
H2 = 1

Bt = min([Wp1 - Wa2, Wp3 - Wp2])

Wc1 = (Wa2 + Wp1)/2
Wc2 = (Wp2 + Wp3)/2

# Delta
delta_a = 10**(-0.05*Aa)
delta_p1 = H1*(10**(0.05*Ap1)-1)/(10**(0.05*Ap1)+1)
delta_p2 = H2*(10**(0.05*Ap2)-1)/(10**(0.05*Ap2)+1)

delta = min([delta_p1, delta_p2, delta_a])

# Alpha_a
if delta < delta_a:
    Aa = -20*np.log10(delta)

# Beta
beta = 0;
if Aa < 21:
    beta = 0
elif Aa <= 50:
    beta = 0.5842*(Aa-21)**0.4 + 0.07886*(Aa-21)
else:
    beta = 0.1102*(Aa-8.7)

# Filter order
if Aa <= 21:
    D = 0.9222
else:
    D = (Aa-7.95)/14.36
M = int(np.ceil(2*np.pi*D/Bt + 1))
N = M - 1
Norig = N 

filterFitsInSpecs = False
first = True
while not filterFitsInSpecs:
    w = signal.kaiser(N+1, beta)

    na = np.arange(-N/2, N/2+1, 1)
    n = np.arange(N+1)

    hID = np.zeros(N+1)
    if N % 2 == 0:
        na = np.delete(na, N//2) # remove index 0
        hID[n != N//2] = H1*(np.sin(Wc2*na)/(na*np.pi) - np.sin(Wc1*na)/(na*np.pi)) + \
                         H2*(np.sin(np.pi*na)/(na*np.pi) - np.sin(Wc2*na)/(na*np.pi))
        hID[N//2] = H1*(Wc2/np.pi - Wc1/np.pi) + H2*(1 - Wc2/np.pi) # lim in na->0
    else:
        hID = H1*(np.sin(Wc2*na)/(na*np.pi) - np.sin(Wc1*na)/(na*np.pi)) + \
              H2*(np.sin(np.pi*na)/(na*np.pi) - np.sin(Wc2*na)/(na*np.pi))

    h = hID * w
    
    numPoints = 1024
    W, H = signal.freqz(h, 1, worN = numPoints)
    Hdb = 20*np.log10(abs(H))
    
    if first:
        first = False
        HdbOriginal = Hdb
        
    DeltaW = np.pi/numPoints
    kSB1 = int(np.ceil(Wa2/DeltaW))
    kPB1 = int(np.floor(Wp1/DeltaW))
    kPB2 = int(np.ceil(Wp2/DeltaW))
    kPB3 = int(np.floor(Wp3/DeltaW))
    
    if  np.all(Hdb[:kSB1] < -AaOrig) and \
        np.all(Hdb[kPB1:kPB2] > 20*np.log10(H1-delta_p1)) and \
        np.all(Hdb[kPB1:kPB2] < 20*np.log10(H1+delta_p1)) and \
        np.all(Hdb[kPB3:] > -Ap2) and \
        np.all(Hdb[kPB3:] < Ap2):
        filterFitsInSpecs = True
    else:
        N += 1

fig = plt.figure(figsize = (12,5))
plt.subplots_adjust(bottom=0.15, 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]$')

# Plot spectra
ax2 = fig.add_subplot(1,2,2)
ax2.plot(W, HdbOriginal, label = f'Original: N = {Norig}')
ax2.plot(W, 20*np.log10(abs(H)), label = f'Modified: N = {N}')
ax2.set_xlabel(r'$\Omega$')
ax2.set_ylabel(r'$|H(e^{j\Omega}))| [dB]$');
ax2.set_xlim([0, np.pi])
ax2.set_ylim([-60, 5])
ax2.legend(loc = 'lower right')

# Stop band lim
ax2.add_patch(mpl.patches.Polygon([[Wa1, -AaOrig ], [Wa2, -AaOrig], [Wa2, -AaOrig+2], [Wa1, -AaOrig+2]], closed=True,
                                  fill=False, hatch='////', color = 'red'))
# First pass band lims
ax2.add_patch(mpl.patches.Polygon([[Wp1, 20*np.log10(H1-delta_p1)-2], [Wp2, 20*np.log10(H1-delta_p1)-2], [Wp2, 20*np.log10(H1-delta_p1)], [Wp1, 20*np.log10(H1-delta_p1)]], closed=True,
                                  fill=False, hatch='////', color = 'red'))
ax2.add_patch(mpl.patches.Polygon([[Wp1, 20*np.log10(H1+delta_p1)], [Wp2, 20*np.log10(H1+delta_p1)], [Wp2, 20*np.log10(H1+delta_p1)+2], [Wp1, 20*np.log10(H1+delta_p1)+2]], closed=True,
                                  fill=False, hatch='////', color = 'red'))
# Second pass band lims
ax2.add_patch(mpl.patches.Polygon([[Wp3, -Ap2-2], [Wp4, -Ap2-2], [Wp4, -Ap2], [Wp3, -Ap2]], closed=True,
                                  fill=False, hatch='////', color = 'red'))
ax2.add_patch(mpl.patches.Polygon([[Wp3, Ap2], [Wp4, Ap2], [Wp4, Ap2+2], [Wp3, Ap2+2]], closed=True,
                                  fill=False, hatch='////', color = 'red'));

Optimizacioni postupak projekovanja FIR filtara

Metod projektovanja FIR filtara ograničavanjem idelnog impulsnog odziva nije optimalan metod ni po kom kriterijumu. Videli smo, vrlo je jednostavno projektovati filtre ovim metodom i to je glavna prednost ovog metoda projektovanja. Zbog Gibsovog fenomena koji nastaje usled odsecanja idealnog impulsnog odziva, najveće varijacije amplitudske karakteristike su na granicama sa prelaznim zonama. Međutim, moguće je ravnomernije rasporediti varijacije amplituda, tako da specifikacije budu ispoštovane, a da se smanji red filtra. Ovo se postiže optimizacionim postupkom.

Cilj optimizacionog postupka je da menja koeficijente filtra sve dok neka optimizaciona funkcija ne dobije minimalnu vrednost. Ova optimizaciona (engl. cost) funkcija može biti različitog oblika. Obično se polazi od greške realne amplitudske karakteristike u odnosu na idealnu amplitudsku karakteristiku:

$$ E(\Omega) = |H_{ID}(e^{j\Omega})| - |H(e^{j\Omega})|. $$

Različite optimizacione funkcije se mogu formirati od navedene greške. Najčešća dva metoda su minimax i metoda najmanjeg kvadrata (least squares). Kod prvog metoda je cilj minimizovati maksimalnu grešku i na taj način ravnomerno rasporediti grešku po celom propusnom i nepropusnom opsegu:

$$ \underset{b_n,\, n \in [0, N-1]}{\min} \left( \max (|E(\Omega)|) \right). $$

Kod drugog metoda je cilj naći koeficijente filtra kod koga će srednja kvadratna greška biti minimalna: $$ \underset{b_n,\, n \in [0, N-1]}{\min} \left(\frac{1}{2\pi} \int_{-\pi}^{\pi} |E(\Omega)|^2 d\Omega \right). $$

Ukoliko se želi postići da varijacije amplitude u različitim opsezima budu različite, onda se funkcija greške može otežiniti, tj. skalirati u odgovarajućim opsezima. Na taj način se menja relativan uticaj greške u tim opsezima i ekvivalentno postiže različito slabljenje. Težinska funkcija za, na primer, NF filtar je:

$$ W(\Omega) = \begin{cases} 1, & 0 \leq |\Omega| \leq \Omega_p \\ \delta_p/\delta_a, & \Omega_a \leq |\Omega| \leq \pi \end{cases}, $$

pa je otežinjena funkcija greške sada:

$$ E(\Omega) = W(\Omega)(|H_{ID}(e^{j\Omega})| - |H(e^{j\Omega})|). $$

Detalji o tome kako rade ovi optimizacioni algoritmi se mogu pronaći u knjizi i materijalima za predavanja. Za sada nam oni nisu neophodni, bitno je da razumemo da se na neki način menjaju koeficijenti filtra sa ciljem minimizacije određene optimizacione funkcije.

Parks-MekLelanov (Remezov) metod projektovanja

Remezov metod koristi minimax kroterijum optimizacije. Poziva se korišćenjem Scipy funkcije signal.remez(). Ovoj funkciji treba da prosledimo red filtra, oblik idealne frekvencijske karakteristike, vektor relativnih graničnih učestanosti i, opciono, vektor težinskih faktora kojim se obezbeđuje relativan uticaj greške u različitim opsezima.

Postoje empirijske formule za određivanje reda filtra. Ove formule će dati samo procenu reda, i najčešće specifikacije neće biti ispoštovane. Zbog toga je, kao i ranije, potrebno da proveramo da li su specifikacije ispoštovane i da povećavamo red filtra sve dok one ne postanu ispoštovane. Empiijske formule se razlikuju za NF/VF funkcije prenosa i za PO/NO funkcije prenosa. Ovde ih navodimo radi kompletnosti.

Za NF/VF filtar: $$ D_{\infty}(\delta_p, \delta_a) = (0.005309\log(\delta_p)^2 + 0.07114\log(\delta_p) - 0.4761)\log(\delta_a) - \\ (0.00266\log(\delta_p)^2 + 0.5941\log(\delta_p) + 0.4278)\\ f(\delta_p, \delta_a) = 11.01217 + 0.51244(\log(\delta_p) - \log(\delta_a)) \\ B_t = 2\pi\frac{|f_p - f_a|}{f_s} \\ N \geq \frac{2\pi D_{\infty}(\delta_p, \delta_a)}{B_t} - \frac{f(\delta_p, \delta_a)B_t}{2\pi} $$

Za PO/NO filtar: $$ D_{\infty}(\delta_p, \delta_a) = (0.01201\log(\delta_p)^2 + 0.09664\log(\delta_p) - 0.51325)\log(\delta_a) + \\ (0.00203\log(\delta_p)^2 - 0.57054\log(\delta_p) - 0.44314)\\ f(\delta_p, \delta_a) = -16.9 + 14.6(\log(\delta_p) - \log(\delta_a)) \\ B_t = \min \left( 2\pi\frac{|f_{p,i} - f_{a,i}|}{f_s}\right) \\ \delta_p = \min(\delta_{p,i}),\, \, \, \delta_a = \min(\delta_{a,i}) \\ N \geq \frac{2\pi D_{\infty}(\delta_p, \delta_a)}{B_t} - \frac{f(\delta_p, \delta_a)B_t}{2\pi} \\ $$

U narednom primeru se koristi Remezov algoritam za projektovanje VF filtra.

np, mpl, plt, signal = importEssentialLibs(USE_WIDGETS)

# Specifikacije
fs = 2000
fa = 500
fp = 700
Ap = 1
Aa = 40
dp = (10**(0.05*Ap)-1)/(10**(0.05*Ap)+1)
da = 10**(-0.05*Aa)

# Red filtra
D = (0.005309*np.log10(dp)**2 + 0.07114*np.log10(dp) - 0.4761)*np.log10(da)
D = D - (0.00266*np.log10(dp)**2 + 0.5941*np.log10(dp) + 0.4278)
f = 11.01217 + 0.51244*(np.log10(dp) - np.log10(da))
Bt = (fp-fa)/fs*2*np.pi
M = int(np.ceil(2*np.pi*D/Bt - f*Bt/(2*np.pi) + 1))

# Relativne ucestanosti
Fp = fp/fs
Fa = fa/fs
# Zeljena idealna amplitudska karakteristika
Hid = [0, 1]
# Vektor relativnih ucestanosti za koje se zadaje Hid(w)
F = [0, Fa, Fp, 0.5]
# Vektor tezinskih faktora
W = [dp/da, 1]
# Projektovanje filtra
h = signal.remez(M, F, Hid, W);

fig, axs = plt.subplots(1, 3, figsize = [12, 4])
plt.subplots_adjust(bottom=0.15, wspace = 0.4)

# Impulsni odziv
n = np.arange(M)
ax = axs[0]
ax.stem(n, h)
ax.set_xlabel(r'$n$')
ax.set_ylabel(r'$h[n]$')

# Spektar
W, H = signal.freqz(h, 1, worN = 1024)

ax = axs[1]
ax.plot(W, abs(H))
ax.set_xlabel(r'$\Omega$')
ax.set_ylabel(r'$|H(e^{j\Omega}))|$');
ax.set_xlim([0, np.pi])

ax.add_patch(mpl.patches.Polygon([[0, da], [2*np.pi*Fa, da], [2*np.pi*Fa, da+0.1], [0, da+0.1]], closed=True,
                                  fill=False, hatch='////', color = 'red'))
ax.add_patch(mpl.patches.Polygon([[2*np.pi*Fp, 1+dp], [np.pi, 1+dp], [np.pi, 1+dp+0.1], [2*np.pi*Fp, 1+dp+0.1]], closed=True,
                                  fill=False, hatch='////', color = 'red'))
ax.add_patch(mpl.patches.Polygon([[2*np.pi*Fp, 1-dp-0.1], [np.pi, 1-dp-0.1], [np.pi, 1-dp], [2*np.pi*Fp, 1-dp]], closed=True,
                                  fill=False, hatch='////', color = 'red'))

ax = axs[2]
ax.plot(W, np.unwrap(np.angle(H)))
ax.set_xlabel(r'$\Omega$')
ax.set_ylabel(r'$H_p(e^{j\Omega}))$');
ax.set_xlim([0, np.pi]);

Primetimo da specifikacije nisu u potpunosti ispoštovane. Dakle, potrebno je iterativno povećavati red filtra sve dok se ne dobije željena amplitudska karakteristika. To ćemo uraditi u narednom primeru višestrukog propusnika opsega kome su zadate sledeće specifikacije:

• Minimalno slabljenje u svim nepropusnim opsezima iznosi $\alpha_a = 30\, \mathrm{dB}$.
• Maksimalno slabljenje u svim propusnim opsezima iznosi $\alpha_p = 0{,}5\, \mathrm{dB}$.
• Granične učestanosti prvog nepropusnog opsega su: $\Omega_{a0} = 0$ i $\Omega_{a1} = 0{,}13\pi$.
• Granične učestanosti prvog propusnog opsega su: $\Omega_{p1} = 0{,}18\pi$ i $\Omega_{p2} = 0{,}38\pi$. Amplitudska karakteristika u ovom opsegu treba da ima vrednost $H_1 = 0{,}45$.
• Granične učestanosti drugog nepropusnog opsega su: $\Omega_{a2} = 0{,}43$ i $\Omega_{a3} = 0{,}64\pi$.
• Granične učestanosti drugog propusnog opsega su: $\Omega_{p3} = 0{,}73\pi$ i $\Omega_{p4} = 0{,}9\pi$. Amplitudska karakteristika u ovom opsegu treba da ima vrednost $H_2 = 1$.
• Granične učestanosti trećeg nepropusnog opsega su: $\Omega_{a4} = 0{,}97$ i $\Omega_{a5} = \pi$.

np, mpl, plt, signal = importEssentialLibs(USE_WIDGETS)

H1 = 0.45
H2 = 1
Wa1 = 0.13*np.pi 
Wa2 = 0.43*np.pi
Wa3 = 0.64*np.pi # razmisli da je promenis
Wa4 = 0.97*np.pi 
Wa5 = np.pi; 
Wp1 = 0.18*np.pi
Wp2 = 0.38*np.pi
Wp3 = 0.73*np.pi
Wp4 = 0.9*np.pi
Ap = 0.5
Aa = 30
dp1 = H1*(10**(0.05*Ap)-1)/(10**(0.05*Ap)+1)
dp2 = H2*(10**(0.05*Ap)-1)/(10**(0.05*Ap)+1)
dp = min([dp1, dp2])
da = 10**(-0.05*Aa)

D = (0.001201*np.log10(dp)**2 + 0.09664*np.log10(dp) - 0.51325)*np.log10(da)
D = D + 0.00203*np.log10(dp)**2 - 0.57054*np.log10(dp) - 0.44314
f = -16.9 - 14.6*(np.log10(dp)-np.log10(da))
Bt1 = abs(Wa1-Wp1)
Bt2 = abs(Wa2-Wp2)
Bt3 = abs(Wa3-Wp3)
Bt4 = abs(Wa4-Wp4)
Bt = min([Bt1, Bt2, Bt3, Bt4])
M = int(np.ceil(2*np.pi*D/Bt - f*Bt/(2*np.pi) + 1))

Nfft = 1024;

Fp_1 = Wp1/2/np.pi 
Fp_2 = Wp2/2/np.pi
Fp_3 = Wp3/2/np.pi
Fp_4 = Wp4/2/np.pi
Fa_1 = Wa1/2/np.pi
Fa_2 = Wa2/2/np.pi
Fa_3 = Wa3/2/np.pi
Fa_4 = Wa4/2/np.pi
Hid = [0, H1, 0, H2, 0]
F = [0, Fa_1, Fp_1, Fp_2, Fa_2, Fa_3, Fp_3, Fp_4, Fa_4, 0.5]
Wstop1 = dp/da
Wstop2 = dp/da
Wstop3 = dp/da
Wpass1 = dp/dp1
Wpass2 = dp/dp2   
Weight = [Wstop1, Wpass1, Wstop2, Wpass2, Wstop3]

fig, axs = plt.subplots(1, 3, figsize = [12, 4])
plt.subplots_adjust(bottom=0.15, wspace = 0.4)

filterFitsInSpecs = False
first = True
while not filterFitsInSpecs:
    h = signal.remez(M, F, Hid, Weight)
    W, H = signal.freqz(h, 1, worN = Nfft)
    Ha = abs(H)
            
    DeltaW = np.pi/Nfft
    kSB1 = int(np.ceil(Wa1/DeltaW))
    kPB1 = int(np.floor(Wp1/DeltaW))
    kPB2 = int(np.ceil(Wp2/DeltaW))
    kSB2 = int(np.ceil(Wa2/DeltaW))
    kSB3 = int(np.ceil(Wa3/DeltaW))
    kPB3 = int(np.floor(Wp3/DeltaW))
    kPB4 = int(np.floor(Wp4/DeltaW))
    kSB4 = int(np.ceil(Wa4/DeltaW))
    
    if  np.all(Ha[:kSB1] < da) and \
        np.all(Ha[kPB1:kPB2] > H1-dp1) and \
        np.all(Ha[kPB1:kPB2] < H1+dp1) and \
        np.all(Ha[kSB2:kSB3] < da) and \
        np.all(Ha[kPB3:kPB4] > H2-dp2) and \
        np.all(Ha[kPB3:kPB4] < H2+dp2) and \
        np.all(Ha[kSB4:] < da):
        filterFitsInSpecs = True
    else:
        M += 1
        print(M)

    # Impulsni odziv
    n = np.arange(len(h))
    ax = axs[0]
    ax.stem(n, h)
    ax.set_xlabel(r'$n$')
    ax.set_ylabel(r'$h[n]$')

    # Spektar
    W, H = signal.freqz(h, 1, worN = 1024)

    ax = axs[1]
    ax.plot(W, abs(H))
    ax.set_xlabel(r'$\Omega$')
    ax.set_ylabel(r'$|H(e^{j\Omega}))|$');
    ax.set_xlim([0, np.pi])

    ax.add_patch(mpl.patches.Polygon([[0, da], [Wa1, da], [Wa1, da+0.1], [0, da+0.1]], closed=True,
                                      fill=False, hatch='////', color = 'red'))
    ax.add_patch(mpl.patches.Polygon([[Wp1, H1+dp1], [Wp2, H1+dp1], [Wp2, H1+dp1+0.1], [Wp1, H1+dp1+0.1]], closed=True,
                                      fill=False, hatch='////', color = 'red'))
    ax.add_patch(mpl.patches.Polygon([[Wp1, H1-dp1-0.1], [Wp2, H1-dp1-0.1], [Wp2, H1-dp1], [Wp1, H1-dp1]], closed=True,
                                      fill=False, hatch='////', color = 'red'))
    ax.add_patch(mpl.patches.Polygon([[Wa2, da], [Wa3, da], [Wa3, da+0.1], [Wa2, da+0.1]], closed=True,
                                      fill=False, hatch='////', color = 'red'))
    ax.add_patch(mpl.patches.Polygon([[Wp3, H2+dp2], [Wp4, H2+dp2], [Wp4, H2+dp2+0.1], [Wp3, H2+dp2+0.1]], closed=True,
                                      fill=False, hatch='////', color = 'red'))
    ax.add_patch(mpl.patches.Polygon([[Wp3, H2-dp2-0.1], [Wp4, H2-dp2-0.1], [Wp4, H2-dp2], [Wp3, H2-dp2]], closed=True,
                                      fill=False, hatch='////', color = 'red'))
    ax.add_patch(mpl.patches.Polygon([[Wa4, da], [np.pi, da], [np.pi, da+0.1], [Wa4, da+0.1]], closed=True,
                                      fill=False, hatch='////', color = 'red'))
    ax = axs[2]
    ax.plot(W, np.unwrap(np.angle(H)))
    ax.set_xlabel(r'$\Omega$')
    ax.set_ylabel(r'$H_p(e^{j\Omega}))$');
    ax.set_xlim([0, np.pi]);

Metod najmanjeg kvadrata

Ako se za optimizacionu funkciju odabere minimalna srednja kvadratna greška, postupak je isti, ali se za pojektovanje filtra koristi funkcija signal.firls().