Analisi di segnali in SciPy

Marco D. Santambrogio – marco.santambrogio@polimi.itLuca Cerina – luca.cerina@mail.polimi.itEmanuele Del Sozzo – emanuele.delsozzo@polimi.itLorenzo Di Tucci – lorenzo.ditucci@mail.polimi.itGiuseppe Natale – giuseppe.natale@polimi.itMarco Rabozzi – marco.rabozzi@polimi.itAlberto Scolari – alberto.scolari@polimi.itVer. aggiornata al 18 Novembre 2016

SciPy

• SciPy è un modulo di Python orientato al calcolo scientifico

• I sotto-moduli di SciPy permettono di agire su diversi applicazioni, come interpolazione, integrazione, ottimizzazione, elaborazione di immagini, statistica, algebra lineare, analisi di segnali, …

• SciPy è disegnato per operare in modo efficiente su array NumPy, quindi conviene importare anche NumPy quando si usa SciPy

• La documentazione di ScyPy è disponibile a questo link:http://docs.scipy.org/doc/scipy/reference/

2

Moduli SciPy

I principali moduli disponibili all’interno di SciPy sono:

3

Sotto-modulo Applicazione

scipy.cluster Quantizzazionevettoriale/Kmeansscipy.constants Costantifisicheematematichescipy.fftpack TrasformatadiFourierscipy.integrate Integrazionescipy.interpolate Interpolazionescipy.io Input&outputdidatiscipy.linalg Algebralinearescipy.ndimage Immaginin-dimensionali

Moduli SciPy

Noi ci concentreremo solamente sull’elaborazione di segnali

4

Sotto-modulo Applicazione

scipy.odr Regressioneortogonalescipy.optimize Ottimizzazionescipy.signal Elaborazionedisegnaliscipy.sparse Matricisparsescipy.spatial Algoritmiestrutturedatispazialiscipy.special Funzionimatematichespecialiscipy.stats Statistica

Filtri a tempo discreto

• I filtri a tempo discreto possono essere classificati in Finite Impulse Response (FIR) e Infinite Impulse Response (IIR)

• SciPy permette di usare entrambi I tipi di filtri tramite varie funzioni, come firwin, firwin2, iirdesign e iirfilter

• SciPy mette a dispozione anche altri tipi di filtri (mediano, Wiener, Hilbert, ecc.)

• Guardate la documentazione o usare help(signal) per vedere quali altri filtri sono a disposizione

5

firwin

• La funzione firwin permette di creare un filtro basato sul metodo delle finestre

• In base agli argomenti dati alla funzione, è possibile ottenere diversi tipi di filtri (passa basso, passa banda, ecc.)

6

output=signal.firwin(numtaps,cutoff,width=None,window=‘hamming’,pass_zero=True,scale=True,nyq=1.0)

Lunghezzadelfiltro(numerodicoefficienti)

Frequenzeditagliofloatolista(banda)

Larghezzadellaregioneditransizione

Finestradausare

SeTrue,ilguadagnoallafrequenza0è1

SeTrue,Icoefficientisonoscalatiperottenere1comerisposta

infrequenzaadunacertafrequenza

FrequenzadiNyquist

NdarrayIcoefficientidel

filtroFIR

freqz

• Una volta calcolati i coefficienti, la funzione freqz, calcola le risposte in frequenza di un filtro digitale

7

w,h=signal.freqz(b,a=1,worN=None,whole=0,plot=None)

Numeratoredelfiltrolineare

Numerodifrequenzedacomputareintornoalcerchio

unitarioDefault=512frequenze

Se0,lefrequenzesonocalcolateda0allafrequenzadiNyquist

SeTrue,da0a2𝜋

E’possibilepassareunafunzionechedisegnidirettamenteweh

w:lefrequenzenormalizzaacuièstatacomputatah(ndarray)h:larispostainfrequenza(ndarray)

Denominatoredelfiltrolineare

Esempio filtro FIR usando la funzione firwin

8

Frequenzeditaglionormalizzate

Applicare filtri sui dati

• Una volta generato un filtro FIR/IIR, è possibile applicarlo su dei dati che devono essere filtrati

• Le principali funzioni che si possono usare sono:• lfilter: filtra i dati su una dimensione applicando un filtro

FIR/IIR• filtfilt: applica un doppio filtro sui dati

9

lfilter e filtfilt

10

output=signal.lfilter(b,a,x,axis=-1,zi=None)

b,a:vettoredeicoefficientidelnumeratore(denominatore) Datidiinput

Asselungoilqualeapplicareilfiltro

Condizioneinizialeperilritardodelfiltro

output=signal.filtfilt(b,a,x,axis=-1,padtype=‘odd’,padlen=None,method=‘pad’,irlen=None)

Datifiltrati

Datifiltrati

Esempio di filtro sui dati

11

Esempio di filtro sui dati

12

Analisi spettrale

