Download Tutorial Matlab - nivel basico and more Study Guides, Projects, Research Matlab skills in PDF only on Docsity!
APPUNTI dall'Esercitazione su MATLAB™
Tenuta dal Prof. Zanzi
Politecnico di Milano, 07 - Aprile - 2000
ESEMPIO # 1 - Introduzione; gli operatori .* .^ ./ ; la FFT; comandi per la grafica ESEMPIO # 2 - Generazione di un rumore bianco ESEMPIO # 3 - Un esempio di filtro a fase lineare ESEMPIO # 4 - Filtri FIR; utilizzo dei comandi fir1, filter, filtfilt ESEMPIO # 5 - Matrici e vettori: il comando reshape; la convoluzione 2D ESEMPIO # 6 - Introduzione al 2D: comando meshgrid e comandi per la grafica ESEMPIO # 7 - FFT2; somma di sinusoidi Autore: Nicola Alagia, [email protected] Nota: Gli appunti non sono stati revisionati da alcun professore.
ESEMPIO # 1
Iniziamo costruendo un vettore di numeri; in Matlab c'e' un modo comodo per fare questo. Supponiamo di voler costruire un vettore 'v' di 10 numeri ordinati da 1 a 10; il comando e': » v=1: La variabile 'v' contiene quindi i seguenti elementi v = 1 2 3 4 5 6 7 8 9 10 Si noti che Matlab usa come passo predefinito il numero 1; se ne vogliamo usare un altro dobbiamo indicarlo esplicitamente. Ecco un esempio con passo - 2 ed un altro con passo 0.5: » v=10:-2: v = 10 8 6 4 2 0 » v=1:0.5:3; Per snellire la scrittura, al posto di 0.5 si puo' indicare semplicemente .5 (v=1:.5:3); Se alla fine dell'istruzione Matlab viene usato il punto e virgola ';' il risultato dell'operazione non viene visualizzato a schermo; questo e' ovviamente fondamentale quando si usano vettori e matrici di grandi dimensioni perche' evita notevoli perdite di tempo. Se si vuole mettere un commento dopo un'istruzione bisogna usare il carattere '%'; Matlab ignora il contenuto di tutta la riga dopo il suddetto carattere: » v=1:0.5:3; % Questo e' un commento » % e questa e' una intera riga di commento Quando si lavora con matrici e vettori il risultato di una moltiplicazione (righe per colonne) dipende anche dall'ordine degli operandi; ad esempio, indicando con A(m,n) una matrice con m righe ed n colonne, con V(1,n) un vettore riga e con V(n,1) un vettore colonna, e' noto che: A(3,2) * B(2,4) = C(3,4) (il risultato e' una matrice di 3 righe e 4 colonne) B(2,4) * A(3,2) = ?? (la moltiplicazione non e' definita) V1(1,n) * V2(n,1) = k (il risultato e' uno scalare) V2(n,1) * V1(1,n) = D(n,n) (il risultato e' una matrice n x n) In Matlab esiste un altro tipo di moltiplicazione (moltiplicazione di array) che, presi due vettori V(1,n) moltiplica elemento per elemento e restituisce come risultato ancora un vettore del tipo V(1,n); per fare questo bisogna usare l'operatore '.' , cioe': V1(1,n) . V2(1,n) = V3(1,n) (il risultato e' un vettore) Un esempio numerico riesce a chiarire subito; digitiamo al prompt di Matlab i seguenti comandi: » v1=[1 2 3]; » v2=[1 4 5]; » v3=v1.*v v3 = 1 8 15
Per visualizzare con dei grafici quello che stiamo facendo, usiamo un po' di nuovi comandi. subplot(m,n,x) serve per dividere lo schermo in una matrice di m righe e n colonne in modo da poter disegnare piu' grafici contemporaneamente; il parametro x serve per indicare quindi in quale parte dello schermo vogliamo disegnare (nota che spesso i comandi sono piu' facili da usare che da tentare di spiegare !!! ). Chi preferisce puo' usare il comando subplot anche senza le virgole, cioe' subplot(mnx)o addirittura anche senza parentesi, cioe' subplot mnx; stem (v) traccia il grafico discreto dei valori contenuti nel vettore v plot (v) simile a stem, ma interpola linearmente i valori tra un campione e l'altro grid serve per disegnare la griglia nell'ultimo grafico visualizzato title('titolo del grafico'); xlabel('etichetta asse delle ascisse'); ylabel('etichetta asse delle ordinate'); Torniamo al nostro esempio e proviamo a dare i seguenti comandi: » subplot(211); % divido lo schermo grafico in 2 parti » plot(p); % disegno i campioni casuali nel tempo » title('Segnale casuale generato da randn(1,16)'); » xlabel('Tempo'); » ylabel('Valore'); » subplot(212); % lavoro nella seconda meta' dello schermo grafico » plot(P); % disegno (o meglio TENTO DI DISEGNARE) la trasformata » title('Spettro del segnale casuale'); » xlabel('Frequenze'); » ylabel('Valore'); Riportiamo qui sotto i grafici generati da Matlab e notiamo che qualcosa non ha funzionato! 0 2 4 6 8 10 12 14 16
- 0 1 2 Segnale casuale generato da randn(1,16) Tempo Valore -3 -2 -1 0 1 2 3 4 5 6
0 5 10 Frequenze Valore Spettro del segnale casuale
Quale errore abbiamo commesso? Ci siamo dimenticati che la trasformata di Fourier di un segnale e', a priori, un vettore di numeri complessi e Matlab lo considera come tale. Per questo motivo il comando plot non e' adatto: ci sarebbe bisogno di un grafico in 3 dimensioni! Solitamente, come ben sappiamo, si preferisce rimanere nelle piu' comode 2 dimensioni e disegnare soltanto una parte delle informazioni contenute nel vettore complesso P. Per fare cio' abbiamo a disposizione i 4 seguenti comandi : abs(P) restituisce il valore assoluto del vettore; angle(P) restituisce la fase; real(P) restituisce la parte reale; imag(P) restituisce la parte immaginaria. Decidiamo quindi di rappresentare su un nuovo grafico (grazie al comando figure) il modulo e la fase della fft del nostro segnale casuale. » figure; % creo una nuova finestra grafica » subplot(211); » plot(abs(P)); % disegno il modulo » title('Modulo della fft del segnale casuale p(n)') » subplot(212); » plot(angle(P)); % disegno la fase » title('Fase della fft del segnale casuale p(n)') » grid % visualizzo la griglia nel grafico della fase 0 2 4 6 8 10 12 14 16
Modulo della fft del segnale casuale p(n) 0 2 4 6 8 10 12 14 16
Fase della fft del segnale casuale p(n)
Frequenze in Hz Fase della fft del segnale casuale p(n) Note sulla definizione dell'asse delle frequenze:
- dobbiamo ricordarci sempre che tanto il segnale nel tempo p(n), quanto la sua trasformata P(k)=FFT(p(n)), sono sottointese essere due sequenze periodiche. Per questo motivo la frequenza a 256 Hz viene esclusa: e' infatti del tutto equivalente alla frequenza 0 (256Hz e' cioe' ancora la componente continua);
- l'aver definito il vettore f a partire da zero non ha nulla a che vedere con il discorso fatto in precedenza quando abbiamo sottolineato che l'indice di un vettore parte sempre da 1. Se vogliamo infatti sapere quanto vale il modulo della fft alla frequenza di Nyquist, qual e' il comando giusto? Possono venire in mente almeno due idee: » P(128) % cioe' il campione a 128 Hz; » P(9) % cioe' il campione numero 9 Quella esatta e' ovviamente la seconda perche' Matlab ragiona sempre in campioni, anche se noi abbiamo definito un asse delle frequenze! La risposta alla nostra richiesta e' quindi: ans =
Quest'ultimo output di Matlab ( ans = 7.6027 ) ci permette di sottolineare un'ultima cosa prima di chiudere il primo esercizio: non avevamo detto che i campioni della fft sono da considerarsi sempre complessi? Non ci dovevamo percio' aspettare una risposta del tipo: ans = 7.6027 + 0.0000i? In realta' noi sappiamo bene (e quindi anche Matlab lo sa) che se il segnale nel tempo e' reale, la sua fft ha i campioni alla continua ed al Nyquist anch'essi reali. Come conferma, digitiamo 'P' al prompt di Matlab in modo da vederne il contenuto, notando che i campioni 1 e 9 sono numeri reali. Prima di incominciare il secondo esercizio puliamo la memoria con il comando clear
ESEMPIO # 2
Vogliamo generare un rumore casuale con spettro costante (cioe' bianco), lavorando nel dominio delle frequenze. Generiamo il modulo e la fase separatamente. Costruiamo un vettore A lungo 100 elementi per il modulo e imponiamolo unitario a tutte le frequenze tranne che alla continua ed al Nyquist. Usiamo poi il comando plot per visualizzarlo e, poiche' l'output grafico ci piace poco, con il comando axis regoliamo gli intervalli di valori sugli assi X e Y. » A=ones(1,100); % genera una matrice 1x100 (cioe' un vettore di tutti 1) » A(1)=0; % cosi' abbiamo valor medio nullo (nel tempo)^1 » A(51)=0; » plot(A); » axis([0 100 0 2]); % asse X da 0 a 100, Y da 0 a 2 » title('Modulo dello spettro di un rumore bianco'); Ricordiamo che in Matlab l'indice di un vettore parte sempre da 1, quindi la continua e' A(1), mentre il Nyquist e' A( N/2 + 1 ) = A(51). Il tutto e' evidenziato dalla figura riportata qui sotto. 0 10 20 30 40 50 60 70 80 90 100
Modulo dello spettro di un rumore bianco Nota che il modulo, essendo costante, soddisfa la simmetria pari (condizione necessaria, ma non sufficiente, per avere il segnale reale nel tempo). Dobbiamo adesso generare la fase (che in questo problema deve essere casuale), stando quindi attenti a soddisfare la simmetria dispari. Un modo semplice e' generare la prima meta' dei campioni in modo casuale, forzare la continua e il Nyquist ad essere nulli ed infine forzare la simmetria dispari per la seconda meta' dei campioni. » a(1)=0; % forzo la continua a fase nulla » a(2:50)=2pirand(1,49)-pi; % genero un vettore di 49 campioni casuali tra - pi e +pi » a(51)=0; % forzo il Nyquist » a(52:100)=-a(50:-1:2); % forzo la simmetria dispari rispetto al Nyquist » subplot(211); % visualizzo il risultato... » plot (a); (^1) Ricorda che, dalla formula di trasformazione della DFT, il campione alla componente continua e'
X(k=0) = ∑
= N n
x n
1
( ). Quindi, in generale, il valor medio nei tempi e'
N
X k
x n
N
x
N n
1
=
» whos Name Size Elements Bytes Density Complex a 1 by 100 100 800 Full No A 1 by 100 100 800 Full No P 1 by 100 100 1600 Full Yes Calcoliamo adesso l'antitrasformata di P ed otteniamo il nostro rumore bianco. » p=IFFT(P); Noi sappiamo che p e' un vettore con tutte le componenti reali (abbiamo 'faticato' tanto per costruirlo!), ma Matlab non lo sa; calcolando la IFFT di una sequenza complessa si aspetta in generale (come e' giusto che sia) un'altra sequenza complessa. Ma e' proprio vero che la parte immaginaria di p=IFFT(P) e' esattamente nulla? Proviamo a disegnarla: » subplot(211); » plot(real(p)); title('Parte reale del rumore bianco'); » subplot(212); » plot(imag(p)); title('Parte immaginaria del rumore bianco'); 0 10 20 30 40 50 60 70 80 90 100
Parte reale del rumore bianco 0 10 20 30 40 50 60 70 80 90 100
x 10 -16 (^) Parte immaginaria del rumore bianco A prima vista la parte immaginaria non e' nulla, ma notiamo che sull'asse delle ordinate e' indicato il numero 10-^16 : vuol dire che c'e' un piccolissimo errore dovuto alla quantizzazione perche' ovviamente Matlab lavora con un numero grande ma finito di bit. L'entita' dell'errore dipende anche dalla capacita' del Personal Computer che si sta usando, comunque deve essere <10-^12. Per risolvere questo problema e' bene estrarre la sola parte reale del vettore p: » p=real(p);
Vediamo infine cosa succede se viene a mancare la simmetria coniugata dello spettro anche per colpa di un solo elemento del vettore. Poniamo quindi, a titolo di esempio, A(100)=0.5 e ripetiamo i comandi necessari fino a visualizzare la parte immaginaria della sequenza nei tempi. » A(100)=0.5; » P=A.exp(ja); » p=IFFT(P); » subplot(211); » plot(real(p)); » title('Parte reale del rumore bianco'); » subplot(212); » plot(imag(p)); » title('Parte immaginaria del rumore bianco'); 0 10 20 30 40 50 60 70 80 90 100
x 10
- (^3) Parte immaginaria del rumore bianco Si vede chiaramente che e' comparsa una sinusoide di ampiezza 5*10-^3 , e quindi non trascurabile.
ESEMPIO # 4
In questo esercizio vediamo come costruire un filtro FIR (con il comando fir1) e come filtrare una sequenza casuale con i comandi filter e filtfilt. Per prima cosa generiamo una sequenza di 100 campioni casuali distribuiti normalmente con media nulla e varianza 9 e la visualizziamo. » s=randn(1,100)*sqrt(9); % segnale casuale » figure; » plot(s); » title('Rumore gaussiano con varianza 9 e media nulla'); 0 10 20 30 40 50 60 70 80 90 100
Rumore gaussiano con varianza 9 e media nulla Prima di proseguire, e' bene sapere che il comando fir1(N,w) costruisce un filtro (di default e' un passa-basso) utilizzando il metodo delle finestre (di default usa la finestra di Hamming). I parametri da specificare sono l'ordine del
filtro N (il filtro sara' lungo quindi N+1 campioni) e la frequenza di taglio w ( con 0 w w=1 e' la pulsazione
di Nyquist) Se w e' un vettore, w=[w1 w2], allora indica un intervallo e Matlab realizza un filtro passa-banda con banda passante da w1 a w2. Con il comando fir1(N,[w1 w2],'stop') realizziamo un arresta-banda (nota che N deve essere pari). Usando invece fir1(N,w,'high') il risultato e' un passa-alto (anche qui N deve essere pari). Esiste anche un comando fir2 per creare un filtro specificando alcune frequenze ed il valore che il modulo deve avere Per maggiori dettagli digitare help fir 1 e help fir2 al prompt di Matlab. Siamo pronti per creare un filtro FIR di ordine 12, passabanda tra 0,3 e 0,6 rispetto alla frequenza di Nyquist. Dopo filtreremo con il comando y=filter(num,den,x), dove num e den sono numeratore e denominatore del filtro, x e' il segnale in ingresso e y quello in uscita. » b=fir1(12,[0.3 0.6]); % costruisco il filtro » p=filter(b,1,s); % filtro la sequenza s(n) ed ottengo p(n) » figure; » subplot(311); » plot(abs(fft(s))); % disegno lo spettro della sequenza originaria » axis([0 100 0 80]); » title('Modulo della FFT di un rumore gaussiano (nota la simmetria pari;... Nyquist=51)'); » subplot(312); » plot(abs(fft(p))); % disegno lo spettro della sequenza filtrata » title('Spettro filtrato con un passa-banda [0.3 0.6] di ordine 12'); » axis([0 100 0 80]);
Se vogliamo un filtraggio a fase zero dobbiamo filtrare due volte (da destra a sinistra e viceversa); fortunatamente esiste il comando filtfilt che fa proprio questo. Attenzione: cosi' facendo abbiamo raddoppiato l'ordine del filtro. » r=filtfilt(b,1,s); » subplot(313); » plot(abs(fft(r))); » title('Spettro filtrato con un passa-banda [0.3 0.6] di ordine 24 (usando ... il comando FiltFilt)'); » axis([0 100 0 80]); 0 10 20 30 40 50 60 70 80 90 100
Modulo della FFT di un rumore gaussiano (nota la simmetria pari; Nyquist=51) 0 10 20 30 40 50 60 70 80 90 100
Spettro filtrato con un passa-banda [0.3 0.6] di ordine 12 0 10 20 30 40 50 60 70 80 90 100
Spettro filtrato con un passa-banda [0.3 0.6] di ordine 24 (usando il comando FiltFilt) Nota che il taglio e' piu' netto (infatti l'ordine del filtro e' doppio)
» V=[1 2 3 4 5 6 7 8 9];
» M=reshape(V,3,3) M = 1 4 7 2 5 8 3 6 9 una possibile applicazione e' la seguente: vogliamo costruire una matrice M quadrata 6x6, di tutti zeri tranne tre elementi a caso che devono avere valore unitario; i passi da compiere sono: » V=zeros(1,36); % parto da un vettore riga di 66 elementi » x=ceil(36rand(1,3)); % scelgo 3 numeri a caso tra 1 e 36 » V(x)=ones(1,3); % utilizzo i tre precedenti numeri casuali come indice per scegliere quali elementi di V devono avere valore unitario » M=reshape(V,6,6) % costruisco la matrice 6x6 a partire dal vettore V M = 0 0 0 0 1 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 Definiamo ora la seguente matrice C a croce » C=[0 2 0 ; 2 2 2 ; 0 2 0] C = 0 2 0 2 2 2 0 2 0 e replichiamola in M centrata proprio nelle posizioni indicate dagli 1 (che sono quindi degli spilli); per farlo convolviamo C ed M con il comando Y=conv2(M,C,'same') dove l'opzione 'same' obbliga Y a mantenere le stesse dimensioni di M » Y=conv2(M,C,'same') Y = 0 0 0 2 2 2 0 2 0 0 2 0 2 2 2 0 0 0 0 2 0 2 0 0 0 0 2 2 2 0 0 0 0 2 0 0
ESEMPIO # 6
Vediamo in questo esercizio alcuni comandi caratteristici dei segnali in due dimensioni. Analogamente al caso 1D, dove definivamo l'asse temporale per poi costruire il segnale da elaborare, nel caso 2D bisogna definire una matrice di punti; per fare questo potremmo ad esempio usare due cicli for: » for x=1: for y=1: a(x,y)=sin(0.2x+0.4y); % sinusoide 2D end; end; e' sicuramente piu' veloce il comando meshgrid, che si usa nel seguente modo: » x=-20:.5:20; % vettore lungo le x » y=-20:.5:20; % vettore lungo le y » [X,Y]=meshgrid(x,y); % costruzione della matrice bidimensionale » A=1./(sqrt(X.^2+Y.^2)+3); % esempio di segnale 2D nota che gli indici del ciclo for devono per forza essere numeri naturali positivi, poiche' servono anche come indici della matrice a(x,y); nel secondo esempio invece, e' stato messo in evidenza che questa restrizione e' superata, quindi: non si usano i cicli for ma il comando meshgrid. Se vogliamo visualizzare su schermo i segnali 2D abbiamo a disposizione diversi comandi, tra i quali imagesc, mesh, contour, surf, plot3; vediamone l'ouput grafico. » figure; » subplot 221; » imagesc(A); » colorbar; % serve per avere la legenda a fianco del grafico » xlabel('asse x'); » ylabel('asse y'); » zlabel('funzione A=...'); » title('comando imagesc'); » subplot 222; » mesh(X,Y,A); » xlabel('asse x'); ylabel('asse y'); zlabel('funzione A=...'); » title('comando mesh'); » subplot 223; » contour(X,Y,A); % per disegnare i livelli » xlabel('asse x'); ylabel('asse y'); zlabel('funzione A=...'); » title('comando contour'); » subplot 224; » surf(X,Y,A); % per disegnare la superficie » xlabel('asse x'); ylabel('asse y'); zlabel('funzione A=...'); » title('comando surf');
ESEMPIO # 7
Generiamo ora due sinusoidi 2D, le sommiamo e ne facciamo la trasformata di Fourier bidimensionale. » N=60; » x=-N:N; % vettore lungo le x » y=-N:N; % vettore lungo le y » [X,Y]=meshgrid(x,y); » a=cos(.2X+.2Y); % coseno con azimuth +45 deg » b=cos(.3X+.5Y); % coseno con azimuth atan(.5/.3)=63 deg » c=a+b; % somma di sinusoidi » C=fft2(c); % DFT in due dimensioni » figure; » subplot 221; imagesc(x,-y,a); title(' a=cos(.2X+.2Y) '); » subplot 222; imagesc(x,-y,b); title(' b=cos(.3X+.6Y) '); » subplot 223; imagesc(x,-y,c); title(' c=a+b '); » subplot 224; imagesc(x/N,-y/N,fftshift(abs(C))); » title('Spettro di c (notare i 4 impulsi)'); » axis([- 30 30 - 30 30]/N); grid; xlabel('f / fNy'); ylabel('f / fNy'); a=cos(.2X+.2Y) -50 0 50
b=cos(.3X+.5Y) -50 0 50
c=a+b -50 0 50
Spettro di c (notare i 4 impulsi) f / fNy f / fNy -0.2 -0.1 0 0.1 0.
Il comando fftshift serve per traslare lo spettro (ricorda che la traslazione e' circolare, data la periodicita') in modo da visualizzare la componente alla frequenza nulla nel centro dello schermo.
