Capitolo 4 Progetto di filtri FIR (II) -...

Appunti diAppunti diElaborazione numerica dei segnaliElaborazione numerica dei segnaliCapitolo 4 Capitolo 4 - - Progetto di filtri FIR (II)(II)

Introduzione ............................................................................................... 1Dettagli analitici ......................................................................................... 4Procedure di ottimizzazione ....................................................................... 6

Algoritmo di Remez............................................................................. 11Procedura di Parks e McClellan..................................................... 13

Procedura Remez di MATLAB................................................................. 22Conclusioni sui metodi di progetto di filtri FIR ........................................ 23

INTRODUZIONE

Così come nel metodo delle finestre, anche il progetto di filtri (numerici) FIR basati sul metododel campionamento in frequenza ha come punto di partenza una funzione di trasferimentoobbiettivo , che indichiamo Hd(f). Per fissare le idee, consideriamo ancora una volta la funzione ditrasferimento del filtro passa-basso ideale. Il passo successivo consiste nell’effettuare uncampionamento di tale Hd(f), ossia sostanzialmente nell’individuare i valori che tale funzione ditrasferimento assume in un set preassegnato di frequenze fk equispaziate (ovviamente semprenell’intervallo non ambiguo [-fC/2,fC/2], dato che la cosa si ripete identica per gli altri periodi):

f (Hz)

Hd(f)

fC/2-fC/2

Campionamento (a frequenze equispaziate) della funzione di trasferimento obbiettivo, che è statascelta coincidente con quella del filtro passa-basso ideale

Il passo ulteriore consiste ora nell’imporre che il filtro H(f) realmente implementato abbia, nellefrequenze fk gli stessi valori ottenuti dal campionamento di Hd(f). Per tutte le altre frequenze, invece,il vincolo è ancora una volta quello per cui la H(f) ottenuta deve rispettare la maschera preassegnata:

Appunti di “Elaborazione numerica dei segnali” - Capitolo 4

Autore: Sandro Petrizzelli2

f (Hz)bandapassante

banda ditransizione

bandaarrestata

Classica maschera per un filtro passa-basso.

Affinché la funzione implementata H(f) rispetti i valori campionati di Hd(f), è evidentementenecessario che la corrispondente funzione di risposta all’impulso presenti gli opportuni terminisinusoidali: infatti, ciascun campione in frequenza corrisponde ad una precisa sinusoide (complessa),per cui la risposta all’impulso dovrà contenere le sinusoidi corrispondenti ai vari campioni. Adesempio, supponiamo che il campione di Hd(f) prelevato a frequenza fk valga 1:

( ) 1fH kd =

Allora, la funzione di risposta all’impulso dovrà contenere un termine tf2j ke π .Sommando tutti i termini sinusoidali del tipo tf2j ke π , otterremo la risposta all’impulso desiderata,

che indichiamo con hT(t). Tuttavia, siamo ancora nel campo tempo-continuo, mentre a noi interessala sequenza hT(n) corrispondente al campionamento di hT(t) a passo T=1/fC. Dobbiamo perciòeffettuare un successivo campionamento nei tempi (a passo T), il quale produce unaperiodicizzazione in frequenza. Tale periodicizzazione non ci dà problemi di aliasing perché, peripotesi, la funzione di trasferimento si estende tra -fC/2 e +fC/2.

Al contrario, abbiamo il solito problema che la funzione di risposta all’impulso non può averelunghezza infinita, per cui dobbiamo ancora una volta finestrare nei tempi. Indichiamo allora conh(n) la sequenza ottenuta finestrando hT(n) su un intervallo di durata NT. Per quanto riguarda lafinestra usata immaginiamo di usare una semplice finestra rettangolare: vedremo che, al contrario delmetodo delle finestre, questa è la scelta migliore che si possa fare.

