Progetto di filtri FIR con il metodo delle finestre

Un filtro ideale passa-basso lascia passare tutto fino a una frequenza di taglio e blocca tutto il resto: ha risposta impulsiva infinita e non causale, quindi non si può realizzare. Il progetto di un filtro FIR (FIR filter design) consiste nel trovare i coefficienti $h[0],\dots,h[N]$ di un filtro causale a risposta impulsiva finita la cui risposta in frequenza (Risposta in frequenza dei sistemi FIRSe all'ingresso di un FIR c'è un esponenziale complesso A e^{jφ} e^{jω̂n} (per ogni n), l'uscita è lo stesso esponenziale moltiplicato per H(ω̂) = Σ b_k e^{-jω̂k}: la frequenza non cambia, ampiezza e fase sono modificate da |H| (guadagno) e ∠H (sfasamento). Per sovrapposizione si trattano somme di sinusoidi. H è periodica di periodo 2π e, per coefficienti reali, hermitiana (|H| pari, fase dispari). La cascata ha H = H1·H2. Esempi: ritardo (fase lineare), differenza prima (passa-alto), {1,2,1} (passa-basso), media mobile di L punti (Dirichlet: |H| = |sin(Lω̂/2)/(L sin(ω̂/2))|, fase lineare -(L-1)ω̂/2).Risposta in frequenza dei sistemi FIR →) approssima quella ideale entro tolleranze assegnate. Qui ci si limita ai filtri a fase lineare, per i quali il progetto riguarda solo il modulo (Sistemi a fase lineare e assenza di distorsioneUn sistema non deforma il segnale se y[n] = K x[n-n0]: modulo costante e fase lineare -n0 w, cioè ritardo di gruppo costante n0. Un sistema reale e causale ha fase lineare (generalizzata) se e solo se è FIR con risposta impulsiva simmetrica h[n] = h[N-n] (ampiezza pari, fase -N w/2) o antisimmetrica h[n] = -h[N-n] (ampiezza dispari, fase -N w/2 + pi/2). Il ritardo di gruppo è N/2: intero se N è pari (vale la condizione di non distorsione), semi-intero se N è dispari (uscita interpolata e ritardata). Un IIR causale non può essere simmetrico, quindi non ha fase lineare esatta.Sistemi a fase lineare e assenza di distorsione →, Filtri FIR a fase lineare - tipi e zeriI FIR a fase lineare di ordine N si dividono in quattro tipi: I (N pari, h simmetrica), II (N dispari, simmetrica), III (N pari, antisimmetrica), IV (N dispari, antisimmetrica). Si scrive H = e^{-jwN/2} e^{jb} A(w), con A somma di coseni (tipi I e II) o di seni (III e IV). Dalla relazione H(z) = ±z^{-N} H(1/z) segue che gli zeri vengono in gruppi (z0, z0*, 1/z0, 1/z0*) e che ci sono zeri forzati: tipo II in z = -1, tipo III in z = 1 e z = -1, tipo IV in z = 1, tipo I nessuno. Quindi II non fa passa-alto, III non fa passa-basso né passa-alto, IV non fa passa-basso. Ogni tipo II, III, IV è un tipo I moltiplicato per (1+z^-1), (1-z^-2), (1-z^-1).Filtri FIR a fase lineare - tipi e zeri →). Tre metodi: finestre (questo è il metodo principale del corso), campionamento in frequenza e minimax.

Le specifiche di un filtro passa-basso

Gli altri tipi di filtro (passa-alto, passa-banda, elimina-banda) si ricavano da un prototipo passa-basso, quindi si descrive solo questo. Le specifiche sono date sul modulo della risposta in frequenza, perché le tolleranze non possono essere nulle e fra la banda passante e quella oscura serve una banda di transizione (transition band).

Definizione (specifiche del passa-basso).

• $\hat\omega_p$: limite della banda passante (pass-band edge); per $|\hat\omega|\le\hat\omega_p$ il segnale passa praticamente inalterato.
• $\hat\omega_s$: limite della banda oscura (stop-band edge); per $\hat\omega_s\le|\hat\omega|\le\pi$ il segnale è fortemente soppresso. Per $\hat\omega_p<|\hat\omega|<\hat\omega_s$ si è in banda di transizione e non si impone nulla.
• $\delta_p$: errore in banda passante, $1-\delta_p\le|H(e^{j\hat\omega})|\le1+\delta_p$ (l'ideale sarebbe $|H|=1$).
• $\delta_s$: errore in banda oscura, $|H(e^{j\hat\omega})|\le\delta_s$ (l'ideale sarebbe $|H|=0$).

In decibel: ondulazione in banda passante (pass-band ripple) $R_p=20\log_{10}\dfrac{1+\delta_p}{1-\delta_p}$ e attenuazione in banda oscura (stop-band attenuation) $A_s=-20\log_{10}\delta_s$.

Esempio. Con $\delta_p=\delta_s=10^{-2.5}=0{,}00316$ si ha $A_s=50$ dB e $R_p=20\log_{10}\frac{1{,}00316}{0{,}99684}=0{,}055$ dB. Con $\delta_p=0{,}05$ e $\delta_s=0{,}01$ si ha $R_p=20\log_{10}\frac{1{,}05}{0{,}95}=0{,}87$ dB e $A_s=40$ dB.

Filtro causale a fase lineare di tipo I

Un FIR causale di ordine $N$ ha $N+1$ coefficienti $h[0],\dots,h[N]$; è a fase lineare se $h[n]=\pm h[N-n]$. Nel tipo I ($N$ pari, simmetria $+$) la risposta in frequenza si scrive

$$H(e^{j\hat\omega})=e^{-j\hat\omega N/2}\,\bar H(\hat\omega),\qquad \bar H(\hat\omega)=\sum_{n=0}^{N/2}p_n\cos(n\hat\omega),$$

con $p_0=h[N/2]$ e $p_n=2h[N/2-n]$ per $n=1,\dots,N/2$. La parte $\bar H$ è reale e si chiama risposta di ampiezza (amplitude response); la fase è esattamente $-\hat\omega N/2$ (ritardo di $N/2$ campioni). Il perché: si porta il filtro in forma simmetrica traslando di $N/2$ campioni, $\bar h[n]=h[n+N/2]$, per $n=-N/2,\dots,N/2$; essendo $\bar h[n]=\bar h[-n]$, le esponenziali $e^{-jn\hat\omega}$ e $e^{jn\hat\omega}$ si accoppiano in $2\cos(n\hat\omega)$ e la somma diventa reale. Negli altri tre tipi di fase lineare cambia solo un fattore fisso (e quindi un vincolo sugli zeri in $z=\pm1$), mentre $\bar H$ resta una combinazione di coseni (o seni).

Esempio. Per $N=4$ e $h=\{1,2,3,2,1\}$ si ha $p_0=h[2]=3$, $p_1=2h[1]=4$, $p_2=2h[0]=2$, quindi $\bar H(\hat\omega)=3+4\cos\hat\omega+2\cos2\hat\omega$ e $\bar H(0)=9=\sum h[n]$.

Obiettivo: la serie di Fourier dell'ampiezza desiderata

Si vuole un filtro la cui risposta approssimi

$$H_d(e^{j\hat\omega})=e^{-j\hat\omega N/2}\,D(\hat\omega),$$

dove $D(\hat\omega)$ è l'ampiezza desiderata: reale, non negativa e periodica di periodo $2\pi$ (per un passa-basso, $1$ per $|\hat\omega|\le\hat\omega_0$ e $0$ altrimenti). Il fattore di fase è il ritardo che garantisce la fase lineare. Poiché $D$ è periodica, ha una serie di Fourier (Serie di FourierUn segnale periodico di periodo $T_p$ si scrive come somma di esponenziali alle frequenze multiple della fondamentale $F=1/T_p$: $s(t)=\sum_nS_ne^{i2\pi nFt}$, con $S_n=\frac1{T_p}\int_{T_p}s(t)e^{-i2\pi nFt}dt$. Si basa sull'ortogonalità degli esponenziali su un periodo. $S_0$ è il valor medio; per segnali reali $S_{-n}=S_n^$ e si passa alla forma con coseni e seni; vale il teorema di Parseval $P=\sum|S_n|^2$. La convoluzione ciclica diventa il prodotto $T_pX_nY_n$ dei coefficienti. Le somme troncate presentano il fenomeno di Gibbs vicino ai salti.Serie di Fourier →) che coincide con una DTFT (Trasformata di Fourier a tempo discretoLa trasformata di Fourier di un segnale discreto è $S(f)=\sum_nT,s(nT),e^{-i2\pi fnT}$: una funzione continua e periodica in $f$ di periodo $F_p=1/T$. L'antitrasformata è l'integrale su un periodo, $s(nT)=\int_0^{F_p}S(f)e^{i2\pi fnT}df$. Le regole sono quelle del caso continuo (traslazione, convoluzione $\leftrightarrow$ prodotto, Parseval $\sum T|s|^2=\int_0^{F_p}|S|^2$), con incremento e somma corrente al posto di derivata e integrale. Per segnali reali $S(f)=S^(-f)$, quindi basta $[0,F_p/2]$. Esempi: $\delta\to1$, $1\to\delta_{F_p}$, rect $\to$ sinc periodico, $a^n\mathbf 1_0\to\frac T{1-ae^{-i2\pi fT}}$.Trasformata di Fourier a tempo discreto →):

$$D(\hat\omega)=\sum_{n=-\infty}^{+\infty}h_d[n]\,e^{-j\hat\omega n},\qquad h_d[n]=\frac1{2\pi}\int_{-\pi}^{\pi}D(\hat\omega)\,e^{j\hat\omega n}\,d\hat\omega .$$

Quindi $h_d[n]$ è la risposta impulsiva (infinita) del filtro ideale con ampiezza $D$ e fase nulla. Siccome $D$ è reale e pari, $h_d$ è reale e pari. Confronto: la risposta di un FIR è una somma finita $\sum_{n=0}^{N}h[n]e^{-j\hat\omega n}$, quella desiderata una somma infinita. Tutta la progettazione sta nel scegliere $h[n]$ in modo che la somma finita approssimi quella infinita. Il filtro ha la fase di $H_d$ per costruzione; resta da approssimare il modulo.

Teorema (i tre metodi).

1. Finestre (windowing, o troncamento della serie di Fourier): si usa una somma parziale della serie di $D$.
2. Campionamento in frequenza (frequency sampling): $h$ è l'IDFT di campioni di $D$, quindi $H$ coincide con $D$ solo su quei campioni.
3. Minimax (algoritmo di Parks-McClellan): si cercano i coefficienti che minimizzano il massimo dell'errore pesato $W(\hat\omega)\,|D(\hat\omega)-\bar H(\hat\omega)|$; il peso $W$ permette errori diversi nelle diverse bande. I filtri risultanti sono equiripple.

Metodo delle finestre

Il troncamento

Per un tipo I si lavora con il filtro non causale simmetrico $\bar h[n]$, $-N/2\le n\le N/2$, che ha risposta $\bar H(\hat\omega)=\sum_{n=-r}^{r}\bar h[n]e^{-j\hat\omega n}$ con $r=N/2$. L'idea è prendere i coefficienti centrali della serie di $D$:

$$\bar h[n]=h_d[n]\,w_R[n],\qquad w_R[n]=\begin{cases}1,&|n|\le N/2\\0,&\text{altrove.}\end{cases}$$

Il filtro causale si ottiene ritardando di $N/2$ campioni:

$$h[n]=\bar h\!\left[n-\tfrac N2\right]=h_d\!\left[n-\tfrac N2\right]w_R\!\left[n-\tfrac N2\right],\qquad n=0,\dots,N.$$

Il ritardo non cambia il modulo, aggiunge soltanto il fattore $e^{-j\hat\omega N/2}$ e rende il filtro realizzabile. $h$ è simmetrica, $h[n]=h[N-n]$, perché lo è $h_d$.

Esempio: il passa-basso ideale

Per $D(\hat\omega)=1$ in $|\hat\omega|\le\hat\omega_0$ e $0$ altrove ($\hat\omega_0=2\pi f_0/F_s$ è la frequenza di taglio normalizzata):

$$h_d[n]=\frac1{2\pi}\int_{-\hat\omega_0}^{\hat\omega_0}e^{j\hat\omega n}d\hat\omega=\frac{\sin(\hat\omega_0 n)}{\pi n}=\frac{\hat\omega_0}{\pi}\operatorname{sinc}\!\Big(\frac{\hat\omega_0 n}{\pi}\Big),\qquad h_d[0]=\frac{\hat\omega_0}{\pi}.$$

(Con $\operatorname{sinc}(x)=\frac{\sin\pi x}{\pi x}$, come nel corso e in numpy.sinc.) Il filtro causale è $h[n]=\frac{\hat\omega_0}{\pi}\operatorname{sinc}\big(\frac{\hat\omega_0}{\pi}(n-\frac N2)\big)w[n]$. Il coefficiente più grande è quello centrale, $h[N/2]=\hat\omega_0/\pi\le1$ qualunque sia il taglio, e più la banda è stretta più tutti i coefficienti sono piccoli (guadagno 1 in continua spalmato su molti campioni).

Esempio. $N=20$, $\hat\omega_0=\pi/3$, finestra rettangolare: $h[10]=\frac13=0{,}3333$ e, per simmetria, $h[10\pm1]=0{,}2757$, $h[10\pm2]=0{,}1378$, $h[10\pm3]=0$ (è $\sin(\pi)=0$), $h[10\pm4]=-0{,}0689$, $h[10\pm5]=-0{,}0551$, $h[10\pm6]=0$, $h[10\pm7]=0{,}0394$, $h[10\pm8]=0{,}0345$, $h[10\pm9]=0$, $h[10\pm10]=-0{,}0276$. La somma dei coefficienti è $1{,}0048\approx1$ (il guadagno in continua).

Perché il troncamento è il migliore possibile, e perché dà Gibbs

Teorema (ottimalità). Fra tutti i filtri FIR a fase lineare di ordine $N=2r$, la somma parziale $\bar H(\hat\omega)=\sum_{n=-r}^{r}h_d[n]e^{-j\hat\omega n}$ minimizza l'errore quadratico medio $E^2=\frac1{2\pi}\int_{-\pi}^{\pi}\big(D(\hat\omega)-\bar H(\hat\omega)\big)^2d\hat\omega$.

Dimostrazione. Per Parseval $E^2=\sum_{n=-\infty}^{\infty}(h_d[n]-\bar h[n])^2=\sum_{|n|\le r}(h_d[n]-\bar h[n])^2+\sum_{|n|>r}h_d[n]^2$, con $\bar h[n]=0$ fuori da $|n|\le r$. L'ultima somma non dipende dalla scelta di $\bar h$; la prima è non negativa e vale zero se e solo se $\bar h[n]=h_d[n]$ per $|n|\le r$. $\square$

Il troncamento è dunque ottimo in media quadratica, ma non è ottimo nel caso peggiore: guardando il massimo dell'errore ci sono filtri migliori (è il minimax). Per capire cosa succede al modulo, si osserva che moltiplicare per $w_R[n]$ in tempo equivale a convolvere in frequenza (Proprietà della trasformata di FourierCon le proprietà (linearità, simmetrie, ritardo $\leftrightarrow e^{-j\omega t_0}$, modulazione $\leftrightarrow$ traslazione in frequenza, scala, dualità, convoluzione $\leftrightarrow$ prodotto, Parseval $E=\frac1{2\pi}\int|X|^2$, derivata $\leftrightarrow j\omega$, moltiplicazione per $t\leftrightarrow j,d/d\omega$, integrazione) quasi tutte le trasformate si ottengono da poche coppie base senza integrare.Proprietà della trasformata di Fourier →):

$$\bar H(\hat\omega)=\frac1{2\pi}\int_{-\pi}^{\pi}D(\theta)\,W_R(\hat\omega-\theta)\,d\theta,\qquad W_R(\hat\omega)=\sum_{n=-r}^{r}e^{-j\hat\omega n}=\frac{\sin(L\hat\omega/2)}{\sin(\hat\omega/2)},$$

con $L=2r+1=N+1$ la lunghezza della risposta impulsiva (è il nucleo di Dirichlet: la serie geometrica di $L$ termini). $W_R$ vale $L$ in $0$ e si annulla in $\hat\omega=\frac{2\pi}{L}k$, $k\ne0$. Il tratto fra i due primi zeri $\pm\frac{2\pi}L$ è il lobo principale (main lobe), largo $\frac{4\pi}{L}$; gli altri sono i lobi laterali (side lobes). Il primo lobo laterale vale circa il 22% del principale (cioè $-13{,}3$ dB) qualunque sia $L$: aumentando $L$ i lobi si restringono ma non si abbassano.

Il risultato dell'integrale è la media pesata di $D$ con il nucleo $W_R$ centrato in $\hat\omega$. Dove $D$ è costante o varia lentamente il nucleo vede un valore quasi costante e $\bar H\simeq D$; vicino a un salto di $D$ i lobi laterali, che hanno segno alterno, entrano e escono dal salto e producono oscillazioni (fenomeno di Gibbs).

Esempio. Con $\hat\omega_0=\pi/3$ il massimo di $\bar H$ è $1{,}0706$ per $L=21$, $1{,}0928$ per $L=61$, $1{,}0908$ per $L=151$ e $1{,}0900$ per $L=401$: il sovraelongo converge a circa il 9% del salto (con minimo $-0{,}089$ nella banda oscura), come si vede nel grafico. L'errore in banda passante e in banda oscura è circa $\delta_p\simeq\delta_s\simeq0{,}09$, cioè $A_s\simeq21$ dB, e non dipende da $L$. Quello che dipende da $L$ è la larghezza della transizione, circa $0{,}9\cdot\frac{2\pi}{L}$, e la rapidità delle oscillazioni.

Altre finestre

Per ridurre i lobi laterali, al posto di $w_R$ si usa una finestra $w[n]$ più regolare, che scende dolcemente verso i bordi. Il prezzo è un lobo principale più largo, cioè una transizione più larga. Le espressioni qui sono quelle non causali di lunghezza $L=N+1$ dispari, $|n|\le N/2$; nel codice dei laboratori le stesse finestre sono scritte con $n=0,\dots,N$, per esempio la Hamming $w[n]=0{,}54-0{,}46\cos\frac{2\pi n}{N}$.

(Hann e Hamming sono "coseni rialzati": un coseno più una costante.)

I loro spettri sono combinazioni di nuclei di Dirichlet traslati: i termini coseno spostano il nucleo in $\pm\frac{2\pi}{N}$ e, se i pesi sono scelti bene, i lobi laterali dei tre nuclei si cancellano a vicenda.

Proprietà delle finestre fisse e valori tipici (tabella delle dispense; i valori di larghezza del lobo principale e di livello del lobo laterale sono stati ricalcolati numericamente per $L=51$):

Tre fatti generali per tutte le finestre fisse:

1. la larghezza della transizione è inversamente proporzionale alla lunghezza: $\dfrac{\hat\omega_s-\hat\omega_p}{2\pi}\simeq\dfrac{\alpha}{L}$, con $\alpha$ che dipende solo dalla finestra (colonna "banda di transizione");
2. gli errori sono quasi uguali nelle due bande, $\delta_p\simeq\delta_s$, e dipendono solo dal tipo di finestra, non da $L$;
3. entrambi non dipendono dalla frequenza di taglio $\hat\omega_0=\frac{\hat\omega_p+\hat\omega_s}2$.

Il punto 2 ha una conseguenza pratica: con le finestre fisse non si possono scegliere $\delta_p$ e $\delta_s$ separatamente. Queste cifre sono regole pratiche (si vedrà nell'esempio sotto che per soddisfare davvero 50 dB serve un ordine un poco più grande di quello dato dalla tabella).

La finestra di Kaiser

La finestra di Kaiser ha un parametro $\beta$ che permette di scegliere l'attenuazione:

$$w[n]=\frac{I_0\!\Big(\beta\sqrt{1-\big(\frac{2n}{N}\big)^2}\Big)}{I_0(\beta)},\quad|n|\le\frac N2,\qquad I_0(x)=1+\sum_{k=1}^{\infty}\Big(\frac{(x/2)^k}{k!}\Big)^2,$$

dove $I_0$ è la funzione di Bessel modificata di ordine zero. Con $\beta=0$ si ottiene la rettangolare ($A_s=21$ dB). Date $A_s$ in dB e la banda di transizione $\hat\omega_s-\hat\omega_p$, si usano le formule empiriche

$$\beta=\begin{cases}0{,}1102\,(A_s-8{,}7),&A_s>50\\0{,}5842\,(A_s-21)^{0{,}4}+0{,}07886\,(A_s-21),&21<A_s\le50\\0,&A_s\le21\end{cases}\qquad N\ge\frac{A_s-8}{14{,}36\,\frac{\hat\omega_s-\hat\omega_p}{2\pi}} .$$

Esempio. Per $A_s=50$ dB e $\hat\omega_s-\hat\omega_p=0{,}125\pi$ (cioè $0{,}0625$ in unità di $2\pi$): $\beta=0{,}5842\cdot29^{0{,}4}+0{,}07886\cdot29=2{,}2466+2{,}2869=4{,}534$ con la seconda formula (con la prima, che vale sopra 50 dB, si avrebbe $4{,}551$: i due valori coincidono quasi sul confine) e $N\ge\frac{42}{14{,}36\cdot0{,}0625}=46{,}8$, quindi $N=48$ (pari, per il tipo I).

Procedura di progetto

1. Si calcola la risposta impulsiva ideale $h_d[n]$ (con taglio $\hat\omega_0=\frac{\hat\omega_p+\hat\omega_s}2$).
2. Si sceglie il tipo di finestra dalla $A_s$ richiesta (la meno costosa che arriva a $A_s$) e l'ordine $N$ dalla larghezza della transizione: $L=N+1\ge\alpha\,\frac{2\pi}{\hat\omega_s-\hat\omega_p}$.
3. Si calcola $\bar h[n]=h_d[n]\,w[n]$ e, per un filtro causale, $h[n]=\bar h[n-N/2]$.
4. Si calcola $H(e^{j\hat\omega})$ e si verifica che le specifiche siano rispettate; se no si aumenta $N$ o si cambia finestra.

Vantaggi del metodo: è semplicissimo, intuitivo e conserva gli zeri di $h_d$ (utile per i filtri di interpolazione, vedi Elaborazione multirate - decimazione e interpolazioneL'elaborazione multirate cambia la frequenza di campionamento di un segnale. Interpolazione per $L$ (da $F_s$ a $LF_s$): l'espansore inserisce $L-1$ zeri fra i campioni, $W(\hat\omega)=X(L\hat\omega)$, lo spettro si ripete $L$ volte nell'asse più ampio; un passabasso con taglio $\pi/L$ e guadagno $L$ elimina le immagini, e con il sinc ideale $h[n]=\operatorname{sinc}(n/L)$ i campioni originali sono conservati ($y[nL]=x[n]$); in pratica FIR a finestre o interpolatore lineare (triangolo). Decimazione per $M$ (da $F_s$ a $F_s/M$): $y[n]=x[nM]$, $Y(\hat\omega)=\frac1M\sum_{k=0}^{M-1}X\big(\frac{\hat\omega-2\pi k}M\big)$, repliche sovrapposte (aliasing) se la banda supera $\pi/M$: si fa precedere da un passabasso con taglio $\pi/M$. Conversione razionale $F_s'=\frac LMF_s$ ($L,M$ coprimi): espansore $L$, un solo passabasso di taglio $\min(\frac\pi L,\frac\pi M)$ e guadagno $L$, decimatore $M$; $y[n]=\sum_kx[k],g[nM-kL]$ non è una convoluzione e il sistema è periodicamente tempo-invariante. Realizzare il tutto in più stadi (fattori piccoli prima per l'interpolazione, grandi prima per la decimazione) e con strutture polifase riduce molto il costo.Elaborazione multirate - decimazione e interpolazione →). Svantaggi: i filtri non sono ottimi (le stesse specifiche si ottengono con ordini minori con il minimax) e gli errori $\delta_p$, $\delta_s$ sono simili, il che non è sempre ciò che serve. Nella pratica le finestre più usate sono Hamming (per la semplicità) e Kaiser (per le prestazioni).

Esempio completo e confronto fra i metodi

Specifiche: $F_s=8$ kHz, banda passante fino a $f_p=1$ kHz, banda oscura da $f_s=1{,}5$ kHz, $A_s=50$ dB ($\delta_p=\delta_s=0{,}00316$). Normalizzando $\hat\omega=2\pi f/F_s$: $\hat\omega_p=0{,}25\pi$, $\hat\omega_s=0{,}375\pi$, $\hat\omega_s-\hat\omega_p=0{,}125\pi$ e taglio $\hat\omega_0=0{,}3125\pi$.

• Rettangolare: non va bene, perché $\delta\simeq0{,}09$ qualunque sia $N$.
• Hamming: $L\simeq3{,}3\cdot\frac{2\pi}{0{,}125\pi}=52{,}8$, quindi $L=53$, $N=52$. Controllo numerico: $\delta_p=0{,}00313$, $\delta_s=0{,}00330$ ($A_s=49{,}6$ dB), appena fuori specifica; con $N=54$ si ottiene $\delta_p=0{,}00258$, $\delta_s=0{,}00220$ e la specifica è soddisfatta. La tabella è una regola pratica, non una garanzia: il passo 4 serve.
• Hann e Blackman: con la ricerca dell'ordine minimo servono $N=76$ (Hann) e $N=74$ (Blackman).
• Kaiser: $\beta=4{,}55$, $N=48$ soddisfa le specifiche ($\delta_p=0{,}00299$, $\delta_s=0{,}00260$); con $N=46$ si ha $0{,}0047$ e non basta.
• Minimax (vedi sotto): $N=42$ ($\delta_p=\delta_s=0{,}00279$); con $N=40$ è $0{,}0043$.

I tre filtri hanno tutti $\delta\le0{,}00316$ ma ordini $54,48,42$: a parità di specifica il minimax è il più economico (e il filtro è più corto, quindi più veloce da calcolare e con minor ritardo, $N/2$ campioni).

Campionamento in frequenza

L'idea è molto semplice: si campiona la risposta desiderata in $L=N+1$ punti equispaziati $\hat\omega_k=\frac{2\pi}{L}k$, $k=0,\dots,N$, si prendono questi $L$ valori come una DFT (Trasformata di Fourier discreta (DFT) e convoluzione circolareLa DFT di N campioni x[0..N-1] è X[k] = sum x[n] e^{-j2pi kn/N}, k = 0..N-1 (IDFT: x[n] = (1/N) sum X[k] e^{j2pi kn/N}); vede il segnale come periodico di periodo N e dà N campioni equispaziati della DTFT, cioè X(z) valutata sulle radici N-esime dell'unità. Se x ha durata M <= N non c'è aliasing temporale e lo zero-padding (N > M) infittisce solo i campioni della stessa DTFT. Proprietà: traslazione circolare, simmetria coniugata X[k] = X*[N-k], prodotto = convoluzione circolare. La convoluzione circolare coincide con la lineare solo se N >= Dx + M - 1; altrimenti la coda si ripiega sull'inizio. La FFT calcola la DFT con (N/2) log2 N moltiplicazioni invece di N^2.Trasformata di Fourier discreta (DFT) e convoluzione circolare →) e la risposta impulsiva è l'IDFT:

$$H[k]=H_d(e^{j\hat\omega_k})=e^{-j\hat\omega_k N/2}D(\hat\omega_k),\qquad h[n]=\mathrm{IDFT}\{H[k]\}.$$

Il filtro ottenuto ha esattamente $H(e^{j\hat\omega_k})=H_d(e^{j\hat\omega_k})$ nei punti campione, ma fra un campione e l'altro non c'è controllo. Equivale a $h[n]=\sum_{m}h_d\!\left[n-\tfrac N2+mL\right]$, $0\le n\le N$, cioè la risposta ideale (infinita) ripiegata su sé stessa con periodo $L$, perché campionare in frequenza significa ripetere in tempo (Trasformata di Fourier discreta (DFT) e convoluzione circolareLa DFT di N campioni x[0..N-1] è X[k] = sum x[n] e^{-j2pi kn/N}, k = 0..N-1 (IDFT: x[n] = (1/N) sum X[k] e^{j2pi kn/N}); vede il segnale come periodico di periodo N e dà N campioni equispaziati della DTFT, cioè X(z) valutata sulle radici N-esime dell'unità. Se x ha durata M <= N non c'è aliasing temporale e lo zero-padding (N > M) infittisce solo i campioni della stessa DTFT. Proprietà: traslazione circolare, simmetria coniugata X[k] = X*[N-k], prodotto = convoluzione circolare. La convoluzione circolare coincide con la lineare solo se N >= Dx + M - 1; altrimenti la coda si ripiega sull'inizio. La FFT calcola la DFT con (N/2) log2 N moltiplicazioni invece di N^2.Trasformata di Fourier discreta (DFT) e convoluzione circolare →): l'aliasing nel tempo della risposta impulsiva ideale è l'errore (verificato numericamente per l'esempio sotto).

Esempio. $N=16$ ($L=17$), $\hat\omega_0=\pi/3$. I campioni cadono in $\hat\omega_k=\frac{2\pi k}{17}$, cioè in $0$, $0{,}118\pi$, $0{,}235\pi$, $0{,}353\pi$, $0{,}471\pi,\dots$: i campioni $k=0,1,2$ sono in banda passante (il taglio è $\pi/3=0{,}333\pi$), già $k=3$ è oltre il taglio. Si pone $D=1$ per $k=0,1,2$ e $k=15,16$ (i simmetrici), $0$ negli altri. Si ottiene $h[0]=0{,}0529$, $h[1]=0{,}0112$, $h[2]=-0{,}0443$, $h[3]=-0{,}0734$, $h[4]=-0{,}0460$, $h[5]=0{,}0404$, $h[6]=0{,}1566$, $h[7]=0{,}2555$, $h[8]=0{,}2941$ (e simmetrici), con somma $1$. L'errore è circa $0{,}10$ in banda passante e $0{,}11$ in banda oscura (dai campioni in poi).

Per attenuare le oscillazioni si ammorbidisce il salto con un campione di transizione: se si pone $D=0{,}5$ nel campione $k=3$ (che sta nella banda di transizione) gli errori scendono a $0{,}030$ (banda passante) e $0{,}032$ (banda oscura), e con $0{,}41$ la banda oscura arriva a $0{,}009$ (con $0{,}043$ in banda passante). Il prezzo è una transizione più larga: la risposta scende da $1$ a $0$ attraverso due intervalli fra campioni ($k=2\to3\to4$) invece che uno. Nel grafico i punti sono i campioni imposti.

Minimax e algoritmo di Parks-McClellan

Il metodo minimax (o Chebyshev approximation) è ottimale nel caso peggiore: dati l'ampiezza desiderata $D(\hat\omega)$ e un peso positivo $W(\hat\omega)$, entrambi continui su un insieme compatto $I\subset[0,\pi]$ (per un passa-basso $I=[0,\hat\omega_p]\cup[\hat\omega_s,\pi]$: la transizione è esclusa), si cerca

$$\min_{h}\ \max_{\hat\omega\in I}\ W(\hat\omega)\,\big|D(\hat\omega)-\bar H(\hat\omega)\big|.$$

Il peso rende più o meno importante l'errore in una banda: per avere $\delta_s$ dieci volte più piccolo di $\delta_p$ si pone $W=\frac1{10}$ in banda passante e $W=1$ in banda oscura.

Teorema (alternanza). Sia $\bar H(\hat\omega)=\sum_{n=0}^{r}p_n\cos(n\hat\omega)$ una combinazione di $r+1$ coseni e $E(\hat\omega)=W(\hat\omega)(D(\hat\omega)-\bar H(\hat\omega))$. $\bar H$ è la soluzione del problema minimax se e solo se $E$ ha almeno $r+2$ alternanze, cioè punti $\hat\omega_0<\dots<\hat\omega_{r+1}$ in $I$ con $E(\hat\omega_i)=-E(\hat\omega_{i+1})$ e $|E(\hat\omega_i)|=\max_I|E|$.

Una combinazione di $r+1$ coseni ha $r+1$ gradi di libertà: serve una alternanza in più. Le alternanze stanno ai bordi delle bande o nei massimi e minimi dell'errore. La conseguenza è che l'errore ha tutti i picchi della stessa altezza: la risposta è equiripple, cioè con ondulazioni di ampiezza uguale. L'algoritmo di Remez parte da una stima dei punti $\hat\omega_i$, risolve il sistema lineare di $r+2$ equazioni $W(\hat\omega_i)\big(D(\hat\omega_i)-\sum_np_n\cos(n\hat\omega_i)\big)=(-1)^i\varepsilon$ nelle incognite $p_0,\dots,p_r,\varepsilon$, aggiorna i punti e ripete, e converge in poche iterazioni.

Esempio. $N=10$ ($r=5$, quindi servono $7$ alternanze), banda passante $[0,0{,}3\pi]$, banda oscura $[0{,}4\pi,\pi]$, pesi uguali. scipy.signal.remez dà $h=\{-0{,}0338,-0{,}1353,-0{,}0002,0{,}1237,0{,}2840,0{,}3493,\dots\}$ (simmetrica) con errore massimo $0{,}1742$ in entrambe le bande. I punti di errore massimo sono $\hat\omega/\pi=0;\,0{,}191;\,0{,}3;\,0{,}4;\,0{,}514;\,0{,}734;\,1$ con segni $-,+,-,+,-,+,-$: sono proprio $7$, di cui $4$ ai bordi delle bande. Con peso $1:10$ (banda oscura più importante) gli errori diventano $0{,}4926$ e $0{,}0494$, in rapporto $9{,}97\simeq10$.

Stima dell'ordine

Per un passa-basso, con $\delta_p$ e $\delta_s$ non necessariamente uguali:

$$N\simeq\frac{-10\log_{10}(\delta_p\delta_s)-13}{14{,}6\,\frac{\hat\omega_s-\hat\omega_p}{2\pi}} .$$

L'ordine è inversamente proporzionale alla larghezza della transizione, non dipende dal taglio e dipende da $\delta_p\delta_s$ (raddoppiare $\delta_p$ e dimezzare $\delta_s$ lascia $N$ invariato). Vale anche per i passa-alto. Per passa-banda ed elimina-banda ci sono due transizioni: si calcolano $N_1$ e $N_2$ e l'ordine sta fra $\max(N_1,N_2)$ e $N_1+N_2$, in pratica più vicino al massimo.

Esempio. Con $\delta_p=\delta_s=0{,}00316$: $-10\log_{10}(10^{-5})-13=37$ e $N\simeq\frac{37}{14{,}6\cdot0{,}0625}=40{,}5$. Il valore effettivo è $N=42$ (con $N=40$ l'errore è $0{,}0043$): la stima è approssimata e va verificata.

Confronto

Domande d'esame

1. Descrivi la tecnica del windowing per il progetto di un filtro FIR e commenta come la scelta della finestra e della sua lunghezza incide sul risultato. Traccia: ampiezza desiderata $D$ periodica e serie di Fourier, coefficienti $h_d[n]$ del passa-basso ideale; troncamento a $|n|\le N/2$ = finestra rettangolare, ritardo di $N/2$ per la causalità e fase lineare; effetto in frequenza = convoluzione con il nucleo di Dirichlet, lobo principale (larghezza $\frac{4\pi}{L}$) e lobi laterali (22%); Gibbs 9% indipendente da $L$; finestre rastremate: lobi laterali più bassi, lobo principale più largo, $\frac{\hat\omega_s-\hat\omega_p}{2\pi}\simeq\frac\alpha L$ e $\delta_p\simeq\delta_s$ fissati dal tipo; stesso compromesso risoluzione/leakage nella STFT (Trasformata di Fourier a tempo breve e spettrogrammaPer un segnale lungo il cui contenuto in frequenza cambia nel tempo si calcolano tante DFT brevi: X[k, l] = sum_{n=0}^{L-1} w[n] x[l+n] e^{-j2pi kn/N}, con finestra di analisi w di lunghezza L, istante l = m R (passo R, di solito L/2) e DFT a N >= L punti. Il modulo (o i dB) di X[k, l] come immagine tempo-frequenza è lo spettrogramma. La finestra impone un compromesso: lunga = ottima risoluzione in frequenza (Delta w circa 8pi/L per Hann) ma transizioni temporali sfocate su circa L/2 campioni; corta = buona localizzazione nel tempo ma righe vicine fuse. Si provano più lunghezze.Trasformata di Fourier a tempo breve e spettrogramma →, Analisi spettrale con la DFT - finestre e leakageAnalizzare lo spettro di un segnale con la DFT vuol dire osservarne solo L campioni (finestra rettangolare) e campionare in frequenza la DTFT. Una sinusoide non dà una riga ma il nucleo di Dirichlet della finestra centrato in w0: lobo principale largo 4pi/L e lobi laterali (il primo a -13.3 dB). Due toni più vicini della mezza larghezza del lobo non si separano (risoluzione, che dipende solo da L), e i lobi laterali di un tono mascherano i toni deboli (leakage). La DFT campiona la DTFT a 2pi k/N: l'ampiezza è corretta solo se il tono cade su un bin (periodi interi nella finestra), altrimenti si perde fino a 3.9 dB; lo zero-padding infittisce i campioni ma non migliora la risoluzione. Le finestre rastremate (Hann, Hamming, Blackman) abbassano i lobi laterali (-31, -42, -58 dB) a prezzo di un lobo principale più largo.Analisi spettrale con la DFT - finestre e leakage →).
2. Confronta i metodi delle finestre, del campionamento in frequenza e minimax. Traccia: tabella sopra; ottimalità e alternanza, $r+2$ alternanze, equiripple, pesi; ordine inferiore a parità di specifiche (esempio $54$ contro $42$); campionamento in frequenza esatto nei campioni, campione di transizione.
3. Dimostra che il troncamento della serie di Fourier minimizza l'errore quadratico fra i filtri FIR di ordine $N$. Traccia: Parseval, separazione dei termini $|n|\le r$ e $|n|>r$, la seconda somma non dipende da $h$; ma non minimizza l'errore massimo (per questo esiste il minimax).