Post on 15-Feb-2019
1
Università degli Studi di Padova - Facoltà di Ingegneria
Corso di Laurea in Ingegneria Biomedica
A.A. 2007-2008
Laboratorio di
Elaborazione di Dati, Segnali e
Immagini Biomediche(Parte 4)
Prof. Giovanni SparacinoDipartimento di Ingegneria dell’Informazione
Università di Padovae-mail: giovanni.sparacino@unipd.it
web: http: www.dei.unipd.it/~gianni
2
Scaricare il file dafiltrare.zip e sompattarlo nella cartella di lavoro
SVOLGERE IN AULA
Esercizi 4.1, 4.2, 4.3
Progetto tramite allocazione zeri e poli e applicazione a dati ECG
Esercizio 4.4
Progetto “assistito dal calcolatore”
FINIRE A CASA PER PREPARAZIONE ALL’ESAME
Esercizi 4.5, 4.6 e 4.7
NB: questo set di diapositive contiene anche materiale di riferimento per l’implementazione in
Matlab delle tecniche di progettazione di filtri viste a lezione (forse utile per le tesine)
ESERCITAZIONE SUL FILTRAGGIO
3
0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 2.2
-0.5
0
0.5
1
1.5
2
2.5
Time in seconds
ECG
FILTRAGGIO NOTCH
4
4.1. Con il comando load ecg60hz, si carica in memora un tracciato ECG
corrotto da rumore di rete a 60 Hz, contenuto nel vettore ecg, riferito al
vettore t. Visualizzare il tracciato in figure(1) usando il comando plot.
Scrivere un M-file che, dopo aver alla prima riga di codice caricato il
tracciato, lo filtri dal rumore mediante un opportuno filtro notch. Nello
specifico, all’interno dell’M-file si deve far
a) costruire (con le tecniche di “allocazione zeri e poli" viste a lezione) un filtro
notch, con due zeri complessi e coniugati e due poli nell’origine, che miri a
ridurre il rumore a 60Hz
b) visualizzare in figure(2) la risposta in frequenza del filtro, in modulo e fase, e
in figure(3) il diagramma zeri-poli
c) applicare il filtro al segnale
d) plottare in figure (4) sullo stesso grafico (o su due riquadri incolonnati) il
segnale originale e quello filtrato
5
4.2.
a) Modificare il filtro dell’esercizio precedente collocando i due poli
(complessi coniugati) all’interno del cerchio con la stessa fase dei due
zeri e ad una distanza dall'origine indicata dall’utente
b) Far eseguire progetto del filtro e plot dei risultati per vari valori (es. 0, 0.2,
0.4, 0.8, 0.9, 0.95, 0.99) della distanza suddetta
6
Soluzione 4.1
Due poli nell’origine. Elimino l’interferenza a 60Hz,
ma anche componenti del segnale alle frequenze
vicine (perdo informazione utile ?)
7
Soluzione 4.2
Due poli vicino agli zeri. Elimino l’interferenza a
60Hz, senza intaccare l’ampiezza delle componenti
del segnale alle frequenze vicine (NB: la fase non è
più lineare)
8
Codice (risolve
simultaneamente 4.1 e 4.2
grazie all’input alla riga 14)
NON GUARDARE LA
SOLUZIONE
9
1 2 3 4 5 6 7 8
-2
-1.5
-1
-0.5
0
0.5
1
1.5
2
2.5
Time in seconds
ECG
FILTRAGGIO RUMORE AD ALTA FREQUENZA
10
4.3. Il file ecghfn.mat contiene nel vettore ecg un tracciato (campionato con
frequenza desumibile dal vettore t) corrotto da rumore ad alta frequenza.
Allocando un solo polo sull’asse reale in modo opportuno, ed uno zero
nell’origine, costruire un sistema del primo ordine che corrisponda ad un
filtro passa-basso a guadagno statico unitario. In particolare:
a) Visualizzare la risposta in frequenza in modulo e fase e il diagramma zero-
poli
b) Applicare il filtro al segnale
c) Plottare sullo stesso grafico il segnale originale e quello filtrato
d) Ripetere per vari valori della posizione del polo, osservando come varia la
risposta in frequenza e l’azione del filtro sul segnale
11
NON GUARDARE LA
SOLUZIONE
12
Soluzione 4.3
Si fa fatica a sagomare la risposta in frequenza
come nel passa-basso ideale.
All’esercizio 4.4 si usa, assistiti dal calcolatore,
l’approccio con filtro di Butterworth
13
4.4
Con riferimento al segnale dell’esercizio 4.3. Costruire (aiutandosi con la
funzione butter del matlab, vd help) dei filtri di Butterworth passa-basso con le
caratteristiche sotto, visualizzarne la risposta in frequenza e il diagramma zero
poli, e applicarli (uno alla volta !), al segnale e confrontare su uno stesso plot, al
variare del filtro, segnale originale e segnale filtrato:
a) ordine 2, cutoff 40 Hz
b) ordine 4, cutoff 40 Hz
c) ordine 2, cutoff 80 Hz
d) ordine 4, cutoff 80 Hz
14
NON GUARDARE LA
SOLUZIONE
15
16
4.5 . Si considerino le seguenti specifiche:
� Frequenza di campionamento 100Hz
� Banda passante 20Hz
� Banda oscura 30Hz
� Ripple doppio in banda passante rispetto banda oscura (pesi=[2 1]).
Costruire in Matlab un filtro FIR con la funzione remez, sintassi b=remez(r,fm,m,pesi), e visualizzarne la risposta in frequenza in
modulo e fase. Studiare la risposta in frequenza del filtro al variare dell’ordine r (es. r =3, 5, 15, 20) e del vettore pesi (es. pesi=[2 1],
[100 1], [1 100])
4.6. Si considerino le seguenti specifiche:
� Frequenza di campionamento 100Hz
� Banda passante 20Hz
� Banda oscura 30Hz
� Rp=0.95 Rs=0.05
Costruire un filtro ellittico (istruzione ellipord e poi ellip) e studiare la risposta in frequenza del filtro al variare delle specifiche Rp
(es.0.99, 0.90) e Rs (es. 0.01, 0.1).
Suggerimento per Soluzione 4.5.
fm=[0 20/50 30/50 1]
m=[1 1 0 0]
r=10
pesi=[ 1 2]
b=remez(r,fm,m,pesi)
Suggerimento per Soluzione 4.6
Rpdb=-20*log10(.95)
Rsdb=-20*log10(.05)
[n,wn]=ellipord(20/50,30/50,Rpdb,Rsdb)
[b,a]=ellip(n,Rpdb,Rsdb,wn)
ESERCIZI DA FARE PER CASA PER
PREPARAZIONE ALL’ESAME
17
4.7 (FILTRAGGIO RUMORE A BASSA FREQUENZA)
Il file ecglfn.mat contiene un tracciato ECG campionato a 1000Hz, corrotto da un
drift a bassa frequenza.
a) Dopo aver studiato bene l’help della funzione ellip, progettare un filtro
ellittico passa-alto di ordine 4 (frequenza di taglio 0.5 Hz, 90% di guadagno
in banda passante, 10% di guadagno in banda oscura) in grado di eliminare
il drift.
b) Visualizzare la risposta in frequenza, in modulo e fase, del filtro e il relativo
diagramma zeri poli
c) Applicare il filtro al segnale
d) Plottare sullo stesso grafico il segnale originale e quello filtrato
18
19
20
MATERIALE DI RIFERIMENTO PER
L’IMPLEMENTAZIONE IN MATLAB DELLE TECNICHE
DI PROGETTAZIONE DI FILTRI VISTE A LEZIONE
21
PROGETTO DI FILTRI: SPECIFICHE
Banda OscuraBanda Passante Banda Transizione
Ripple in banda passante Rp
Ripple in banda oscura
Modulo risposta in frequenza (dB)
Wp Ws f
|H(f)|
22
Nell’analogico: esistono tecniche analitiche per determinare H(s)=B(s)/A(s), a
partire da specifiche su Wp, Ws, e sui ripple tollerati in banda passante ed in
banda oscura. Ad esempio:
Butterworth
•monotono sia in banda passante che in banda oscura (no ripple)
Chebyshev
•ripple (uniforme) in banda passante e monotono in banda oscura (tipo 1)
•monotono in banda passante e ripple (uniforme) in banda oscura (tipo 2)
Ellittici
•ripple sia in banda passante che in banda oscura
PROGETTO DI FILTRI IIR (SENZA VINCOLI DI FASE)
23
Approccio nel numerico: si sfruttano i risultati dell’analogico,
determinando da H(s)=B(s)/A(s) una funzione di trasferimento
H(z)=B(z)/A(z) che la “equivale” attraverso vari approcci (invarianza
all’impulso, soluzione numerica dell’equazione differenziale,
trasformazione bilineare, …)
24
Esempio
1000 200 400
Frequenza di campionamento 2000 Hz
Wp=200/1000=0.2
Ws=400/1000=0.4
Rp=3Rs=15
Ripple in banda passante Rp
Attenuazione
in banda
oscura Rs
25
[b,a] = butter(n,Wn)
0 50 100 150 200 250 300 350 400 450 5000
0.5
1
1.5modulo
frequenza effettiva
0 50 100 150 200 250 300 350 400 450 500-4
-2
0
2
4fase
frequenza effettiva
Esempio
Fc=1000;
Ftaglio=200;
Wn=Ftaglio/(Fc/2);
[b,a]=butter(8,Wn);
…
[b,a]=butter(4,Wn);
Filtro di Butterworth passa-basso
26
[b,a] = butter(n,Wn,‘high')
0 50 100 150 200 250 300 350 400 450 5000
0.5
1
1.5modulo
frequenza effettiva
0 50 100 150 200 250 300 350 400 450 500-4
-2
0
2
4fase
frequenza effettiva
Filtro di Butterworth passa-alto
Esempio
Fc=1000;
Ftaglio=100;
Wn=Ftaglio/(Fc/2);
[b,a]=butter(8,Wn,'high');
…
[b,a]=butter(4,Wn,'high');
27
[b,a] = butter(n,[Wn1 Wn2])
Filtro di Butterworth passa-banda
Fc=1000;
Ftaglio1=100;
Wn1=Ftaglio1/(Fc/2);
Ftaglio2=400;
Wn2=Ftaglio2/(Fc/2);
Wn=[Wn1 Wn2]
[b,a]=butter(8,Wn);
…
[b,a]=butter(4,Wn);
0 50 100 150 200 250 300 350 400 450 5000
0.5
1
1.5modulo
frequenza effettiva
0 50 100 150 200 250 300 350 400 450 500-4
-2
0
2
4fase
frequenza effettiva
28
29
Filtro Ellittico passa-basso
0 50 100 150 200 250 300 350 400 450 5000
0.2
0.4
0.6
0.8
1modulo
frequenza effettiva
0 50 100 150 200 250 300 350 400 450 500-4
-2
0
2
4fase
frequenza effettiva
Fc=1000;
Ftaglio=100;
Wn=Ftaglio/(Fc/2);
Rp=0.5;
Rs=20;
[b,a]=ellip(4,Rp,Rs, Wn);
[b,a] = ellip(n,Rp,Rs,Wn)
30
Filtro Ellittico passa-basso
Fc=1000;
Ftaglio=100;
Wn=Ftaglio/(Fc/2);
Rp=0.5;
Rs=20;
[b,a]=ellip(6,Rp,Rs, Wn);
[b,a] = ellip(n,Rp,Rs,Wn)
0 50 100 150 200 250 300 350 400 450 5000
0.2
0.4
0.6
0.8
1modulo
frequenza effettiva
0 50 100 150 200 250 300 350 400 450 500-4
-2
0
2
4fase
frequenza effettiva
31
ELLIPORD Elliptic filter order selection.
[N, Wn] = ELLIPORD(Wp, Ws, Rp, Rs) returns the order N of the lowest
order digital elliptic filter that loses no more than Rp dB in the
passband and has at least Rs dB of attenuation in the stopband.
Wp and Ws are the passband and stopband edge frequencies, normalized
from 0 to 1 (where 1 corresponds to pi radians/sample).
ELLIPORD also returnsWn, the elliptic natural frequency to
use with ELLIP to achieve the specifications
Come trovo l’ordine “ottimo” ?
>> ellipord(100/500,120/500,0.5,20)
ans =
4
Funzioni simili per gli altri filtri:
BUTTORD Butterworth filter order selection.
REMEZORD FIR order estimator (lowpass, highpass, bandpass, multiband)
32
PROGETTO DI FILTRI FIR
Filtri FIR:
•Possono avere fase lineare
•Stabili
•I metodi di progetto sono lineari
•Realizzabilita’ HW (solo parte MA)
•Transiente limitato (finite impulse response ! )
•Rispetto ai filtri IIR servono ordini più alti
33
Fc=1000;
Ftaglio=100;
Wn=Ftaglio/(Fc/2);
[b]=fir1(10,Wn);
a=1;
N=1000;
[H,f]=freqz(b,a,N,Fc);
figure(1)
clf
subplot(3,1,1)
plot(f,abs(H))
title(['modulo'])
xlabel('frequenza effettiva')
subplot(3,1,2)
plot(f,angle(H))
title(['fase'])
xlabel('frequenza effettiva')
subplot(3,1,3)
stem(b)
title('risposta impulsiva')
orient tall
0 50 100 150 200 250 300 350 400 450 5000
0.2
0.4
0.6
0.8
1modulo
frequenza effettiva
0 50 100 150 200 250 300 350 400 450 500-4
-2
0
2
4fase
frequenza effettiva
1 2 3 4 5 6 7 8 90
0.05
0.1
0.15
0.2
0.25
0.3
0.35risposta impulsiva
Filtro FIR passa-basso
[b] = fir1(n, Wn)
34
Idea: fittare dei campioni del modulo risposta in frequenza (f,M) con
un filtro FIR di ordine pre-fissato
Il ripple dei filtri ottenuti è uniforme in banda passante ed in banda
oscura
istruzione remez
Metodo di Parks-McClellan
35
REMEZ Parks-McClellan optimal equiripple FIR filter design.
B=REMEZ(N,F,A) returns a length N+1 linear phase (real, symmetric
coefficients) FIR filter which has the best approximation to the
desired frequency response described by F and A in the minimax
sense. F is a vector of frequency band edges in pairs, in ascending
order between 0 and 1. 1 corresponds to the Nyquist frequency or half
the sampling frequency. A is a real vector the same size as F
which specifies the desired amplitude of the frequency response of the
resultant filter B. The desired response is the line connecting the
points (F(k),A(k)) and (F(k+1),A(k+1)) for odd k; REMEZ treats the
bands between F(k+1) and F(k+2) for odd k as "transition bands" or
"don't care" regions. Thus the desired amplitude is piecewise linear
with transition bands. The maximum error is minimized.
B=REMEZ(N,F,A,W) uses the weights in W to weight the error. W has one
entry per band (so it is half the length of F and A) which tells
REMEZ how much emphasis to put on minimizing the error in each band
relative to the other bands.
36
Fc=1000;
f=[0 200 300 500]; %campioni asse frequenze
M=[1 1 0 0]; % campioni del modulo da fittare
Wn=f/(Fc/2);
ord=12; %ordine del filtro (restituisce ord+1 coefficienti)
[b]=remez(ord,Wn,M);
a=1;
N=1000;
[H,w]=freqz(b,a,N);
subplot(3,1,1)
plot(Fc*w/(2*pi),abs(H),f,M,'r*')
title(['modulo'])
xlabel('frequenza effettiva')
subplot(3,1,2)
plot(Fc*w/(2*pi),angle(H))
title(['fase'])
xlabel('frequenza effettiva')
subplot(3,1,3)
stem(b)
title('risposta impulsiva')
0 50 100 150 200 250 300 350 400 450 5000
0.2
0.4
0.6
0.8
1
1.2
1.4modulo
frequenza effettiva
0 50 100 150 200 250 300 350 400 450 500-4
-2
0
2
4fase
frequenza effettiva
1 2 3 4 5 6 7 8 9 10 11-0.1
0
0.1
0.2
0.3
0.4
0.5
0.6risposta impulsiva
b = remez(n,f,M)
Filtro FIR a fase lineare passa-basso (Parks-McClellan)
37
Fc=1000;
f=[0 200 300 500]; %campioni asse frequenze
M=[1 1 0 0]; % campioni del modulo da fittare
Wn=f/(Fc/2);
ord=30; %ordine del filtro (restituisce ord+1 coefficienti)
[b]=remez(ord,Wn,M);
a=1;
N=1000;
[H,w]=freqz(b,a,N);
subplot(3,1,1)
plot(Fc*w/(2*pi),abs(H),f,M,'r*')
title(['modulo'])
xlabel('frequenza effettiva')
subplot(3,1,2)
plot(Fc*w/(2*pi),angle(H))
title(['fase'])
xlabel('frequenza effettiva')
subplot(3,1,3)
stem(b)
title('risposta impulsiva')
0 50 100 150 200 250 300 350 400 450 5000
0.2
0.4
0.6
0.8
1
1.2
1.4modulo
frequenza effettiva
0 50 100 150 200 250 300 350 400 450 500-4
-2
0
2
4fase
frequenza effettiva
0 5 10 15 20 25 30 35-0.1
0
0.1
0.2
0.3
0.4
0.5
0.6risposta impulsiva
b = remez(n,f,M)
Filtro FIR a fase lineare passa-basso (Parks-McClellan)
38
Fc=2000;
f=[0 100 150 250 350 500 550 750 850 1000]; %campioni asse frequenze
M=[0 0 1 1 0 0 1 1 0 0]; % campioni del modulo da fittare
Wn=f/(Fc/2);
ord=30; %ordine del filtro (restituisce ord+1 coefficienti)
[b]=remez(ord,Wn,M);
a=1;
N=1000;
[H,F]=freqz(b,a,N,Fc);
subplot(3,1,1)
plot(F,abs(H),f,M,'r*',f,M,'g--')
title(['modulo'])
xlabel('frequenza effettiva')
subplot(3,1,2)
plot(F,angle(H))
title(['fase'])
xlabel('frequenza effettiva')
subplot(3,1,3)
stem(b)
title('risposta impulsiva')
orient tall
Filtro FIR a fase lineare multi-banda (Parks-McClellan)
0 100 200 300 400 500 600 700 800 900 10000
0.2
0.4
0.6
0.8
1
1.2
1.4modulo
frequenza effettiva
0 100 200 300 400 500 600 700 800 900 1000-4
-2
0
2
4fase
frequenza effettiva
0 5 10 15 20 25 30 35-0.2
-0.1
0
0.1
0.2
0.3
0.4
0.5risposta impulsiva
39
Filtro FIR a fase lineare multi-banda (Parks-McClellan)
Fc=2000;
f=[0 100 150 250 350 500 550 750 850 1000]; %campioni asse frequenze
M=[0 0 1 1 0 0 1 1 0 0]; % campioni del modulo da fittare
Wn=f/(Fc/2);
ord=60; %ordine del filtro (restituisce ord+1 coefficienti)
[b]=remez(ord,Wn,M);
a=1;
N=1000;
[H,F]=freqz(b,a,N,Fc);
subplot(3,1,1)
plot(F,abs(H),f,M,'r*',f,M,'g--')
title(['modulo'])
xlabel('frequenza effettiva')
subplot(3,1,2)
plot(F,angle(H))
title(['fase'])
xlabel('frequenza effettiva')
subplot(3,1,3)
stem(b)
title('risposta impulsiva')
orient tall
0 100 200 300 400 500 600 700 800 900 10000
0.2
0.4
0.6
0.8
1
1.2
1.4modulo
frequenza effettiva
0 100 200 300 400 500 600 700 800 900 1000-4
-2
0
2
4fase
frequenza effettiva
0 10 20 30 40 50 60 70-0.2
-0.1
0
0.1
0.2
0.3
0.4
0.5risposta impulsiva
SERVE UN’ORDINE ALTO
(problemi on-line, transienti…)
IIR
40
Filtro IIR a fase non lineare multi-banda
0 100 200 300 400 500 600 700 800 900 10000
0.2
0.4
0.6
0.8
1
1.2
1.4modulo
frequenza effettiva
0 100 200 300 400 500 600 700 800 900 1000-4
-2
0
2
4fase
frequenza effettiva
1 2 3 4 5 6 7 8 9 10 11-0.4
-0.3
-0.2
-0.1
0
0.1
0.2
0.3risposta impulsiva
% filtro IIR
Fc=2000;
f=[0 100 150 250 350 500 550 750 850 1000]; %campioni asse frequenze
M=[0 0 1 1 0 0 1 1 0 0]; % campioni del modulo da fittare
Wn=f/(Fc/2);
ord=10; %ordine del filtro (restituisce ord+1 coefficienti)
[b,a]=yulewalk(ord,Wn,M);
N=1000;
[H,F]=freqz(b,a,N,Fc);
subplot(3,1,1)
plot(F,abs(H),f,M,'r*',f,M,'g--')
title(['modulo'])
xlabel('frequenza effettiva')
subplot(3,1,2)
plot(F,angle(H))
title(['fase'])
xlabel('frequenza effettiva')