Avendo finestrato nei tempi, abbiamo una convoluzione in frequenza, tra lo spettro G(f) dellafinestra e lo spettro HT(f) corrispondente ad hT(n):

)f(G*)f(H)f(H)n(g)n(h)n(h TT =→⋅=

E’ la stessa cosa vista nel metodo delle finestre, con la fondamentale differenza che, mentre inquel caso tale convoluzione ci creava dei problemi (in quanto smussava la funzione di trasferimentodesiderata), in questo caso risulta molto opportuna. Vediamo di capire perché.

Avendo campionato in frequenza la funzione obbiettivo Hd(f), è come se ci fossimo posti, comeobbiettivo, una funzione di trasferimento fatta nel modo seguente:

Progetto di filtri FIR (parte II): metodo del campionamento in frequenza


f (Hz)

HT(f)

Conosciamo perciò i valori che la funzione di trasferimento finale H(f) dovrà avere in determinatefrequenze fk, mentre non sappiamo niente per i valori delle altre frequenze. Abbiamo perciò unproblema di interpolazione dei valori mancanti (che dovranno sempre rientrare nella maschera),interpolazione che deve portarci ad un andamento quanto più possibile simile ad un rettangolo.Allora, nel momento in cui convolviamo G(f) e HT(f), di fatto effettuiamo una interpolazione (infrequenza), dove la funzione interpolante è proprio G(f). Se abbiamo scelto una finestra rettangolare,la funzione interpolante è la migliore che ci sia, ossia il sin(f)/f. Questo dimostra che, con questometodo, la finestratura giunge quanto mai opportuna.

Si pone un ulteriore aspetto, legato sempre all’interpolazione: abbiamo osservato che il nostroobbiettivo è quello di ottenere, dalla convoluzione tra G(f) e HT(f), una funzione H(f) quanto piùpossibile simile ad un rettangolo; in realtà, a noi interessa solo che la H(f) ottenuta rispetti i vincoliimposti dalla maschera; ancora meglio, mentre desideriamo che la H(f), nella banda passante e nellabanda arrestata somigli il più possibile al rettangolo, ci interessa poco quello che accade all’internodella banda di transizione, purché tale banda venga rispettata. Possiamo allora pensare diricercare delle procedure di ottimizzazione che, a spese di quelloche accade nella banda di transizione, ottimizzino l’andamento diH(f) nella banda passante ed in quella arrestata. Ovviamente, perché questeprocedure siano “libere di agire” nella banda di transizione, non dobbiamo fissare i campioni relativia tale banda, per cui, nel caso considerato prima, il set di campioni da considerare sarà il seguente:

f (Hz)fC/2-fC/2

Nel caso da noi considerato, c’era un solo campione per la banda di transizione (sia a frequenzepositive sia a frequenze negative), per cui lo eliminiamo, in modo che le procedure di ottimizzazione nonabbiano vincoli sulla banda di transizione, salvo ovviamente quelli imposti dalla maschera. In generale, i

campioni nella banda di transizione sono più di uno e vanno tutti eliminati



DETTAGLI ANALITICI

Facciamo una serie di passaggi analitici a supporto delle considerazioni essenzialmente qualitativefatti fino ad ora.

Abbiamo detto che il primo passo, nel metodo di campionamento in frequenza, è quello dicampionare la funzione obbiettivo Hd(f) in un set di frequenze (o pulsazioni) equispaziate. Seindichiamo con M il numero di campioni che vogliamo prendere (comprendendo sia quelli afrequenza positiva sia quelli a frequenza negativa), le pulsazioni da considerare sono le seguenti:

=π

=ωpari Mper 1-

2

M0,1,....,

dispari Mper 2

1-M0,1,....,

k M

2kk

(abbiamo fatto riferimento, per semplicità di notazioni, alle pulsazioni). E’ ovvio che l’ultimapulsazione differisce a seconda che M venga scelto pari o dispari.

L’obbiettivo è di determinare una funzione di risposta all’impulso h(n) cui corrisponda unafunzione di trasferimento H(ω) che “onori” i campioni appena definiti, il che significa cheH(ωk)=Hd(ωk), e che rispetti i vincoli della maschera.

Per definizione, possiamo esprimere H(ω) come trasformata di h(n):

∑−

=

ω−=ω1M

0n

nTje)n(h)(H

In questa espressione, abbiamo evidentemente considerato una sequenza h(n) di lunghezza M ecausale (cioè con il primo campione non nullo posizionato in n=0). Il nostro obbiettivo è calcolareproprio i campioni che costituiscono h(n).

Calcoliamo adesso H(ω) nelle pulsazioni ωk:

∑∑−

=

π−−

=

ω− ==ω1M

0n

nTM

2jk1M

0n

nTjk e)n(he)n(h)(H k

C’è da fissare la variabilità dell’indice k. Per quanto detto prima, ci interessano M/2 campioni afrequenza positiva e altrettanti a frequenza negativa. Ma se noi immaginiamo i campioni H(ωk) comequelli ottenuti tramite una DFT (basata su M campioni) della sequenza h(n), sappiamo bene che ilvettore ottenuto dalla DFT contiene prima gli M/2 campioni a frequenza positiva e poi quelli afrequenza negativa. Di conseguenza, in quella relazione possiamo tranquillamente prenderek=0,1,...,M-1, rispettando quanto detto prima circa i valori delle ωk:

1-M0,1,....,k e)n(h)(H1M

0n

nTM

2jk

k ==ω ∑−

=

π−

Adesso, dobbiamo imporre che H(ωk)=Hd(ωk):

1-M0,1,....,k e)n(h)(H1M

0n

nTM

2jk

kd ==ω ∑−

=

π−



Il passo successivo è quello di invertire questa relazione, esprimendo h(n) in funzione deicampioni Hd(ωk) presi dalla funzione obbiettivo. Non è difficile effettuare questa inversione. Per

prima cosa, moltiplichiamo ambo i membri di quella espressione per l’esponenziale T

M

mk2j

eπ

, conm=0,1,...,M-1:

1-M0,1,....,m

1-M0,1,....,k ee)n(he)(H

1M

0n

TM

mk2jnT

M

2jkT

M

mk2j

kd ==

=ω ∑−

=

ππ

−π

Adesso effettuiamo, sempre su entrambi i membri, una sommatoria su indice k:

1-M0,1,....,m ee)n(he)(H1M

0k

1M

0n

TM

mk2jnT

M

2jk1M

0k

TM

mk2j

kd ==ω ∑∑∑−

=

−

=

ππ

−−

=

π

In base ad una nota proprietà degli esponenziali, la doppia sommatoria a secondo membro siriduce ad un unico termine e precisamente h(n), per cui concludiamo che

1-M0,1,....,n e)(H)n(h1M

0k

TM

mk2j

kd =ω= ∑−

=

π

Questa è dunque la relazione che ci permette di calcolare i valori di h(n) a partire dai campioni infrequenza presi dalla funzione obbiettivo Hd(ω). E’ evidente che si tratta semplicemente di calcolarela IDFT di Hd(ω).

Mettiamoci adesso in una ipotesi particolare e precisamente quella di h(n) reale (come noivogliamo che sia). Sappiamo allora che la corrispondente H(ω) soddisfa la seguente proprietà disimmetria hilbertiana:

π

−=

π

M

2)kM(H

M

2kH *

Questa condizione consente evidentemente di ridurre il numero di campioni in frequenza daprendere: anziché considerare M campioni, ne potremo prendere (M+1)/2 se M è dispari oppure M/2se M è pari. Sotto queste ipotesi, le equazioni lineari da risolvere per passare da Hd(ωk) ad h(n)risultano senz’altro semplificate.

Come si vedrà nel prossimo paragrafo, ulteriori semplificazioni si ottengono se h(n), oltre adessere reale, presenta anche una simmetria pari rispetto ad un campione centrale (per M dispari,figura a sinistra) oppure rispetto ad un istante centrale (per M pari, figura a destra):

n

h(n)

n

h(n)



Come ben sappiamo, funzioni di risposta all’impulso di questo tipo corrispondono a filtri a faselineare, del tipo e-jωkT, dove k è il numero di ritardi temporali (ciascuno di 1 passo dicampionamento) necessari affinché h(n) risulti reale e pari.

PROCEDURE DI OTTIMIZZAZIONE

Concentriamoci adesso sulle procedure di ottimizzazione: come già anticipato, esse hanno loscopo di ottenere la migliore somiglianza possibile tra H(ω) e Hd(ω)nella banda passante e nellabanda arrestata, a spese di quello che accade nella banda di transizione (bisogna comunque accertarsiche H(ω) non vada al di fuori della maschera prefissata).

Cominciamo con l’esprimere H(ω) come trasformata di Fourier (DTFT) della corrispondentefunzione di risposta all’impulso h(n):

∑−

=

ω−=ω1N

0n

nTje)n(h)(H

Dobbiamo decidere come effettuare l’ottimizzazione di H(ω). Il modo più classico è quello diminimizzare lo scarto esistente tra H(ω) e la funzione obbiettivo Hd(ω). Tale scarto è

)(H)(H)(e d ω−ω=ω

Per come abbiamo posto fino ad ora il problema, non ha tuttavia molto senso minimizzare questaquantità, mentre è più sensato riferirsi allo scarto (o errore) quadratico:

[ ]2

d2 )(H)(H)(e ω−ω=ω

Anche qui, però, c’è da considerare che e2(ω) dipende dalla frequenza, per cui risulteràverosimilmente diverso frequenza per frequenza. Si preferisce perciò calcolare la media e2(ω) sullabanda di interesse, ossia l’errore quadratico medio1:

[ ]∫∫ ωω−ω=ωω=B

2

d

B

2 d)(H)(Hd)(eE

Essendo H(ω) la trasformata di h(n), è evidente che, al variare di h(n) varia il valore di E.

∫ ∑ ω

ω−=

−

=

ω−

B

2

d

1N

0n

nTj d)(He)n(hE

Allora, la configurazione degli h(n) che minimizza il valore di E è la soluzione del nostroproblema.

1 Per comprendere a pieno questi discorsi, ci basta considerare un caso del tutto analogo, relativo però al dominio dei tempi: se

vogliamo approssimare una data funzione del tempo f(t), incognita, con una stima g(t) da noi stessi calcolata di f(t), noicerchiamo di determinare g(t) minimizzando l’errore quadratico medio tra g(t) ed f(t); non ha invece senso minimizzare nél’errore istantaneo e(t)=g(t)-f(t) né quello quadratico istantaneo, e2(t)=(g(t)-f(t))2: infatti, a noi non interessa che istante peristante lo scarto sia minimo, ma interessa che, nel complesso, cioè per tutta la durata temporale di nostro interesse, le duefunzioni risultino il più simili possibile. Questi discorsi saranno ampiamente ripresi nel capitolo sulla stima spettrale.



Fin qui, il discorso va bene per qualsiasi filtro. Noi, però, stiamo considerando filtri FIR e, inparticolare, ci fa comodo avere filtri FIR a fase lineare. Vediamo allora se questo ci può daredelle informazioni in più.

Sappiamo bene che, se il filtro FIR è a fase lineare, la funzione di risposta all’impulso h(n)presenta una simmetria pari rispetto ad un campione centrale (per N dispari) oppure rispetto ad unistante centrale (per N pari). In generale, possiamo esprimere questa simmetria con l’uguaglianza

)1nN(h)n(h −−=

Ad esempio, se N=5 (numero dispari), il campione h(0) coincide con il campione h(4) ed ilcampione h(1) coincide con h(3):

n

h(n)

Questa proprietà ci consente di semplificare l’espressione di H(ω) come trasformata di h(n):

∑−

=

ω−=ω1N

0n

nTje)n(h)(H

Infatti, consideriamo la somma dei due termini generici, in quella sommatoria, legati ai campionih(n) e h(N-n-1) simmetrici e quindi uguali:

[ ][ ] ( )nTcos)n(h2ee)n(h

eee)n(he)n(he)n(he)1nN(he)n(hnTjnTj

nTjT)1N(jnTjT)1nN(jnTjT)1nN(jnTj

ω=+=

=+=+=−−+ωω−

ω−ω−ω−−−ω−ω−−−ω−ω−

dove abbiamo tenuto conto che il termine T)1N(je −ω− determina la fase del filtro:

T2

1Nj

R e)(H)(H−

ω−ω=ω

Tenendo conto di questo risultato, possiamo evidentemente riscrivere HR(ω) (modulo dellafunzione di trasferimento)2 come sommatoria basata non più su N campioni, ma su un numeroinferiore di campioni (praticamente metà). In particolare, se per semplicità supponiamo N dispari,possiamo scrivere che

( )∑−

=

ω=ω2/)1N(

0nR nTcos)n(h2)(H

2 Tutti i discorsi seguenti sono relativi al modulo di H(ω)



E’ chiaro che il fattore 2 può essere eliminato ponendo banalmente b(n)=2h(n), dove però bisognastare attenti al fatto che le sequenze b(n) ed h(n) non hanno la stessa lunghezza (mentre h(n) è lungaN, b(n) è lunga (N-1)/2):

( )∑−

=

ω=ω2/)1N(

0n

nTcos)n(b)(H

Di conseguenza, per un filtro FIR a fase lineare, il problema dideterminare i coefficienti h(n) dell’espansione di H(ω) inesponenziali si è ricondotto al problema di determinare icoefficienti b(n) dell’espansione di H(ω) in coseni. Potremmo a questopunto ripetere i discorsi fatti prima a proposito dell’errore E, in modo da scrivere che

∫ ∑∫ ω

ω−ω=ωω=

−

=B

2

d

2/)1N(

0nB

2 d)(H)nTcos()n(bd)(eE

Si può però perfezionare il metodo, con la seguente osservazione: così come abbiamo definito E,noi attribuiamo uguale importanza agli errori commessi alle varie frequenza, ossia quindi cidisinteressiamo della posizione (in frequenza) in cui ciascun errore H(ω)-Hd(ω) si verifica; questo,però, non è molto sensato, in base alle stesse considerazioni fatte in precedenza: un errore in bandapassante o in banda arrestata è per noi molto più grave di un errore in banda di transizione.

Di conseguenza, per attribuire un peso diverso agli errori a seconda di dove vengono commessi, siusa una opportuna funzione peso W(ω): questa consente di definire un errore pesato, definitocome

[ ] )(W)(H)(H)(e d ωω−ω=ω

Questa funzione peso conferisce sostanzialmente un grado di libertà in più al nostro progetto3.Si pone ora nuovamente il problema della minimizzazione. Potremmo procedere come fatto

prima, considerando l’integrale di e2(ω) sulla banda considerata, ma in realtà non conviene. Alcontrario, conviene a questo punto sfruttare un particolare teorema, noto come teoremadell’alternanza di Chebyshev, del quale andiamo perciò a descrivere l’enunciato.

In primo luogo, questo teorema si applica ai casi in cui si vuole approssimare una data funzionemediante una funzione polinomiale trigonometrica, ossia sostanzialmente una somma pesata di Senio Coseni. E’ chiaro che si tratta proprio del nostro caso, in quanto noi vogliamo approssimare lafunzione Hd(ω) mediante una funzione H(ω) che è proprio una somma pesata di coseni:

( )∑−

=

ω=ω2/)1N(

0n

nTcos)n(b)(H

Premesso questo, il teorema afferma che una approssimazione polinomialetrigonometrica è ottima se, nella banda di interesse, la funzione di

errore pesato e(ω) presenta r+1 frequenze estremali, dove 2

1N1r

−=− .

Cominciamo a capire cosa sono queste frequenze estremali. Intanto, dobbiamo considerare chee(ω) è una funzione di ω (e quindi di f), per cui assumerà valori verosimilmente diversi per i vari ω.In particolare, ci saranno pulsazioni alle quali e(ω) assume picchi positivi (cioè massimi locali),pulsazioni in cui e(ω) assume picchi negativi (cioè minimi locali), pulsazioni in cui e(ω)=0 e

3 E’ ovvio che la scelta di una W(f)=1 su tutte le frequenze elimina la pesatura, per cui si torna sostanzialmente ai discorsi di prima



pulsazioni in cui e(ω) assume valori intermedi. Allora, le frequenze estremali di e(ω) sono quellefrequenze in cui e(ω) assume tutti i propri massimi ed i propri minimi locali e, in particolare, talimassimi e minimi corrispondono ad uguali spostamenti (a segno alternato) rispetto al valore medio0.

Tradotto in formule, se indichiamo con δ il massimo spostamento dell’errore rispetto a zero, lefrequenze estremali fk, per k=1,2,....,r+1, sono tali che

( ) 1r1,2,.....,k )1(fe kk +=δ−=

Osserviamo subito una cosa: le frequenze estremali fk sono semplicemente r+1 particolarifrequenze nell’intervallo in cui è stata definita la funzione errore; in generale, queste frequenze nonsono equispaziate tra loro, al contrario delle frequenze usate per il campionamento di Hd(ω). E’ benedunque distinguere le frequenze estremali di e(ω) dalle frequenze a cui è stata campionata Hd(ω).

Tornando adesso all’enunciato del teorema, esso afferma che, se la funzione e(ω) presenta r+1

frequenze estremali, allora il polinomio ( )∑−

=

ω=ω2/)1N(

0n

nTcos)n(b)(H rappresenta la migliore

approssimazione della funzione obbiettivo Hd(ω). Non si tratta, ovviamente, della miglioreapprossimazione in assoluto, ma solo con riferimento alla funzione W(ω) che si è scelta.

Notiamo nuovamente che stiamo scegliendo una ottimizzazione diversa da quella vista prima:prima, infatti, consideravamo come funzione errore (cioè come parametro dello scostamento tra H(ω)e Hd(ω)), una misura integrata dell’errore quadratico; in questo caso, invece, stiamo considerando difatto il valore massimo dello scostamento tra H(ω) e Hd(ω). La differenza fondamentale tra i duecriteri è che con la minimizzazione dell’errore quadratico medio non si ha alcun controllo sugliscostamenti alle varie frequenze, mentre con il nuovo metodo visto otteniamo scostamenti pesaticostanti:

Fatte queste considerazioni, si capisce che non abbiamo ancora trovato il modo di calcolare icoefficienti h(n) che minimizzino lo scarto tra H(ω) e Hd(ω). Abbiamo semplicemente appurato qualisiano le condizioni sotto le quali questo problema ammette una soluzione: nella banda di interesse, lafunzione errore deve presentare r+1 frequenze estremali. Adesso dobbiamo quindi capire comedeterminare h(n): sostanzialmente, dobbiamo combinare la definizione della funzione errore col fattoche ci devono essere r+1 frequenze estremali.

Andiamo allora semplicemente ad applicare la definizione di frequenze estremali: sono lefrequenze alle quali risulta ( ) δ−= k

k )1(fe ; sostituendo l’espressione di e(f), otteniamo

( ) ( ) ( )[ ] 1r1,2,....,k )1(fHfHfW kkkdk +=δ−=−

Essendo r+1 le frequenze estremali, possiamo scrivere r+1 equazioni di questo tipo:



( ) ( ) ( )[ ]( ) ( ) ( )[ ]

( ) ( ) ( )[ ]

δ−=−

δ=−δ−=−

++++

1r1r1rd1r

22d2

11d1

)1(fHfHfW

...

fHfHfW

fHfHfW

Le incognite, in questo sistema sono evidentemente δ ed i coefficienti b(n) in funzione dei qualibisogna esprimere H(f). Sono invece noti i valori Hd(fk), in quanto la funzione obbiettivo si supponenota.

Conviene portare i valori W(fk) e H(fk) nei secondi membri delle varie equazioni:

( ) ( ) ( )

( ) ( ) ( )

( ) ( ) ( )

δ−+=

δ+=

δ−=

+

+++

1r

1r1r1rd

222d

111d

fW)1(fHfH

...

fWfHfH

fWfHfH

Sostituendo adesso al posto dei vari Hd(fk) le rispettive espressioni in funzione dei b(n), otteniamo

( ) ( ) ( )

( ) ( ) ( )

( ) ( ) ( )

=δ

−+π

=δ

+π

=δ

−π

++

+−

=+

−

=

−

=

∑

∑

∑

1rd1r

1r2/)1N(

0n1r

2d2

2/)1N(

0n2

1d1

2/)1N(

0n1

fHfW

)1(nTf2cos)n(b

...

fHfW

nTf2cos)n(b

fHfW

nTf2cos)n(b

Avendo detto che 2

1N1r

−=− , i coefficienti b(n) sono in numero r. Di conseguenza, il sistema è

di r+1 equazioni in r+1 incognite (oltre ai b(n) bisogna infatti considerare δ). Lo scriviamo allora informa matriciale:

( ) ( ) ( ) ( )( ) ( ) ( ) ( )( ) ( ) ( ) ( )

( ) ( ) ( ) ( )

( )( )

( )

=

δ

−

−π−ππ

−π−ππ

π−ππ

−π−ππ

+

+

+

+++

1rd

2d

1d

1r

1r

1r1r1r

3333

2222

1111

fH

...

fH

fH

2

1Nb

....

)1(b

)0(b

fW

)1(Tf)1N(cos....nTf4cosTf2cos1

........................fW

1Tf)1N(cos....nTf4cosTf2cos1

fW

1Tf)1N(cos...Tf4cosTf2cos1

fW

1Tf)1N(cos...nTf4cosTf2cos1



Questo sistema è facilmente risolvibile se le sue dimensioni sono sufficientemente piccole, ossiasostanzialmente se la lunghezza N del filtro è piccola. In caso contrario, si usa una proceduraiterativa basata sull’algoritmo di Remez, che sarà oggetto del prossimo paragrafo.

Algoritmo di Remez

Per comprendere come funzioni l’algoritmo iterativo di Remez, procediamo con un esempioconcreto.

Vogliamo progettare un filtro equiripple (cioè con δ costante) la cui risposta all’impulso h(n)sia costituita da N=13 campioni. In base ai discorsi fatti prima, questo significa che dobbiamodeterminare 7 coefficienti b(n), i quali definiscono la funzione di trasferimento H(ω):

( )∑=

ω=ω6

0n

nTcos)n(b)(H

Avendo detto che 2

1N1r

−=− , è chiaro che le frequenze estremali dovranno essere 8. Il punto

di partenza dell’algoritmo è allora quello di fissare a caso 8valori da attribuire alle frequenze estremali. Dato che le frequenze estremalisono relative a tutta la banda di interesse, dobbiamo distribuirle tra la banda passante e la bandaarrestata. Per esempio, ne prendiamo 3 in banda passante e 5 in banda arrestata:

f (Hz)bandapassante

bandaarrestata

Scelta (casuale) delle frequenze estremali in banda passante ed in banda arrestata. Sono evidenziatein figura sia le frequenze scelte sia i corrispondenti valori che si vuole che il polinomio H(ω) assuma in

loro corrispondenza

Da notare che nessuno ci impedisce di prendere le frequenze estremali equispaziate, sia in bandapassante sia in banda arrestata (anzi questa è la scelta che si fa generalmente), ma è bene ricordareancora una volta che tali frequenze, in generale, non sono equispaziate. Dato che la prima scelta ècasuale, possiamo scegliere come ci pare.

A questo punto, imponiamo che il polinomio ( )∑=

ω=ω6

0n

nTcos)n(b)(H assuma dei valori ben

precisi in corrispondenza di queste frequenze: in particolare, esso dovrà valere 1 nelle frequenzeestremali in banda passante, mentre dovrà valere 0 nelle frequenze estremali in banda arrestata:

( )

∈ω∈ω

=ω=ω ∑= arrestata banda se 0

passante banda se 1nTcos)n(b)(H

k

k6

0nkk



Imponendo queste condizioni, otteniamo un sistema di 8 equazioni in 7 incognite.......(?)...... Lasoluzione del sistema corrisponde ai valori dei coefficienti b(n) tali che il polinomio passi per i puntifissati (quelli cioè indicati nell’ultima figura).

Il polinomio così determinato non soddisfa però necessariamente le condizioni desiderate, inquanto i massimi ed i minimi non necessariamente (a meno di un colpo di fortuna) assumono ugualevalore in modulo. Allora, il passo successivo consiste nel determinare i massimi ed i minimi delpolinomio appena determinato e soprattutto le frequenze fk corrispondenti a tali massimi e minimi.Queste frequenze determinano univocamente la matrice dei coefficienti del sistema trovato prima equi di seguito riproposto:

( ) ( ) ( ) ( )( ) ( ) ( ) ( )( ) ( ) ( ) ( )

( ) ( ) ( ) ( )

( )( )

( )

=

δ

−

−π−ππ

−π−ππ

π−ππ

−π−ππ

+

+

+

+++

1rd

2d

1d

1r

1r

1r1r1r

3333

2222

1111

fH

...

fH

fH

2

1Nb

....

)1(b

)0(b

fW


........................fW


fW


fW


Come si può ben vedere, una volta note le frequenze fk, la matrice dei coefficienti diventa nota, inquanto la funzione finestra W(f) è stata fissata a priori e se ne conoscono perciò tutti i valori.

Il sistema può dunque essere risolto, il che significa determinare i valori dei coefficienti b(n) delfiltro, nonché il valore di δ. In effetti, però, ciò che più interessa in questo momento è il valore di δ.Di conseguenza, invece di risolvere l’intero sistema, si procedere a determinare solo il valore di δ,tramite la seguente formula analitica:

)f(W

a)1(....

)f(W

a

)f(W

a

)f(Ha....)f(Ha)f(Ha

1r

rr

2

1

1

0

1rdr2d11d0

+

+

−++−

+++=δ

dove bisogna porre

∏+

≠= ω−ω

=1r

kn1n nk

k coscos

1a

Questa espressione analitica scaturisce direttamente dall’espressione matriciale del sistema.Al primo passo, dunque, fissate a caso le r+1 frequenze estremali, si può calcolare δ.

Successivamente, anziché calcolare esattamente il polinomio tramite la risoluzione del sistema, sipuò effettuare una interpolazione. Si pone allora il problema della scelta della funzione interpolante,ma la soluzione è semplice: sicuramente, non è possibile usare la funzione sin(x)/x, per il semplicemotivo che i punti fissi non saranno necessariamente equidistanti; è lecito allora usare i polinomi diLagrange (oppure anche le splines) Ciò significa esprimere il polinomio H(ω) nel modo seguente:

∑

∑+

=

+

=

ω−ωβ

ω−ωβ

ω=ω

1r

1k k

k

1r

1k k

kk

coscos

coscos)(H

)(H



Per applicare questa formula servono evidentemente i valori di H(ω) nelle frequenze estremalifissate prima. Ma questi valori si calcolano facilmente tramite la stessa relazione usata per costruireil sistema:

( ) ( ) ( ) 1r1,2,....,k fW

)1(fHfH

k

k

kdk +=δ−

−=

Una volta calcolato δ, quindi, siamo in grado di determinare i coefficienti b(n) del polinomioH(ω) e quindi anche i coefficienti dell’errore pesato:

[ ] )(W)(H)(H)(e d ωω−ω=ω

Siamo perciò in grado di calcolare e(ω) delle frequenze estremali ωk selezionate inizialmente (acaso). Si tratta ora di valutare se queste sono le vere frequenze estremali (o una loro accettabileapprossimazione) oppure no. I casi possibili sono allora due:

a) se risulta |e(ω)|≤δ per tutte le frequenze ωk, significa che abbiamo trovato l’approssimazioneottima, per cui l’algoritmo può già essere arrestato;

