Esercitazioni MATLAB

Esercitazioni con MATLAB per il corso di TRATTAMENTO DEI SEGNALI Laurea Magistrale in Fisica

Dott. Silvia Maria Alessio - Dipartimento di Fisica Universit` di Torino - Anno Accademico 2011/2012 a 18 aprile 2012

IntroduzioneIn questo documento sono raccolte le esercitazioni che accompagnano il corso di Trattamento dei segnali. Esse si dividono in esercitazioni guidate, intese per essere svolte in classe sotto la diretta supervisione della docente, ed esercitazioni libere, che lo studente svolger` per proprio conto secondo le indicazioni fora nite. In seguito preparer` una breve relazione sui risultati ottenuti, da portare a allesame.

ii

Indice

Introduzione 1 Esercitazioni guidate 1.1 Un primo approccio a MATLAB . . . . . . . . . . . . . . . . . . 1.1.1 1.1.2 1.1.3 1.1.4 1.2 Premessa . . . . . . . . . . . . . . . . . . . . . . . . . . . Generalit` su MATLAB . . . . . . . . . . . . . . . . . . . a Lambiente di lavoro MATLAB . . . . . . . . . . . . . . . Generazione di dati sintetici . . . . . . . . . . . . . . . . .

ii 4 4 4 5 7 20 23 23 24 25 26 28 29 30 30 34 35

Trasformata discreta di Fourier . . . . . . . . . . . . . . . . . . . 1.2.1 1.2.2 1.2.3 1.2.4 1.2.5 1.2.6 Premessa . . . . . . . . . . . . . . . . . . . . . . . . . . . DFT e spettro di ampiezza: caratteristiche . . . . . . . . La function t di MATLAB . . . . . . . . . . . . . . . . . Prove con dati sintetici . . . . . . . . . . . . . . . . . . . Prova coi dati di pressione . . . . . . . . . . . . . . . . . . Aliasing . . . . . . . . . . . . . . . . . . . . . . . . . . . .

1.3

Filtraggio nel dominio del tempo . . . . . . . . . . . . . . . . . . 1.3.1 1.3.2 1.3.3 Premessa . . . . . . . . . . . . . . . . . . . . . . . . . . . Filtraggio via equazione alle dierenze . . . . . . . . . . . Filtraggio mediante convoluzione lineare . . . . . . . . . . 1

2 1.3.4 1.3.5 1.3.6 1.4

INDICE Filtraggio nel dominio della frequenza . . . . . . . . . . . Filtraggio a fase nulla . . . . . . . . . . . . . . . . . . . . Filtraggio con decimazione . . . . . . . . . . . . . . . . . 36 37 37 40 40 41 43 43 44 44 45 46 48 48 48 49 59 61 61 61 62 63 63 65 65

Caratteristiche di un ltro numerico . . . . . . . . . . . . . . . . 1.4.1 1.4.2 1.4.3 1.4.4 Premessa . . . . . . . . . . . . . . . . . . . . . . . . . . . Risposta in frequenza . . . . . . . . . . . . . . . . . . . . Risposta allimpulso . . . . . . . . . . . . . . . . . . . . . Deconvoluzione . . . . . . . . . . . . . . . . . . . . . . . .

1.5

Progetto di ltri FIR equiripple . . . . . . . . . . . . . . . . . . . 1.5.1 1.5.2 1.5.3 1.5.4 Premessa . . . . . . . . . . . . . . . . . . . . . . . . . . . Procedura di progetto: un passa-banda . . . . . . . . . . Stima dellordine minimo . . . . . . . . . . . . . . . . . . Progetto di ltri passa-basso e passa-alto . . . . . . . . .

1.6

Metodi di stima dello spettro di potenza . . . . . . . . . . . . . . 1.6.1 1.6.2 1.6.3 Premessa . . . . . . . . . . . . . . . . . . . . . . . . . . . Metodi non-parametrici . . . . . . . . . . . . . . . . . . . Metodi parametrici . . . . . . . . . . . . . . . . . . . . . .

2 Esercitazioni libere 2.1 Analisi di una serie temporale . . . . . . . . . . . . . . . . . . . . 2.1.1 2.1.2 2.1.3 2.1.4 2.1.5 2.1.6 Premessa . . . . . . . . . . . . . . . . . . . . . . . . . . . Analisi statistica elementare . . . . . . . . . . . . . . . . . Autocorrelazione . . . . . . . . . . . . . . . . . . . . . . . Spettro di potenza . . . . . . . . . . . . . . . . . . . . . . Signicativit` dei picchi . . . . . . . . . . . . . . . . . . . a Whitening . . . . . . . . . . . . . . . . . . . . . . . . . . .

INDICE 2.1.7 2.1.8 2.2 Rumore rosso . . . . . . . . . . . . . . . . . . . . . . . . . Stazionariet` . . . . . . . . . . . . . . . . . . . . . . . . . a

3 66 67 68 68 69 73 77 78 79 80

Analisi spettrale evolutiva: la CWT . . . . . . . . . . . . . . . . 2.2.1 2.2.2 2.2.3 2.2.4 2.2.5 2.2.6 2.2.7 Premessa . . . . . . . . . . . . . . . . . . . . . . . . . . . Le function . . . . . . . . . . . . . . . . . . . . . . . . . . Il programma . . . . . . . . . . . . . . . . . . . . . . . . . Interpretazione dello scalogramma . . . . . . . . . . . . . Prove con wavelets diverse . . . . . . . . . . . . . . . . . . Spettro globale di wavelet . . . . . . . . . . . . . . . . . . Scelta delle scale . . . . . . . . . . . . . . . . . . . . . . .

Capitolo 1

Esercitazioni guidate1.11.1.1

Un primo approccio a MATLABPremessa

Lo scopo dellesercitazione ` quello di imparare ad usare MATLAB 1 ad un e livello di base, utilizzando questo breve testo per provare al computer alcune operazioni elementari, seguendo passo per passo le istruzioni date.

Attenzione: il testo che segue ` volto esclusivamente a fornire agli studenti e del corso di Trattamento dei Segnali le nozioni indispensabili per eseguire le esercitazioni. Non va pertanto considerato una introduzione a MATLAB con caratteristiche di generalit` e completezza. a

Nella prossima sottosezione viene data unidea generale della struttura di MATLAB e delle parti che lo compongono; nel terzo si accede a MATLAB, si introducono al suo interno alcune variabili e si impara a maneggiarle e a gracarle in modo semplice. Segue una sottosezione in cui si parla di generazione di segnali sintetici (cio` dati simulati matematicamente, non frutto di vere mie sure). Nellelaborazione numerica dei segnali si presenta spesso la necessit` di a generare tali segnali, per utilizzarli al ne di testare algoritmi e sistemi. poich, come si vedr`, dati possono essere importati in ambiente MATLAB e poi e a esportati, si consiglia di creare nella propria area di lavoro una directory o cartella apposita, in cui poi risiederanno tutti i les relativi agli scambi tra MATLAB e lambiente di lavoro esterno (les da importare e les esportati). Durante le esercitazioni si useranno alcuni les (contenenti set di dati, ltri numerici e pro1 Ci

si riferisce alla versione 7.8.0, R2009a di MATLAB.

4

1.1. UN PRIMO APPROCCIO A MATLAB

5

grammi) che lo studente dovr` copiare nella propria area di lavoro, prelevandoli a dalla pagina web del Corso su Campusnet (http://sica.campusnet.unito.it) 2 .

1.1.2

Generalit` su MATLAB a