• SciPy fornisce funzioni per calcolare la densità spettrale

• Ad esempio, si può usare il metodo del periodogramma o il metodo di Welch

13

f,Pxx=signal.periodogram(x,fs=1.0,window=None,nfft=None,detrend=‘constant’,return_onesided=True,scaling=‘density’,axis=-1)

f,Pxx=signal.welch(x,fs=1.0,window=‘hanning’,nperseg=256,noverlap=None,nfft=None,detrend=‘constant’,return_onesided=True,scaling=‘density’,axis=-1)

x =serietemporaledeivalorimisurati

fs =requenzadicampionamento

f =ArraydifrequenzecampionatePxx =densitàspettraledix

Esempio con metodo del periodogramma

14

Plotinscalalogaritmica

Punticasualisecondounadistribuzionenormalescale=deviazionestandardsize=shapedioutput

Esempio con metodo di Welch

15

Plotinscalalogaritmica

Punticasualisecondounadistribuzionenormalescale=deviazionestandardsize=shapedioutput

Altri filtri

• SciPy fornisce altri tipi di filtri che possono essere applicati ad immagini

• E’ possibile applicare filtri di covoluzione e correlazione:• Edge detection• Blur• …

16

Esempio filtro su immagine

17

Ilmodulomiscforniscevariefunzionalitàcomeimmaginisucuiapplicarefiltri

ApplicaunaconvoluzioneaidueinputusandolatrasformatadiFourierveloce

Accesso a file con NumPy e SciPy

• Anche NumPy e SciPy forniscono funzioni per agire su filehttp://docs.scipy.org/doc/numpy/reference/routines.io.html

• Le funzioni:np.load(nomeFile, …)

np.save(nomeFile, array, …)

servono per caricare e salvare un array NumPy su file .npy, .npz

18

Accesso a file con NumPy e SciPy

• La funzione:np.savez(nomeFile, *args, **kwds)

permette di salvare più array NumPy sotto forma di un dizionario (file .npz)

• Le funzioni:np.loadtxt(nomeFile, …)

np.savetxt(nomeFile, array, …)

servono per caricare e salvare un array NumPy su file di testo

19

Esempio accesso file con NumPy

20

Esercizio 1

Dato il segnale ecg presente nel file ecg.txt, scrivere un programma Python che applica un filtro FIR sul segnale per eliminare i 50HzIl segnale ecg è stato campionato a 1KHzRappresentare sia un grafico della risposta in frequenza che del segnale prima e dopo il filtro

21

Soluzione 1

22

MATERIALE AGGIUNTIVO SUI FILTRI

23

firwin2

• La funzione firwin2 permette di definire risposte in frequenza arbitrarie specificando una lista di frequenze e i loro relativi guadagni

24

output=signal.firwin2(numtaps,freq,gain,nfreqs=None,window=‘hamming’,nyq=1.0,antisymmetric=False)

Lunghezzadelfiltro(numerodicoefficienti)

Listadifrequenzecampione

ListadiguadagniFinestradausare Seilrisultanteimpulsoèsimmetrico

oantisimmetrico

FrequenzadiNyquist

NdarrayIcoefficientidel

filtroFIR

Dimensionedellameshdiinterpolazione

Esempio filtro FIR usando la funzione firwin2

25

Filtri IIR - iirfilter• SciPy fornisce due funzioni per I filtri IIR: iirdesign e iirfilter

• E’ possibile impostare il tipo di filtro da applicare

26

output=signal.iirfilter(N,Wn,rp=None,rs=None,btype=‘band’,analog=False,ftype=‘butter’,output=‘ba’)

Ordinedelfiltro

PerIfiltridiChebyshevoellittici,èl’oscillazione(ripple)massimanellabandapassante

PerIfiltridiChebyshevoellittici,èl’attenuazioneminima

nellabandaeliminata

TipodifiltroIIRdaprogettare:• Butterworth• ChebyshevI• ChebyshevII• Causer/elliptic• Bessel/Thomson

Tipodioutput:• Numeratore/denominatore(‘ba’)• Polo-zero(‘zpk’)• Sezionidisecondoordine(‘sos’)

Frequenzaangolare(perfiltrianalogici)Normalizzatatra0a1(perfiltridigitali)

Tipodifiltro:- bandpass- lowpass- highpass- bandstop

SeTrue,restituisceunfiltroanalogico,digitalealtrimenti

Filtri IIR - iirdesign

27

output=signal.iirdesign(wp,ws,gpass,gstop,analog=False,ftype=‘ellip’,output=‘ba’)

wp ews sono,rispettivamente,lefrequenzepassaedeliminabanda

gpass egstop sono,rispettivamente,laperditamassimainbandapassanteelaminimaattenuazioneineliminabanda

Stessiargomentiiirfilter

Esempio filtro IIR digitale

28

Esempio filtro IIR analogico

29

freqs calcolalarispostainfrequenzadiunfiltroanalogico