b) se invece risulta |e(ω)|>δ per almeno una frequenza ωk, allora dobbiamo scegliere un ulterioreinsieme di frequenze. Il modo più sensato di fare la nuova scelta è quello di individuare lefrequenze alle quali la funzione e(ω) risulta avere i propri massimi e minimi locali. Puòcapitare che queste frequenze siano più di r+1: in questo caso, bisogna selezionare quelle in cui|e(ω)| risulta maggiore.

Se si è verificata la situazione (b), per cui è stato scelto un nuovo set di frequenze estremali,l’algoritmo va ripetuto come prima

L’ultima condizione imposta fa si, quindi, che le nuove frequenze fk siano in corrispondenza deimassimi e dei minimi della funzione errore pesato e(ω) calcolata al primo passo. La successivaiterazione porta dunque a ricalcolare δ (con la formula analitica) e successivamente i coefficienti delpolinomio H(ω) (tramite l’interpolazione di Lagrange). Si può quindi ricalcolare e(ω) e fare le steseverifiche di prima circa la condizione |e(ω)|≤δ. Non appena questa condizione risulta verificata perqualsiasi frequenza del set individuato, l’algoritmo si può arrestare.

Quando l’algoritmo è stato arrestato (il che accade, sostanzialmente, quando i valori di fk e di δnon variano più in modo apprezzabile), l’ultimo insieme calcolato dei coefficienti di H(ω) è quelloche fornisce la migliore approssimazione. Da esso si può risalire alla funzione di risposta all’impulsoh(n) del filtro.