MATLAB ` un sistema software interattivo per e calcolo, visualizzazione e graca, analisi ed elaborazione di dati, simulazione e modellazione, programmazione in un proprio linguaggio. Lelemento di base ` una matrice che non richiede dimensionamento: MATLAB e sta in eetti per Matrix Laboratory ed in esso tutte le variabili sono matrici. Il sistema software MATLAB ` costituito da un pacchetto fondamentale (MAe TLAB Toolbox), pi` numerosi Toolboxes per applicazioni speciche: ad esempio u elaborazione numerica dei segnali, wavelets, statistica, ecc... Il pacchetto-base include cinque parti principali. 1. Lambiente di lavoro MATLAB (Development Environment): set di strumenti che consentono di usare les e function di MATLAB (le function sono programmi che implementano i vari algoritmi; si veda pi` avanu ti). Tramite tali strumenti si manipolano, nel cosiddetto workspace di MATLAB, che ospita linsieme delle variabili denite in un dato istante di una sessione di lavoro, le variabili stesse; si importano ed esportano dati nel e dal workspace, ecc... Il desktop di MATLAB include la Command Window (nestra dove, al prompt di MATLAB, si possono digitare direttamente i comandi); unaltra nestra col diario dei comandi dati in precedenza (Command History); un elenco dei componenti installati 3 ; browser per accedere allhelp, per elencare le variabili contenute nel workspace, i les presenti nelle directories/cartelle... Oltre a lavorare da linea comandi, lutente pu` accedere a o vari strumenti specialistici (Tools, come ad es. sptool relativo al Toolbox di Signal Processing) che si presentano come Graphical User Interfaces (GUIs) 4 .2 Il lettore che non disponga delle credenziali necessarie pu` richiedere il materiale o allautrice, tramite e-mail allindirizzo alessioph.unito.it. 3 I Toolboxes specialistici si acquistano separatamente dal pacchetto di base. 4 In queste esercitazioni si lavorer` prevalentemente da linea comandi e/o tramite semplici a programmi: si ` ritenuto pi` istruttivo per lo studente accedere cos` piuttosto che dalle GUIs, e u , ai vari algoritmi. Tuttavia le GUIs sono molto comode e lo studente ` caldamente incoraggiato e ad esplorarle.

6

CAPITOLO 1. ESERCITAZIONI GUIDATE 2. La libreria di funzioni matematiche (function): algoritmi di calcolo che vanno dalle funzioni elementari, come seno e coseno, a funzioni pi` u complicate, come ricerca degli autovalori di una matrice o FFT. Le function sono insiemi di istruzioni, cio` programmi in linguaggio MAe TLAB, indicati come M-les perch caratterizzati da una estensione .m. e Ogni Toolbox non ` nientaltro che una raccolta di M-les. Tra le categorie e di M-les disponibili nel pacchetto-base citiamo: matrici e loro manipolazione;

funzioni matematiche elementari; analisi dati e FFT;

funzioni matematiche specialistiche; interpolazione, polinomi;

soluzione di equazioni dierenziali, ecc... Gli M-les corrispondenti a funzioni relative ad un certo argomento si trovano raggruppate in una singola directory/cartella di MATLAB; per esempio, funzioni della libreria matematica si trovano in elmat manipolazione elementare di matrici, funzioni matematiche elementari,

datafun analisi dati e FFT,

elfun

polyfun interpolazione e polinomi.

3. Il sistema graco (Handle Graphics): comprende comandi high level, per visualizzare dati in due e tre dimensioni (e altro, p.es. animazioni) e comandi low level per personalizzare i graci ottenuti. Esiste anche la possibilit` di costruire GUIs (Graphical User Interfaces), con menu per a gestire le applicazioni personali, sviluppate dallutente usando il linguaggio di programmazione MATLAB. 4. Il linguaggio di programmazione MATLAB, che serve allutente per sviluppare i propri algoritmi di calcolo, i quali vengono posti in nuovi M-les (script e function). Le espressioni matematiche coinvolgono intere matrici, per cui si ha una grande concisione nel linguaggio, con conseguente rapidit` di programmazione. a Uno script le ` un M-le contenente una sequenza di istruzioni in line guaggio MATLAB. Supponiamo di aver creato, con un editor qualsiasi o con quello disponibile in MATLAB (v. descrizione del desktop nella prossima sottosezione), un nuovo M-le detto nuovo.m. Dando il comando >>nuovo le istruzioni verranno eseguite; si tratta dunque dellanalogo di una subroutine del Fortran. Un function le ` anchesso un M-le, salvo che in testa ha una denition e line in cui vengono deniti esplicitamente linput e loutput: per esempio function h = circle(x,y,r) %CIRCLE=function to draw a circle of radius r centered in x,y

1.1. UN PRIMO APPROCCIO A MATLAB %Call syntax: circle(x,y,r); %Input: (x,y) = center coordinates; r = radius %Output: none th = 0:pi/50:2*pi; xunit = r * cos(th) + x; yunit = r * sin(th) + y; h = plot(xunit, yunit); ` E quindi analogo ad una function Fortran.

7

5. La Application Program Interface (API): una libreria che consente di scrivere programmi in FORTRAN e C che interagiscono in vari modi con MATLAB. Le modalit` duso di MATLAB pi` complicate non verranno prese in considea u razione in queste pagine. Infatti per potere svolgere le esercitazioni gli studenti hanno soltanto bisogno di sapere: come introdurre dati nel workspace e come estrarli; come maneggiare le variabili nel workspace, anche mediante semplici espressioni matematiche che le contengano; tali espressioni spesso includeranno function di MATLAB; come produrre e stampare semplici graci bidimensionali, con assi dalle caratteristiche a piacere, delle label per gli assi stessi ed eventualmente un titolo; come scrivere e far girare semplici programmi che svolgano i compiti elencati ai punti precedenti. Pertanto nel seguito parleremo esclusivamente di quanto concerne queste operazioni di base.

1.1.3

Lambiente di lavoro MATLAB

Questa sottosezione descrive come maneggiare le variabili nel workspace e come importare ed esportare dati. 1. Per entrare in MATLAB battere al prompt del sistema 5 : matlab Appare il desktop di MATLAB, che presenta, nelle condizioni di default :5 Dora

in poi, si consiglia di dare eettivamente i comandi citati!

8

CAPITOLO 1. ESERCITAZIONI GUIDATE in alto, la barra coi vari menu (File, Edit, ..., Desktop, Window, Help), il pulsante per invocare lhelp (?), leditor per M-les (icona col foglio bianco), la nestra per vedere o cambiare la directory/cartella corrente ecc... al centro, la Command Window col prompt (>>); a sinistra, lelenco dei les contenuti nella Current Directory; a destra in alto, il contenuto del workspace; a destra in basso, la Command History (storia dei comandi dati, anche nelle sessioni precedenti). Le varie nestre si possono chiudere per far pi` posto ad altre; per variare u laspetto del desktop andare su: Desktop/Desktop Layout. Nella Command Window si danno i comandi, si fanno girare gli M-les 6 , si inseriscono dati... Alcuni comandi utili: >>type nome_file mostra le istruzioni in linguaggio MATLAB contenute nel M-le il cui nome ` specicato nella stringa nome_file. Il comando non funziona per le e funzioni built-in di MATLAB, funzioni di base come seno, coseno, t ecc... >>what matlab/elfun elenca i les contiene una data directory/cartella di MATLAB (nellesempio, quella delle funzioni matematiche elementari). Inoltre funzionano sia il comando Unix ls, sia lanalogo MATLAB, che ` dir. e In ogni istante della sessione di lavoro si pu` accedere alla documentao zione su Web e allhelp per vedere quali funzioni sono disponibili e come utilizzarle. Esiste, come si ` detto, un apposito browser per lhelp; qui e facciamo solo alcuni esempi da linea comandi: >>helpbrowser richiama il browser dellhelp; >>help nome_function d` sintassi e comportamento della funzione il cui nome ` dato nella stringa a e nome_function; analogamente, >>help nome_comando d` informazioni sul comando stesso; a >>doc nome_function apre la pagina della documentazione relativa alla funzione il cui nome ` e dato nella stringa nome_function; analogamente, >>doc nome_comando apre la pagina del comando stesso; >>help elfun elenca e descrive brevemente le funzioni in una particolare directory/cartella,

6 Quando MATLAB cerca un le, lo fa secondo un certo path pressato; ogni le che si vuole chiamare deve essere o nella directory/cartella corrente, o sul path di ricerca. Per vedere quale sia il path o apportarvi cambiamenti: File/Set Path.


9

qui elfun 7 ; >>help produce una lista delle directories di MATLAB con descrizione della categoria di funzioni contenute in ciascuna; >>lookfor parola_chiave cerca, tra le funzioni disponibili, quella che ha qualche cosa a che fare con la parola_chiave data. Infatti ogni function ` corredata di alcune righe e di intestazione, in cui viene spiegato lo scopo e luso della function stessa; il comando lookfor esamina la prima riga di ogni intestazione (header ) alla ricerca della data parola_chiave. Per inciso, si noti che quando si prepara un nuovo M-le, conviene sempre mettere in testa (in uno script) o dopo la linea di denizione (in una function) un header, ossia una linea commentata che spieghi concisamente lo scopo del le. In tal modo il comando lookfor potr` operare anche a sugli M-les preparati dellutente. Attenzione: quando, come accadr` nel seguito delle esercitazioni, si deve a ripetere pi` volte una certa sequenza di comandi, ` sempre consigliabile u e memorizzarla in un nuovo M-le; cos` si risparmia molto lavoro. Notiamo anche che volendo modicare una function di MATLAB (sia ad es. cwt.m), si pu` copiarla nella propria directory/cartella di lavoro ed edio tarla. MATLAB utilizzer` automaticamente la cwt.m che si trova nellarea a di lavoro dellutente: infatti quando viene dato il comando cwt, MATLAB cerca cwt.m prima nellarea dellutente e solo dopo, non trovandola, cerca nellarea di installazione. Per richiamare comandi gi` dati (e poi eventualmente modicarli) si usano a le frecce sulla tastiera; per interrompere lesecuzione di un comando, >>Ctrl c 2. Per inserire a mano una matrice 8 , sia A: si digita ad esempio >>A=[1 2 3; 4 5 6; 7 8 10] In risposta, MATLAB visualizza la matrice cos` introdotta nel workspace: A= 1 2 3 4 5 6 7 8 10 Cliccando poi sul nome della matrice nel workspace browser, si apre lArray Editor (editor di matrici 9 ) che consente di modicarla. Per indicare ( 1) nei numeri complessi si usa i oppure j, indierentemente. Bisogna per`, naturalmente, fare attenzione a non ridenire o incautamente queste variabili nei propri script e function, usandoli come indici in cicli di for (che sono gli analoghi MATLAB dei cicli di do del FORTRAN), ecc... poich inoltre alcune function di MATLAB (come, per e esempio, cwt.m del Wavelet Toolbox Wavelet) usano i o j come indice, ` e7 Vi sono contenute le funzioni matematiche elementari; tuttavia le funzioni pi` basilari, u come seno e coseno, sono built-in, cio` parte del core di MATLAB, molto ecienti ma non e visibili come M-les. 8 Attenzione: MATLAB distingue le maiuscole dalle minuscole (come si dice, ` case e sensitive), per cui la variabile A ` diversa dalla variabile a. e 9 A tutti gli eetti pratici possiamo usare come sinonimi le parole matrix e array.

10

CAPITOLO 1. ESERCITAZIONI GUIDATE bene controllare durante la sessione di lavoro il valore di i e j, prima di o usarli come 1. In alternativa, un numero complesso pu` essere ottenuto dalle sue parti reale ed immaginaria mediante la function complex. Se non si desidera la visualizzazione del risultato di una operazione, si mette dopo il comando un punto e virgola: ad esempio >>A=[1 2 3; 4 5 6; 7 8 10]; Si possono anche dare due o pi` comandi in una stessa linea, separandoli u con una virgola: ad esempio, per assegnare due scalari scal1e scal2 si pu` digitare o >>scal1=1,scal2=2 Per continuare a nuova riga una linea comandi troppo lunga si usano tre punti: ad esempio, supponendo di aver denito degli scalari da a a g e dei vettori di ugual lunghezza da h a l, potremo scrivere >>s=a+b+c+d+e+f+g+10*log10(abs(fft(... h+i+l))) 3. Uso delle matrici e relativi indici: in MATLAB tutte le variabili sono matrici. Uno scalare ` una matrice 1x1; un vettore ` una matrice 1xN e e (vettore-riga) oppure Nx1 (vettore-colonna). Esistono non solo matrici di numeri, ma anche di stringhe character. Le dimensioni (M,N) di una matrice A contenuta nel workspace possono essere rese note e memorizzate col comando >>[M,N]=size(A) Infatti quando si desidera assegnare pi` variabili in output, le si include u in parentesi quadre, separandole con virgole. Quando non si assegna esplicitamente alcuna variabile di output, il risultato di una operazione viene di default assegnato alla variabile ans (da answer, risposta), il cui valore viene visualizzato. Se si considera, anzich` una matrice vera e propria, un semplice vettore, e ad esempio >>v=[1 2 3 4 5] per conoscerne le dimensioni si pu` usare anche la function length: o >>L=length(v) Una matrice reale (come A o v) pu` venire trasposta digitando o >>A=A, v=v; In questo modo un vettore-riga diviene un vettore-colonna, e viceversa. Se A ` complessa, A ` la trasposta complesso-coniugata, mentre la seme e plice trasposta ` A.. e Luso degli indici per maneggiare gli elementi di matrice obbedisce alle usuali regole: detta A la matrice, A(i,j) ` lelemento nella riga i e coe lonna j. Lultima riga e lultima colonna di una matrice possono essere indicate con end: ad esempio, A(1,end) indica lelemento nella prima riga e ultima colonna. Se si tenta di usare indici fuori dallintervallo denito in precedenza, si ha un errore, ma se si deniscono nuovi elementi fuori range, la matrice originaria viene allargata no a contenere i nuovi elementi. Gli elementi non deniti, posti tra quelli pre-esistenti e quelli deniti ex-novo, vengono uguagliati a zero. Ad esempio: >>A

1.1. UN PRIMO APPROCCIO A MATLAB (matrice con tre righe e tre colonne denita in precedenza...) >>A(4,4)=25 (matrice allargata).

11

Quanto al numero massimo di dimensioni di una matrice, le uniche limitazioni sono dovute alla memoria disponibile ed al massimo numero di elementi non dimensioni in una matrice. La memoria disponibile dipende dal computer in uso; il numero massimo di elementi permesso dipende dalla versione di Matlab installata. Per conoscere il numero massimo di elementi (maxsize) si usa il comando computer che fornisce informazioni sul sistema in uso: [C,maxsize]=computer Ogni matrice denita dovr` dunque avere dimensioni tali che il loro proa dotto (numero di righe numero di colonne ...) sia inferiore a (maxsize). Invece per conoscere i limti di memoria sul proprio sistema si usa il comando memory. ` E possibile maneggiare parti di matrici: ad esempio A(1:k,j) indica gli elementi nella colonna j, dalla prima alla k-esima riga; A(:,j) indica tutti gli elementi della colonna j-esima. Si possono, inne, cancellare righe e colonne: ad es. >>A(:,2)=[] cancella la seconda colonna di A. 4. Operazioni elementari con matrici (quindi con qualsiasi variabile): esempi. Nota: le operazioni qui citate sono fatte elemento per elemento. Esse richiedono che gli operandi abbiano le stesse dimensioni; se un operando ` e uno scalare, viene applicato a tutti gli elementi dellaltro operando (scalar expansion). Quando utilizzate insieme, in una determinata espressione, tutte queste operazioni hanno uguale precedenza e vengono valutate da destra verso sinistra; si pu` cambiare questordine di precedenza con luso di parentesi. o Addizione elemento per elemento: >>a=[0 1 0 1 0 1] >>b=[2 2 2 2 2 2] >>a+b Sottrazione elemento per elemento: >>a-b Cambiamento di segno: >>-a Moltiplicazione elemento per elemento >>a.*b Divisione elemento per elemento: >>a./b10 Lo stesso comando senza il punto implementa invece il prodotto matriciale, righe per colonne.

10

:

12

