Guida alla progettazione di filtri digitali in MATLAB per l'esame di Teoria ed elaborazione dei segnali

MATLAB

ons, in particolare: e

Signal Processing Toolbox Control System Toolbox.

Specifiche di progetto e premesse

Per progettare un filtro stabile (BIBO) e causale si devono avere i poli all’interno del cerchio

unitario del piano complesso.

Le tipologie di filtri digitali si dividono in:

• FIR: sono filtri sempre stabili (poli in z=0), possono essere causali e a fase lineare

(hanno riposta in frequenza simmetrica), la classe di simmetria della risposta impulsiva

dipende dalla tipologia di filtro scelta (vedi p.44 "Teoria ed elaborazione dei segnali");

• IIR: possono presentare fase lineare solo se il sistema è non causale; la stabilità BIBO

richiede che la ROC comprenda la circonferenza unitaria, cioè che tutti i poli siano

all’interno di essa per un sistema causale (attenzione in fase di progetto).

In generale i filtri FIR sono più semplici da implementare, i filtri IIR hanno a loro vantaggio

un ordine tipicamente più basso dei FIR per pari selettività. I FIR hanno sempre fase lineare,

gli IIR no.

Esempi di filtri digitali in MATLAB

BPF con traslazione del LPF (FIR)

Si progetti un filtro passa banda, il segnale è campionato a 10.000.000 Hz, il filtro deve essere

centrato a 1Mhz. La larghezza di banda deve essere 100kHz, la larghezza delle bande laterali

10kHz. Attenuazione fuori banda almeno 40dB, guadagno in banda 1.

Risoluzione

clc

clear all

close all

Dopo l’esecuzione di molti script, la Command Window e la memoria di MATLAB possono

contenere variabili residue o figure aperte. I tre comandi: - ‘clc‘ cancella il contenuto della

4

Simone Randino Teoria ed elaborazione dei segnali Progettazione

Command Window, - ‘clear all‘ rimuove tutte le variabili dal workspace, - ‘close all‘ chiude

tutte le finestre grafiche aperte. In questo modo l’esecuzione inizia sempre in un ambiente

“pulito”.

T = 1*10^(-7); % periodo di campionamento

fs = 1/T; % frequenza di campionamento (10 MHz)

N = 11101; % numero di campioni

df = 1/(N*T); % risoluzione in frequenza

Viene definito un sistema con frequenza di campionamento MHz (T ns). Il

f = 10 = 100

s

numero di campioni stabilisce la lunghezza delle FFT e, con essa, la risoluzione

N = 11101

.

in frequenza 1

∆f = N T

B = 5*(10^4); % banda passante monolatera (50 kHz)

Bt = 1.1*B; % banda con margine di transizione

Si desidera un filtro passa-basso con banda utile kHz e fascia di transizione estesa

B = 50

fino a (55 kHz). La transizione viene poi utilizzata per generare un profilo lineare

B = 1,1 B

t

del guadagno.

H = zeros(1, N);

H conterrà i valori discreti della risposta in frequenza del filtro (definiti da 0 a con

−∆f

f s

passo ). Tutti gli elementi sono inizialmente posti a zero.

∆f

for k = 1:N

f(k) = (k-1)*df;

if f(k) <= B

H(k) = 1; % guadagno banda passante

elseif f(k) <= Bt

H(k) = 1 - (f(k)-B)/(Bt-B); % rampa di discesa

elseif (f(k) > (fs - Bt)) && (f(k) <= (fs - B))

H(k) = (f(k)-(fs - Bt))/(Bt-B);% rampa di salita (simmetria)

elseif f(k) > (fs - B)

H(k) = 1; % replica banda passante

end

end

Il ciclo scorre tutte le frequenze discrete f(k) ed assegna il valore corrente di secondo

H(k)

queste regole:

• Banda passante guadagno unitario.

[0, B]: 5

Simone Randino Teoria ed elaborazione dei segnali Progettazione

• Fascia di transizione discesa lineare da 1 a 0.

[B, B ]:

t

• Simmetria attorno a : viene creata la “copia” a elevata frequenza per assicurare che

f s

la risposta sia reale (proprietà di un FIR reale).

• Banda passante alta nuovo guadagno unitario.

−

[f B, f ]:

s s

h = ifft(H);

L’ifft restituisce la risposta impulsiva periodica corrispondente alla definita. Poiché

h[n] H(f )

è reale e coniugata simmetrica, risulterà reale.

H(f ) h[n]

M = (N-1)/2; % indice centrale

h1 = zeros(1, N);

h1(1:M) = h(M+2:N);

h1(M+1:N) = h(1:M+1);

L’ifft produce una sequenza circolare con la parte causale in coda; le tre righe sopra spostano

la sequenza così che l’indice (campione causale) corrisponda alla posizione M+1. Si

n = 0

ottiene dunque la forma classica “simmetrica” di un FIR lineare di fase.

fc = 1*(10^6); % frequenza di traslazione (1 MHz)

n = 0:(N-1); % vettore degli indici tempo

h_bp = 2 * h1 .* cos(2*pi*fc*(n - M)*T);