Si è quindi progettato un filtro equiripple a fase lineare, a 13 campioni, che approssima al meglioun andamento di tipo passa-basso secondo l’approssimazione di Chebyshev.

Procedura di Parks e McClellan

Nel paragrafo precedente abbiamo spiegato per via essenzialmente qualitativa come si svolge ilprogetto di un filtro FIR a fase lineare tramite l’applicazione dell’algoritmo di Remez. Vediamoadesso, in modo più operativo e rigoroso dal punto di vista matematico, come bisogna procedere. Laprocedura che illustriamo è quella proposta da Parks e McClellan.

I dati di cui bisogna disporre inizialmente sono 3:



• funzione di trasferimento obbiettivo, Hd(ω);

• funzione di pesatura dell’errore, W(ω);

• lunghezza del filtro da implementare, N.

La funzione di trasferimento obbiettivo viene ovviamente specificata tramite il vettore deicampioni che si vuole siano “onorati” dalla funzione di trasferimento H(ω) ottenuta dal progetto.

La prima considerazione da fare, per comprendere l’applicazione del metodo, riguarda proprio laH(ω) che si desidera ottenere. Il nostro obbiettivo è quello di ottenere un filtro FIR a fase lineare; inbase ai discorsi fatti già in precedenza (soprattutto nel capitolo sulla trasformata zeta), questosignifica che H(ω) è esprimibile come prodotto di un termine HR(ω) reale e di un termine di fasepuramente rettilineo; in particolare, se N è la lunghezza del filtro (supposta per semplicità dispari),