CAPITOLO 1. ESERCITAZIONI GUIDATE Elevamento a potenza elemento per elemento: >>b.^2 Radice quadrata elemento per elemento: >>sqrt(b) che ` come e >>b.^(1/2) ma pu` risultare pi` accurato. o u Addizione di uno scalare scal a tutti gli elementi di una matrice a: >>a+scal (si ricordi la scalar expansion). Analogamente si procede per la sottrazione. Logaritmo elemento per elemento, in base e, in base 2 o in base 10: >>log(b) oppure >>log2(b) oppure ancora >>log10(b) 5. Importazione di dati in MATLAB: per farla esistono diverse possibilit`. Limitiamoci alle pi` semplici, menzionando solo il fatto che a u in MATLAB esiste anche un Import Wizard, richiamabile col comando uiimport. Fornire i dati come lista esplicita di elementi, come gi` visto; ad a esempio: >>B=[1 1 1 0] Usare leditor preferito oppure quello di MATLAB per creare un M-le che faccia la stessa cosa. Ad esempio si pu` editare un le o contenente listruzione B=[1 1 1 0], salvarlo come B.m, poi dare il comando >>B Usare lArray Editor. Caricare i dati da un le ascii esterno contenente solo la variabile da caricare. Questa diverr` una matrice MxN, se il le conteneva M righe a di N dati ciascuna. Il nome della variabile sar` quello del le, privato a dellestensione. Per provare il procedimento: copiare nella propria directory/cartella di lavoro il le pd.dat, che contiene una serie temporale di 25921 dati di pressione differenziale in Pascal 11 ; lintervallo di campionamento ` Tc = e 5 s; copiare anche tr.dat, che contiene i contemporanei dati di trasparenza atmosferica registrati con lo stesso Tc ; digitare >>load pd.dat >>load tr.dat e controllare come si presentano le nuove variabili cos` caricate.

11 Si tratta di microuttuazioni di pressione atmosferica, misurate con un microbarometro, registrate in inverno a Trino Vercellese (Torino) durante un esperimento volto a rivelare linuenza di onde di gravit` in atmosfera sulla visibilit`, durante episodi di nebbia. a a


13

Allo scopo, da linea comandi si pu` usare o >>who oppure >>whos che per ogni variabile d` anche le dimensioni ed il tipo di array: a double array= matrice normale, char array=matrice character... 6. Workspace e uso della memoria: le variabili create o importate nel workspace si possono rimuovere; ad es., per cancellare la variabile A si digita >>clear A; mentre se si vuole ripulire tutto il workspace, >>clear Le function utilizzano, per le loro variabili, unarea di memoria distinta dal workspace principale (tant` vero che le variabili denite allinterno e di una function non sono visibili col comando who). Per ripulire anche questarea di memoria, si d` il comando a >>clear function 7. Funzioni per generare matrici speciali: >>z=zeros(2,4) z= 0 0 0 0 0 0 0 0 >>q=ones(2,3) q= 1 1 1 1 1 1 >>r=rand(3,2) (matrice con tre righe e due colonne di numeri pseudocasuali con distribuzione uniforme, compresi tra zero e uno). >>r=randn(10,1) (matrice con dieci righe e una colonna di numeri pseudocasuali con distribuzione gaussiana a media nulla e varianza unitaria). La routine rand produce numeri pseudocasuali; la sequenza di numeri eettivamente generati ` determinata dallo stato del generatore. poich e e questo ha memoria, ogni volta che lo si chiama esso genera una sequenza diversa. Se vogliamo fargli generare due volte la stessa sequenza, dobbiamo leggere lo stato interno del generatore: rand_state=rand(state); generare la prima sequenza: x1=rand(10,1); resettare lo stato interno del generatore: rand(state,rand_state); generare la seconda sequenza: x2=rand(10,1); Daltra parte, poich lo stato interno del generatore viene resettato ogni e volta che MATLAB viene lanciato, la sequenza di numeri generata la prima

14