Moltiplicare per equivale a traslare la risposta in frequenza cen-

−

h [n] 2 cos 2πf (n M )T

1 c

trandola su MHz (passaggio da LPF a BPF ). Il fattore 2 compensa

f = 1 band-pass filter

c

l’energia suddivisa tra le due bande laterali.

fvtool(h_bp, 1, ’Fs’, fs);

fvtool è lo strumento interattivo di MATLAB per ispezionare ampiezza, fase, risposta im-

pulsiva e stabilità del filtro. 6

Simone Randino Teoria ed elaborazione dei segnali Progettazione

Visualizzazione in fvtool

LPF con tecnica delle finestre (FIR)

Progettare un filtro passa-basso, il sistema ha una frequenza di campionamento pari a 1Mhz,

il filtro deve avere una frequenza di taglio 100kHz, una banda di transizione 10kHz, attenu-

azione in banda arrestata almeno 40dB.

Risoluzione

clc

clear all

close all

Come già spiegato, azzera Command Window, variabili e figure.

fc = 10^6; % frequenza centrale (1 MHz)

B = 100e3; % banda (100 kHz)

Bt = 10e3; % zona transizione (10 kHz)

N = 4001; % semi-lunghezza h[n]

T = 1/fc; % periodo di campionamento

Gli stessi simboli del primo script: qui fissiamo un ordine per il FIR simmet-

2N +1 = 8003

rico. 7

Simone Randino Teoria ed elaborazione dei segnali Progettazione

w = zeros(1, 2*N+1); % prealloc

w(1:N+1) = 0:1/N:1; % rampa ascendente

w(N+2:2*N+1) = 1-1/N:-1/N:0; % rampa disc.

Costruiamo manualmente la finestra di Bartlett (triangolare), lunga campioni, da 0

2N +1

a 1 e ritorno a 0.

h(1:2*N+1) = 2*B*T * sinc( 2*B*T * (-N:N) );

Versione discreta di sinc(2BT troncata a ∈

h (n) = 2BT n), n [−N, N ].

ideale

q = h .* w; % finestra triangolare -> riduce lobi laterali

Il prodotto puntuale attenua i lobi laterali (minor ripple in banda-arrestata) a costo di

transizione più larga.

fvtool(q);

fvtool mostra modulo, fase e risposta impulsiva del FIR rifinito q.

Considerazioni sulle finestre

Finestra rettangolare

% MANUALE (M punti)

M = 801; % lunghezza finestra

w_rect = ones(1, M); % tutti 1

Finestra semplice, ma brusca, ripple elevato in banda di transizione. Adatta quando serve

risoluzione in frequenza estrema e si può tollerare leakage.

Hamming

% MANUALE (M punti)

n = 0:M-1;

w_hamm = 0.54 - 0.46*cos(2*pi*n/(M-1));

% MATLAB built-in

w_hamm = hamming(M).’;

Compromesso classico: lobo principale simile a quello della rettangolare, ma lobi laterali

attenuati.

Hann

w_hann = hann(M).’; % built-in 8

Simone Randino Teoria ed elaborazione dei segnali Progettazione

Formula: .

2πn

−

w[n] = 0.5 0.5 cos −1

M

Blackman

w_black = blackman(M).’; % built-in

Manuale: Usata quando lo in banda arrestata il ripple

− −

0.42 0.5 cos+0.08 cos(4πn/(M 1)).

deve essere molto basso.

Bartlett (triangolare)

w_bart = bartlett(M).’; % built-in

È la finestra già implementata manualmente nello script, impiego tipico nei filtri FIR che

richiedono compromesso “soft” e calcolo leggero.

ampiezza del lobo principale attenuazione dei lobi laterali, finestre

Trade-off principale: vs.

più complesse hanno maggior costo computazione, ma consentono di ottenere risultati migliori

e minor oscillazioni.

LPF FIR con campionamento in frequenza (N pari)

clear all; close all; clc % vedi sez. iniziali

T = 1e-6; % periodo di campionamento

N = 10000; % numero di campioni (PARI)

df = 1/(N*T); % risoluzione di freq.

B = 1e5; % banda passante 100 kHz

Bt = 1.1e5; % banda transizione 110 kHz

Parametri identici a quelli già visti; l’unica novità è che qui è pari, aspetto che richiederà

N

un piccolo “riallineamento” dell’impulso più avanti.

H = zeros(1,N); % risposta in frequenza

for k = 1:N

f(k) = (k-1)*df; % asse delle frequenze

t(k) = (k-1)*T; % asse temporale (solo debug)

if f(k) <= B

H(k) = exp(-j*pi*f(k)*T); % guadagno 1, fase lineare

elseif f(k) <= Bt

H(k) = (1-(f(k)-B)/(Bt-B)) * exp(-j*pi*f(k)*T);

elseif f(k) >= 1/T-Bt && f(k) <= 1/T-B

H(k) = ((f(k)-1/T+Bt)/(Bt-B))* exp(-j*pi*(f(k)-1/T)*T);

elseif f(k) > 1/T-B

H(k) = exp(-j*pi*(f(k)-1/T)*T);

9