sappiamo che il termine di fase è 2

1Nj

e−

ω− (corrispondente al ritardo temporale da imporre alla

funzione di risposta all’impulso perché essa, da essere reale e pari, diventi reale e causale:

n

h(n)

n

Se la funzione di risposta all’impulso è reale e pari, di lunghezza N dispari, sono necessari (N-1)/2ritardi temporali (ciascuno di un passo di campionamento) per farla diventare causale e quindi

fisicamente realizzabile. Nel caso di figura, N=5, per cui bisogna ritardare di 2 passi di campionamentola sequenza perché diventi causale. Questo comporta che la funzione di trasferimento H(ω) abbia una fase

data da e-jω2T

Quindi, possiamo scrivere in generale che

2

1Nj

R e)(H)(H−

ω−ω=ω

L’utilità di questo passaggio è evidente: nel progetto del filtro FIR tramite l’applicazionedell’algoritmo di Remez, noi ragioniamo solo sul modulo della funzione di trasferimento da otteneree non sulla fase, che è chiaramente fissata una volta stabilita la lunghezza N del filtro. Diconseguenza, nell’evoluzione dell’algoritmo ci si porta dietro solo il termine HR(ω) reale.

Abbiamo inoltre trovato prima che

( )∑−

=

ω=ω2/)1N(

0n

nTcos)n(b)(H

In effetti, questa espressione è valida solo se H(ω) è reale (cioè a fase nulla), il che significa chequella è rigorosamente l’espressione di HR(ω):

( )∑−

=

ω=ω2/)1N(

0nR nTcos)n(b)(H



Sulla base di queste osservazioni, è HR(ω) la funzione che deve rispettare i vincoli imposti dallamaschera rappresentata nella figura seguente:

Caratteristiche del modulo della funzione di trasferimento di un filtro fisicamente realizzabile. Ci sonole oscillazioni in banda passante, dove la funzione varia in un intervallo di ampiezza 2δ1, e le oscillazioni

in banda arrestata, dove la funzione varia in un intervallo di ampiezza δ2>δ1

I vincoli imposti da questa maschera possono essere sinteticamente rappresentati nel modoseguente:

S2R2

P1R1

)(H- arrestata banda

1)(H-1 passante banda

ω>ωδ≤ω≤δ→

ω≤ωδ+≤ω≤δ→

Tornando adesso alla funzione HR(ω), i discorsi fatti erano validi per una lunghezza N del filtrodispari. In effetti, se la lunghezza del filtro è pari, le cose non cambiano molto: si può infattidimostrare che, in questo caso, HR(ω) assume l’espressione

( )∑−

=

ωω

=ω1)2/N(

0nR nTcos)n(a

2cos)(H

dove i coefficienti a(n) sono linearmente legati ai coefficienti b(n) dalle seguenti relazioni:

=

−

−=−−=

=

2

Nb21

2

Na

22

N1,2,....,k )1k(a)k(b2)k(a

)1(b2

1)0(a

Entrambi questi casi (N dispari o pari) valgono nell’ipotesi che h(n) sia simmetrica (rispetto ad uncampione centrale se N è dispari o rispetto ad un istante centrale se N è pari). Ci sono poi altri duecasi da considerare: il caso in cui h(n) è antisimmetrica (cioè a simmetria dispari) con N pari o con N



dispari. E’ possibile allora rappresentare contemporaneamente tutti e 4 i casi possibili tramite leseguenti posizioni:

)(P)(Q)(HR ω⋅ω=ω

( )

−

−

=ωα=ω ∑=

pari N ed ricaantisimmet h(n) 1-2

N

dispari N ed ricaantisimmet h(n) 2

3N

pari N ed simmetrica h(n) 1-2

N

dispari N ed simmetrica h(n) 2

1N

L dove nTcos)n()(P

dove

L

0n

ωω

ω

=ω

pari N ed ricaantisimmet h(n) 2

sin

dispari N ed ricaantisimmet h(n) sin

pari N ed simmetrica h(n) 2

cos

dispari N ed simmetrica h(n) 1

)(Q

In queste espressioni, la sequenza { })k(α è sostanzialmente quella dei parametri del filtro. Come

visto in precedenza, questi coefficienti { })k(α sono linearmente legati ai campioni della rispostaall’impulso h(n).

Queste precisazioni sono opportune in quanto, nell’eventualità di ideare una proceduramatematica (ad esempio con Matlab) da implementare sul calcolatore, è possibile selezionare unoqualsiasi dei 4 casi. Si nota, inoltre, che, per h(n) simmetrica e di lunghezza N dispari, ritroviamoesattamente gli stessi risultati trovati prima.

A questo punto, possiamo supporre di aver fissato la lunghezza N (pari o dispari) del filtro, il tipodi risposta all’impulso (simmetrica o antisimmetrica) e la funzione obbiettivo Hd(ω). In particolare,su quest’ultima possiamo fare ragionamenti analoghi a quelli fatti per H(ω) e HR(ω): abbiamo infattidetto che l’algoritmo necessita in ingresso dei campioni (equispaziati) di Hd(ω) che si vuole venganoonorati dalla H(ω) finale; dobbiamo fornire manualmente questi campioni, il che significa chedaremo dei campioni reali; ad esempio, nel caso di un filtro passa basso ideale, imporremo che talicampioni valgano 1 in banda passante e 0 in banda arrestata. Allora, possiamo esprimerci dicendoche la funzione obbiettivo non è Hd(ω), ma la sua versione reale HdR(ω), ossia fondamentalmente ilsuo modulo. Nel caso del filtro passa-basso ideale, avremo quanto segue:

f (Hz)

HdR(f)

fC/2-fC/2



L’ultima informazione da dare in ingresso all’algoritmo è la funzione W(ω) con la qualeintendiamo pesare gli errori del filtro HR(ω) ottenuto rispetto al filtro desiderato HdR(ω):

[ ] )(W)(H)(H)(e RdR ωω−ω=ω

In particolare, è conveniente normalizzare W(ω) all’unità nella banda arrestata e porre invece

12 /)(W δδ=ω nella banda passante:

ωωδδ

=ωarrestata bandain 1

passante bandain /)(W 12

Se fosse δ1=δ2 (filtro equiripple), la funzione peso non avrebbe ovviamente più senso, in quantovarrebbe 1 per tutte le frequenze: è questo il caso in cui si vuole attribuire ugual peso a tutti gli erroridel filtro implementato rispetto al filtro ideale desiderato. Se invece δ1≠δ2, allora la funzione pesointerviene a dare maggiore importanza all’errore nell’intervallo di frequenze che richiede scarti piùpiccoli: per esempio, se δ1<δ2 la funzione peso attribuisce peso maggiore agli errori in bandapassante che non a quelli in banda arrestata.

Tornando all’espressione di e(ω), possiamo sostituire l’espressione di HR(ω) e fare qualchepassaggio:

[ ] )(Q)(W)(P)(Q

)(H)(W)(P)(Q)(H)(e dR

dR ωω

ω−

ωω

=ωωω−ω=ω

Per questioni di comodità matematica, ci conviene definire una funzione peso modificata W’(ω)ed una funzione di trasferimento obbiettivo modificata H’dR(ω):

[ ] )('W)(P)('H)(e

)(Q)(W)('W

)(Q

)(H)('H

dR

dRdR ωω−ω=ω→

ωω=ωωω

=ω

E’ evidente che, fissate a priori HdR(ω) e W(ω) nonché la lunghezza N del filtro ed il tipo dirisposta all’impulso, anche le due nuove funzioni sono automaticamente fissate.

Sostituendo infine l’espressione di P(ω), che contiene le incognite del nostro problema, otteniamo

( ) )('WnTcos)n()('H)(eL

0ndR ω

ωα−ω=ω ∑

=

A questo punto, il problema di fondo è determinare i parametri { })k(α del filtro che minimizzano

il massimo valore assoluto di e(ω) su tutta la banda di frequenza che stiamo considerando: in terminimatematici, possiamo dire che noi cerchiamo la soluzione del problema rappresentato da

{ }[ ])(emaxmin

B)k(ω

∈ωα

dove ovviamente ω varia solo nella banda B in cui intendiamo approssimare il filtro idealedesiderato.



Sostituendo l’espressione dell’errore trovata poco fa, otteniamo

{ }( )

ω

ωα−ω ∑

=ωα

)('WnTcos)n()('HmaxminL

0ndR

)k(

Come già osservato in precedenza, la soluzione a questo problema è indicata dal teoremadell’alternanza di Chebyshev, in base al quale la condizione necessaria e sufficiente affinché ilpolinomio P(ω) sia la unica e migliore approssimazione di H’dR(ω) nella banda B è che la funzionee(ω) abbia almeno L+2 frequenze estremali nella banda B: se 1L10 .... +ω<<ω<ω sono tali frequenze

estremali, deve risultare che( ) ( )( ) ( )ω=ω

ω−=ω

∈ω

+

emaxe

ee

Bi

1ii

A questo punto dobbiamo ripetere gli stessi ragionamenti fatti in precedenza circa il modo con cuiimporre le condizioni del teorema dell’alternanza. Riportiamo perciò direttamente la formamatriciale del sistema da risolvere, adottando ovviamente la nuova simbologia introdotta in questoparagrafo:

( ) ( ) ( ) ( )( ) ( ) ( ) ( )( ) ( ) ( ) ( )

( ) ( ) ( ) ( )

( )

( )( )( )

( )

=

δα

αα

−π−ππ

π−ππ

−π−ππ

π−ππ

+

+

+

+++

1LdR

2dR

1dR

0dR

1L

1L

1L1L1L

2222

1111

0000

f'H

...

f'H

f'H

f'H

L

....

)1(

)0(

fW


........................fW


fW


fW


Le incognite di questo sistema sono gli L+1 coefficienti α(k) del polinomio P(ω) e il δ.Per quanto detto in precedenza, al primo passo dell’interazione non conviene

risolvere il sistema, ma solo calcolare δ, tramite la seguente formula analitica:

)f('W

a)1(....

)f('W

a

)f('W

a

)f('Ha....)f('Ha)f('Ha

1L

1L1L

1

1

0

0

1LdR1L1dR10dR0

+

++

++

−++−

+++=δ

dove bisogna porre

∏+

≠= ω−ω

=1L

kn0n nk

k coscos

1a

Al primo passo, dunque, fissate a caso le L+2 frequenze estremali, si può calcolare δ.Possiamo successivamente calcolare il polinomio P(ω) mediante una interpolazione di Lagrange:

ciò significa che intendiamo costruire il polinomio



∑

∑

=

=

ω−ωβ

ω−ωβ

ω=ω

L

0k k

k

L

0k k

kk

coscos

coscos)(P

)(P

dove i campioni di P(ω) si ottengono dalla stessa relazione utilizzata per costruire il sistema, vale adire

( ) ( ) ( ) 1L0,2,....,k f'W

)1(f'HfP

k

k

kdRk +=δ−

−=

A questo punto, disponiamo sia di δ sia dei coefficienti α(k) del polinomio P(ω). Possiamo perciòcalcolare la funzione errore nelle frequenze estremali prescelte o, comunque, in un set abbastanzaalto di frequenza4:

( ) 1L0,2,....,k )('WnTcos)n()('H)(e k

L

0nkkdRk +=ω

ωα−ω=ω ∑

=

I casi possibili sono adesso due:

a) se risulta |e(ω)|≤δ per tutte le frequenze ωk, significa che abbiamo già trovatol’approssimazione ottima, per cui l’algoritmo può già essere arrestato;