CAPITOLO 1. ESERCITAZIONI GUIDATE volta che si chiama rand in ogni sessione di lavoro sar` la stessa, a meno di a cambiare lo stato del generatore stesso. Considerazioni analoghe valgono per la function randn. 8. Funzioni speciali che danno i valori di certe costanti: pi, ;

i oppure j, ( 1) (rammentare il pericolo di ridenirli usandoli come indici);

realmin, il pi` piccolo numero reale ammesso; u realmax, il pi` grande numero reale ammesso; u Inf (o inf), innito;

eps, precisione relativa oating;

NaN (o nan), not-a-number: pu` essere, ad esempio, usato per rapo presentare dati mancanti in una serie temporale; la sua posizione nel vettore o matrice pu` poi essere rintracciata con luso di operatori o logici (v. oltre). Questi valori si possono ridenire, ad es. digitando >>eps=1.e-6, e poi riportare ai valori di default: >>clear eps 9. Espressioni matematiche: da linea comandi si possono usare le varie funzioni. Ne citiamo alcune: (a) funzioni trigonometriche: sin, asin, cos, acos, tan, atan, cot, acot ecc...; (b) funzioni esponenziali: exp, log, ,log2, log10, sqrt; (c) funzioni complesse: abs, angle, conj, real, imag, complex; (d) funzioni di arrotondamento: fix (arrotondamento verso lo zero), floor (arrotondamento verso lintero pi` vicino nella direzione di u meno innito), ceil (arrotondamento verso lintero pi` vicino nella u direzione di pi` innito), round (arrotondamento verso lintero pi` u u vicino), rem (resto dopo divisione), sign (funzione segno); (e) funzioni di statistica ed analisi dati elementare, quali la media e la deviazione standard 12 . Ad esempio: >>mupd=mean(pd) trova la media del vettore pd; >>mupd=median(pd) trova la mediana; >>sigmapd=std(pd)12 In MATLAB, le funzioni relative a questo argomento risiedono nella directory/cartella datafun e lavorano anche su vere matrici, non soltanto su vettori. A questo proposito, notiamo che per lanalisi statistica, in generale multivariata, MATLAB ` column oriented: assume che e la matrice dei dati abbia per colonne le singole sequenze (p.es. serie temporali) di cui i dati sono composti. Pertanto le varie operazioni (di media, deviazione standard ecc...) vengono condotte per colonne e il risultato ` un vettore riga. Quali vettori di dati, MATLAB accetta e per altro sia vettori riga, sia vettori colonna.


15

trova la deviazione standard; inoltre vi sono min (minimo assoluto), max (massimo assoluto), sum (somma degli elementi), diff (dierenze, prima approssimazione della derivata), trapz (integrazione numerica trapezoidale), ecc... Quando si cerca il massimo o minimo assoluto, si pu` nel contempo o trovare lindice del vettore a cui questo si verica: >>[mapd,indmapd]=max(pd) e analogamente >>[mipd,indmipd]=min(pd) Operando solo su parte del vettore, si possono limitare le operazioni precedenti a un dato intervallo di indici. Per togliere la media agli elementi di un vettore di dati (cio` per e centrare la sequenza), si operer` con la scalar expansion: a >>pd_c=pd-mupd; Invece per togliere la media ad una matrice vera e propria di dati bisogna matricizzare loperazione, come nellesempio che segue. >>A=[1 4 7; 2 5 8; 3 6 9] >>muA=mean(A) (vettore riga) >>[M,N]=size(A) >>e=ones(N,1) (vettore colonna) >>A_c=A-e*muA (A meno il prodotto matriciale del vettore-colonna e, fatto di N valori unitari, con il vettore-riga muA). 10. Formato dei valori numerici che appaiono sullo schermo (i calcoli vengono comunque sempre fatti in doppia precisione): ` possibile sceglierlo e col comando format. Ad esempio >>x=[4/3 1.2345e-6]; >>format short >>x 1.3333 0.0000 >>format long >>x 1.33333333333333 0.00000123450000 >>format short e >>x 1.3333e+00 1.2345e-06 ecc... Il default ` il formato short. e 11. Espressioni con operatori relazionali e logici: in esse si confrontano elementi analoghi di matrici equidimensionali (salvo scalar expansion). Gli operatori relazionali sono: < > >= >A=[1 2 3 4] >>B=[1 0 0 4] >>A==B ans= 1001 Sebbene questi operatori vengano normalmente usati nella programmazione in linguaggio MATLAB, in conditional statements quali if - then - else per ramicare algoritmi, essi possono essere utilizzati per eseguire manipolazioni di matrici: per esempio >>v=[1.2 0.8660 0.3] >>u=v(v>=sin(pi/3)) trova tutti gli elementi vi del vettore v tali che vi sin( ) e li memorizza 3 nel vettore u. Gli operatori relazionali vengono valutati dopo gli operatori aritmetici. Gli operatori logici sono: & and | or ~ not ed operano come quelli relazionali. Esistono anche function logiche, su cui non ci soermiamo; citiamo solo la function find, che determina gli indici degli elementi di una matrice che soddisfano ad una data condizione logica: ad esempio >>i=find(A>2) i= 3 4 e la function isfinite(x) che esamina la variabile x alla ricerca di valori Inf o NaN: essa diverr` true (ossia 1) in corrispondenza dei valori numerici a e false (ossia 0) in corrispondenza di quelli Inf o NaN.

12. Esportazione di dati da MATLAB: il contenuto del workspace pu` o essere salvato in un le binario, detto MAT-le, con estensione .mat. Per esempio, da linea comandi, >>save 211212 salver` il contenuto del workspace nel le 211212.mat. In seguito, riena trando in MATLAB dopo esserne usciti, tutte le variabili salvate potranno essere ricollocate nel workspace con il comando load. Proviamo questa procedura: >>quit per uscire da MATLAB (attenzione a non farlo mai senza aver prima salvato il risultato delle proprie fatiche!!!). Poi al prompt del sistema matlab per rientrare; >>who per vericare che il workspace ora ` vuoto; e >>load 211212 per restaurare il precedente contenuto del workspace; >>who per vericarlo.


17

Attenzione: variabili omonime pre-esistenti, denite tra il momento del riingresso in MATLAB ed il momento del caricamento, verranno ricoperte. Si pu` anche salvare solo certe variabili, ad esempio o >>save file1 v A (dove v, A sono le variabili da salvare), oppure decidere di salvare una variabile alla volta in un le ascii, ad esempio >>save A.dat A -ascii (8 digits) o >>save A.dat A -double (16 digits). In seguito la variabile potr` tornare nel workspace con a >>load A.dat 13. Graci in due dimensioni, nozioni essenziali: si utilizza la function plot. Si voglia gracare la variabile pd in funzione del tempo, in scala linearelineare. Allora occorre prima di tutto costruire un vettore che rappresenti la variabile indipendente (in questo caso, un vettore di tempi equispaziati, espressi, p.es., in secondi). Per questo si possono usare due comandi, di seguito esemplicati. (a) >>M=25921 >>T_c=5. >>t_i=0. >>t_f=(M-1)*T_c >>t=t_i:T_c:t_f; che crea un vettore-riga t contenente i valori da t_i a t_f, con passo T_c (se non specicato vale 1.); (b) >>M=25921 >>T_c=5. >>t_i=0. >>t_f=(M-1)*T_c >>t=linspace(t_i,t_f,M); che crea un vettore-riga con valori da t_i a t_f, contenente M valori. Ora per plottare il vettore pd in funzione di t si digita >>plot(t,pd) ed il graco appare, in una nestra separata. Al graco si pu` poi assegnare un titolo e delle label per gli assi con o lindicazione dellunit` di misura: questa ed altre operazioni sono rese a agevoli dai comandi e dai menu della barra che appare nella nestra della gura. Citiamo ad esempio Edit per FigureProperties e AxesProperties, Insert per le label degli assi, la lente per zoomare, ecc... Naturalmente le stesse operazioni sono possibili anche da linea comandi: ad esempio >>xlabel(Tempo [s]); ylabel(Pressione diff. [Pa]) >>title(Dati esperimento Trino 21.12.89) 1313 In

questi comandi il testo tra apici viene interpretato secondo le regole di Te.

18

CAPITOLO 1. ESERCITAZIONI GUIDATE Se si desidera un grigliato di fondo, >>grid on mentre per disattivare il grigliato, >>grid off Si possono anche tracciare pi` linee su un graco, in stili o colori diversi; ad u esempio possiamo fare il graco di pd e tr sovrapposti. Lunica avvertenza ` che trattandosi in questo caso di grandezze siche diverse, dovremo prima e (centrarle e) standardizzarle, dividendole ciascuna per la loro deviazione standard. In questo modo diverranno variabili adimensionali confrontabili. Supponendo di averlo fatto, daremo poi il comando plot(t,pd,k-,t,tr,b--) graca le variabili rispettivamente in nero continuo e blu tratteggiato. Il graco pu` riguardare anche solo una parte dei dati: ad esempio o MM=fix(M/3); plot(t(1:MM),pd(1:MM),k-,t(1:MM),tr(1:MM),b--) che graca insieme il primo terzo dei valori di pressione dierenziale e trasparenza. In questo contesto pu` essere utile una legenda: o >>legend(Pressione differenziale,Trasparenza,2) (il 2 specica la posizione della legenda nel graco, in alto a sinistra). Altri possibili stili per le linee nei graci sono: : (puntinata) e -. (tratto e punto). Si possono pure fare graci per punti, usando i seguenti markers: +, o, *, x, ., s (square=quadrato), d (diamond=rombo), ^ (triangolo), ecc... Per avere markers sovrapposti a linee, occorre specicare entrambi: ad es. :^ dar` una linea puntinata con triangoli in corrispondenza dei valori. a Si pu` scegliere il colore delle linee che si tracciano; ad esempio o >>plot(t(1:200),pd(1:200),g) dar` una linea verde (su uno schermo a colori, linee multiple in un graco a appaiono automaticamente di colori diversi). Per aggiungere un ulteriore elemento (per esempio una curva) ad un graco gi` esistente, battere, prima di dare il secondo comando di plot, a >>hold on (per disattivare, >>hold off). Per stabilire esplicitamente lintervallo di valori da plottare sugli assi, si usa >>axis([xmin xmax ymin ymax]) dove ovviamente xmin ecc.. sono da sostituire coi valori desiderati. Se per uno di questi estremi si desidera un valore scelto automaticamente, lo si mette uguale a -Inf (meno innito) se ` un limite inferiore, +Inf se ` un e e limite superiore. Queste manipolazioni degli assi consentono di zoomare solo su parte della serie, esattamente come si farebbe limitando il comando plot a una parte del vettore pd o con lo strumento zoom. Anche laspect ratio degli assi si pu` scegliere: ad esempio o >>axis equal produce intervalli di uguale ampiezza sugli assi.


19

Quanto ai tick, essi vengono automaticamente equispaziati, sebbene sia possibile controllarne la posizione con il menu Edit/CurrentObjectProperties nella nestra della gura. Si apre un Property Editor che consente di variare a piacere le propriet` delloggetto graco (sia esso gura, linea, punto, a asse, ecc...) su cui si ` cliccato per ultimo 14 . e Si pu` intervenire sui tick anche da linea comandi: ad esempio, dopo aver o dato un comando di plot (gura n.1), digitando >>h=gca(1) >>get(h) >>set(h,Ytick,[0 40 80 120 160 200]) >>set(h,YMinorTick,on) si metteranno i major ticks nelle posizioni indicate in parentesi quadre e si abiliter` il tracciamento dei minor ticks. a Per avere graci in scale diverse da quella lineare-lineare, si usano, al posto di plot, le funzioni loglog, semilogx, semilogy, con ovvio signicato. Citiamo inne il fatto che si possono fare pi` graci in una sola gura u usando listruzione subplot, creare graci con doppio asse delle ordinate usando yyplot, ecc... 14. Salvataggio di graci: una volta creato un graco, si pu` salvarlo in un o formato graco a scelta. Si esporta il graco in un le, ad es. postscript, (sia figura.ps), col comando >>print -dps2 figura.ps (b/n, postscript livello 2), oppure >>print -dpsc2 figura.ps (colore, postscript livello 2). Naturalmente sono possibili altri formati di output, quali EPS (encapsulated postscript ), TIFF, JPEG ecc... Supponiamo ora che si voglia intervenire successivamente su un certo graco che ` stato creato. Distinguiamo fra due diverse esigenze: e (a) in una stessa sessione di lavoro si desidera avere pi` nestre-gura, u richiamabili e modicabili. Ci` ` possibile col comando figure: si o e creano in successione diverse nestre-gura, ciascuna caratterizzata da un proprio handle (un numero intero), e si plottano dierenti set di dati nelle singole nestre-gura. Per esempio: >>...... >>h1=figure crea una gura vuota, caratterizzata da un handle che viene memorizzato in h1; >>...... >>plot...... riempie la gura con un graco; >>...... >>h2=figure14 Ogni singolo oggetto graco riceve un identicatore univoco (handle) che pu` essere reo cuperato (cio` si pu` venire a conoscenza del suo valore) e poi usato per cambiare le propriet` e o a delloggetto stesso con opportuni comandi. Questa ` la base del sistema Handle graphics di e MATLAB.

20

CAPITOLO 1. ESERCITAZIONI GUIDATE crea una seconda gura vuota con un handle diverso, che viene memorizzato in h2; >>...... >>plot...... riempie la seconda gura, ecc... Ogni volta si crea una gura diversa ed il suo handle viene memorizzato. Per usare nuovamente una gura gi` creata (per aggiungervi un nuovo plot, o semplicemente rivederla a o modicarla in qualche modo) si clicca sulla gura scelta o si digita >>figure(h) dove h ` lo handle della gura voluta. Si procede poi alle modiche e desiderate. Le gure cos` create si posso chiudere, anche digitando da linea comandi >>close(h) ecc... Lo handle serve anche per modicare le propriet` della gura (o di a qualsiasi altro oggetto graco) da linea comandi, come gi` accennato: a ad es., per una gura con un handle h, digitando >>get(h) si ha una lista delle propriet`, che si possono cambiare con comandi a del tipo >>set(h,PropertyName,PropertyValue) Si pu` scegliere cos` ad esempio, il valore A3 per il formato della o , carta (mentre il default ` A4): e >>set(h,PaperType,A3) Inoltre si possono fare graci con orientamento landscape anzich e portrait, ecc... (b) Si vuole trasmettere un graco da una sessione di lavoro ad unal` tra. E in eetti possibile salvare la gura corrente (quella su cui si ` cliccato per ultima) per poi rigenerarla quando lo si desidera; e allo scopo si fa s` che MATLAB crei un M-le con i comandi utili a ricreare la gura (che pu` poi essere rilanciato a piacere) e, se neo cessario, anche un le .mat avente lo stesso nome del precedente e contenente i dati che servono per ricreare la gura medesima. Ad esempio, digitando >>print -dmfile mio_grafico dopo aver creato una gura, si otterranno un mio_grafico.m e un mio_grafico.mat; in ogni momento, poi, digitando >>mio_grafico si ricreer` la gura. a

1.1.4

Generazione di dati sintetici

Spesso, per studiare le caratteristiche e le prestazioni di algoritmi e sistemi, si eseguono simulazioni numeriche mediante segnali costruiti (sintetizzati) appositamente. Nel caso pi` comune, essi vengono composti con un numero arbitrario u di sinusoidi con certe ampiezze, frequenze e fasi iniziali. Se lo scopo ` quello e

1.1. UN PRIMO APPROCCIO A MATLAB Frequenza f0 f1 f2 f3 f4 f5 Valore 6.0 102 3.5 102 1.5 102 8.0 103 4.0 103 2.0 103

21

Hz Hz Hz Hz Hz Hz

Tabella 1.1: Frequenze analogiche delle sei sinusoidi che andranno a formare i dati sintetici. di simulare dati random da misure siche, ` bene che la fase iniziale di ogni e componente sia scelta a caso tra 0 e 2. Per la stessa ragione, alle sinusoidi si aggiunge rumore bianco gaussiano 15 con media e varianza opportuna. Nonostante nelle tecniche di elaborazione di dati a tempo discreto lintervallo di campionamento dei dati non giochi alcun ruolo, qui e nel seguito si desidera mantenere costantemente un aggancio con un ipotetico contesto sico, in cui si suppone sia esistito un segnale a tempo continuo da cui i dati a tempo discreto deriverebbero, per campionamento a passo Tc ; quindi assegneremo anche ai dati sintetici che ora genereremo un valore di Tc da noi arbitrariamente scelto. Per ssare le idee, decidiamo di generare un segnale contenente N = 8192 dati con Tc = 5 s, includendovi: sei sinusoidi, tutte di ampiezza unitaria, fase iniziale random e frequenze date nella Tabella 1.1. rumore bianco a media nulla e varianza 0.25 (r = 0.5). A tal ne creeremo un M-le, sia16 dat_sint.m, che contenga le seguenti istruzioni:15 Campioni di rumore bianco gaussiano possono venire generati partendo da un generatore di numeri pseudo-casuali e modicandone luscita per avere la gaussianit`. a Supponiamo di scegliere una serie di numeri reali secondo le regole seguenti:

tutti i numeri devono essere compresi tra zero e uno; non si possono fare discriminazioni in favore di (o contro) alcun numero o intervallo di numeri; tutti devono avere la stessa probabilit` di essere scelti; a ogni numero deve essere scelto indipendentemente dagli altri; non si debbono ricordare i numeri gi` scelti e cambiare le proprie decisioni in base ad essi. a Cos` facendo si genera rumore bianco con distribuzione uniforme. Diverso ` il caso gaussiano. Supponiamo di captare, con unantenna, energia da tutte le e sorgenti di segnali radio dellUniverso, sommando tali energie. La tensione risultante sar` a casuale, ma non uniforme: la maggior parte del tempo questi segnali tenderanno a cancellarsi lun laltro e genereranno tensioni prossime allo zero; una volta ogni tanto invece interferiranno costruttivamente e genereranno termini elevati. La distribuzione delle tensioni sar` gaussiana; a per generare rumore bianco gaussiano, le regole di scelta sono le seguenti: i numeri devono essere estratti da una popolazione con distribuzione gaussiana; devono essere scelti indipendentemente luno dallaltro.16 Non

si pu` assegnare lo stesso nome ad un M-le e ad una variabile. o

22

CAPITOLO 1. ESERCITAZIONI GUIDATE

N=8192 T_c=5 Nsin=6 A=[1 1 1 1 1 1] f_sint=[6e-2 3.5e-2 1.5e-2 8e-3 4e-3 2e-3] t_ds=0:T_c:(N-1)*T_c; sigmar=0.5 ds=zeros(1,N); rand(state,0) (porta il generatore rand nel suo stato iniziale) for i=1:Nsin pharand=rand(1,1)*2*pi; ds=ds+A(i)*sin(pharand+2*pi*f_sint(i)*t_ds); end randn(state,0) (porta il generatore randn nel suo stato iniziale) rumb=sigmar*randn(1,N); ds=ds+rumb; Listruzione for produce un loop sullindice i, essendo, come gi` accennato, a la versione di MATLAB di un do loop del FORTRAN.

Generati cos` i dati, sar` utile eseguire un graco complessivo ed uno espanso, a in una scala opportuna per vedere bene le uttuazioni. Si noti che con Tc = 5 s, le sinusoidi hanno periodi di 16.7, 28.6, 66.7, 125., 250. e 500. s rispettivamente. Riessione: rumore e segnale utile possono prevalere luno sullaltro, a seconda delle situazioni; pi` avanti, nelle esercitazioni sullanalisi spettrale, si parler` u a di rapporto segnale-rumore. Per ora potrebbe essere istruttivo tenere sse le ampiezze delle sinusoidi e aumentare (diminuire) la varianza del rumore, osservando no a quando ` possibile riconoscere le sinusoidi (il rumore). e

Seguendo la stessa traccia, lo studente prepari ora altri tre programmi, che serviranno per le prove successive: 1. dat_rumb.m che generi N campioni di rumore bianco gaussiano a media nulla e con varianza r ; 2. dat_sint_2sin.m che generi N campioni di rumore bianco gaussiano a media nulla e con varianza r sovrapposto una sinusoide di ampiezza A, e frequenza ssata, con fase iniziale casuale; 3. dat_sint_2sin.m che generi N campioni di rumore bianco gaussiano a media nulla e con varianza r sovrapposto a due sinusoidi di ampiezze A0 , A1 e frequenze diverse, con fasi iniziali casuali. I valori dei parametri r , A ecc... possono per ora essere ssati arbitrariamente; in seguito editeremo questi prototipi di le per inserire i valori opportuni, a

1.2. TRASFORMATA DISCRETA DI FOURIER seconda della necessit`. a

23

Concluso il lavoro, usciamo da MATLAB, con quit o selezionando Exit MATLAB dal menu File. Se al termine di una sessione si desiderasse salvare delle variabili, useremmo save e al login seguente potremmo ricaricare le variabili nel workspace con load. Nel caso di questa sessione tuttavia non serve salvare nulla.

1.21.2.1

Trasformata discreta di FourierPremessa

Con questa sezione delle esercitazioni guidate si entra nel campo specico delle applicazioni di elaborazione numerica dei segnali, che in MATLAB corrispondono a funzioni del Signal Processing (SP) Toolbox 17 . Lo scopo dellesercitazione ` il calcolo della trasformata di Fourier discreta di una sequenza di dati, una e quantit` complessa di cui si esaminer` il modulo (spettro di ampiezza). Per le a a prove pratiche useremo i dati sintetici ds e i dati veri pd. Prima di proseguire ` bene chiarire un punto concettualmente importante. La e nozione di spettro di ampiezza si applica ai segnali deterministici, mentre per i segnali casuali si adotta un approccio statistico, che porta alla nozione di spettro di potenza del processo casuale (innitamente persistente e di energia innita) sotteso alla generazione dei dati stessi. Per tale processo casuale non ` denito e alcun spettro di ampiezza, perch la DFT non esiste; esiste invece la DFT della e sequenza di autocovarianza del processo, cio`, appunto, lo spettro di potenza. e I dati utilizzati nelle presenti esercitazioni sono da un lato misure vere, da un altro lato dati sintetici, randomizzati con linserimento di rumore bianco per farli assomigliare il pi` possibile a dati veri (che sono quelli che ci interessano). A u tutti gli eetti, la loro rappresentazione in frequenza va dunque data in termini di spettro di potenza. Tuttavia, a ni esercitativi pratici, per poter ripartire adeguatamente nel tempo le esercitazioni pratiche in rapporto alla trattazione teorica, parleremo inizialmente di spettro di ampiezza anche per dati di questo tipo; semplicemente, lo faremo consapevoli del fatto che allora stiamo guardando i dati in questione non come segmenti di durata nita di realizzazioni di processi casuali ma come singole sequenze deterministiche di durata nita. Al termine dellesercitazione di questa sezione si eseguono anche alcune prove per visualizzare, nel dominio della frequenza ed in quello del tempo, il fenomeno dellaliasing. Allo scopo si generano sequenze di dati sintetici discreti (sempre pensati come derivanti per campionamento, a passo Tc arbitrariamente scelto,17 In realt` lalgoritmo FFT ` contenuto nel Toolbox di base, ma altre function pi` a e u specialistiche sono nel Toolbox SP.

24


da un segnale a tempo continuo) 18 , contenenti ciascuna una sinusoide la cui frequenza in Hz sia esterna allintervallo (0, fN y ). Successivamente viene stimato il relativo spettro di ampiezza; in questo modo ` possibile constatare che le e frequenze esterne allintervallo (0, fN y ) vengono tutte riesse al suo interno.

1.2.2

DFT e spettro di ampiezza: caratteristiche

1. La DFT X[k] di una sequenza, x[n] (n = 0, N 1), viene calcolata mediante la function che implementa lalgoritmo FFT. Per lecienza dellalgoritmo, ` bene lavorare su un numero di dati uguale ad una potenza e di due, sia M = 2 . Pertanto la sequenza di dati deve: (a) essere presa com `, se M = N ; e (b) essere allungata, aggiungendo al fondo degli zeri (padding di zeri), se si decide per M > N . Leventuale allungamento con zeri si fa sempre dopo aver ridotto i dati disponibili a media nulla, sottraendo da essi il valore medio. Questa operazione di centratura dei dati si fa sempre e comunque prima di eseguire la trasformata, anche nel caso non sia necessario alcun padding, al ne di evitare che la media (trasformata a k = 0) inuenzi i valori spettrali a bassa frequenza. La trasformata X[k] ` complessa in generale e di essa prendiamo il modulo e |X[k]|, k = 0, M 1, comprendente M valori. Per conoscere le frequenze a cui tali valori corrispondono (siano fk [Hz]) si deve tenere conto delle frequenze angolari discrete che corrispondono alle stime spettrali, k = 2k , M

dove M ` la lunghezza della DFT, e del periodo di campionamento Tc dei e dati. Infatti la relazione tra k e fk ` e fk = da cui fk = k . M Tc k 2Tc

e Lo spettro ` dunque esteso su frequenze discrete da 0 a M1 , cio` da 0 e MTc 1 alla frequenza analogica che precede la frequenza di campionamento Tc . In frequenze angolari discrete si va da 0 incluso a 2 escluso.qui ancora una volta che nelle operazioni descritte si ` voluto mantenere il e riferimento al campionamento e al relativo passo temporale Tc , cio` alla variabile frequenza f e espressa in Hz, anche se, trattandosi di segnali a tempo discreto, sarebbe stato possibile luso della sola variabile (oppure ), che ` del tutto svincolata dallidea di campionamento e dal e Tc .18 Ricordiamo

1.2. TRASFORMATA DISCRETA DI FOURIER

25

Il controllo del calcolo di |X[k]| va fatto vericando che sia soddisfatta la relazione di Parseval per la DFT: 1 MM1 N

k=1

|X[k]| =

2

x2 [n],n=0

dove N ` la lunghezza della sequenza utilizzata. Riassumendo, si dovr`: e a (a) costruire il vettore di dati a media nulla da trasformare, prendendo tal quale e centrando oppure centrando e allungando con zeri la sequenza x[n]; (b) fare il calcolo del secondo membro della relazione di Parseval (quello che riguarda i dati e ne rappresenta lenergia); (c) eseguire la DFT; (d) prenderne il modulo; (e) con questultimo calcolare il valore del primo membro della relazione di Parseval ed accertarsi che sia uguale al secondo (salvo errori di arrotondamento); (f) calcolare le frequenze in Hz corrispondenti alle stime eettuate; (g) gracare il risultato. 2. La X[k] ` simmetrica, essendo la trasformata di una sequenza reale; di e conseguenza anche lo spettro di ampiezza lo ` e, se M ` pari come di solito e e avviene, la seconda parte di esso (da f M +1 a fM1 ) ripete la prima parte 2 (da f1 a f M 1 ) in modo speculare:2

|X[k]| = |X [M k]|

k = 0, M 1.

Alla frequenza zero (k = 0) e alla frequenza di Nyquist (k = M/2) la trasformata ` reale; il valore a k = 0 rappresenta la media dei dati, che e sono stati centrati, e quindi deve risultare piccolissimo (idealmente, nullo). Pertanto linformazione utile ` contenuta nella prima parte dello spettro, e 1 1 da k = 1 a k = M , cio` da f1 = MTc a f M = fN y = 2Tc , e quindi per dati e 2 2 reali centrati e M pari ` suciente gracare lo spettro su questultimo e intervallo 19 . Il graco spettrale si fa normalmente in scala log-log, ma anche scale lineari o semilineari sono talvolta utili (solo se lasse f ` lineare e si include la frequenza zero).

1.2.3

La function t di MATLAB

Attenzione: dora in poi si consideri che nella teoria si fa uso di indici che partono dal valore 0 mentre in MATLAB gli indici partono da 1.+1 1 M dispari la situazione ` simile, per` fNy cade tra due campioni (tra M2 e M2 ) e o M 1 M +1 e i campioni da 2 a M 1 replicano in modo speculare quelli da 1 a 2 . Dora in poi assumeremo che la lunghezza della DFT sia pari; se non lo fosse, si potr` sempre aggiungere a uno zero al fondo della sequenza. 19 Per

26


Per la DFT dei dati via Fast Fourier Transform si usa la function fft; essa pu` accettare in input un numero qualsiasi di dati, ma lalgoritmo veloce viene o eettivamente applicato solo se la lunghezza della trasformata da calcolare (M = nfft) ` uguale ad una potenza di due; in caso contrario viene fatto il calcolo e diretto della DFT, assai pi` lento. Per la scelta di nfft pu` essere daiuto la u o function nextpow2, che, dato un intero, calcola lesponente della pi` vicina pou tenza di due, per eccesso: ad esempio, nel caso dei dati ds, il digitare >>N=8192 >>enfft=nextpow2(N) >>nfft=2^enfft d` semplice conferma del fatto che era gi` stato scelto un N uguale ad una a a potenza di due. Invece per N = 25921, che ` la lunghezza dei dati pd, si ottiene e il valore 215 = 32768. Per usare fft, ad esempio sui dati ds generati da dat_sint, se si vuole una FFT di lunghezza uguale alla lunghezza dei dati si scrive 20 DS=fft(ds); Se la si vuole di lunghezza diversa, specicata dal valore di nfft, si scrive invece DS=fft(ds,nfft); Ove necessario, il padding di zeri viene fatto automaticamente.

1.2.4

Prove con dati sintetici

Iniziamo coi i dati sintetici ds (N = 8192), scegliendo, per la lunghezza della FFT, M = N = 8192 = 213 . Si riduce a media nulla con ds=ds-mean(ds); Si calcola poi la DFT con nfft=N; DS=fft(ds,nfft); e di essa si prende il modulo con modDS=abs(DS); Questo ` lo spettro di ampiezza. In seguito si graca il risultato in scala semie logaritmica su tutto lintervallo di frequenze con i comandi: f_Ny=1/(2*T_c) deltf=1/(length(DS)*T_c) f=0:deltf:(2*f_Ny-deltf); semilogy(f,modDS) Si noti la simmetria del graco. Poi si pu` zoomare sui particolari dello o spettro, o limitando lintervallo di indici delle variabili, o dando specici limiti per i valori sugli assi, o usando gli strumenti della nestra graca.20 Preparando uno script sia spe amp.m con le istruzioni date in questa sottosezione per il calcolo e la visualizzazione dello spettro di ampiezza dei dati ds, si verr` a disporre di un a prototipo modicabile, utile in seguito. Questo ` conveniente anche per tutte le prove che e seguono, e non verr` ripetuto. a


27

Per un graco in scala log-log, da eseguire in unaltra gura, si scrive invece loglog(f,modDS) (la frequenza zero viene automaticamente esclusa), mentre per un graco linearelineare, da fare ancora in unaltra gura 21 , si scrive plot(f,modDS) Inne notiamo che talvolta si usa anche rappresentare lo spettro in dB. In tal caso il comando ` e plot(f,10*log10(modDS))

Riessione: come mai i graci in scala semilogaritmica (fatti con semilogy o in bB) appaiono cos` schiacciati? Come mai questo non succede col graco fatto con loglog? La posizione dei picchi spettrali, i quali sono superiori a 1. 103 , pu` essere o trovata coi comandi i=find(modDS>10^3); format long e fi=f(i) che danno le frequenze dei picchi che valicano il limite specicato, in forma ` di vettore-colonna facilmente leggibile. E cos` possibile vericare la simmetria attesa per la DFT e il fatto che i picchi cadono (approssimativamente) alle frequenze scelte per generare le sinusoidi 22 .

Se si desidera vedere anche lo spettro di fase, si opera cos` per calcolare la fase : angleDS=angle(DS); poi per dispiegarla, cio` per rimuovere i salti di fase maggiori di , riducendoli e al loro complemento a 2 in modo da togliere le discontinuit` ttizie introdotte a dalla funzione arcotangente usata per il calcolo della fase, uangleDS=unwrap(angleDS); e inne per gracarla in scala lineare-lineare (in una gura a parte) plot(f,uangleDS)

Ricordiamo che talvolta pu` essere utile ribaltare la parte della trasformata che o sta oltre fN y sulle frequenze negative, tra 0 e fN y . In questo modo i nuovi intervalli sono: in frequenze analogiche, da fN y inclusa a +fN y esclusa; in numeri darmonica, M k M 1. 2 2 La procedura ` la seguente: e DSshift=fftshift(DS); modDSshift=abs(DSshift); Per il graco di modDSshift occorre generare un vettore frequenze apposito: fnegpos=(-f_Ny:deltf:f_Ny-deltf);21 Dora in avanti ` considerato implicito che sequenze diverse vadano gracate in gure e diverse se non condividono almeno in parte i valori della variabile indipendente. 22 Alcuni picchi appariranno due volte nellelenco, a causa della loro larghezza nita.

28


plot(fnegpos,modDSshift) (oppure semilogy; non loglog, a causa della presenza di frequenze negative). Si ricordi di aprire una nestra separata anche questa volta.

Si osservi di nuovo la simmetria del graco; per rendersi conto quantitativamente dellaccaduto si potrebbe anche eettuare la ricerca delle frequenze di massimo pure per modDSshift e confrontarle con quelle per modDS.

Osservazioni Guardando il graco in scala log-log, si nota bene uno zoccolo di valori spettrali rapidamente variabili, specialmente evidente alle frequenze pi` u alte. Tale caratteristica ` dovuta al rumore. Lo si pu` constatare provando e o a generare un set di dati sintetici uguale al precedente ma senza rumore (con le sole sinusoidi) e vedendo che cosa cambia nello spettro. Riessione: ` la linearit` della DFT che ci consente di estendere il modello e a del segnale (segnale utile+rumore) dal dominio del tempo al dominio della frequenza. I picchi spettrali corrispondenti alle sinusoidi discrete introdotte nei dati sintetici non appaiono innitamente sottili, ma hanno piuttosto una certa larghezza nita (e infatti quando ne abbiamo determinato le posizioni, alcuni picchi sono apparsi due volte nella lista). Il fatto non sorprende, visto quanto appreso riguardo agli eetti del windowing e al campionamento spettrale e; inoltre nei dati vi ` anche del rumore. Leetto del e campionamento in frequenza pu` essere studiato generando nuovi dati o sintetici contenenti sinusoidi di frequenze esattamente uguali ad alcune delle frequenze discrete a cui lo spettro viene campionato e osservando come laspetto dello spettro si modica. Riessione: nel primo caso, con frequenze scelte come in Tabella 1.1, la sequenza generata ` periodica? Nel secondo caso ` periodica? Perch? e e e

1.2.5

Prova coi dati di pressione

Eseguiamo ora una prova coi dati pd, la cui lunghezza ` N=25921 e il cui Tc ` e e di 5 s (lo stesso valore scelto arbitrariamente per i dati sintetici), facendo una fft su un numero di punti pari alla potenza di due immediatamente superiore a N . Prima di procedere sottrarremo ai dati la media. Calcoleremo lo spettro dampiezza e faremo il graco, ad esempio con loglog. Osservando lo spettro ottenuto, come semilogy o loglog, lo studente vedr` che a esso ha una forma (tipica dei dati geosici: atmosferici, marini, sismici....) che assegna molta potenza alle basse frequenze e poi sempre meno andando verso quelle alte. Come sappiamo, un segnale con uno spettro di questo tipo ` detto e rumore di colore caldo (ad esempio rumore rosso). Esso ` tipico dellambiente e


29

atmosferico e oceanico e la sua forma denuncia la memoria del sistema, ossia la persistenza del fenomeno sico sotteso alla generazione dei dati (inerzia del sistema o processo che prosegue nel tempo). In altre parole, il rumore ` il e risultato dellinterazione tra una forzante di rumore bianco e i componenti, a risposta lenta, del sistema sico esaminato; ad es., in una serie climatica di temperature oceaniche, la persistenza ` dovuta al fatto che linerzia termica e delloceano integra nel tempo il forcing dovuto al tempo atmosferico (che qui viene visto come rumore). Sovrapposti a uno spettro di questo tipo possono poi esserci, nei dati geosici, picchi associati a speciche periodicit`. a

1.2.6

Aliasing

Lo scopo dellesercitazione ` la visualizzazione del fenomeno dellaliasing. Ci` si e o ottiene generando dati sintetici discreti (pensati come derivanti, tramite campionamento con un certo passo Tc arbitrario, da un segnale continuo), contenenti sinusoidi di frequenza esterna allintervallo (0, fN y ). Sia f1 una frequenza in (0, fN y ). Si considerino poi le innite frequenze f2 = k f1 Tc

con k = (0, ) intero, che invece sono esterne allintervallo in questione. Si pensi ora ad una sequenza contenente una sinusoide con una delle frequenze f2 , campionata a passo Tc : calcolando il relativo spettro di ampiezza ` possibile e constatare che f2 viene riessa in f1 , nel senso che nello spettro il picco compare a frequenza f1 .

Procedimento: 1. Generare quattro diverse sequenze sintetiche, tutte arbitrariamente riferite a Tc = 5 s, fN y = 101 Hz, contenenti ciascuna: del rumore bianco gaussiano, con valor medio nullo e deviazione standard pari a 0.05, e una sinusoide con frequenza f2 tale da riettersi, secondo la formula data sopra, in f1 = 4. 102 Hz. Ad esempio, applicare la formula per f2 , due volte con k = 1 e due volte con k = 2, una volta col segno meno ed una col segno pi` nei due casi, u ottenendo quattro valori di f2 da usare per costruire quattro sequenze diverse; 2. calcolare lo spettro di ampiezza 23 e vericare la riessione nei quattro casi, con graci in scala lineare (o anche semilogy, loglog oppure in dB). Poich la simmetria della DFT completa ` gi` stata osservata nei graci e e a23 Attenzione: da questo momento in poi non verr` pi` ricordata ogni volta lopportunit` a u a del togliere la media ai dati prima di calcolarne lo spettro!

30

CAPITOLO 1. ESERCITAZIONI GUIDATE fatti in precedenza, dora in poi ci si pu` limitare a riportare nei graci o spettrali solo i valori da k = 0 a k = M ; in termini di indici in Matlab, ci` o 2 signica gracare gli elementi del vettore spettro di ampiezza dallindice 1 allindice M + 1 = nfft/2+1. Anzi, se usiamo semilogy o la scala di dB, 2 visto che prima di eettuare la trasformata abbiamo tolto la media ai dati e quindi il modulo della trasformata alla prima frequenza, la frequenza zero, ` piccolissimo, ` meglio partire dallindice 2: in questo modo si evita che la e e parte interessante del graco appaia schiacciata. Il problema ovviamente non si presenta in scala lineare e nemmeno con loglog, perch con esso e la frequenza zero viene automaticamente esclusa; 3. vericare laliasing anche mediante graci temporali, in una scala di tempo adatta a visualizzare qualche ciclo della sinusoide aliased ; 4. generare una quinta sequenza derivante da un segnale continuo sinusoidale di frequenza pari alla quarta precedente (f2 = 4.4 101 Hz), ma con un passo di campionamento che sia dieci volte pi` piccolo di quello usato u prima, cio` 0.5 s. Ovviamente cambia, rispetto ai casi precedenti, la free quenza di Nyquist (come diventa?). Perci` una sinusoide di tale frequenza o pu` ora essere rivelata correttamente, perch inferiore alla nuova fN y (` o e e vero?...). Vericare questo con un graco temporale ed uno spettrale.

1.31.3.1

Filtraggio nel dominio del tempoPremessa

Lo scopo dellesercitazione ` eettuare, nel dominio del tempo, il ltraggio di una e serie temporale di dati con un ltro FIR di cui si conosce la risposta allimpulso h[n]. Verranno impiegate, nei prossimi due paragra, tre funzioni tra quelle disponibili in MATLAB per operazioni in questo campo. La prima ` filter, che pu` essere e o vista come una implementazione eciente dellequazione alle dierenze; quindi essa pu` servire per applicare sia ltri IIR, sia ltri FIR di cui si conoscano i o coecienti costanti dellequazione alle dierenze. Nel seguito chiameremo di tipo a i coecienti del primo membro dellequazione alle dierenze (quelli che moltiplicano la sequenza di output, ritardata o no) e di tipo b i coecienti del secondo membro (quelli che moltiplicano la sequenza di input, ritardata o no). In particolare, nel caso FIR, quello cui noi ci limiteremo, i coecienti di tipo b sono i campioni della risposta allimpulso del ltro h[n] mentre quelli di tipo a sono rappresentati da un unico valore unitario, quello che moltiplica y[n]: gli altri sono nulli. La seconda function utilizzata ` conv, che invece implementa la convoluzione e lineare della serie di dati con h[n].

1.3. FILTRAGGIO NEL DOMINIO DEL TEMPO

31

Inne, la terza function usata ` fftfilt, che implementa un eciente algoritmo e di ltraggio nel dominio della frequenza, basato sulla FFT e chiamato overlapadd. Sar` fatto anche un cenno ad un metodo di ltraggio caratterizzato da ritardo a nullo. Le operazioni di ltraggio verranno esemplicate utilizzando, come dati da ltrare, quelli sintetici ds con rumore e sei sinusoidi. Per quanto concerne il ltro, sono disponibili le risposte allimpulso di tre diversi ltri FIR passa-banda, tutte di lunghezza Mf = 55 (ordine Mf 1 = 54); i nomi dei les sono bpf_h.mat, bpf_m.mat+,bpf_l.mat+;le corrispondenti variabili si chiamano hh, hm, hl. Nelle Fig. 1.1, 1.2 e 1.3 sono mostrati i moduli delle rispettive risposte in frequenza, in unit` logaritmiche (moduli quadri) e lineari. a10 0 10 20 30 40 50 60 70 80 90 100 110 120 130 0 0.1 0.2 1.1 1 0.9 0.8 0.7

20 log10 |H |

|H |0.3 0.4 0.5

0.6 0.5 0.4 0.3 0.2 0.1

0 0

0.1

0.2

0.3

0.4

0.5

Figura 1.1: Filtro bpf h Risposta in frequenza: modulo quadro in dB (a sinistra) e modulo in scala lineare (a destra), in funzione della frequenza adimensionale . Le linee verdi indicano lattenuazione minima ed il ripple in banda passante. = 27 campioni: il senso La fase ` lineare e comporta un ritardo di f2 e di questo ritardo ` che se si devono gracare insieme la serie originale e quella e ltrata, in funzione del tempo, la seconda deve iniziare con un ritardo di 27 campioni rispetto al primo campione della serie originale. La serie di output, senza transitori, conta Nf = N Mf + 1 campioni, per un ingresso composto da N campioni 24 . Le bande passanti dei tre ltri passa-banda, espresse in termini di frequenze adimensionali , sono indicate in Tabella 1.2.24 Siano x[l ], con l = (0, N 1), i dati di input. Lindice l e i successivi indici con lo 0 0 0 stesso pedice hanno valori che partono da zero.

M 1

32


10 0 10 20 30 40 50 60 70 80 90 100 110 120 130 0 0.1 0.2 1.1 1 0.9 0.8 0.7

20 log10 |H |

|H |0.3 0.4 0.5

0.6 0.5 0.4 0.3 0.2 0.1

0 0

0.1

0.2

0.3

0.4

0.5

Figura 1.2: Filtro bpf m Risposta in frequenza: modulo quadro in dB (a sinistra) e modulo in scala lineare (a destra), in funzione della frequenza adimensionale .10 0 10 20 30 40 50 60 70 80 90 100 110 120 130 0 0.1 0.2 1.1 1 0.9 0.8 0.7

20 log10 |H |

|H |0.3 0.4 0.5

0.6 0.5 0.4 0.3 0.2 0.1

0 0

0.1

0.2

0.3

0.4

0.5

Figura 1.3: Filtro bpf m Risposta in frequenza: modulo quadro in dB (a sinistra) e modulo in scala lineare (a destra), in funzione della frequenza adimensionale .

Se si scartano i transitori, la convoluzione con cui il ltro viene applicato,Mf 1

y[j0 ] =k0 =0

h[k0 ]x[j0 k0 ],

deve avere gli indici delloutput y[j0 ] compresi tra j0 = Mf 1 e j0 = N 1, in modo da avere j0 k0 0; j0 k0 N 1.

1.3. FILTRAGGIO NEL DOMINIO DEL TEMPO Filtro bpf l bpf m bpf h Intervallo passante di 0.05 0.1075 0.1075 0.232 0.232 0.4

33

Tabella 1.2: Bande passanti dei ltri citati nel testo, espresse in termini di frequenza adimensionale . Ricordando che per i dati ds si ` assunto un passo di campionamento di 5 s, e ` immediato calcolare le frequenze in Hz che limitano le bande passanti dei e tre ltri, quando applicati a dati con quellintervallo di campionamento. Quali frequenze, fra quelle contenute nei dati ds (Tabella 1.1), cadono nella banda passante del primo ltro? ...del secondo? ...del terzo? Lo studente, rispondendo a queste domande, si far` unidea dellaspetto che dovr` avere lo spettro dei a a dati ds ltrati con ciascuno dei tre ltri; poi ltrer` eettivamente, ad esempio aPer riferirsi, come in MATLAB, ad indici che corrono da 1 anzich` da 0, linput deve per` e o essere riscritto come x[l], con l = (1, N ) e la convoluzione deve essere riscritta comeMf

y[j] =k=1

h[k]x[j k + 1].

Escludendo i transitori, lindice delloutput varier` da Mf a N . Per riportare allora lindice a delloutput allintervallo partente da 1, si deve passare a un indice i = j Mf + 1, col quale si scriver` la convoluzione come aMf

y[i] =k=1

h[k]x[j k + 1],

dove i = 1, N Mf + 1 senza i transitori. Ma come collocheremo temporalmente il primo campione delloutput rispetto a quello dellinput? Consideriamo che un dato y[i] ` composto da una combinazione lineare di valori dellingresso e x[l] aventi indici l compresi fra l1 = j k + 1 = (i + Mf 1) Mf + 1 = i per k = Mf e l2 = j k + 1 = (i + Mf 1) 1 + 1 = i + Mf 1 per k = 1. Da ci` si ricava o l2 l1 t pertanto si assegna y[i] al tempo ti = tl + dove tl = (l 1)Tc ` il tempo di x[l]. e campioni, un Si ha dunque, come previsto, un ritardo in tempo reale corrispondente a f 2 numero che nel nostro caso ` intero perch Mf 1 = 54 ` pari; rispetto allinput, loutput e e e della convoluzione priva di transitori presenta in meno in coda.Mf 1 2 M 1

= =

Mf 1, (l2 l1 )Tc = (Mf 1)Tc ;

Mf 1 t = tl + Tc , 2 2

campioni in meno in testa e altrettanti

34


con bpf_m, e vericher` lavvenuto ltraggio mediante lesame dello spettro di a ampiezza dei dati ltrati e il suo confronto con quello dei dati ds originali, calcolato nellesercitazione della sottosezione 1.2.4. Inne, nellultima parte della presente sezione, si ltra nuovamente via convoluzione lineare, ma prima di ltrare si adotta una tecnica di decimazione dei dati, al ne di far corrispondere la banda passante del ltro usato, originariamente data in termini di frequenze adimensionali (Tabella 1.2), con un intervallo di frequenze in Hz diverso da quello su cui si cadrebbe mantenendo il passo di campionamento originale: si ricordi infatti che f = /Tc . La decimazione sposta la frequenza di Nyquist, riducendola di un fattore pari a quello di decimazione. Pertanto prima di decimare si eettua il necessario pre-ltraggio passabasso anti-aliasing, atto ad eliminare dai dati tutte le componenti di frequenza superiore alla nuova frequenza di Nyquist.

1.3.2

Filtraggio via equazione alle dierenze

Supponiamo di voler ltrare una sequenza col ltro bpf_m. Il ltraggio si si esegue fornendo in input alla function filter i coecienti hm, lunico coeciente di tipo a uguale a 1 ed i dati da ltrare: dsfed=filter(hm,1,ds); Il nome suggerito per il vettore dei dati ltrati richiama il fatto che provengono da dati ds e sono filtrati via equazione alle dierenze. Il vettore di dati ltrati, come si verica con length, ` lungo come linput, cio` ha una lunghezza e e di 8192 campioni. Per comprendere che cosa questo comporti, conviene innanzi tutto gracare i dati ltrati e specialmente la testa della sequenza: plot(dsfed(1:100)) Si nota chiaramente che ` presente un transitorio iniziale e questo ` in accordo e e col fatto che filter calcola loutput partendo dal primo passo temporale e assumendo condizioni iniziali nulle: i valori passati, cio` antecedenti a n = 0, e della sequenza di ingresso al ltro FIR sono presi uguali a zero; pertanto un numero di campioni delloutput pari allordine del ltro ` aetto dal transitorio e iniziale. Al fondo della sequenza ltrata non si nota invece niente che abbia laspetto di un transitorio: lo studente verichi questo con un graco degli ultimi 100 punti (indici end-99:end). Ci` signica che filter si ferma allultimo o campione della sequenza di ingresso. La percezione dellavvenuto ltraggio si ha calcolando e gracando lo spettro di ampiezza della sequenza di dati ltrati e confrontandolo con quello dei dati ds originali. Sceglieremo sempre nfft=N; calcoleremo come fatto nella sottosezione 1.2.4 lo spettro dei dati ds; poi faremo lo stesso con dsfed: DSFED=fft(dsfed,nfft); modDSFED=abs(DSFED);

1.3. FILTRAGGIO NEL DOMINIO DEL TEMPO loglog(f,modDSFED)

35

Il vettore frequenze ` lo stesso utilizzato per plottare modDS perch le due fft, e e di ds e dsfed, hanno la stessa lunghezza, che ` pari al numero di dati, ossia e 8192. Si vede dal graco fatto che i picchi, corrispondenti a frequenze non incluse nella banda passante del ltro usato, vengono considerevolmente attenuati: la sola componente frequenziale trasmessa ` quella a frequenza f = 3.5 102 e Hz (periodo 28.6 s) e infatti la banda passante copre, per il ltro bpf_m e con Tc = 5 s, lintervallo da 2.15 102 a 4.64 102 Hz (periodi da 21.6 a 46.5 s). Lo spettro permette di osservare periodi da 10 s a 11.4 ore, corrispondenti rispettivamente alla frequenza di Nyquist ed alla prima frequenza non nulla in assenza di padding di zeri, N1 c (frequenza di Rayleigh). T Suggerimento: gracare insieme modDS e modDSFED in colori diversi. Inoltre fare un graco del modulo quadro dello spettro in dB (gracare 20 log 10(modDS), ecc...), per poter ragionare sullattenuazione.

Scartando i primi Mf 1 = 54 campioni delloutput (cio` conservandone 8192e 54=8138) e ricalcolando lo spettro su 8192 valori di frequenza, si pu` notare che o dal punto di vista spettrale quei 54 dati iniziali cambiano poco o nulla il risultato. Ecco il procedimento: dsfed_ok=dsfed(55:8192); DSFED_OK=fft(dsfed_ok,nfft); (DFT su 8192 punti; padding automatico) modDSFED_OK=abs(DSFED_OK); loglog(f,modDSFED_OK)

1.3.3

Filtraggio mediante convoluzione lineare

Abbiamo detto che escludendo dal risultato del ltraggio i transitori iniziale e nale, con Mf = 55 campioni per la risposta allimpulso del ltro e N = 8192 dati, la convoluzione lineare contiene N Mf + 1 = 8192 55 + 1 = 8138 valori. Accettando invece i transitori iniziale e nale la convoluzione ` pi` lunga: vie u sto che ogni coda di transitorio dura quanto lordine del ltro, si avranno (N Mf + 1) + (Mf 1) = N = 8192 valori includendo il solo transitorio iniziale (` il caso corrispondente alluso di filter) e (N Mf + 1) + 2(Mf e 1) = N + Mf 1 = 8192 + 55 1 = 8246 valori includendo sia il transitorio iniziale, sia quello nale. Questultimo caso rappresenta appunto il comportamento della function conv che implementa la convoluzione lineare: questo lo si

36 vede dando i comandi25


dsfconv=conv(hm,ds); length(dsfconv) il cui risultato ` proprio 8246. Il nome dsfconv ricorda che sono dati ds filtrati e via convoluzione). Ora grachiamo i primi 100 dati, notando il transitorio iniziale: plot(dsfconv(1:100)) Poi grachiamo gli ultimi 100 dati, notando il transitorio nale: plot(dsfconv(end-99:end)) e inne i valori di dsfconv corrispondenti agli ultimi 100 di dsfed, sovrapposti a quelli, notando che sono uguali: plot(dsfconv(N-99:N),m,LineWidth,2) hold on plot(dsfed(end-99:end),k) hold off Riassumendo, per tenere in considerazione solo la parte pulita delloutput, si dovr` a usando conv, scartare Mf 1 = 54 dati in testa e in coda alla serie ltrata, conservando cos` solo 8246-108 = 8138 dati (dal 55o al 8192o); usando filter, scartare Mf 1 = 54 dati in testa, conservandone 8192-54 = 8138 (dal 55o al 8192o). Facoltativamente, lo studente ripeta ora il ltraggio con filter e/o conv ed il calcolo dello spettro di ampiezza dei dati ltrati utilizzando gli altri due ltri, bpf_h e bpf_l, e commenti i risultati.

1.3.4

Filtraggio nel dominio della frequenza

Teoricamente, questa via prevedrebbe di fare il prodotto, elemento per elemento, della trasformata dei dati con quella della risposta allimpulso del ltro (cio` con e la risposta in frequenza); in seguito, antitrasformando, si otterrebbero i dati ltrati. Le singole trasformate, per evitare laliasing nel dominio del tempo, dovrebbero essere calcolate su un numero di campioni suciente a contenere tutta la lunghezza della convoluzione. Esiste per` in MATLAB una function che implementa questo metodo: si chiao ma fftfilt e ltra una sequenza usando un eciente algoritmo, basato sulla FFT, detto overlap-add (una tecnica adatta solo a ltri FIR). Supponendo di voler ltrare i dati sintetici ds col ltro bpf_m in questo modo e controllare poi25 Conviene, come argomenti di input a questa function, dare prima il vettore pi` corto e u poi il pi` lungo: in questo modo essa lavora pi` velocemente. u u


37

il risultato ottenuto con un graco temporale e uno spettrale, si proceder` cos` a : dsff=fftfilt(hm,ds,length(ds)); length(dsff) (risulta uguale a 8192: si tratta di un caso analogo a quello di filter, con transitorio iniziale); plot(dsff(1:100)) (visualizza il transitorio iniziale, che non ci preoccuperemo di scartare). Poi calcoleremo lo spettro di ampiezza come al solito DSFF=fft(dsff,nfft); modDSFF=abs(DSFF); e lo gracheremo, sovrapposto a quello dei dati ds, in colori diversi. Il vettore frequenze ` sempre lo stesso dei casi precedenti, con 8192 frequenze. Si e pu` vedere, dal graco che si ottiene, che il risultato ` quello atteso, in base alle o e frequenze contenute nei dati sintetici ed alla banda passante del ltro usato.

1.3.5

Filtraggio a fase nulla

Accenniamo al fatto che in MATLAB ` implementato anche un metodo di le traggio detto zero-phase forward and reverse (digital) ltering, tramite la function filtfilt (sintassi: y=filtfilt(b,a,x) con ovvio signicato dei simboli). Dopo aver ltrato i dati in avanti, la sequenza ltrata viene rovesciata e ltrata nuovamente; loutput della function (y) ` il risultato del secondo ltraggio, e rovesciato temporalmente. Il risultato ` uno sfasamento nullo, ossia nessun rie tardo; lo spettro di ampiezza dellinput viene modicato secondo il quadrato del modulo della risposta in frequenza del ltro. Nella function vengono pure applicati accorgimenti volti a minimizzare i transitori iniziale e nale. Lingresso deve essere pi` lungo di tre volte lordine del ltro, che pu` essere sia FIR, sia u o IIR. Con questo metodo non verranno eseguite prove.

1.3.6

Filtraggio con decimazione

Si voglia ora ltrare gli stessi dati sintetici ds, caratterizzati da Tc = 5 s, fN y = 101 Hz, in modo tale da preservare la componente di frequenza f = 4. 103 Hz; tale frequenza dovr` allora cadere entro la banda passante del ltro che a si user`. Si supponga inoltre di voler riutilizzare il ltro bpf_m, la cui banda a passante, espressa in termini di frequenza adimensionale , va da 0.1075 a 0.232. Con Tc = 5 s ci` implicherebbe frequenze da 2.15 102 Hz a 4.64 102 Hz, o cio` un intervallo da cui f ` esclusa. e e

38


Se per` si decimano i dati con un fattore K = 8, col che il passo di campionao mento passa a Tc = KTc = 40 s, lo stesso intervallo di frequenze adimensionali viene riportato ad un intervallo di frequenze in Hz contenente f .

Vericare col calcolo: qual ` il nuovo intervallo? In particolare, quanto vale la e nuova frequenza di Nyquist? Risposte: il nuovo intervallo va da 2.69 103 a 5.80 103 Hz; la nuova frequenza di Nyquist ` 1.25 102 Hz. e Per il necessario pre-ltraggio anti-aliasing occorrerebbe un ltro ideale passabasso con frequenza di taglio uguale alla nuova frequenza di Nyquist dopo la 1 1 th decimazione, ossia, in termini adimensionali, c = 2K = 16 = 0.0625. Non essendo realizzabile un tale ltro ideale, e dovendosi accontentare di un passa basso con una banda di transizione di larghezza non nulla, si progetta un ltro avente un limite di banda passante inferiore al cuto teorico, in modo tale th che alla frequenza c il modulo della risposta in frequenza del ltro anti-aliasing realizzabile sia gi` scesa ad un valore accettabilmente basso. a La risposta allimpulso di un ltro passa-basso di ordine 54 con le caratteristiche desiderate ` contenuta nel le lpf_dec.mat. e Nella Fig. 1.4 ` riportato il modulo della corrispondente risposta in frequenza: e il limite di banda ` c = 0.05 (fc = 1. 102 Hz con Tc = 5 s, contro e 1.25 102 Hz del cuto teorico) e ad una frequenza adimensionale di 0.0625 lattenuazione ` gi` di circa 11 dB (arriva a 30 dB a c = 0.07). Esso rimuover` e a a dunque ecientemente tutte le componenti frequenziali di frequenza superiore alla nuova frequenza di Nyquist. Dopo aver copiato nella propria directory/cartella di lavoro il le e aver caricato nel workspace, con load, la corrispondente risposta allimpulso che si chiama haa, si dovr`: a eseguire la convoluzione dei dati ds col ltro anti-aliasing, scartando i transitori: Mf= length(haa)% K=8% dsfaa=conv(haa,ds);% length(dsfaa) dsfaa_ok=dsfaa(Mf:N);% length(dsfaa_ok) i filtri hannola stessa lunghezza fattore di decimazione dati filtrati con filtro anti-aliasing transitori scartati = 8,

eettuare la decimazione, prendendo soltanto un dato ogni K partendo dal primo ed andando no a...?

Considerando che N = 8192 e Mf = 55 (ordine Mf 1 = 54), avendo scar-


39

10 0 10 20 30 40 50 60 70 80 90 100 110 120 130 0 0.1 0.2 1.1 1 0.9 0.8 0.7

20 log10 |H |

|H |0.3 0.4 0.5

0.6 0.5 0.4 0.3 0.2 0.1

0 0

0.1

0.2

0.3

0.4

0.5

Figura 1.4: Filtro lpf dec Risposta in frequenza: modulo quadro in dB (a sinistra) e modulo in scala lineare (a destra), in funzione della frequenza adimensionale . tato i transitori disponiano di N Mf campioni su cui eettuare la decimaN Mf +1 = 819254 zione con un fattore K = 8. La parte intera, sia NN, di K 8 ` 1017, quindi lultimo campione da prendere nella decimazione ha indice e NN*K = 10178: NN=fix((N-Mf+1)/K) dsfdec=dsfaa_ok(1:K:NN*K);%decimazione length(dsfdec) Questi dati ds ltrati con ltro anti-aliasing e decimati saranno 1017; eseguire la convoluzione dei dati decimati col ltro passa-banda: dsfdecf=conv(hm,dsfdec); length(dsfdecf) Questi dati ds ltrati con ltro anti-aliasing, decimati e nuovamente ltrati saranno di lunghezza 1017+54 = 1071 con i transitori; poi i transitori verranno scartati: dsfdecf_ok=dsfdecf(Mf:NN);%transitori scartati length(dsfdecf_ok) ottenendo una lunghezza di 1071-54*2 = 963 senza transitori; calcolarne lo spettro di ampiezza e gracarlo nel modo usuale. Attenzione al fatto che il set di frequenze discrete non ` lo stesso: decimando cambia e infatti la frequenza di Nyquist.

40


Lo spettro che si otterr`, confrontato con quello dei dati di partenza, mostrer` a a con chiarezza il tipo di ltraggio eettuato. Si noti come la decimazione com1 porti, tra laltro, un fattore K = 1 sullampiezza della trasformata. 8 Riessione: Nella pratica, leliminazione dei transitori in genere non ` strettae mente necessaria, perch come gi` detto dal punto di vista spettrale i transitori e a hanno un peso non determinante. Pertanto pi` semplicemente si potrebbe procedere come segue: u ...... Mf= length(haa) K=8 ......%dati di partenza e loro spettro di ampiezza su nfft=N frequenze dsfaa=conv(haa,ds); dsfdec=dsfaa(1:K:end); dsfdecf=conv(hm,dsfdec); dsfdecf=dsfdecf-mean(dsfdecf); nfft1=2048 SPE=fft(dsfdecf,nfft1); modSPE=abs(SPE); ff_Ny=1/(2*T_c*K); deltff=1/(nfft1*T_c*K); ff=0:deltff:2*ff_Ny-deltff; loglog(ff(1:nfft1/2+1),modSPE(1:nfft1/2+1),r) Si ottiene in questo modo un graco da zero alla nuova frequenza di Nyquist e su 1025 frequenze, a cui sovrapporre un analogo graco di modDS, da zero alla vecchia frequenza di Nyquist su 4097 frequenze. In questo caso i transitori non eliminati sono due perch conv viene usata due e volte. Se si confronta lo spettro ottenuto per ultimo con quello precedente in cui i transitori venivano eliminati, si nota che il loro eetto ` visibile, ma non e tale da alterare radicalmente il risultato.

Facoltativamente, lo studente ripeta ora il ltraggio con decimazione applicando gli altri due ltri.

1.41.4.1

Caratteristiche di un ltro numericoPremessa

Lo scopo dellesercitazione ` lo studio, nei domini del tempo e della frequenza, e delle caratteristiche di un ltro FIR di cui inizialmente supporremo di conoscere

1.4. CARATTERISTICHE DI UN FILTRO NUMERICO la risposta allimpulso h[n], di lunghezza Mf .

41

Si inizier` col graco (stem plot ) della h[n]. Seguir` il calcolo della risposta in a a frequenza del ltro, data dalla DTFT di h[n], secondo la formulaMf 1

H(ej ) =n=0

h[n]ejn .

H(ej ) verr` valutato su un insieme sucientemente tto di frequenze angolari a discrete in 0 < , in modo da poterla gracare in modulo e fase (il che in pratica corrisponde a sostituire la DTFT con una DFT, calcolata via FFT). Supponendo poi, al contrario, di conoscere H(ej ) ma non h[n], si vedr` come ` a e possibile trovare la risposta allimpulso di un ltro FIR di cui sia nota la risposta in frequenza. La procedura verr` applicata ai quattro ltri FIR a fase lineare a utilizzati per lesercitazione precedente. Come ultima fase dellesercitazione, si studier` come operare per dedurre la risposta allimpulso di un ltro, conoscendo a linput e loutput (deconvoluzione).

1.4.2

Risposta in frequenza

Prima di tutto si visualizzi landamento delle risposte allimpulso dei quattro ltri con degli stem plots: ad esempio stem(hm) Si vede bene la simmetria pari dei coecienti attorno al punto n = (di che tipo (I, II, II o IV) ` questo ltro FIR a fase lineare?). e Lo studente operi analogamente sulle altre tre risposte allimpulso. Per il calcolo della risposta in frequenza sulla met` superiore del cerchio unitaa e rio nel piano z, ad esempio su NF2F T = 512 punti, listruzione ` per il ltro bpf_m [Hm,wm]=freqz(hm,1,512); dove largomento 1 sta per lunico coeciente di tipo a, nel caso di un ltro FIR: infatti la routine pu` calcolare la risposta in frequenza di un qualsiasi o ltro, sia IIR sia FIR. Il secondo output, wm, ` il vettore di frequenze angolari discrete nellintervallo e 0 < , su cui Hm viene campionata. La separazione fra i valori di ` e 2 2 = 1024 = 0.006134. NF F T Se si desiderasse invece in output il corrispondente vettore di frequenze adimensionali nellintervallo da 0 < 0.5, si assegnerebbe una frequenza di campionamento unitaria:M1 2

= 27

42


f_c=1 [Hm,num]=freqz(hm,1,512,f_c);

Ora si grachi il modulo della risposta in frequenza: ad esempio, per un graco del modulo quadro in dB in funzione di , una forma di rappresentazione molto usata, si scrive plot(wm,20*log10(abs(Hm))) oppure per un graco del modulo in scala lineare-lineare plot(wm,abs(Hm))

Per la fase si scrive plot(wm,unwrap(angle(Hm))) dove il ruolo della function unwrap ` quello di dispiegare la fase, eliminando e i salti di fase superiori a 2, che sono discontinuit` ttizie introdotte dalluso a della funzione arcotangente nel calcolo della fase stessa: la fase infatti ` calcoe lata proprio come arcotangente del rapporto tra la parte immaginaria e quella reale di Hm.

Lo studente ripeta via via la procedura anche per gli altri tre ltri, esaminando ogni volta i graci ottenuti.

Riessioni Le gure per i moduli delle risposte in frequenza concordano con quelle date nella sezione 1.3? Con quale metodo presumibilmente sono stati progettati i ltri? Quale minima attenuazione si ha tra banda passante e banda oscura (si osservi, per ciascun ltro, il graco in dB)? Nello stesso graco, quanti dB sono ampie allincirca le oscillazioni in banda passante? Il ripple (si veda il graco in scala lineare) ` maggiore in banda passante e (p ) o attenuata (s )? Qual `, approssimativamente, il rapporto ep s ?

Qual`, approssimativamente, la larghezza delle bande di transizione? e A proposito della fase: si pu` aermare che essa ` lineare? o e Con quale limitazione?

1.4. CARATTERISTICHE DI UN FILTRO NUMERICO

43

Come si dedurrebbe dal graco il ritardo di fase in numero di campioni, introdotto dal ltro sul segnale in uscita?

1.4.3

Risposta allimpulso

Intendiamo qui ricavare, nota che sia la risposta in frequenza di un ltro, la corrispondente risposta allimpulso. A questo scopo occorre, anzich` la risposta e in frequenza su frequenze angolari discrete

Esercitazioni MATLAB

Documents

Transcript of Esercitazioni MATLAB