b) se invece risulta |e(ω)|>δ per almeno una frequenza ωk, allora dobbiamo scegliere un ulterioreinsieme di frequenze. Scegliamo le frequenze alle quali la funzione e(ω) risulta avere i proprimassimi e minimi locali. Può capitare che queste frequenze siano più di L+2: in questo caso,bisogna selezionare quelle in cui |e(ω)| risulta maggiore.

Se si è verificata la situazione (b), per cui è stato scelto un nuovo set di frequenze estremali,l’algoritmo va ripetuto come prima. Dobbiamo perciò ricalcolare δ (con la formula analitica) esuccessivamente i coefficienti α(k) del polinomio P(ω) (tramite l’interpolazione di Lagrange). Si puòquindi ricalcolare e(ω) e fare le stesse verifiche di prima circa la condizione |e(ω)|≤δ. Non appenaquesta condizione risulta verificata per qualsiasi frequenza del set individuato, l’algoritmo si puòarrestare.

In pratica, l’algoritmo fa in modo da incrementare il valore di δ fin quando esso raggiunge il suolimite superiore, che rappresenta la soluzione migliore.

Quando l’algoritmo è stato arrestato (il che accade, sostanzialmente, quando i valori di fk e di δnon variano più in modo apprezzabile), l’ultimo insieme calcolato dei coefficienti di P(ω) è quelloche fornisce la migliore approssimazione.

Noto P(ω), risulta automaticamente determinata la funzione di trasferimento HR(ω): sappiamoinfatti che

)(P)(Q)(HR ω⋅ω=ω

e che Q(ω) è un dato del problema.

4 A questo proposito, ricordiamo che il numero di frequenze in cui calcolare e(ω), al fine di verificare i valori del suo modulo, deve

essere abbastanza elevato; in particolare, si è verificato che è sufficiente un numero 16N di punti, dove N è sempre la lunghezzadel filtro



Da HR(ω) si risale alla funzione di trasferimento del filtro: ad esempio, abbiamo visto che, per Ndispari e per una h(n) a simmetria pari, si ha che

2

1Nj

R e)(H)(H−

ω−ω=ω

Infine, con le opportune formule di antitrasformazione si ottiene h(n) a partire da H(ω).Il programma al computer che implementa tutto questo meccanismo è stato proposto da Parks e

McClellan. In particolare, l’algoritmo è fatto in modo che l’uscita finale non siano i coefficienti diP(ω), ma direttamente i campioni della risposta all’impulso h(n) del filtro.

Facendo un riepilogo generale della procedura appena illustrata, possiamo inquadrarla in 4 passifondamentali:

1) introduzione dei dati di ingresso: bisogna fissare la funzione di trasferimento Hd(ω)desiderata, la funzione W(ω) di pesatura dell’errore e la lunghezza N del filtro;

2) formulazione del problema equivalente secondo l’approssimazione diChebyshev: significa individuare, a partire dai dati in ingresso, le funzioni H’dR(ω), W’(ω) eP(ω);

3) applicazione dell’algoritmo di Remez: significa eseguire la procedura iterativa checonduce a determinare i coefficienti α(k) del polinomio P(ω) che offre la miglioreapprossimazione possibile di H’dR(ω); si parte da una scelta iniziale casuale delle frequenzeestremali;

4) calcolo della risposta all’impulso: bisogna risalire ai campioni della sequenza h(n) apartire dai coefficienti del polinomio P(ω) ottenuti al passo precedente.

Nella figura seguente riportiamo un diagramma di flusso di come evolve l’algoritmo di Remez:



Parametri in ingresso

Scelta iniziale casuale delleL+2 frequenze estremali

Calcolo del valore di delta sullabase delle freq. estremali scelte

Determinazione di P(w) su L+2punti (freq. estremali)

Valutazione dell'errore e(w) e delsuo modulo

δ≤ω)(eω∀

STOP(migliore approssimazione

ottenuta)

SI

NO

Determinazione dei max e minlocali della funzione e(w) e delle

corrispondenti frequenze

Trovatepiù di L+2frequenze

Selezione delle L+2 frequenze incui |e(w)| assume i valori maggiori

NO

SI

Diagramma di flusso dell’algoritmo di Remez

Per concludere, facciamo una serie di considerazioni sulla procedura appena illustrata e, inparticolare, sull’importanza di W(ω):

• in primo luogo, l’uso di una funzione W(ω) di pesatura dell’errore, scelta con δ1=δ2 (filtroequiripple), garantisce che l’escursione della funzione H(ω) ottenuta, rispetto a quelladesiderata Hd(ω), è costante, dopo la pesatura, su tutto l’intervallo di frequenze considerato;

• in secondo luogo, nella procedura illustrata è comunque rimasto insoluto un punto: è vero chel’uso della funzione W(ω) consente di ottenere errore pesato equiripple, ma è anche vero che ilvalore del ripple δ non è sotto il nostro controllo, in quanto è quello fornito dall’algoritmoall’atto della sua convergenza; il parametro da cui dipende δ è la lunghezza N del filtro; diconseguenza, il valore di N non è fisso a priori: noi lo fissiamo arbitrariamente all’avviodell’algoritmo, dopo di che verifichiamo il valore che l’algoritmo ci da per δ; se tale valore farientrare il filtro nella maschera, abbiamo finito; in caso contrario, invece, dobbiamonecessariamente incrementare N (lasciando invariato tutto il resto) e ripetere la procedura.



Sulla base di queste considerazioni, si comprende l’importanza di scegliere una funzione W(ω)opportuna, pena l’uso di un numero N di campioni maggiore rispetto a quello effettivamentenecessario.

Per quanto riguarda, invece, la scelta di N pari o dispari, le considerazioni sono altre e peraltromolto semplice: se N dispari, abbiamo visto che il ritardo temporale (necessario per rendere h(n)causale) è di un numero intero di passi di campionamento, mentre invece, se N è pari, il ritardo è diun numero intero di mezzi passi di campionamento. Nel primo caso avremo una fase, nella funzionedi trasferimento H(ω), del tipo e-jωkT, con k intero, mentre nell’altro caso avremo una fase e-jωkT/2,sempre con k intero. La scelta dipende perciò dal contesto.

PROCEDURA REMEZ DI MATLABIl Matlab contiene una procedura standard, denominata Remez.m, che implementa la procedura

di Parks-McLellan prima descritta.Nella versione più semplice, la procedura richiede in ingresso 3 parametri:

• la lunghezza N (iniziale) della risposta all’impulso del filtro;

• il vettore F contenente le frequenze in cui si vuole che la funzione di trasferimento ottenutaonori i campioni della funzione di trasferimento desiderata; si tratta, in particolare, di digitarefrequenze normalizzate alla frequenza di Nyquist fC/2: si richiede infatti che i valori inclusi inquesto vettore siano compresi nell’intervallo (0,1), dove 1 corrisponde proprio alla frequenza diNyquist; si richiede inoltre che il numero di valori digitati sia pari, in quanto (???)

• il vettore M contenente proprio i campioni della funzione di trasferimento desiderata: sirichiede, in particolare, che questi valori siano uguali nella banda passante (ad esempio tutti 1)e nella banda arrestata (ad esempio tutti 0).

Una tipica chiamata di questa procedura è dunque la seguente:

F=[0.05,0.1,0.15,0.2,0.25,0.3,0.35,0.4,0.69,0.7,0.75,0.8,0.85,0.9,0.95,0.99]M=[1,1,1,1,1,1,1,1,0,0,0,0,0,0,0,0]h=remez (13,F,M)

Il vettore F contiene 16 frequenze, non equispaziate, in cui vengono imposti precisi valori(contenuti nel successivo vettore M) che la funzione di trasferimento ottenuta dovrà rispettare:vengono dichiarati 8 valori unitari per la banda passante e 8 valori nulli per la banda arrestata. Lasuccessiva istruzione fa eseguire l’algoritmo, richiedendo una lunghezza del filtro pari a 13.

Mandando in esecuzione, il programma restituisce i seguenti risultati:

h = [-0.0150,-0.0012,0.0416,-0.0204,-0.0992,0.1086,0.4773,0.4773,0.1086,-0.0992,

-0.0204,0.0416,-0.0012,-0.0150]

I campioni ottenuti sono 14 (a conferma che non è detto che la lunghezza scelta inizialmente siaquella sufficiente), sono tutti reali e sono inoltre simmetrici rispetto ad un istante centrale collocatotra il settimo e l’ottavo campione (a conferma che il filtro è a fase lineare).

Una chiamata leggermente più complessa della procedura è quella in cui viene indicata anche lafunzione di pesatura dell’errore.



CONCLUSIONI SUI METODI DI PROGETTO DI FILTRI FIRAbbiamo esaminato due diversi metodi di progetto di filtri FIR: il metodo delle finestre ed il

metodo del campionamento in frequenza (basato sull’algoritmo di Remez). Disegnando in scalalogaritmica i filtri con lo stesso numero di coefficienti, ma progettati con i vari metodi visti, si notacome il metodo basato su Remez è quello che riesce a garantire la migliore attenuazione.

Autore: SANDRO PETRIZZELLIe-mail: [email protected]

sito personale: http://users.iol.it/sandrysuccursale: http://digilander.iol.it/sandry1

Capitolo 4 Progetto di filtri FIR (II) -...

Documents

Transcript of Capitolo 4 Progetto di filtri FIR (II) -...