Teoria e Tecniche di Elaborazione...

Università degli Studi di Udine Facoltà di Scienze Matematiche, Fisiche e Naturali

Laurea Specialistica in Informatica A.A. 2004 – 2005

Relazione per il Corso di

Teoria e Tecniche di Elaborazione dell’Immagine

prof. Vito Roberto Stefano Maraspin – Matr. 74625 Stefano Valle – Matr. 69163

Relazione per il Corso di Teoria e Tecniche di Elaborazione dell’Immagine

2

Stefano Maraspin – Stefano Valle

3

Indice – Scheda 1 Esercizio 1 – Scheda 1 .......................................................................................................................10

1.1 Traccia................................................................................................................................10 1.2 Svolgimento .......................................................................................................................10

1.2.1 Referenza alle primitive di Matlab utilizzate.............................................................10 1.3 Codice Matlab ....................................................................................................................10

1.3.1 Codice onda Sinusoidale............................................................................................10 1.3.2 Codice onda Rettangolare ..........................................................................................11 1.3.3 Codice onda Triangolare............................................................................................12

1.4 Risultati ..............................................................................................................................14 1.5 Commenti...........................................................................................................................27

Esercizio 2 – Scheda 1 .......................................................................................................................28 2.1 Traccia................................................................................................................................28 2.2 Svolgimento .......................................................................................................................28

2.2.1 Referenza alle primitive di Matlab utilizzate.............................................................30 2.2.2 Chiamate a funzioni da noi definite in precedenza ....................................................30

2.3 Codice Matlab ....................................................................................................................30 2.3.1 Codice dell’Autocorrelazione fra le tre forme d’onda ...............................................30 2.3.2 Codice della Correlazione incrociata .........................................................................32 2.3.3 Codice della Convoluzione ........................................................................................34

2.4 Risultati ..............................................................................................................................36 2.4.1 Autocorrelazione........................................................................................................36 2.4.2 Correlazione Incrociata ..............................................................................................39 2.4.3 Convoluzione .............................................................................................................42

2.5 Commenti...........................................................................................................................45 Esercizio 3 – Scheda 1 .......................................................................................................................47



3.3 Codice Matlab ....................................................................................................................48 3.3.1 Codice funzione di supporto interv ............................................................................48 3.3.2 Codice funzione gaussiana.........................................................................................48 3.3.3 Codice funzione esponenziale....................................................................................49 3.3.4 Codice funzione impulso ...........................................................................................49 3.3.5 Codice funzione principale ........................................................................................50

3.4 Risultati ..............................................................................................................................51 3.5 Commenti...........................................................................................................................65


4

Esercizio 4 – Scheda 1 .......................................................................................................................72 4.1 Traccia................................................................................................................................72 4.2 Svolgimento .......................................................................................................................72

4.2.1 Referenza alle primitive di Matlab utilizzate.............................................................72 4.3 Codice Matlab ....................................................................................................................73

4.3.1 Codice Rumore Uniforme..........................................................................................73 4.3.1 Codice Rumore Gaussiano.........................................................................................74

4.4 Risultati ..............................................................................................................................74 4.4.1 Rumore a distribuzione uniforme ..............................................................................75 4.4.1 Rumore a distribuzione gaussiana .............................................................................80

4.5 Commenti...........................................................................................................................86 Esercizio 5 – Scheda 1 .......................................................................................................................87



5.3 Codice Matlab ....................................................................................................................88 5.3.1 Codice rumore gaussiano (variante) ..........................................................................88 5.3.2 Codice funzione principale ........................................................................................88

5.4 Risultati ..............................................................................................................................90 5.5 Commenti.........................................................................................................................100


5

Indice – Scheda 2

Esercizio 1 – Scheda 2 .....................................................................................................................103

6.1 Traccia..............................................................................................................................103 6.2 Svolgimento .....................................................................................................................103

6.2.1 Referenza alle primitive di Matlab utilizzate...........................................................103 6.2.2 Chiamate a funzioni da noi definite in precedenza ..................................................103

6.3 Codice Matlab ..................................................................................................................104 6.4 Risultati ............................................................................................................................105 6.5 Commenti.........................................................................................................................112

Esercizio 2 – Scheda 2 .....................................................................................................................114 7.1 Traccia..............................................................................................................................114 7.2 Svolgimento .....................................................................................................................114

7.2.1 Chiamate a funzioni da noi definite in precedenza ..................................................114 7.3 Codice Matlab ..................................................................................................................115

7.3.1 Codice funzione di somma dei segnali ....................................................................115 7.4 Risultati ............................................................................................................................115

7.4.1 DFT di un segnale periodico....................................................................................116 7.4.2 DFT della somma di due segnali sinusoidali ...........................................................123 7.4.3 DFT della somma di un segnale periodico con rumore gaussiano ..........................126

7.5 Commenti.........................................................................................................................132 Esercizio 3 – Scheda 2 .....................................................................................................................134



8.3 Codice Matlab ..................................................................................................................135 8.4 Risultati ............................................................................................................................138 8.5 Commenti.........................................................................................................................151



9.3 Codice Matlab ..................................................................................................................154 9.3.1 Codice Funzione Esercizio.......................................................................................154 9.3.2 Codice Funzione a Gradino......................................................................................156 9.3.3 Codice Derivatore ....................................................................................................157

9.4 Risultati ............................................................................................................................159


6

9.5 Commenti.........................................................................................................................168 Esercizio 5 – Scheda 2 .....................................................................................................................170



10.3 Codice Matlab ..................................................................................................................171 10.3.1 Codice Funzione Esercizio.......................................................................................171

10.4 Risultati ............................................................................................................................173 10.5 Commenti.........................................................................................................................184


7

Indice – Scheda 3

Esercizio 1 – Scheda 3 .....................................................................................................................187


11.2.1 Referenza alle primitive di Matlab utilizzate...........................................................190 11.3 Codice Matlab ..................................................................................................................190

11.3.1 Codice Programma Principale .................................................................................190 11.3.2 Codice Generatore Immagini Sintetiche 1 ...............................................................192 11.3.3 Codice Generatore Immagini Sintetiche 2 ...............................................................193 11.3.4 Codice Loader Immagini Pre-Esistenti ....................................................................194




12.3 Codice Matlab ..................................................................................................................214 12.3.1 Codice Programma Principale .................................................................................214 12.3.2 Codice Generatore Rumore Gaussiano a Media Nulla ............................................215



13.2.1 Referenza alle primitive di Matlab utilizzate..........................................................231 13.2.2 Chiamate a funzioni da noi definite in precedenza ..................................................231

13.3 Codice Matlab ..................................................................................................................232 13.3.1 Codice Programma Principale .................................................................................232 13.3.2 Codice Generatore Filtro Gaussiano 2D ..................................................................235




14.3 Codice Matlab ..................................................................................................................268


8

14.3.1 Codice Programma Principale .................................................................................268 14.4 Risultati ............................................................................................................................270 14.5 Commenti.........................................................................................................................285


9

Scheda 1


10

Esercizio 1 – Scheda 1

1.1 Traccia Si crei e visualizzi un segnale periodico con ampiezza, frequenza e fase date dall’utente. Si trattino i casi delle:

• Onde sinusoidali • Onde rettangolari • Onde triangolari

1.2 Svolgimento Vengono create e visualizzate le tre forme d’onda richieste: Segnale sinusoidale, utilizzando la funzione di libreria Matlab sin, per la rappresentazione

dell’onda sinusoidale, sulla base dei parametri forniti dall’utente (ampiezza, frequenza, fase, intervallo e numero dei campioni da considerare).

Segnale rettangolare, ricavato dal calcolo del segno dell’onda sinusoidale vista in precedenza e moltiplicato per l’ampiezza data.

Segnale triangolare, calcolandolo in tre fasi distinte (ripetute per ogni periodo); dapprima si calcolano i valori nel primo quarto di periodo, poi nel secondo e così via fino ad ottenere un’onda triangolare rappresentata in un intero periodo. In particolare viene calcolato, nel primo intervallo (trattandosi di onda periodica), il coefficiente angolare delle semirette che costituiranno l’onda triangolare:

o Prima fase: da 0 a T/4, semiretta crescente o Seconda fase: da T/4 a 3T/4, semiretta descrescente o Terza fase: da 3T/4 a T, semiretta crescente

1.2.1 Referenza alle primitive di Matlab utilizzate sin(x): restituisce l’onda sinusoidale calcolata sui valori di input abs(x): calcola il valore assoluto del parametro x sign(x): calcola il segno del parametro x length(x): calcola la lunghezza del vettore x mod(a,b): calcola il modulo fra i due parametri a e b

1.3 Codice Matlab

1.3.1 Codice onda Sinusoidale % Funzione sinus (amp, freq, fase, range, ncomp)


11

% Argomenti:

% - amp: ampiezza onda

% - freq: frequenza onda

% - fase: fase onda

% - range: intervallo di valori nella rappresentazione (formato [inizio,fine])

% - ncomp: numero di campioni considerati

%

% Esempio di utilizzo:

% sinus (5, 2, pi, [1,3], 1000)

%

% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

% Autori: Stefano Maraspin & Stefano Valle

% Data: 14/10/2004

% Versione: 1.0

function y = sinus(amp, freq, fase, range, ncomp)

% Viene ricavato il valore del periodo T

T = 1./freq;

% Viene ricavato i valori di ascissa iniziali e finali

xmin = range(1,1);

xmax = range(1,2);

% Viene calcolato l'intervallo tra valori successivi delle ascisse

passo = abs(xmax - xmin) / (ncomp - 1);

% Viene generato un vettore contenente i punti sull'ascissa che saranno

% considerati

x = [xmin : passo : xmax]';

% Per ciascun punto x sull'asse delle ascisse, viene calcolata la relativa

% immagine, applicando la formula per la costruzione di un'onda sinusoidale

y = amp * sin(x.*2*pi./T+fase);

plot(x,y)

1.3.2 Codice onda Rettangolare % Funzione rectan (amp, freq, fase, range, ncomp)

% Argomenti:



% - fase: fase onda



%


12


% rectan (5, 2, pi, [1,3], 1000)

%

% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -


% Data: 14/10/2004

% Versione: 1.0

function y = rectan(amp, freq, fase, range, ncomp)


T = 1./freq


xmin = range(1,1);

xmax = range(1,2);




% considerati

x = [xmin : passo : xmax]'

% Per ciascun punto x sull'asse delle ascisse, viene calcolata la relativa

% immagine, valutando il segno di un'onda sinusoidale di pari caratteristiche

y = amp * sign(sinus(amp,freq,fase,range,ncomp));

plot(x,y)

1.3.3 Codice onda Triangolare % Funzione trian (amp, freq, fase, range, ncomp)

% Argomenti:



% - fase: fase onda



%


% trian (5, 2, pi, [1,3], 1000)

%

% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -


% Data: 14/10/2004

% Versione: 3.0


13

function y = trian(amp, freq, fase, range, ncomp)


T = 1./freq


xmin = range(1,1);

xmax = range(1,2);




% considerati

x = [xmin : passo : xmax]';

% Calcoliamo il valore della fase (inserito in randianti) nel sistema di

% riferimento del periodo T utilizzato

diff_fase = (T * fase) / (2 * pi)

% Per effettuare i calcoli successivi è necessario traslare l'asse di

riferimento

% riportandolo alle coordinate d'origine

base = x - xmin + diff_fase;

% Iteriamo su ciascun elemento dell'ascissa

for i = 1:length(x)

% ne calcoliamo la posizione relativa al periodo considerato

rel_pos = mod(base(i),T);

% e, a seconda del quarto di periodo in cui ci troviamo valutiamo

% i valori di ordinata corrispondenti

if (rel_pos < (T/4))

y(i) = (rel_pos) * amp / (T/4);

elseif (rel_pos > (T/4) & rel_pos < (3*T/4))

y(i) = ((T/2 - rel_pos) * amp / (T/4));

elseif (rel_pos > (3*T/4) & rel_pos < T)

y(i) = ((T/2 - rel_pos ) * (- amp) / (T/4)) - 2*amp;

% (condizione mai verificata, inserita solo per pulizia logica

% in questo modo si è in grado di identificare eventuali errori nei dati)

else

y(i) = 0;

end

end


14

plot(x,y)

1.4 Risultati Si riportano di seguito alcuni dei risultati più significativi, selezionati fra i test effettuati in laboratorio.


15

1. Caso elementare di onda sinusoidale di frequenza unitaria, visualizzata in una finestra di dimensione 2.

Ampiezza 1 Frequenza 1

Fase 0 Range 0 ÷ 1

N. Campioni 1000

Chiamata sinus(1,1,0,[0,1],1000)


16

2. Si ripete la situazione precedente, in questo caso con un numero di campioni considerevolmente inferiore.


Fase 0 Range 0 ÷ 1

N. Campioni 10



17

3. Vengono in questo caso modificati i valori di ampiezza e frequenza.


Fase 0 Range 0 ÷ 1

N. Campioni 1000



18

4. Mantenendo i valori di ampiezza e frequenza visti nel caso 3, si opera ora a livello di fase, con l’intenzione di traslare la forma d’onda di π/2.


Fase π/2 Range 0 ÷ 1

N. Campioni 1000

Chiamata sinus(5,2,pi/2,[0,1],1000)


19

5. Caso elementare di onda rettangolare di frequenza unitaria, visualizzata in una finestra di dimensione 2.


Fase 0 Range 0 ÷ 1

N. Campioni 1000

Chiamata rectan(1,1,0,[0,1],1000)


20



Fase 0 Range 0 ÷ 1

N. Campioni 10



21



Fase 0 Range 0 ÷ 1

N. Campioni 1000



22




N. Campioni 1000

Chiamata rectan(5,2,pi/2,[0,1],1000)


23

9. Caso elementare di onda triangolare di frequenza unitaria, visualizzata in una finestra di dimensione 2.


Fase 0 Range 0 ÷ 1

N. Campioni 1000

Chiamata trian(1,1,0,[0,1],1000)


24



Fase 0 Range 0 ÷ 1

N. Campioni 10



25



Fase 0 Range 0 ÷ 1

N. Campioni 1000



26




N. Campioni 1000

Chiamata trian(5,2,pi/2,[0,1],1000)


27

1.5 Commenti Non abbiamo riscontrato particolari anomalie, o comunque comportamenti degni di nota; come poteva essere deducibile, notiamo infatti che:

• Il parametro fornito come ampiezza determina il valore massimo (in modulo) che l’onda assumerà in ordinata.

• La frequenza fornita determina il periodo della forma d’onda e quindi il numero di ripetizioni che si avranno nella finestra temporale visualizzata. Nel risultato 1 (e analogamente nei risultati 5 e 9) è infatti possibile osservare questo tipo di comportamento: la frequenza data è pari ad 1, la finestra di visualizzazione è anch’essa di dimensione 1 e quindi si ottiene forma d’onda rappresentata in un unico periodo. Nel risultato 3 (e analogamente nei risultati 7 e 11) invece si può notare come raddoppiando la frequenza, a parità di finestra di visualizzazione, raddoppiano i periodi che vengono mostrati.

• La fase fornita determina la traslazione dell’onda verso sinistra o destra (rispettivamente se positiva o negativa). Nel risultato 4 (e analogamente nei risultati 8 e 12) è possibile osservare l’effetto che si ottiene nella rappresentazione impostando un valore di fase positivo (precisamente π/2): l’intera forma d’onda subisce una traslazione verso sinistra di ¼ di periodo.

• Il numero di campioni fornito determina in buona parte la precisione con cui verranno effettuati tutti i calcoli per la visualizzazione delle varie onde. Nel risultato 2 (e analogamente nei risultati 6 e 10) è possibile osservare questo tipo di comportamento: riducendo il numero di campioni si ottiene una forma d’onda più sgranata e meno precisa. La riduzione del numero di campioni permette a Matlab di essere più veloce computazionalmente nell’effettuare la rappresentazione (o eventuali calcoli inerenti la forma d’onda), ma pregiudica la visualizzazione dei risultati.


28


2.1 Traccia Si calcoli e visualizzi l’autocorrelazione delle tre forme d’onda dell’esercizio 1. Si calcolino le correlazioni incrociate e le convoluzioni.

2.2 Svolgimento L’esercizio è stato svolto in tre passi sequenziali, ognuno dei quali prendeva in considerazione le tre forme d’onda base (segnale sinusoidale, rettangolare e triangolare) e ne applicava le funzioni seguenti:

• Autocorrelazione attraverso la funzione di libreria xcorr (nella versione con un solo parametro) che restituisce un vettore contenente la rappresentazione dell’autocorrelazione della funzione data in input.

• Correlazione incrociata fra le varie forme d’onda utilizzando ancora la funzione di libreria xcorr (nella versione a due parametri).

• Convoluzione fra le forme d’onda, utilizzando questa volta la funzione di libreria conv. Da notare che le funzioni utilizzate sia per la correlazione che per la convoluzione forniscono in output un risultato (un vettore) di dimensioni diverse rispetto a quelle del parametro fornito in input. In particolare, se si forniscono in input le rappresentazioni vettoriali di due forme d’onda, la prima corrispondente ad un vettore di dimensione n e la seconda corrispondente ad un vettore di dimensione m, l’output fornito attraverso la convoluzione o la correlazione delle due funzioni sarà un vettore di dimensione m+n+1. Riprendendo alcuni concetti teorici, definiamo di seguito alcuni importanti processi sui segnali utilizzati in questo esercizio. Traslazione Preso un reale τ, definiamo:

)()(' τ+≡ tftf Rt ∈∀

In cui )(' tf rappresenta la traslazione nel dominio dei tempi del segnale originale )(tf .

Correlazione Dati due segnali analogici e periodici dello stesso periodo )(1 tf e )(2 tf , definiamo l’integrale di

correlazione come il numero complesso ottenuto dalla seguente:


29

dttftfT

R

T

T

⋅+⋅= ∫−

2

2

2*

112 )()(1 τ

In cui abbiamo che:

• T1 rappresenta il fattore di normalizzazione della funzione

• )(*1 tf è il complesso coniugato della funzione 1f

• )()( 2*

1 τ+⋅ tftf è il prodotto termine a termine (in t) delle due funzioni considerate (in cui

per la prima prendiamo il complesso coniugato) L’integrale di correlazione quindi si definisce attraverso i seguenti passi:

• Traslazione di 2f di un valore τ

• Moltiplicazione per il complesso coniugato di 1f

• Calcolo dell’integrale della funzione ottenuta Tale integrale, quando esiste, definisce una funzione )(12 τR , chiamata funzione di correlazione

incrociata. Questa funzione effettua un confronto fra le due forme d’onda per cercare di quantificare in qualche modo quanto queste si assomiglino. Questo confronto viene effettuato facendo un raffronto fra la funzione 2f traslata di un valore τ variabile e la funzione 1f che invece

viene mantenuta fissa. Evidentemente il momento di correlazione massima si avrà nel momento in cui le due forme d’onda sono allineate, ovvero quando τ = 0. Generalizzando il concetto nel caso di forme d’onda non periodiche, si ottiene la seguente:

dttftfT

R ⋅+⋅= ∫∞

∞−)()(1

2*

112 τ

Nel caso in cui 1f = 2f si parla di autocorrelazione. Posso utilizzarla nella realtà per cercare di

capire se ho a che fare con un segnale periodico, in cui noterei quindi delle ripetizioni nella visualizzazione dell’autocorrelazione. Convoluzione Dati due segnali analogici e periodici dello stesso periodo )(1 tf e )(2 tf , definiamo l’integrale di

convoluzione come il numero complesso ottenuto dalla seguente:

dttftfT

T

T

⋅−⋅= ∫−

2

2

2112 )()(1 τρ

In cui abbiamo che:


30

• T1 rappresenta il fattore di normalizzazione della funzione

• )()( 21 tftf −⋅ τ è il prodotto termine a termine (in t) delle due funzioni considerate, in cui

nel caso della 2f è effettuato un cambio di segno della variabile (rotazione di 180° dell’asse

x – nota come operazione di folding) L’integrale di convoluzione quindi si definisce attraverso i seguenti passi:

• Traslazione di 2f di un valore τ

• Moltiplicazione per 1f

• Calcolo dell’integrale della funzione ottenuta Anche in questo caso si richiede l’esistenza di tale integrale per poter definire una funzione di convoluzione.

2.2.1 Referenza alle primitive di Matlab utilizzate Chiamate di primitive Matlab utilizzate: xcorr(x): calcola l’autocorrelazione del segnale x xcorr(x,y): calcola la correlazione tra il segnale x e il segnale y abs(x): calcola il valore assoluto del parametro x conv(x,y): calcola la convoluzione tra il segnale x e il segnale y

2.2.2 Chiamate a funzioni da noi definite in precedenza Chiamate a funzioni esterne utilizzate: sinus: rif. Es 1 rectan: rif. Es 1 trian: rif. Es 1

2.3 Codice Matlab

2.3.1 Codice dell’Autocorrelazione fra le tre forme d’onda % Funzione autocorr(amp, freq, range, ncomp)

% Argomenti:






31

%


% autocorr(5, 2, [-1,1], 1000)

%

% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -


% Data: 14/10/2004

% Versione: 1.0

function y = autocorr(amp, freq, range, ncomp)

% Vettore che rappresenta un'onda sinusoidale

valori_sinus = sinus(amp, freq, 0, range, ncomp);

% Vettore che rappresenta un'onda rettangolare

valori_rectan = rectan(amp, freq, 0, range, ncomp);

% Vettore che rappresenta un'onda triangolare

valori_trian = trian(amp, freq, 0, range, ncomp);

% Vettore che rappresenta l'autocorrelazione dell'onda sinusoidale,

% normalizzata secondo il numero di campioni considerato

autocorr_sin = xcorr(valori_sinus) ./ (ncomp);

% Vettore che rappresenta l'autocorrelazione dell'onda rettangolare,


autocorr_rect = xcorr(valori_rectan) ./ (ncomp);

% Vettore che rappresenta l'autocorrelazione dell'onda triangolare,


autocorr_tri = xcorr(valori_trian) ./ (ncomp);

% Vengono ricavati i valori di ascissa iniziali e finali per la

% rappresentazione dei segnali base

xmin_base = range(1,1);

xmax_base = range(1,2);


% rappresentazione dell'autocorrelazione fra i segnali

xmin_corr = - abs(xmin_base - xmax_base);

xmax_corr = abs(xmin_base - xmax_base);


% della funzione originale

passo = abs(xmax_base - xmin_base) / (ncomp - 1);



32

% considerati nella rappresentazione dei segnali di base

x = [xmin_base : passo : xmax_base]';


% considerati nella rappresentazione dei segnali autocorrelati

x_autocorr = [xmin_corr : passo : xmax_corr]';

% Rappresentazione grafica in una matrice di figure 3 x 2

% in cui la prima è l'onda di riferimento e la seconda è

% l'onda autocorrelata

subplot(3,2,1)

plot(x,valori_sinus)

subplot(3,2,2)

plot(x_autocorr,autocorr_sin)

subplot(3,2,3)

plot(x,valori_rectan)

subplot(3,2,4)

plot(x_autocorr,autocorr_rect)

subplot(3,2,5)

plot(x,valori_trian)

subplot(3,2,6)

plot(x_autocorr,autocorr_tri)

2.3.2 Codice della Correlazione incrociata % Funzione corr(amp, freq, range, ncomp)

% Argomenti:





%


% autocorr(5, 2, [-1,1], 1000)

%

% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -


% Data: 14/10/2004

% Versione: 1.0

function y = autocorr(amp, freq, range, ncomp)






33



% Vettore che rappresenta la correlazione tra l'onda sinusoidale e

% l'onda rettangolare

corr_sin_rect = xcorr(valori_sinus,valori_rectan) ./ (ncomp);

% Vettore che rappresenta la correlazione tra l'onda triangolare e


corr_tri_rect = xcorr(valori_trian,valori_rectan) ./ (ncomp);

% Vettore che rappresenta la correlazione tra l'onda triangolare e

% l'onda sinusoidale

corr_tri_sin = xcorr(valori_trian,valori_sinus) ./ (ncomp);

% Vengono ricavati i valori di ascissa iniziali e finali





xmin_corr = - abs(xmin_base - xmax_base);

xmax_corr = abs(xmin_base - xmax_base);








% considerati nella rappresentazione dei segnali autocorrelati

x_corr = [xmin_corr : passo : xmax_corr]';


% in cui vengono visualizzate in ordine: la prima onda base, la seconda

% onda base e la correlazione delle prime due

subplot(3,3,1)


subplot(3,3,2)


subplot(3,3,3)

plot(x_corr,corr_sin_rect)

subplot(3,3,4)



34

subplot(3,3,5)


subplot(3,3,6)

plot(x_corr,corr_tri_rect)

subplot(3,3,7)


subplot(3,3,8)


subplot(3,3,9)

plot(x_corr,corr_tri_sin)

2.3.3 Codice della Convoluzione % Funzione myconv(amp, freq, range, ncomp)

% Argomenti:


% - freq: frequenza ondaù



%


% myconv(5, 2, [-1,1], 1000)

%

% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -


% Data: 14/10/2004

% Versione: 1.0

function y = myconv(amp, freq, range, ncomp)







% Vettore che rappresenta la convoluzione tra l'onda sinusoidale e


conv_sin_rect = conv(valori_sinus,valori_rectan) ./ (ncomp);

% Vettore che rappresenta la convoluzione tra l'onda triangolare e


conv_tri_rect = conv(valori_trian,valori_rectan) ./ (ncomp);


35

% Vettore che rappresenta la convoluzione tra l'onda triangolare e

% l'onda sinusoidale

conv_tri_sin = conv(valori_trian,valori_sinus) ./ (ncomp);






xmin_conv = - abs(xmin_base - xmax_base);

xmax_conv = abs(xmin_base - xmax_base);








% considerati nella rappresentazione della convoluzione fra i segnali

x_conv = [xmin_conv : passo : xmax_conv]';


% in cui vengono visualizzate in ordine: la prima onda base, la seconda

% onda base e l'onda che rappresenta la convoluzione delle prime due

subplot(3,3,1)


subplot(3,3,2)


subplot(3,3,3)

plot(x_conv,conv_sin_rect)

subplot(3,3,4)


subplot(3,3,5)


subplot(3,3,6)

plot(x_conv,conv_tri_rect)

subplot(3,3,7)


subplot(3,3,8)



36

subplot(3,3,9)

plot(x_conv,conv_tri_sin)

2.4 Risultati Si riportano di seguito alcuni dei risultati più significativi, selezionati fra i test effettuati in laboratorio. Il codice della funzione riportato in precedenza, in alcuni casi è stato modificato per permettere di far stampare a schermo solo gli elementi di interesse della particolare sezione.

2.4.1 Autocorrelazione 1. Visualizzazione dell’autocorrelazione dell’onda sinusoidale.


Range 0 ÷ 2 N. Campioni 1000

Chiamata autocorr(5,2,[0,2],1000)


37

2. Visualizzazione dell’autocorrelazione dell’onda rettangolare.





38

3. Visualizzazione dell’autocorrelazione dell’onda triangolare.





39

2.4.2 Correlazione Incrociata 1. Visualizzazione dell’autocorrelazione tra l’onda sinusoidale e l’onda rettangolare.



Chiamata corr(5,2,[0,2],1000)


40

2. Visualizzazione dell’autocorrelazione tra l’onda sinusoidale e l’onda triangolare.





41

3. Visualizzazione dell’autocorrelazione tra l’onda rettangolare e l’onda triangolare.





42

2.4.3 Convoluzione 1. Visualizzazione della convoluzione tra l’onda sinusoidale e l’onda rettangolare.



Chiamata conv(5,2,[0,2],1000)


43

2. Visualizzazione della convoluzione tra l’onda sinusoidale e l’onda triangolare.





44

3. Visualizzazione della convoluzione tra l’onda rettangolare e l’onda triangolare.





45

2.5 Commenti Dai cenni teorici presentati inizialmente possiamo osservare come le due operazioni utilizzate (correlazione e convoluzione) siano in realtà piuttosto simili nella loro definizione e nei risultati che poi forniscono. Analizzando infatti l'integrale di correlazione, notiamo come questo con i segnali discreti corrisponda ad una traslazione delle due funzioni una sull'altra (in particolare una rimane fissa, mentre la seconda trasla) e al calcolo del valore ottenuto dalla loro moltiplicazione punto a punto. Questo comportamento è simile a quello dell'integrale di convoluzione, tranne che per due particolari: la traslazione avviene in senso opposto e la funzione che prima veniva considerata nella rappresentazione complessa coniugata, ora viene utilizzata direttamente nel calcolo nella sua forma naturale. Osservando tutti i risultati ottenuti, sia di correlazione che di convoluzione, vediamo come gli estremi dei vari grafici (ovvero i bordi destro e sinistro della finestra di rappresentazione) tendono sempre ad un valore nullo. Questo comportamento si riscontra per il fatto che le varie onde da noi utilizzate nell'esperimento non corrispondono effettivamente a dei segnali periodici infiniti, come idealmente dovrebbero essere, ma sono racchiusi in una finestra di rappresentazione (definita dal range impostato in input). Questo fatto implica che nella sovrapposizione e relativa traslazione dei due segnali, in ogni momento (tranne inizialmente, quando τ = 0) uno dei due segnali viene moltiplicato anche per valori nulli, corrispondenti al fatto che si è giunti oltre il limite della finestra dell'altro segnale. Nei limiti della finestra in cui si rappresentano la correlazione e la convoluzione si ha l'apice di questo comportamento: si ha la moltiplicazione di un segnale con valori tutti nulli, questo ha come risultato un integrale anch'esso nullo. Quindi ecco spiegato come il segnale di output tenda a 0 ai limiti della finestra. Nell'ipotesi ideale di utilizzare segnali infiniti, questo comportamento non si avrebbe ed anzi il segnale che si otterrebbe avrebbe tutti i valori massimi di eguale valore (ovvero dei picchi massimi ripetuti lungo tutto il grafico). Per fornire ulteriori osservazioni riguardo i risultati ottenuti, procediamo ad analizzare separatamente i tre passi dell'esercizio. Autocorrelazione Si osserva come il risultato dell’autocorrelazione presenti un comportamento simile nel caso di onde sinusoidali e triangolari (nella situazione ideale l’output corrisponderebbe ad un’onda sinusoidale),mentre nel caso dell’onda rettangolare si ha un comportamento diverso (nella situazione ideale questa volta l’output corrisponderebbe ad un’onda triangolare). Ciò è probabilmente dovuto al fatto che le prime due classi di segnali hanno un andamento graduale crescente o decrescente, mentre l’onda rettangolare è composta da un andamento caratterizzato da periodi costanti e ampie discontinuità.


46

Correlazione La correlazione fra le varie forme d’onda presenta ancora un comportamento simile fra i vari casi analizzati (a meno di differenze nei valori di ampiezza). Questo si deve al fatto che la moltiplicazione fra le diverse forme d’onda e il successivo integrale, presenta un andamento continuo paragonabile nel caso ideale ancora ad un’onda sinusoidale. Convoluzione Il comportamento è analogo al caso della correlazione, tranne per un andamento simmetrico del risultato rispetto all’asse delle ordinate. Questo è dovuto al fatto che l’integrale di correlazione calcola il complesso coniugato di una delle due funzioni, mentre l’integrale di convoluzione lascia le due funzioni inalterate: questo causa un cambiamento di segno del risultato.


47


3.1 Traccia Si calcoli e visualizzi la convoluzione tra una forma d’onda periodica ed i seguenti segnali: esponenziale decrescente, impulso quadro, gaussiana.

3.2 Svolgimento Essendo numerose le funzioni che entrano in gioco in questo esercizio si è cercato di modulare il più possibile le funzioni create in modo da semplificare notevolmente l’utilizzo delle stesse. Inizialmente quindi sono stati creati i segnali definiti nella traccia: esponenziale decrescente (funzione espon), impulso quadro (funzione impulso) e gaussiana (funzione gaussiana), di cui si riporta di seguito la rappresentazione.

Il codice per il calcolo della convoluzione è stato creato in precedenza, quindi si è proceduto realizzando la funzione conv_segn che, riutilizzando codice e concetti noti, permettesse di fare la convoluzione di due segnali forniti in input. La stessa funzione poi si occupa di visualizzare a schermo sia i sengnali originali che la convoluzione degli stessi. Si è poi richiamata tale procedura in cui comparivano come parametri sia le tre funzioni periodiche di base creata in precedenza (sinusoidale, rettangolare, triangolare) che i tre segnali definiti dalla traccia. Viene inoltre utilizzata la funzione interv come supporto alla procedura principale.

3.2.1 Referenza alle primitive di Matlab utilizzate Chiamate di primitive Matlab utilizzate: exp(x): calcola l’esponenziale del parametro x zeros(x,1): restituisce un vettore di dimensione x riempito con elementi pari a zero length(x): calcola la lunghezza del vettore x abs(x): calcola il valore assoluto del parametro x conv(x,y): calcola la convoluzione tra il segnale x e il segnale y


48


3.3 Codice Matlab

3.3.1 Codice funzione di supporto interv % Funzione interv(range, ncomp)

% Argomenti:



%


% interv([1,3], 1000)

%

% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -


% Data: 19/10/2004

% Versione: 1.0

function y = interv(range, ncomp)


xmin = range(1,1);

xmax = range(1,2);




% considerati

y = [xmin : passo : xmax]';

3.3.2 Codice funzione gaussiana % Funzione gaussiana(range, ncomp, varianza)

% Argomenti:



% - varianza: varianza della distribuzione gaussiana

%


% gaussiana([0,2], 1000, 10)


49

%

% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -


% Data: 26/10/2004

% Versione: 1.0

function y = gaussiana(range, ncomp, varianza)

% Calcolo dei valori delle ascisse mediante la chiamata della

% funzione 'interv'

x = interv(range, ncomp);

% Calcolo della funzione gaussiana, secondo un valore di varianza

% dato

y = exp(-(x).^2 / 2 * varianza^2);

3.3.3 Codice funzione esponenziale % Funzione espon(range, ncomp, amp, alfa)

% Argomenti:



% - amp: ampiezza segnale

% - alfa: fattore di allungamento della forma del segnale

%


% espon([0,2], 1000, 1, 1)

%

% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -


% Data: 26/10/2004

% Versione: 1.0

function y = espon(range, ncomp, amp, alfa)

% Calcolo dei valori delle ascisse mediante la chiamata della

% funzione 'interv'


% Calcolo della funzione esponenziale, secondo un parametro alfa fornito

% in input

y = amp * exp(-(alfa*x));

3.3.4 Codice funzione impulso % Funzione impulso(amp, ncomp, impda, lungh)

% Argomenti:

% - amp: ampiezza segnale


50


% - impda: distanza dell'impulso dall'origine, espressa in numero di campioni

% - lungh: lunghezza dell'impulso, espressa in numero di campioni

%


% impulso(1, 1000, 100, 50)

%

% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -


% Data: 26/10/2004

% Versione: 1.0

function y = impulso(amp, ncomp, impda, lungh)

% Generazione di un vettore di dimensione ncomp, riempito con una serie

% di valori pari a zero

y = zeros(ncomp,1);

% Assegnazione di un valore pari a amp ai soli elementi che rientrano

% nell'intervallo definito dai nei parametri forniti in input

y(impda+1:impda+lungh+1) = amp;

3.3.5 Codice funzione principale % Funzione conv_segn(range, segnale1, segnale2)

% Argomenti:


% - segnale1: primo segnale di cui calcolare la convoluzione

% - segnale2: secondo segnale di cui calcolare la convoluzione

%


% conv_segn([0,2], sinus (5, 2, 0, [0,2], 1000), gaussiana([0,2],1000,10))

%

% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -


% Data: 26/10/2004

% Versione: 1.0

function y = conv_segn(range, segnale1, segnale2)

% Viene recuperato il numero di componenti

ncomp = length(segnale1);

% Vengono definiti i valori dell'asse delle ascisse mediante una chiamata

% alla funzione interv




51

xmin = range(1,1);

xmax = range(1,2);

% Vengono ricavati i valori di ascissa iniziali e finali per la rappresentazione

% della convoluzione

xmin_conv = - abs(xmin - xmax);

xmax_conv = abs(xmin - xmax);


% considerati per la rappresentazione della convoluzione

x_conv = interv([xmin_conv,xmax_conv],ncomp*2-1);

% Convoluzione dei due segnali, normalizzata rispetto al numero di campioni

% considerato

y = conv(segnale1, segnale2) / ncomp;

% Visualizzazione dei segnali e dei risultati

subplot(2,2,1)

plot(x,segnale1)

subplot(2,2,2)

plot(x,segnale2)

subplot(2,2,[3,4])

plot(x_conv,y)



52

1. Esempio di convoluzione tra i tre segnali periodici di base e il segnale ad impulso

Range -2 ÷ 2 Segnale1 sinus, trian, rectan Segnale2 impulso(1,1000,100,100)

Chiamate conv_segn([-2,2],sinus(5,2,0,[-2,2],1000),impulso(1,1000,100,100)) conv_segn([-2,2],rectan(5,2,0,[-2,2],1000),impulso(1,1000,100,100)) conv_segn([-2,2],trian(5,2,0,[-2,2],1000),impulso(1,1000,100,100))


53


54

2. Esempio di convoluzione tra i tre segnali periodici di base e il segnale gaussiano

Range -2 ÷ 2 Segnale1 sinus, trian, rectan Segnale2 gaussiana([-2,2],1000,5)

Chiamate conv_segn([-2,2],sinus(5,2,0,[-2,2],1000),gaussiana([-2,2],1000,5)) conv_segn([-2,2],rectan(5,2,0,[-2,2],1000),gaussiana([-2,2],1000,5)) conv_segn([-2,2],trian(5,2,0,[-2,2],1000),gaussiana([-2,2],1000,5))


55


56

3. Esempio di convoluzione tra i tre segnali periodici di base e il segnale esponenziale

Range 0 ÷ 4 Segnale1 sinus, trian, rectan Segnale2 espon([0,4],1000,1,1)

Chiamate conv_segn([0,4],sinus(5,2,0,[0,4],1000),espon([0,4],1000,1,1)) conv_segn([0,4],rectan(5,2,0,[0,4],1000),espon([0,4],1000,1,1)) conv_segn([0,4],trian(5,2,0,[0,4],1000),espon([0,4],1000,1,1))


57


58

4. Gli esempi successivi verranno effettuati considerando esclusivamente l’onda sinusoidale, questo perché dagli esempi precedenti si è potuto notare come il comportamento dell’operatore di convoluzione sia simile per tutti e tre i segnali periodici. Esempio di convoluzione dell’onda sinusoidale con un segnale rettangolare traslato verso destra.

Range -2 ÷ 2

Segnale1 sinus(5,2,0,[-2,2],1000) Segnale2 impulso(1,1000,600,100)

Chiamate conv_segn([-2,2],sinus(5,2,0,[-2,2],1000),impulso(1,1000,600,100))


59

5. Esempio di convoluzione dell’onda sinusoidale con un segnale rettangolare più ampio.

Range -2 ÷ 2 Segnale1 sinus(5,2,0,[-2,2],1000) Segnale2 impulso(1,1000,300,400)



60

6. Esempio di convoluzione dell’onda sinusoidale sfasata di 180° con un segnale rettangolare identico al precedente.

Range -2 ÷ 2

Segnale1 sinus(5,2,0,[-2,2],1000) Segnale2 impulso(1,1000,300,400)



61

7. Esempio di convoluzione dell’onda sinusoidale con un segnale gaussiano ampio.

Range -2 ÷ 2 Segnale1 sinus(5,2,0,[-2,2],1000) Segnale2 gaussiana([-2,2],1000,2)

Chiamate conv_segn([-2,2],sinus(5,2,0,[-2,2],1000),gaussiana([-2,2],1000,2))


62

8. Esempio di convoluzione dell’onda sinusoidale con un segnale gaussiano stretto.

Range -2 ÷ 2 Segnale1 sinus(5,2,0,[-2,2],1000) Segnale2 gaussiana([-2,2],1000,30)

Chiamate conv_segn([-2,2],sinus(5,2,0,[-2,2],1000),gaussiana([-2,2],1000,30))


63

9. Esempio di convoluzione dell’onda sinusoidale con un segnale esponenziale a discesa lenta.

Range 0 ÷ 4 Segnale1 sinus(5,2,0,[0,4],1000) Segnale2 espon([0,4],1000,1,0.5)

Chiamate conv_segn([0,4],sinus(5,2,0,[0,4],1000),espon([0,4],1000,1,0.5))


64

10. Esempio di convoluzione dell’onda sinusoidale con un segnale esponenziale a discesa veloce.

Range 0 ÷ 4

Segnale1 sinus(5,2,0,[0,4],1000) Segnale2 espon([0,4],1000,1,20)

Chiamate conv_segn([0,4],sinus(5,2,0,[0,4],1000),espon([0,4],1000,1,20))


65

3.5 Commenti Dai risultati ottenuti si può subito notare come le forme d’onda fornite in output dall’operatore di convoluzione siano piuttosto strane o comunque diverse rispetto a quelle viste negli esercizi precedenti. Il maggiore fattore di differenza riguarda senz’altro la non-periodicità del segnale di output, dovuta al fatto che la convoluzione avveniva tra un segnale periodico (onde quadre, triangolari, rettangolari) ed uno non periodico (esponenziale decrescente, impulso quadro e gaussiana). Passiamo ad analizzare singolarmente i risultati ottenuti, partendo dall’ultimo gruppo visto: la convoluzione tra il segnale sinusoidale e le forme d’onda impulsive. Il primo gruppo di risultati infatti evidenzia un comportamento molto simile rispetto al secondo gruppo. Iniziando dal risultato 5, cerchiamo di capire il significato della forma ottenuta simulando il comportamento dell’operatore di convoluzione nelle varie fasi della sua applicazione alle due forme d’onda in input. Per fare questo è necessario tener presente la definizione di convoluzione riportata nell’esercizio 2 – scheda 1. Possiamo identificare diverse fasi, corrispondenti a diversi comportamenti di input, simulando lo scorrimento da sinistra a destra del segnale impulsivo rispetto al segnale sinusoidale:

1. Prima fase: la moltiplicazione avviene durante lo scorrimento iniziale, che include la prima sezione del segnale rettangolare, quella con valore nullo (le aree evidenziata in rosso nell’immagine sottostante identificano le zone interessate rispettivamente del segnale impulsivo e poi del corrispondente output generato).


66

Il risultato della moltiplicazione dei due segnali è quindi nullo, come ovviamente l’integrale di convoluzione, fatto che rende nulla anche la prima parte del segnale in output (risultato 5).

2. Seconda fase: la moltiplicazione avviene nello scorrimento che include la sezione del segnale rettangolare con valore massimo, come visibile nell’immagine seguente.


67

Questo permette di ottenere un risultato non nullo dalla moltiplicazione dei due segnali e quindi anche dall’integrale di convoluzione, che corrisponde ad una forma d’onda che ripete il segnale sinusoidale iniziale di una lunghezza pari alla durata dell’ampiezza massima dell’onda rettangolare.

3. Terza fase: la moltiplicazione avviene nello scorrimento del segnale rettangolare attraverso tutto il segnale sinusoidale, come visibile nell’immagine seguente.


68

Il risultato che si ottiene dalla moltiplicazione dei due segnali non è altro che un’onda sinusoidale di durata pari al range dell’onda periodica a cui si sottrae la durata dell’impulso quadro.

4. Quarta fase: la moltiplicazione avviene nello scorrimento del segnale rettangolare nella fase finale dell’onda sinusoidale, in maniera opposta rispetto alla fase 2 e visibile nell’immagine seguente.


69

Il risultato è analogo a quanto visto per la fase 2, ma di segno opposto. 5. Quinta fase: la moltiplicazione avviene durante lo scorrimento finale, nel momento in cui il

segnale rettangolare ha valore nullo, come visibile nell’immagine seguente.


70

Il risultato è analogo a quanto visto per la fase iniziale. È stato analizzato il risultato 5 entrando minuziosamente nei dettagli, questo non si farà con i restanti risultati, in quanto si nota un comportamento paragonabile a quanto visto anche negli altri esempi riportati. Si osservino a riguardo i risultati 4 e 7. Il risultato 6 presenta invece un comportamento inverso rispetto a quello analizzato in dettaglio; infatti invertendo la fase dell’onda periodica si ottiene un inversione di segno nell’intero output. Il risultato 8 invece fa notare un’importante proprietà dell’operatore di convoluzione: si osserva infatti come l’output tra la convoluzione di un segnale periodico ed un segnale gaussiano (molto stresso, che approssima quindi l’impulso di dirac) sia il realtà il segnale periodico stesso, che viene riportato in posizione centrale nel risultato finale. Da questo si può dedurre come l’impulso di dirac sia in realtà l’elemento neutro per l’operatore di convoluzione. Lo stesso comportamento si nota anche nel risultato 10, in cui si osserva un picco iniziale molto breve (che simula ancora l’impulso di dirac) e che permette all’operatore di convoluzione di riportare il segnale periodico inalterato nel risultato.


71

Per quanto riguarda il risultato 9, sebbene cambi la forma del risultato fornito in output, le considerazioni fatte in precedenza sono ancora valide. La particolare forma del risultato è infatti da attribuirsi unicamente al segnale esponenziale utilizzato.


72


4.1 Traccia Si crei una sequenza di rumore uniforme in un intervallo [a,b]. Si verifichi che la sequenza creata è a distribuzione approssimativamente uniforme. Si ripetano analoghe valutazioni per una sequenza di rumore gaussiano in cui media e varianza siano parametri forniti dall’utente.

4.2 Svolgimento Per la creazione di una sequenza di rumore uniforme si è sfruttata la funzione di libreria Matlab rand che restituisce una sequenza di numeri pseudo-casuali con distribuzione nell’intervallo [0,1]. Effettuando una chiamata alla funzione rand, si otteneva immediatamente un vettore di dimensione data dall’utente valorizzato con numeri pseudo-casuali. I singoli valori sono poi stati moltiplicati per il valore massimo voluto di ampiezza, sempre fornito in input dall’utente. Per la verifica dell’uniformità della sequenza di rumore creata, si è calcolata la media della stessa attraverso la funzione di libreria mean. Analoga procedura è stata seguita per la realizzazione della seconda parte dell’esercizio. Si è utilizzata ora la funzione di libreria Matlab randn che genera una sequenza di numeri pseudo-casuali con distribuzione gaussiana. La chiamata a questa funzione restituiva un vettore di dimensione fornita dall’utente valorizzato con numeri pseudo-casuali. Questo vettore è stato poi moltiplicato per opportuni parametri, sempre forniti dall’utente, che permettessero di impostare la media e la varianza della distribuzione ottenuta. In particolare in corrispondenza di ogni elemento del vettore di numeri casuali sono state effettuate le seguenti operazioni matematiche:

• È stato sommato il valore di media fornito dall’utente • È stato moltiplicato per la deviazione standard, ovvero la radice quadrata della varianza,

anch’essa fornita in input dall’utente Per effettuare inoltre i controlli di correttezza dei risultati ottenuti sono stati confrontati i valori di media e varianza del vettore generato (utilizzando le funzioni di libreria mean e std rispettivamente) con quelli teorici forniti dall’utente.

4.2.1 Referenza alle primitive di Matlab utilizzate rand: restituisce un valore (o un vettore di valori) pseudo-casuale a distribuzione uniforme

nell’intervallo [0,1] randn: restituisce un valore (o un vettore di valori) pseudo-casuale a distribuzione gaussiana

nell’intervallo [0,1] sqrt(x): calcola la radice quadrata del parametro x


73

mean(x): calcola la media del vettore x std(x): calcola la deviazione standard del vettore x abs(x): calcola il valore assoluto del parametro x length(x): calcola la lunghezza del vettore x zeros(x,1): restituisce un vettore di dimensione x riempito con elementi pari a zero

4.3 Codice Matlab

4.3.1 Codice Rumore Uniforme % Funzione noise(range, max, ncomp) % Argomenti: % - range: intervallo di valori nella rappresentazione (formato [inizio,fine]) % - max: valore massimo della sequenza pseudo-casuale creata % - ncomp: numero di campioni considerati % % Esempio di utilizzo: % noise([1,3], 10, 1000) % % - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - % Autori: Stefano Maraspin & Stefano Valle % Data: 19/10/2004 % Versione: 1.0 function y = noise(range, max, ncomp) % Vengono ricavati i valori di ascissa iniziali e finali xmin = range(1,1); xmax = range(1,2); % Viene calcolato l'intervallo tra valori successivi delle ascisse passo = abs(xmax - xmin) / (ncomp - 1); % Viene generato un vettore contenente i punti sull'ascissa che saranno % considerati x = [xmin : passo : xmax]'; % Viene inizializzato un vettore di dimensione fornita in input counter = zeros(1,length(x)); % Chiamata della funzione rand ed estensione del valore del numero % pseudo-casuale restituito al valore massimo che il vettore deve avere y = (rand(ncomp,1)*max); % Calcolo della media del vettore appena creato media = mean(y) % Funzioni di stampa in una matrice 2 x 1 in cui la prima funzione rappresenta % la distribuzione uniforme creata e la seconda rappresenta l'istogramma % della distribuzione subplot(2,1,1) plot(x,y) subplot(2,1,2) hist(y,max)


74

4.3.1 Codice Rumore Gaussiano % Funzione noisegauss(range, media, varianza, ncomp) % Argomenti: % - range: intervallo di valori nella rappresentazione (formato [inizio,fine]) % - media: media della distribuzione gaussiana % - varianza: varianza della distribuzione gaussiana % - ncomp: numero di campioni considerati % % Esempio di utilizzo: % noisegauss([0,2], 10, 5, 1000) % % - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - % Autori: Stefano Maraspin & Stefano Valle % Data: 20/10/2004 % Versione: 1.0 function y = noisegauss(range, media, varianza, ncomp) % Vengono ricavati i valori di ascissa iniziali e finali xmin = range(1,1); xmax = range(1,2); % Viene calcolato l'intervallo tra valori successivi delle ascisse passo = abs(xmax - xmin) / (ncomp - 1); % Viene generato un vettore contenente i punti sull'ascissa che saranno % considerati x = [xmin : passo : xmax]'; % Generazione dei vettori di numeri pseudocasuali con distribuzione gaussiana y = media + randn(ncomp,1) * sqrt(varianza); % calcolo di media e varianza per il confronto con i valori forniti dall'utente media = mean(y) var = std(y)*std(y) % Funzioni di stampa in una matrice 2 x 1 in cui la prima funzione rappresenta % la distribuzione gaussiana creata e la seconda rappresenta l'istogramma % della distribuzione subplot(2,1,1) plot(x,y) subplot(2,1,2) hist(y)



75

4.4.1 Rumore a distribuzione uniforme 11. Esempio di rumore a distribuzione uniforme

Range 0 ÷ 2

Massimo 10 N. Campioni 1000

Chiamata noise([0,2],10,1000)

Risultati:

Teorica Calcolata Media 5 5,0293


76

2. Esempio di rumore a distribuzione uniforme creata considerando un numero di campioni ridotto.

Range 0 ÷ 1 Massimo 10

N. Campioni 50


Risultati:



77

3. Esempio di rumore a distribuzione uniforme creata considerando un numero di campioni molto alto.


N. Campioni 10000


Risultati:



78

4. Esempio di rumore a distribuzione uniforme creata considerando un numero di campioni medio ed un valore massimo ridotto rispetto al caso precedente.


N. Campioni 1000


Risultati:

Teorica Calcolata Media 2,5 2,4704


79

5. Esempio di rumore a distribuzione uniforme creata considerando un numero di campioni medio ed un valore massimo molto più alto rispetto al caso precedente.


N. Campioni 1000


Risultati:



80

4.4.1 Rumore a distribuzione gaussiana 6. Esempio di rumore a distribuzione gaussiana di valore medio nullo e varianza ridotta.

Range 0 ÷ 2 Media 0

Varianza 1 N. Campioni 1000

Chiamata noisegauss([0,2],0,1,1000)

Risultati:


Varianza 1 0,9470


81

7. Esempio di rumore a distribuzione gaussiana di valore medio positivo e varianza ridotta.




Risultati:


Varianza 1 1,0338


82

8. Esempio di rumore a distribuzione gaussiana di valore medio negativo e varianza ridotta.

Range 0 ÷ 2 Media -5


Chiamata noisegauss([0,2],-5,1,1000)

Risultati:

Teorica Calcolata Media -5

Varianza 1


83

9. Esempio di rumore a distribuzione gaussiana di valore medio nullo e varianza aumentata rispetto ai casi precedenti.




Risultati:


Varianza 10 10,9486


84

10. Esempio di rumore a distribuzione gaussiana di valore medio nullo e varianza ancora più grande rispetto ai casi precedenti.




Risultati:


Varianza 1000 1055


85

11. Esempio di rumore a distribuzione gaussiana di valore medio nullo e varianza come al punto (3), ma con un numero di campioni considerevolmente più alto.




Risultati:


Varianza 100 99,9810


86

4.5 Commenti Osservando i risultati ottenuti si può sicuramente affermare come la funzione di librearia rand utilizzata nella prima parte dell’esercizio permetta di creare effettivamente una sequenza di numeri pseudo-casuali a distribuzione uniforme. Questo è riscontrabile sia visivamente, osservando le immagini riportate, sia analiticamente confrontando i valori teorici con quelli reali, calcolati sul singolo esempio riportato. Analoga considerazione va fatta per la funzione di libreria randn la quale permette di creare numeri pseudo-casuali con distribuzione gaussiana molto buoni. Gli esempi riportati sono stati selezionati cercando di evidenziare i comportamenti significativi delle due funzioni create, in relazione a determinate classi di parametri forniti in input. Entrando nei dettagli dei singoli esempi (primi 5 riguardanti la distribuzione uniforme e successivi riguardanti la distribuzione gaussiana), è possibile effettuare alcune osservazioni, valide in modo eguale per entrambe le sezioni di questo esercizio:

• Le due funzioni utilizzare restituiscono valori casuali molto buoni, di media e varianza pressoché identiche ai valori teorici, la cui precisione dipende comunque dal numero di campioni e dal valore massimo concesso agli stessi.

• Riducendo il numero di campioni si osserva che, com’è facile dedurre, l’istogramma generato (risultato 2) è piuttosto discontinuo.

• Aumentando il numero di campioni invece (risultati 3 e 11) si osserva come l’istogramma generato sia piuttosto piatto nel caso di distribuzione uniforme ed invece presenti perfettamente la tradizionale forma “a campana” nel caso della distribuzione gaussiana. Questo risultato, assieme al precedente, evidenzia come la bontà della distribuzione ottenuta dipenda dal numero di campioni considerato. In realtà questo è un comportamento corretto, infatti ipoteticamente solo con un numero di campioni infinito riusciremmo ad ottenere distribuzioni che presentano un grafico ideale (piatto nel primo caso, perfettamente “a campana” nel secondo).

• Riducendo inoltre il valore massimo della distruzione (dist. uniforme – risultato 4), si ritrova a parità di numero di campioni una distribuzione più precisa.

Osservazioni riguardanti prettamente la distribuzione gaussiana: • Variando positivamente la media della distribuzione (risultato 7), si osserva come

l’istogramma rappresentante i risultati ottenuti si sposti correttamente verso destra, concentrando il centro della distribuzione sul valore positivo impostato dai dati in ingresso.

• Variando negativamente la media della distribuzione (risultato 8), si osserva un risultato opposto rispetto al caso precedente, ovvero lo spostamento verso sinistra dell’istogramma ed un posizionamento della parte centrale della distribuzione in corrispondenza del valore medio fornito.

• Modificando il valore di varianza dato in input si osserva invece un progressivo allargamento della tipica forma “a campana” della distribuzione. É possibile verificare questo comportamento osservando i valori di ascissa presenti nei risultati 1, 9, 10.


87


5.1 Traccia Si creino e visualizzino forme d’onda a cui venga sommato rumore gaussiano di media e varianza fornite dall’utente. Si calcolino autocorrelazione, correlazione e convoluzione e se ne visualizzino i risultati.

5.2 Svolgimento Le funzioni necessarie allo svolgimento erano in parte già state realizzate nelgli esercizi precedenti. Inoltre è stata creata una variante della funzione noisegauss, denominata noisegauss2, che ha il solo compito di restituire un vettore di valori pseudo-casuali, senza l’utilizzo di funzioni di supporto per la visualizzazione. Si è quindi creata una funzione modulare noise_per per effettuare le seguenti operazioni:

• Creare un vettore contenente un segnale periodico a cui venisse sommato del rumore gaussiano.

• Applicare a questo vettore e al segnale periodico pulito le operazioni di autocorrelazione, correlazione e convoluzione.

• Effettuare la stampa a schermo dei risultati ottenuti al fine di facilitare il confronto degli stessi

5.2.1 Referenza alle primitive di Matlab utilizzate Chiamate di primitive Matlab utilizzate: randn: restituisce un valore (o un vettore di valori) pseudo-casuale a distribuzione gaussiana

nell’intervallo [0,1] sqrt(x): calcola la radice quadrata del parametro x abs(x): calcola il valore assoluto del parametro x conv(x,y): calcola la convoluzione tra il segnale x e il segnale y xcorr(x): calcola l’autocorrelazione del segnale x xcorr(x,y): calcola la correlazione tra il segnale x e il segnale y



88

5.3 Codice Matlab

5.3.1 Codice rumore gaussiano (variante) % Funzione noisegauss2(media, varianza, ncomp)

% Argomenti:

% - media: media della distribuzione gaussiana

% - varianza: varianza della distribuzione gaussiana


%


% noisegauss2(10, 5, 1000)

%

% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -


% Data: 21/10/2004

% Versione: 1.0

function y = noisegauss(media, varianza, ncomp)

% Generazione dei vettori di numeri pseudocasuali con distribuzione gaussiana

y = media + randn(ncomp,1) * sqrt(varianza);

5.3.2 Codice funzione principale % Funzione noise_per(segnale, range, ncomp, media, varianza)

% Argomenti:

% - segnale: vettore contenente la rappresentazione del segnale periodico



% - media: media rumore gaussiano

% - varianza: varianza

%


% noise_per(sinus (5, 2, 0, [1,3], 1000), [1,3], 1000, 0, 2)

%

% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -


% Data: 21/10/2004

% Versione: 1.0

function y = noise_per(segnale, range, ncomp, media, varianza)

% Vengono ricavati i valori di ascissa iniziali e finali per rappresentare

% il segnale

xmin_segnale = range(1,1);

xmax_segnale = range(1,2);


89

% Vengono ricavati i valori di ascissa iniziali e finali per rappresentare

% correlazione e convoluzione

xmin_corrconv = - abs(xmin_segnale - xmax_segnale);

xmax_corrconv = abs(xmin_segnale - xmax_segnale);

% Viene calcolato l'intervallo tra valori successivi delle ascisse per

% la rappresentazione del segnale

passo_segnale = abs(xmax_segnale - xmin_segnale) / (ncomp - 1);

% Viene calcolato l'intervallo tra valori successivi delle ascisse per

% la rappresentazione di correlazione e convoluzione

passo_corrconv = abs(xmax_corrconv - xmin_corrconv) / (2*ncomp - 2);


% considerati per la rappresentazione del segnale

x_segnale = [xmin_segnale : passo_segnale : xmax_segnale]';


% considerati per la rappresentazione di correlazione e convoluzione

x_corrconv = [xmin_corrconv : passo_corrconv : xmax_corrconv]';

% Generazione del vettore di rumore gaussiano

rumore = noisegauss2(media, varianza, ncomp);

% Generazione del vettore contenente il segnale a cui viene sommato

% rumore gaussiano

y = segnale + rumore;

% Calcolo dell'autocorrelazione del segnale con rumore

autocorr = xcorr(y) ./ ncomp;

% Calcolo della correlazione fra segnale e segnale con rumore

correlazione = xcorr(segnale,y) ./ ncomp;

% Calcolo della convoluzione fra segnale e segnale con rumore

convoluzione = conv(segnale,y) ./ ncomp;

% Funzioni di stampa

subplot(4,2,1)

plot(x_segnale,segnale)

title('Segnale')

subplot(4,2,2)

plot(x_segnale,y)

title('Segnale con rumore')

subplot(4,2,[3,4])


90

plot(x_corrconv,autocorr)

title('Autocorrelazione')

subplot(4,2,[5,6])

plot(x_corrconv,correlazione)

title('Correlazione')

subplot(4,2,[7,8])

plot(x_corrconv,convoluzione)

title('Convoluzione')



91

12. Segnale sinusoidale

Segnale sinus(5, 2, 0, [0,2], 1000) Range 0 ÷ 2

N. Campioni 1000 Media 0

Varianza 1

Chiamata noise_per(sinus(5, 2, 0, [0,2], 1000), [0,2], 1000, 0, 1)


92

13. Segnale rettangolare

Segnale rectan(5, 2, 0, [0,2], 1000) Range 0 ÷ 2


Varianza 1

Chiamata noise_per(rectan(5, 2, 0, [0,2], 1000), [0,2], 1000, 0, 1)


93

14. Segnale triangolare

Segnale trian(5, 2, 0, [0,2], 1000) Range 0 ÷ 2


Varianza 1

Chiamata noise_per(trian(5, 2, 0, [0,2], 1000), [0,2], 1000, 0, 1)


94

15. Segnale sinusoidale a cui si somma rumore gaussiano con varianza maggiore rispetto ai casi precedenti



Varianza 20



95

16. Segnale rettangolare a cui si somma rumore gaussiano con varianza maggiore rispetto ai casi precedenti



Varianza 20



96

17. Segnale triangolare a cui si somma rumore gaussiano con varianza maggiore rispetto ai casi precedenti



Varianza 20



97

18. Segnale sinusoidale a cui si somma rumore gaussiano con media non nulla



Varianza 10



98

19. Segnale rettangolare a cui si somma rumore gaussiano con media non nulla



Varianza 10



99

20. Segnale triangolare a cui si somma rumore gaussiano con media non nulla



Varianza 10



100

5.5 Commenti La presentazione dei risultati è stata suddivisa nel paragrafo precedente in tre gruppi, in base alla media e alla varianza del rumore che veniva sommato al segnale periodico. Si utilizzerà ora la stessa suddivisione per la presentazione delle osservazioni. Nel primo gruppo di risultati, ovvero quello in cui il rumore che viene sommato al segnale periodico è piuttosto ridotto, si può notare come i grafici dell'autocorrelazione del segnale con rumore e della correlazione tra segnale con rumore e segnale pulito siano piuttosto simili. Questo è dovuto indubbiamente al fatto che la quantità di rumore è così bassa da non modificare il comportamento dell'operatore di correlazione in maniera significativa. L'unica discontinuità rispetto ai due grafici si può osservare nel punto centrale, quello in cui τ = 0. Questa discontinuità è ancor più visibile nel secondo gruppo di risultati. Quando infatti la quantità di rumore aumenta (a parità di ampiezza del segnale) si osserva che la correlazione massima si ha nel punto in cui τ = 0 (come ci si aspetta evidentemente), invece spostandosi seppur di poco da questa posizione, la correlazione risulta molto ridotta. Il rumore infatti è così alto che un movimento lieve nella posizione del segnale fa diminuire notevolmente il valore calcolato dall'integrale di correlazione, che tende quindi a diminuire più velocemente. Nei punti più distanti da τ = 0 invece la forma del segnale di autocorrelazione è paragonabile al caso iniziale. In entrambi questi gruppi di risultati la convoluzione tra i due segnali invece non presenta particolari di rilievo. Diverso è il comportamento nel terzo gruppo di risultati, ovvero quelli in cui il valore medio del segnale con rumore è non nullo. Il grafico rappresentante l'autocorrelazione infatti ha una forma piuttosto diversa rispetto ai casi fin'ora visti: si ha un picco centrale molto elevato, che diminuisce in maniera analoga per valori di τ crescenti in valore assoluto. Questo comportamento è dovuto al fatto che nel momento in cui il grafico del segnale con rumore è presente con ogni suo valore nel calcolo dell'autocorrelazione si ha un valore molto alto (in quanto il rumore risulta preponderante quantitativamente rispetto al solo segnale periodico), quando invece il calcolo inizia a traslarlo considerando sempre meno valori di esso, si ha una progressiva discesa del valore dell'autocorrelazione. Questo esercizio ha permesso anche di evidenziare un utilizzo concreto dell'operatore di autocorrelazione: osservando infatti i risultati 4, 5 e 6 si nota come, nonostante il segnale con rumore appaia notevolmente compromesso, si riesca comunque ad evidenziare una periodicità dello stesso osservando la sua autocorrelazione. Tale periodicità è sì visibile guardando il segnale, ma non in modo così marcato come invece viene evidenziato dall’operatore.


101


102

Scheda 2


103


6.1 Traccia Eseguire e visualizzare la Trasformata di Fourier discreta (DFT) di un impulso rettangolare (boxcar) al variare della lunghezza T dell’impulso stesso. Ricavare e visualizzare anche gli spettri di potenza e fase.

6.2 Svolgimento Il codice per la rappresentazione di un impulso rettangolare (boxcar) è già stato definito precedentemente, quindi ci si è potuti concentrare da subito sulla creazione di una funzione generica che permettesse di calcolare la Trasformata di Fourier discreta (DFT) di un segnale fornito in input. Tale funzione, denominata nel seguito dft_segnale, utilizza chiamate successive alle funzioni di libreria Matlab fft e fftshift per calcolare la DFT del segnale fornito. La funzione fftshift permette in particolare di rappresentare al centro dello spetto la componente di frequenza nulla. La DFT così calcolata presenta entrambe le componenti reali ed immaginarie, per ricavare i vari spettri, si è svolto nel modo seguente:

• spettro di ampiezza: si è calcolato il modulo (utilizzando la funzione abs) della DFT; • spettro di fase: si è utilizzata la funzione di libreria angle, applicata alla DFT del segnale • spettro di potenza: si è calcolato il quadrato dello spettro di ampiezza.

6.2.1 Referenza alle primitive di Matlab utilizzate Chiamate di primitive Matlab utilizzate: abs(x): calcola il valore assoluto del parametro x max(x): calcola il valore massimo del vettore x fft(x): calcola la DFT del vettore x fftshift(x): posiziona al centro dello spettro il risultato della fft angle(x): calcola l’angolo in radianti della DFT di un segnale

6.2.2 Chiamate a funzioni da noi definite in precedenza Chiamate a funzioni esterne utilizzate: impulso: rif. Es 3 interv: rif. Es 3


104

6.3 Codice Matlab % Funzione dft_segnale(segnale, range, ncomp)

% Argomenti:

% - segnale: segnale di cui calcolare la DFT



%


% dft_segnale(impulso(3,1000,100,100), [0,2], 1000)

%

% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -


% Data: 4/11/2004

% Versione: 1.1

function y = dft_segnale(segnale, range, ncomp)

% Vengono definiti i valori dell'asse delle ascisse mediante una chiamata

% alla funzione interv


% Calcolo della frequenza di campionamento del segnale

fcamp = ncomp / (abs(range(1,1)-range(1,2)))

% Viene fissata la scala delle ascisse nelle rappresentazioni degli spettri,

% nel dominio delle frequenze

xs = [-fcamp/2:fcamp/(ncomp-1):fcamp/2];

% Viene calcolata la DFT del segnale fornito in input

fft_segnale = fftshift(fft(segnale)) / length(xs);

% Viene calcolato lo spettro di ampiezza del segnale

ampiezza_fft = abs(fft_segnale);

% Viene calcolato lo spettro di potenza del segnale

potenza_fft = ampiezza_fft .* ampiezza_fft;

% Viene calcolato lo spettro di fase del segnale

fase_fft = angle(fft_segnale);

% Calcolo della frequenza di Nyquist

nyquist = fcamp/2;

% Procedure di stampa a video dei risultati

subplot(4,2,[1,2])

plot(x,segnale)


105

axis([-range(1,1),range(1,2),min(segnale)-

(10*abs(max(segnale))/100),max(segnale)+(10*abs(max(segnale)/100))])

title('Impulso')

grid on

subplot(4,2,[3,4])

plot(xs,fft_segnale)

axis([-nyquist,nyquist,min(fft_segnale)-

(10*abs(max(fft_segnale))/100),max(fft_segnale)+(10*abs(max(fft_segnale)/1

00))])

title('Spettro completo del segnale')

grid on

subplot(4,2,5)

plot(xs,ampiezza_fft)

title('Spettro di ampiezza')

axis([-nyquist,nyquist,min(ampiezza_fft)-

(10*abs(max(ampiezza_fft))/100),max(ampiezza_fft)+(10*abs(max(ampiezza_fft

)/100))])

grid on

subplot(4,2,6)

plot(xs,fase_fft)

title('Spettro di fase')

axis([-nyquist,nyquist,min(fase_fft)-

(10*abs(max(fase_fft))/100),max(fase_fft)+(10*abs(max(fase_fft)/100))])

grid on

subplot(4,2,[7,8])

plot(xs,potenza_fft)

title('Spettro di potenza')

axis([-nyquist,nyquist,min(potenza_fft)-

(10*abs(max(potenza_fft))/100),max(potenza_fft)+(10*abs(max(potenza_fft)/1

00))])

grid on



106

1. Esempio di Trasformata di Fourier di un segnale boxcar.

Segnale impulso(5,1000,10,50) Range 0 ÷ 2

N. Campioni 1000

Chiamata dft_segnale(impulso(5,1000,10,50),[0,2],1000)


107

2. DFT di un boxcar di lunghezza molto ridotta rispetto al caso precedente.


N. Campioni 1000



108

3. DFT di un boxcar di lunghezza ridotta ulteriormente, tale da approssimare il segnale di dirac.


N. Campioni 1000



109

4. DFT di un boxcar in cui viene aumentata la lunghezza dell’impulso stesso.


N. Campioni 1000



110

5. DFT di un boxcar identico a quello riportato nel risultato (2), ma traslato in avanti.


N. Campioni 1000



111

6. DFT di un boxcar identico a quello riportato nel risultato (2), ma con un numero di campioni ridotto.


N. Campioni 300



112

6.5 Commenti Come si può notare dai risultati riportati nel paragrafo precedente, la trasformata di Fourier di un impulso quadro (boxcar) è caratterizzata da una forma che si ripete in tutti i casi visti: vi è un picco alla basse frequenze molto ampio, nel cui intorno sono presenti vari picchi sia d’intensità che di frequenza decrescente. Si è cercato di capire meglio la struttura di questo spettro andando a variare tutti i parametri dell’impulso quadro. Riportiamo di seguito le osservazioni emerse:

• Diminuzione della lunghezza dell’impulso (2) L’armonica in bassa frequenza subisce un calo di ampiezza e diventano più consistenti le armoniche man mano che la frequenza aumenta. Il comportamento è corretto in quanto il restringimento dell’impulso comporta un effetto sali-scendi rapido che in termini di frequenza nella DFT viene tradotto in un aumento dell’importanza delle frequenze più alte, questo perchè la discontinuità del segnale è aumentata.

• Aumento della lunghezza dell’impulso (4)

Il comportamento è contrapposto al caso precedente (come le motivazioni): le armoniche si restringono attorno alle basse frequenze ed infatti aumenta in modo consistente l’ampiezza dell’armonica fondamentale.

• Approssimazione dell’impulso di dirac (3)

Comportamento diverso rispetto agli altri casi presenta questa situazione, in cui si è cercato di approssimare l’impulso di dirac. Per rappresentare la DFT della delta di dirac (anche se quella da noi rappresentata ne è solo un’approssimazione) infatti sono necessarie tutte le frequenze dello spettro con ampiezze che decrescono all’aumentare della frequenza.

• Spostamento dell’impulso (5)

Lo spostamento dell’impulso nel grafico temporale, non genera nessun effetto per quanto riguarda gli spettri di ampiezza e fase, rispetto al caso 2 (si noti infatti l’analogia tra gli spettri nei due casi). L’unica variazione che lo spostamento implica, riguarda lo spettro di fase, come ci si poteva attendere.

• Riduzione del numero di campioni del segnale (6)

Riducendo il numero di campioni del segnale, questo non comporta variazioni nei vari spettri (si noti l’analogia con il caso 2). L’unico fattore che viene modificato riguarda la frequenza di campionamento che subisce una riduzione. Tale riduzione è proporzionale alla


113

riduzione che subisce il numero di campioni, infatti la frequenza di campionamento viene calcolata attraverso la formula seguente:

ervallοn

f campcamp int

=


114


7.1 Traccia Eseguire la DFT di un segnale periodico al variare della frequenza; ricavare anche lo spettro di potenza e fase. Eseguire la DFT della somma di due segnali sinusoidali a frequenza variabile e relativi spettri. Eseguire la DFT di un segnale periodico di frequenza variabile a cui è stato sommato del rumore gaussiano a media nulla e varianza variabile.

7.2 Svolgimento Lo svolgimento dell’esercizio è stato suddiviso in tre fasi distinte, che andremo ad analizzare separatamente. 1. DFT di un segnale periodico al variare della frequenza

Si è utilizzata una funzione definita in precedenza (dft_segnale) che permettesse il calcolo e la rappresentazione dei vari spettri richiesti.

2. DFT della somma di due segnali sinusoidali a frequenza variabile Si è creata la funzione somma_segnali, che permettesse di sommare due segnali distinti (a parità di numero di campioni e intervallo). L’output di tale funzione è stato utilizzato come parametro per la procedura dft_segnale creata in precedenza.

3. DFT di un segnale periodico a cui è stato sommato rumore gaussiano Utilizzando una funzione definita in precedenza (noisegauss2) che restituisse un vettore rappresentante del rumore gaussiano con media e varianza definibili dall’utente, si è creata la sequenza di rumore voluta. Tale sequenza e un vettore rappresentante il segnale periodico da analizzare sono stati sommati utilizzando la procedura creata al punto precedente. L’output di tale funzione poi è stato utilizzato come parametro per la procedura dft_segnale.

7.2.1 Chiamate a funzioni da noi definite in precedenza Chiamate a funzioni esterne utilizzate: dft_segnale: rif. Es 2 – Scheda 2 noisegauss2: rif. Es 5 – Scheda 1


115

7.3 Codice Matlab

7.3.1 Codice funzione di somma dei segnali % Funzione somma_segnali(segnale1, segnale2)

% Argomenti:

% - segnale1: primo segnale da sommare

% - segnale2: secondo segnale da sommare

%


% somma_segnali(sinus(5,10,0,[0,2],1000),sinus(5,30,0,[0,2],1000))

%

% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -


% Data: 4/11/2004

% Versione: 1.0

function y = somma_segnali(segnale1,segnale2)

% Calcolo della somma dei due segnali

y = segnale1 + segnale2;



116

7.4.1 DFT di un segnale periodico 1. DFT di un segnale sinusoidale

Segnale sinus(5,50,0,[0,2],1000) Range 0 ÷ 2

N. Campioni 1000

Chiamata dft_segnale(sinus(5,50,0,[0,2],1000),[0,2],1000)

campf 500


117

2. DFT di un segnale sinusoidale di frequenza maggiore rispetto al precedente.


N. Campioni 1000


campf 500


118

3. DFT di un segnale sinusoidale oltre la frequenza di Nyquist.


N. Campioni 1000


campf 500


119

4. DFT di un segnale rettangolare.

Segnale rectan(5,10,0,[0,2],1000) Range 0 ÷ 2

N. Campioni 1000

Chiamata dft_segnale(rectan(5,10,0,[0,2],1000),[0,2],1000)

campf 500


120

5. DFT di un segnale rettangolare di frequenza maggiore rispetto al precedente.

Segnale rectan(5,50,0,[0,2],1000) Range 0 ÷ 2

N. Campioni 1000

Chiamata dft_segnale(rectan(5,50,0,[0,2],1000),[0,2],1000)

campf 500


121

6. DFT di un segnale triangolare.

Segnale trian(5,10,0,[0,2],1000) Range 0 ÷ 2

N. Campioni 1000

Chiamata dft_segnale(trian(5,10,0,[0,2],1000),[0,2],1000)

campf 500


122

7. DFT di un segnale triangolare di frequenza maggiore rispetto al precedente.

Segnale trian(5,50,0,[0,2],1000) Range 0 ÷ 2

N. Campioni 1000

Chiamata dft_segnale(trian(5,50,0,[0,2],1000),[0,2],1000)

campf 500


123

7.4.2 DFT della somma di due segnali sinusoidali 1. DFT della somma di due segnali sinusoidali con frequenze diverse l’uno dall’altro.

Segnale 1 sinus(5,10,0,[0,2],200) Segnale 2 sinus(5,20,0,[0,2],200)


campf 100


124

2. DFT della somma di due segnali sinusoidali con frequenze diverse rispetto al caso precedente (una delle due al limite dello spetto di Nyquist).

Segnale 1 sinus(5,20,0,[0,2],200) Segnale 2 sinus(5,50,0,[0,2],200)


campf 100


125

3. DFT della somma di due segnali sinusoidali identici a quelli utilizzati nel punto 1, di cui uno sfasato di 180°.

Segnale 1 sinus(5,10,0,[0,2],200) Segnale 2 sinus(5,20,pi,[0,2],200)


campf 100


126

7.4.3 DFT della somma di un segnale periodico con rumore gaussiano 1. DFT della somma di un segnale sinusoidale con rumore gaussiano a media nulla.

Segnale 1 sinus(5,10,0,[0,2],200) Segnale 2 noisegauss2(0,l,200)


campf 100


127

2. DFT della somma di un segnale sinusoidale con rumore gaussiano a media nulla e varianza maggiore rispetto al caso precedente.

Segnale 1 sinus(5,10,0,[0,2],200) Segnale 2 noisegauss2(0,l0,200)


campf 100


128

3. DFT della somma di un segnale rettangolare con rumore gaussiano a media nulla. Segnale 1 rectan(5,10,0,[0,2],200) Segnale 2 noisegauss2(0,l,200)


campf 100


129

4. DFT della somma di un segnale rettangolare con rumore gaussiano a media nulla e varianza maggiore rispetto al caso precedente.

Segnale 1 rectan(5,10,0,[0,2],200) Segnale 2 noisegauss2(0,l0,200)


campf 100


130

5. DFT della somma di un segnale triangolare con rumore gaussiano a media nulla. Segnale 1 trian(5,10,0,[0,2],200) Segnale 2 noisegauss2(0,l,200)


campf 100


131

6. DFT della somma di un segnale triangolare con rumore gaussiano a media nulla e varianza maggiore rispetto al caso precedente.

Segnale 1 trian(5,10,0,[0,2],200) Segnale 2 noisegauss2(0,l0,200)


campf 100


132

7.5 Commenti La DFT di un segnale periodico ha evidenziato un comportamento distinto a seconda del tipo di segnale che si andava ad analizzare:

• Segnale sinusoidale: com’era possibile prevedere la DFT di un segnale sinusoidale corrisponde nello spettro ad una sola armonica, di frequenza corrispondente a quella del segnale sinusoidale stesso.

• Segnale rettangolare: la DFT corrisponde nello spettro ad una serie (infinita) di armoniche che riempiono l’intero grafico a partire dalle frequenza del segnale stesso. È ben visibile la situazione nello spettro di ampiezza (oltre che in quello di potenza, anche se in modo minore): l’armonica di base è quella di ampiezza più elevata, man mano che ci si allontana da essa, l’ampiezza diminuisce gradualmente. Si è potuto notare quindi come sia la somma di un numero infinito di segnali sinusoidali che permette di creare un segnale rettangolare.

• Segnale triangolare: la DFT visibile nello spettro di ampiezza è simile al caso precedente, ma nello spettro completo (in cui per la visualizzazione concorrono in qualche modo anche le componenti immaginarie) vediamo una serie di armoniche che alternano valori di ampiezza positivi e negativi, dovuti a valori di fase contrapposti al caso precedente, situazione che poi viene chiaramente appiattita negli spettri di ampiezza e potenza (in quest’ultimo è però talmente elevata la prima armonica da non permettere alle successive di essere visibili correttamente).

Si può quindi identificare un “pattern” di comportamento distinto che ogni segnale presenta nella rappresentazione spettrale. La DFT di due segnali sinusoidali sommati tra loro, corrisponde alla somma di due spettri di segnali sinusoidali presi singolarmente: si ha quindi la presenza di due armoniche posizionate alle rispettive frequenze. Nella visualizzazione dello spettro della DFT di un segnale periodico a cui viene sommato del rumore gaussiano, si nota invece come lo spettro visto inizialmente venga disturbato da tutta una serie di armoniche di entità molto lieve, ma che si posizionano in tutto lo spettro stesso. Tali armoniche costituiscono la rappresentazione spettrale del rumore. È interessante notare come, in condizioni in cui il rumore sia talmente elevato da non rendere più distinguibile il segnale originale (risultati 2, 4, 6), osservando lo spettro di ampiezza si vedano invece predominare le frequenze corrispondenti al segnale periodico di partenza. Questo comportamento è utile per evidenziare delle periodicità all’interno del segnale rumoroso. Inoltre in alcuni casi si riescono anche a ritrovare i pattern propri dei segnali periodici analizzati (esempio, risultato 4 – spettro di ampiezza). Osservazione finale riguarda la composizione stessa dello spettro rappresentante la DFT. In ogni risultato riportato si osserva come in presenza di una sola armonica, si abbia una visualizzazione


133

sotto forma di 2 valori distinti nello spettro. In realtà ciò che è veramente importante e significativo ai fini della rappresentazione nel dominio delle frequenze è la sola parte destra dello spettro. Quello che compare a sinistra invece è ridondante in quanto corrisponde a componenti complesse che non aggiungono nessuna informazione utile alla rappresentazione.


134


8.1 Traccia Si esegua la convoluzione nel dominio dei tempi tra due segnali periodici di frequenza variabile. Si confronti il risultato con quello ottenuto moltiplicando termine a termine i rispettivi spettri.

8.2 Svolgimento Per la risoluzione dell'esercizio è stata in primo luogo considerata la generazione di due segnali,

sulla base dei parametri forniti dall'utente (tipo di segnale, ampiezza, frequenza, fase). A tale scopo sono state utilizzate le funzioni precedentemente definite per la costruzioni di segnali

periodici (sinus, trian, rectan) e qui duplicemente chiamate.

Successivamente sono state utilizzate le funzioni matlab fft, fftshift e abs, applicate ad entrambi i segnali generati in precedenza, così da ottenere i due rispettivi spettri.

È stata fatta in seguito una convoluzione tra i due segnali e si è ricavato lo spettro della stessa, per

mezzo delle funzioni matlab conv prima e fft, fftshift e abs poi, come visto in precedenza, sui segnali originali.

Si è così ottenuto il primo termine di paragone richiesto dall'esercizio, ovvero lo spettro della

convoluzione tra i due segnali. Per ottenere il secondo termine di paragone, come suggerito dalla traccia fornita, è bastato moltiplicare termine a termine i valori degli spettri nell'intervallo definito.

8.2.1 Referenza alle primitive di Matlab utilizzate Chiamate di primitive Matlab utilizzate: abs(x): calcola il valore assoluto del parametro x fft(x): calcola la DFT del vettore x fftshift(x): posiziona al centro dello spettro il risultato della fft conv(x,y): esegue una convoluzione tra i segnali x e y.

8.2.2 Chiamate a funzioni da noi definite in precedenza Chiamate a funzioni utilizzate: sinus: rif. Es 1 – Scheda 1 trian: rif. Es 1 – Scheda 1 rectan: rif. Es 1 – Scheda 1


135

8.3 Codice Matlab % Funzione es23(tipo, amp, freq, fase, range, ncomp) % % Esegue la convoluzione tra due segnali e ne deriva la FFT. % Vengono anche ricavati gli spettri dei due segnali e moltiplicati % tra loro. I risultati vengono poi confrontati. % % Argomenti: % - tipo: tipo di segnale da considerare per il confronto % - amp: ampiezza dei segnali periodici % - freq: frequenza dei segnali periodici % - fase: fase dei due segnali % - range: intervallo da considerare (nel dominio spaziale) % - ncomp: numero di campioni considerati % % Esempio di utilizzo: % es23(['s','s'], [3,2], [2,4],[0,0], [-4,4], 4000) % % - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - % Autori: Stefano Maraspin & Stefano Valle % Data: 11/11/2004 % Versione: 1.0 % % Rev. % 02/01/2005 - Normalizzazione degli output prodotti da FFT function y = es23(tipo, amp, freq, fase, range, ncomp) % Migliorata la precisione di macchina format long % Qui generiamo diversi tipi di segnali periodici, in base alla % preferenza dell'utente % Primo Segnale if (tipo(1,1)=='t') valori_onda1 = trian(amp(1,1), freq(1,1), fase(1,1), range, ncomp); elseif (tipo(1,1)=='r') valori_onda1 = rectan(amp(1,1), freq(1,1), fase(1,1), range, ncomp); else % Default (seno) valori_onda1 = sinus(amp(1,1), freq(1,1), fase(1,1), range, ncomp); end % 0 Padding al segnale temp_onda1 = valori_onda1; temp_onda1(ncomp*2-1) = 0; % Secondo Segnale if (tipo(1,2)=='t') valori_onda2 = trian(amp(1,2), freq(1,2), fase(1,2), range, ncomp); elseif (tipo(1,2)=='r') valori_onda2 = rectan(amp(1,2), freq(1,2), fase(1,2), range, ncomp); else


136

valori_onda2 = sinus(amp(1,2), freq(1,2), fase(1,2), range, ncomp); end % 0 Padding temp_onda2 = valori_onda2; temp_onda2(ncomp*2-1) = 0; % FFT dei due segnali e successiva normalizzazione spettro_onda1 = (fftshift(abs(fft(temp_onda1)))); spettro_onda2 = (fftshift(abs(fft(temp_onda2)))); spettro_onda1 = (spettro_onda1*2)/length(spettro_onda1); spettro_onda2 = (spettro_onda2*2)/length(spettro_onda2); % Moltiplicazione dei due spettri spettro_m = (spettro_onda1 .* spettro_onda2); % Convoluzione tra i due segnali (prima dello 0 padding) convoluzione = double(conv(valori_onda1,valori_onda2))/(ncomp); % Viene ricavato lo spettro della convoluzione e successivamente % normalizzato spettro_conv = abs(fftshift(fft(convoluzione))); spettro_conv = (spettro_conv*2)/length(spettro_conv); % Ricavati i valori in ascissa xmin = range(1,1); xmax = range(1,2); % Viene calcolato l'intervallo tra valori successivi delle ascisse passo = abs(xmax - xmin) / ((ncomp*2-1) - 1); % Viene generato un vettore contenente i punti sull'ascissa che saranno % considerati nella rappresentazione dei segnali di base x = [xmin : passo : xmax]'; % Generazione Intervallo per la rappresentazione della Convoluzione % Procedimento utilizzato stesso che in precedenza, per i segnali. xmax_conv = abs(xmin - xmax); xmin_conv = - abs(xmin - xmax); passo2 = abs(xmax - xmin) / ((ncomp) - 1); x_conv = [xmin_conv : passo2 : xmax_conv]'; % Generazione Intervallo per la rappresentazione degli Spettri delle Immagini fcamp = round(ncomp / (abs(range(1,1)-range(1,2)))); x_spettro = [-fcamp/2:fcamp/(ncomp*2-1):fcamp/2]; % Generazione Intervallo per la rappresentazione dello Spettro della Convoluzione fconv = length(x_conv) / (abs(xmin_conv-xmax_conv)); x_spconv = [-fconv/2:fconv/(length(x_conv)-1):fconv/2]; % Adattamento per la rappresentazione spettro_onda1(ncomp*2)=0; spettro_onda2(ncomp*2)=0; % Rapporto usato per la verifica di correttezza: rapporto = (spettro_m ./ spettro_conv) ; % Predisposizione alla Visualizzazione set(gcf, 'color', 'white'); set(gcf, 'InvertHardCopy', 'off'); nyquist = fcamp/2;


137

% Rappresentazione subplot(4,2,1) plot(x,valori_onda1) title('Primo Segnale') grid on subplot(4,2,2) plot(x_spettro,spettro_onda1) axis([-nyquist,nyquist,min(spettro_onda1),max(spettro_onda1)]) grid on title('Spettro Primo Segnale') subplot(4,2,3) plot(x,valori_onda2) grid on title('Secondo Segnale') subplot(4,2,4) plot(x_spettro,spettro_onda2) axis([-nyquist,nyquist,min(spettro_onda2),max(spettro_onda2)]) grid on title('Spettro Secondo Segnale') subplot(4,2,[5,6]) plot(x_conv,convoluzione) grid on title('Convoluzione tra i due Segnali') subplot(4,2,7) plot(x_spconv,spettro_conv) axis([-nyquist,nyquist,min(spettro_conv),max(spettro_conv)]) grid on title('Spettro della Convoluzione') subplot(4,2,8) plot(x_spconv,spettro_m) axis([-nyquist,nyquist,min(spettro_m),max(spettro_m)]) grid on title('Moltiplicazione dei due Spettri')


138

8.4 Risultati Si riportano di seguito alcuni dei risultati più significativi, selezionati fra i test effettuati in laboratorio. Nota: Adottiamo negli esempi e relative indicazioni le frequenze comprese tra 0 e 9 come basse frequenze, le frequenze comprese tra 10 e 99 come medie frequenze e quelle uguali o maggiori di 100 come alte. La fase sarà da intendersi nulla, salvo dove espressamente specificato. 1. Analisi di due segnali seno (entrambi a “media” frequenza)

Segnale 1 Segnale 2

sinus(2,10,0,[-1,1],1000) sinus(2,50,0,[-1,1],1000)

Range -1 ÷ 1 N. Campioni 1000

Frequenze Segnali 10, 50

Chiamata es23(['s','s'], [2,2], [10,50],[0,0], [-1,1], 1000)


139

Estratto dell’output del rapporto ottenuto tra la moltiplicazione degli spettri e la trasformata della convoluzione:

1.00050024965194 1.00050025023322 1.00050024957151 1.00050025007927 1.00050025116307 1.00050025008773 1.00050025062835 1.00050025012510 1.00050025009817 1.00050025014064 1.00050025008638 1.00050025013648 1.00050025280779 1.00050025012867 1.00050025118994 1.00050025016869 1.00050025136317 1.00050025008739 1.00050025003967 1.00050025008382 1.00050025025553 1.00050025011678 1.00050025161287 1.00050025012670 1.00050025016563 1.00050025014834 1.00050025060872 1.00050025011781 1.00050024935169 1.00050025012840 1.00050024918001 1.00050025009452 1.00050024983812 1.00050025010984 1.00050025005605 1.00050025012451 1.00050024963289 1.00050025012481 1.00050025053439 1.00050025012288 1.00050024988262 1.00050025013434 1.00050025002125 1.00050025015233 1.00050025081042 1.00050025009553


140

2. Analisi di due segnali seno (entrambi a “media” frequenza) con numero maggiore di campioni:

Segnale 1 Segnale 2

sinus(2,10,0,[-1,1],10000) sinus(2,50,0,[-1,1],10000)


Frequenze Segnali 10,50

Chiamata es23(['s','s'], [2,2], [10,50],[0,0], [-1,1], 10000)


141

Estratto dell’output del rapporto ottenuto tra la moltiplicazione degli spettri e la trasformata della convoluzione: 1.00003040089673 1.00006421822579 1.00008062853994 1.00006068719216 1.00002610001019 1.00002314088993 1.00003928632984 1.00000382261068 1.00005275623477 1.00002833760374 1.00009504588775 1.00006719717977 1.00003983318427 1.00005944062831 1.00004183392241 1.00011546748635 1.00005097966319 1.00004874407369 1.00008884056895 1.00000819329971 1.00003166164637 1.00004719012799 1.00007728737544 1.00010377386153 1.00005851024020 1.00011616289572 1.00012996108348 1.00003041670578 0.99998558761900 1.00002479834265 1.00005162308479 1.00004992250612 1.00004273194199 1.00008255003255 1.00001091632832 1.00003758780681 1.00004436836269 0.99999501691386


142

3. Analisi di due segnali seno (entrambi a “bassa” frequenza), con numero inferiore di campioni

Segnale 1 Segnale 2

sinus(2,2,0,[-1,1],100) sinus(2,4,0,[-1,1],100)



Chiamata es23([‚s’,’s’], [2,2], [2,4],[0,0], [-1,1], 100)


143

4. Analisi di un segnale seno e un’onda triangolare (entrambi a bassa frequenza)

Segnale 1 Segnale 2

sinus(2,2,0,[-1,1],100) trian (2,4,0,[-1,1],100)



Chiamata es23([‚s’,’t’], [2,2], [2,4],[0,0], [-1,1], 100)


144

5. Analisi di un segnale seno e un’onda rettangolare (entrambi a “bassa” frequenza)

Segnale 1 Segnale 2

sinus(2,4,0,[-1,1],100) rectan(2,2,0,[-1,1],100)



Chiamata es23([‚s’,’r’], [2,2], [4,2],[0,0], [-1,1], 100)


145

6. Analisi di due onde rettangolari (entrambi a “bassa” frequenza)

Segnale 1 Segnale 2

rectan(2,1,0,[-1,1],100) rectan(2,3,0,[-1,1],100)



Chiamata es23([‚r’,’r’], [2,2], [1,3],[0,0], [-1,1], 100)


146

7. Analisi di un’onda rettangolare e un’onda triangolare (entrambi a “bassa” frequenza)

Segnale 1 Segnale 2

rectan(2,1,0,[-1,1],100) trian(2,4,0,[-1,1],100)



Chiamata es23([‚r’,’t’], [2,2], [1,4],[0,0], [-1,1], 100)


147

8. Analisi di due onde triangolari (entrambi a “bassa” frequenza)

Segnale 1 Segnale 2

trian(2,1,0,[-1,1],100) trian(2,4,0,[-1,1],100)



Chiamata es23([’t’,’t’], [2,2], [1,4],[0,0], [-1,1], 100)


148

9. Analisi di un segnale seno e un’onda rettangolare (entrambi ad “alta” frequenza)

Segnale 1 Segnale 2

sinus(2,1000,0,[-1,1],8000) rectan(2,250,0,[-1,1],8000)



Chiamata es23(['s','r'], [2,2], [1000,250],[0,0], [-1,1], 8000)


149

10. Analisi di un’onda triangolare e un’onda rettangolare (entrambi ad “alta” frequenza)

Segnale 1 Segnale 2

sinus(2,400,0,[-1,1],8000) rectan(2,250,0,[-1,1],8000)



Chiamata es23(['s','r'], [5,2], [400,250],[0,0], [-1,1], 8000)


150

11. Analisi di due segnali seno (entrambi a “bassa” frequenza), di cui il secondo con una fase di ½*pi

Segnale 1 Segnale 2

sinus(2,4,0,[-1,1],100) sinus(2,2,0,[-1,1],100)


Frequenze Segnali 4,2 Fase ½*pi

Chiamata es23(['s','s'], [4,2], [1,2],[0,(1/2)*pi], [-1,1], 100)


151

8.5 Commenti La prima cosa da notare nell’ambito di questo esercizio è senza dubbio il fatto che in questa sede è emerso un importante concetto nell'ambito dell'analisi dei segnali; questo esercizio ha permesso di evidenziare infatti quelli che sono gli effetti del cosiddetto "Teorema di Convoluzione", il quale asserisce che la convoluzione tra due segnali che ammettono rappresentazione spettrale, ammette anch'essa rappresentazione spettrale e lo spettro risultante è l'insieme dei prodotti termine a termine degli elementi dei due spettri di partenza. Tale risultato è ovviamente emerso anche durante le nostre esercitazioni e, come controprova, abbiamo utilizzato il rapporto punto a punto tra i valori dello spettro della convoluzione e quelli dello spettro ottenuto moltiplicando i due spettri iniziali. Tale rapporto si è sempre mantenuto pressoché costante in un intorno di 1 e, quando la precisione di macchina non era sufficientemente accurata da permettere di cogliere le discrepanze tra valori che si discostavano tra loro per quantità molto piccole, il rapporto si manteneva costante a 1. Per evidenziare così le differenze tra i due spettri, causate anche qui da errori di precisione, abbiamo aumentato la precisione di macchina utilizzata, per mezzo della chiamata format long, così da lavorare sempre su numeri a virgola fissa, ma a 14 cifre decimali, anziché solamente 4. In questo caso siamo stati in grado di cogliere gli errori dovuti alla precisione di calcolo, sebbene abbiamo comunque potuto notare che lo scostamento da 1 nel rapporto tra i valori dei due spettri si manteneva sempre e comunque molto basso, soprattutto quando si andava ad utilizzare un numero di campioni sufficientemente elevato (es. ε < 0.0001 con ncamp >= 10000). Abbiamo inoltre notato come in genere, all’aumentare di campioni, corrispondesse una maggiore omogeneità tra i due risultati. Tale fatto ha sicuramente confermato le nostre ipotesi iniziali e ha provato che eventuali discostamenti da tale aspettativa derivavano da errori di precisione. Abbiamo quindi scoperto che la convoluzione (procedimento complicato nel dominio spaziale, dove vi è la necessità di operare attraverso traslazione, folding e integrazione) si riduce ad un'operazione notevolmente più semplice nella rappresentazione spettrale; senza dubbio è da riconoscere che tale proprietà è quindi particolarmente importante nell’analisi dei segnali, dato che ci apre la possibilità di operare sullo spettro del segnale, ad esempio per eliminare intervalli di frequenze (filtraggio in frequenza), anziché dover ricorrere ad un’operazione complessa (e quindi computazionalmente più dispendiosa) quale la convoluzione. Dal punto di vista più tecnico (o per lo meno pratico, nell’ambito di questo esercizio), abbiamo comunque deciso di utilizzare per lo più un numero di campioni relativamente basso (100 elementi). Infatti, sebbene la scelta di un numero maggiore di campioni avrebbe portato ad una maggior precisione nei risultati (cosa peraltro non notabile dai grafici e quindi non particolarmente rilevante), l’aumento della frequenza di campionamento avrebbe avuto come ovvia conseguenza anche l’aumento delle frequenze considerabili (veniva in pratica esteso il limite di Nyquist); lo svantaggio di tale aspetto (che è in altri ambiti è da considerarsi positivo senza dubbio) è che la


152

“dilatazione” dell’intervallo di considerazione dello spettro, rendeva meno evidenti i risultati ottenibili. Abbiamo quindi fatto ricorso a frequenze relativamente basse, così da concentrare l’intervallo di spettro considerato ed evidenziarne i risultati. Per completezza, abbiamo fatto ricorso anche a frequenze di campionamento maggiori, ma solo dove ve ne fosse la necessità (in pratica dove volevamo prendere in considerazione segnali a frequenze più elevate), così da verificare la valenza del teorema in casi diversi.


153


9.1 Traccia Si calcoli la convoluzione tra due funzioni a gradino; una con n discontinuità poste rispettivamente in t1, t2, … , tn; l’altra un derivatore con lunghezza variabile N. Si calcolino e si visualizzino anche gli spettri delle due funzioni e della loro convoluzione.

9.2 Svolgimento Come prima cosa abbiamo creato due funzioni di supporto all’esercizio che andassero a generare due funzioni a gradino, secondo le specifiche richieste dall’esercizio. Per quanto riguarda la prima funzione, è stato scritto del codice che andasse a generare un segnale a gradino, finestrato in un intervallo [0,T], all’interno del quale la funzione risulta crescente a gradini, sino a T/2, decrescente di seguito. In secondo luogo, abbiamo generato il derivatore, semplicemente definendo in un intervallo [0,t] una funzione la cui ampiezza fosse uguale ad uno sino a t/2, fosse uguale a meno uno nella parte restante dell’intervallo. Abbiamo poi ricavato gli spettri delle due funzioni, eseguito la convoluzione sulle stesse e ottenuto anche lo spettro di quest’ultima.

9.2.1 Referenza alle primitive di Matlab utilizzate Chiamate di primitive Matlab utilizzate: abs(x): calcola il valore assoluto del parametro x sign(x): restituisce il segno del parametro x rand(m): genera un array di lunghezza m di numeri casuali uniformemente distribuiti

nell’intervallo [0,1] length(x): restituisce la lunghezza del vettore x ones(n,m): crea una matrice n x m riempita con valori 1 round(x): approssima il numero x all’intero più vicino fft(x): calcola la DFT del vettore x fftshift(x): posiziona al centro dello spettro il risultato della fot conv(s1,s2) esegue la convoluzione tra i segnali s1 e s2

9.2.2 Chiamate a funzioni da noi definite in precedenza Chiamate a funzioni esterne utilizzate: sinus: rif. Es 1 – Scheda 1


154

9.3 Codice Matlab

9.3.1 Codice Funzione Esercizio % Funzione es24(amp, gradini, xinterv, ncomp,

lunghezza_derivatore)

% Argomenti:

% - amp: limite superiore all'ampiezza massima dei gradini

% - gradini: numero di gradini

% - xinterv: intervallo della funzione gradino

% - ncomp: il numero di punti nell'intervallo

% - lunghezza_derivatore: la lunghezza del derivatore

%


% es24(5,4,[-1,1] ,1000, 5)

%

% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -


% Data: 20/11/2004

% Versione: 1.0

function y = es24(amp, gradini, xinterv, ncomp,

lunghezza_derivatore)

gradini = round(gradini/2);

xmin = xinterv(1,1);

xmax = xinterv(1,2);

% Funzione a Gradino

xgradino = gradino(amp/gradini, gradini, xinterv, ncomp);

intervallo = interv(xinterv,ncomp);

% Derivatore

xderivatore = derivatore(lunghezza_derivatore);

intd = interv(xinterv,lunghezza_derivatore);

% Viene generato l'intervallo da considerare per la convoluzione

convoluzione = conv(xgradino, xderivatore);

x_conv = interv([xmin,xmax],length(xgradino)+length(xderivatore)-

1);

% Generazione Spettro della Funzione Gradino


155

spettro_g = abs(fftshift(fft(xgradino)));

spettro_g = (spettro_g*2)/length(spettro_g);

% Generazione Spettro del Derivatore

spettro_d = abs(fftshift(fft(xderivatore)));

spettro_d = (spettro_d*2)/length(spettro_d);

% Generazione Spettro della Convoluzione

spettro_conv = abs(fftshift(fft(convoluzione)));

spettro_conv = (spettro_conv*2)/length(spettro_conv);

% Predisposizione alla Visualizzazione

set(gcf, 'color', 'white');

set(gcf, 'InvertHardCopy', 'off');

% SEGNALI

subplot(3,2,1)

plot(intervallo, xgradino);

title('Funzione a Gradino')

subplot(3,2,3)

plot(intd, xderivatore);

title('Derivatore')

subplot(3,2,5)

plot(x_conv, convoluzione);

title('Convoluzione tra i Segnali')

% SPETTRI

subplot(3,2,2)

plot(intervallo, spettro_g);

title('Spettro della Funzione Gradino')

subplot(3,2,4)

plot(intd, spettro_d);

title('Spettro del Derivatore')

subplot(3,2,6)

plot(x_conv, spettro_conv);


156

title('Spettro della Convoluzione')

9.3.2 Codice Funzione a Gradino % Funzione gradino (amp, freq, fase, intervallo)

% Argomenti:

% - amp: ampiezza della funzione base seno

% - freq: frequenza della stessa

% - fase: fase, sempre della stessa funzione seno

% - intervallo: intervallo della rappresentazione

%


% gradino (5, 2, 0, interv([-1,1], 100))

%

% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -


% Data: 16/11/2004

% Versione: 2.0

function y = gradino(amp, freq, fase, intervallo)


T = 1./freq

% Ricavo il segno di una funzione seno definita nell'intervallo

% considerato. Ad ogni semiperiodo (cambio di segno) corrisponderà

un

% gradino nella nostra funzione finale. Per permettere all'utente

di

% poter influire sul numero di gradini presenti, infatti abbiamo

% sfruttato le caratteristiche di periodicità di una funzione

seno.

range(1,1) = min(intervallo)

range(1,2) = max(intervallo)

ncomp = length(intervallo)

q = sign(sinus(amp,freq,fase+0.001,range,ncomp))

% Ricaviamo il segno iniziale

qsign = q(1);


157

% Ricaviamo un valore casuale k

k = rand(1) * amp;

for i = 1:length(intervallo)

% Se il segno non è cambiato, siamo sempre sullo stesso gradino

if (q(i)==qsign)

y(i)=k;

% Il segno cambia, vi è un nuovo gradino

else

y(i) = k;

qsign=q(i);

% Qui determiniamo le fasi crescente e descrescente della nostra

funzione

if (i<length(intervallo)/2)

k = k + rand(1);

else

k = k -rand(1);

end

end

end

% Le tre righe seguenti ci sono servite durante lo sviluppo per

meglio

% mettere in corrispondenza il rapporto esistente tra la funzione

seno

% e la nostra funzione a gradino

% subplot(2,1,1)

% plot(x,sinus(amp,freq,fase+0.001,range,ncomp));

% subplot(2,1,2)

plot(intervallo,y);

9.3.3 Codice Derivatore % Funzione derivatore (intervallo)

% Argomenti:

% - lunghezza: lunghezza n del derivatore


158

%


% derivatore(4)

%

% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -


% Data: 16/11/2004

% Versione: 1.0

function y = derivatore(lunghezza)

y = ones(lunghezza,1);

y(ceil(lunghezza/2)+1:lunghezza) = -1


159


1. Convoluzione tra una funzione con 4 gradini posti a t distanza l’uno dall’altro e un derivatore di 40 elementi

N° Gradini 4

N. Campioni 1000 Dim. Derivatore 40

Chiamata es24(5,4,[-1,1],1000,40)


160


N° Gradini 8


Chiamata es24(5,8,[-1,1],1000,40)


161


N° Gradini 16


Chiamata es24(5,16,[-1,1],1000,40)


162


N° Gradini 8


Chiamata es24(5,8,[-1,1],1000,80)


163


N° Gradini 8


Chiamata es24(5,8,[-1,1],1000,20)


164

6. Convoluzione tra una funzione con 8 gradini posti a t distanza l’uno dall’altro e un derivatore di 20 elementi.

N° Gradini 8


Chiamata es24(5,8,[-1,1],100,20)


165


N° Gradini 8


Chiamata es24(5,8,[-1,1],100,2)


166


N° Gradini 8


Chiamata es24(5,8,[-1,1],100,60)


167


N° Gradini 8


Chiamata es24(5,8,[-1,1],10000,60)


168

9.5 Commenti Questo esercizio ci ha permesso di capire qual è l’effetto dell’operazione di convoluzione tra una funzione a gradino e un cosiddetto derivatore. Abbiamo notato anche altri particolari interessanti, come ad esempio qual è l’effetto che deriva dalla modifica del numero di campioni del derivatore stesso. Abbiamo iniziato la nostra analisi con un caso (esempio 1) nel quale venivano individuati quattro gradini, per mezzo della convoluzione di un segnale di 1000 elementi con un derivatore di lunghezza 40, capendo così che il derivatore è in grado di individuare le discontinuità del segnale (si osserva in modo particolare che il risultato della convoluzione riportava dei picchi in prossimità di quelle che sono le discontinuità del segnale). Nell’esempio 2 abbiamo aumentato ad otto il numero dei gradini del segnale iniziale, confermando quanto visto nel caso precedente e verificando anche che il derivatore era in grado di mostrare l’entità delle discontinuità; l’esempio 3 non ha fatto altro che confermare tale affermazione e accentuarne gli effetti, anche graficamente, soprattutto negli ultimi 5 gradini della fase crescente del segnale (dove si vede che a discontinuità più o meno accentuate corrispondono simili risultati nella convoluzione). Nell’esempio 4 siamo ritornati ad otto gradini, per il segnale considerato, ma abbiamo aumentato ad 80 i campioni considerati nel derivatore, scendendo a 20 nel caso successivo (esempio 6); tale breve excursus basato sulla modifica del numero di campioni considerati per il derivatore, ci ha permesso di verificare che al decrescere dei campioni del derivatore, maggiore diventava la precisione con la quale venivano individuate le discontinuità. Per capire tali risultati, è opportuno e probabilmente sufficiente pensare a come viene eseguita l’operazione di convoluzione, tra i due segnali considerati (per ulteriori referenze, si veda es. 2 – scheda 1); rappresentiamo di seguito un ingrandimento approssimato di due situazioni, in cui il derivatore ha rispettivamente un numero di campioni elevato (o comunque maggiore del numero di campioni che costituiscono un singolo gradino della prima funzione) e un numero di campioni pari a due.

Notiamo in particolare, come nel primo caso vengano erroneamente considerate nella convoluzione anche discontinuità estranee a quella che si vorrebbe individuare e quindi come vi sia una perturbazione dell’approssimazione di derivata che vogliamo effettuare; nel secondo caso, invece


169

essendo solo due gli elementi del segnale iniziale considerati, le discontinuità saranno individuate precisamente nei punti dove si verificano e non saranno influenzate da altre discontinuità, collateralmente prese in considerazione. Per tale motivo, possiamo affermare che 2 rappresenta senza dubbio il numero di campioni più adatto per costruire un derivatore; d’altro canto, un peggioramento dell’approssimazione dell’operazione di derivata che stiamo facendo dipenderà non solo dal numero di campioni considerati per il derivatore, ma anche dal numero di campioni del segnale originale e conseguentemente dal numero di campioni che compongono ciascun gradino. Vale la pena di notare infatti che nell’esempio 6, benchè il numero di campioni del derivatore sia analogo al caso precedente, ad un abbassamento del numero di campioni del segnale corrisponde una peggior qualità dell’approssimazione dell’operazione di derivazione, approssimata dalla convoluzione. Tale affermazione è confermata nei due casi seguenti; in entrambi abbiamo infatti usato un derivatore da 60 campioni, sempre in un intervallo [–1,1], tuttavia mentre il primo segnale constava di 100 elementi, il secondo era formato da 10000 campioni. In questo caso, si possono fare considerazioni e confronti interessanti anche con i casi precedenti; infatti vediamo che l’esempio 8 presenta risultati peggiori rispetto al caso 4, in cui avevamo 80 campioni per il derivatore, ma anche 1000 per il segnale, mentre il caso 9 si presenta migliore dei casi 6, ma anche 5, dove avevamo solo 20 campioni per il derivatore. Ovviamente, una diminuzione drastica del numero di campioni considerati per il derivatore, implica anche la presenza di un numero di campioni troppo basso per poter generare e quindi visualizzare uno spettro significativo per tale segnale, ma anche per visualizzare il segnale stesso (caso esemplificativo di ciò è ovviamente il 7). Al di là di questa osservazione, che vuole solo giustificare la poca significatività di alcune rappresentazioni, concludiamo comunque dicendo che il numero di campioni per un derivatore che meglio approssima l’operazione di derivazione è 2, mentre il risultato peggiorerà al crescere del rapporto (campioni del derivatore/campioni del segnale).


170


10.1 Traccia Si calcoli convoluzione tra la prima delle due funzioni dell’esercizio precedente ed il boxcar dell’esercizio 1, di lunghezza variabile N. Si calcolino e si visualizzino anche gli spettri delle due funzioni e della loro convoluzione.

10.2 Svolgimento Per la risoluzione dell’esercizio abbiamo definito una funzione es25 che, accettati in input i parametri ampg (ovvero l’ampiezza massima della nostra funzione a gradino), ampi (ampiezza dell’impulso), starti (inizio dell'impulso), lunghi (lunghezza dell'impulso), xinterv (l’intervallo considerato per la funzione gradino), ncomp (il numero di campioni del primo segnale) e gradini (il numero di gradini, sempre del primo segnale), andasse a generare le due funzioni richieste e ad eseguirne la convoluzione, visualizzandone i relativi spettri. La costruzione del segnale a gradino è stata effettuata per mezzo della funzione gradino, impiegata nell’esercizio precedente (Es 2.4, Scheda 2), mentre per la costruzione del Boxcar, si è utilizzata la funzione impulso (Es 1.3, Scheda 1). Dopo aver generato gli spettri dei due segnali (usando la primitiva Matlab fft) e la convoluzione tra essi, abbiamo generato anche lo spettro di questa ed eseguito le diverse rappresentazioni; per la visualizzazione degli spettri, abbiamo usato una scala logaritmica, poiché ci è sembrato questa fosse più adatta a rappresentare in modo significativo i segnali analizzati.

10.2.1 Referenza alle primitive di Matlab utilizzate Chiamate di primitive Matlab utilizzate: abs(x): calcola il valore assoluto del parametro x length(x): restituisce la lunghezza del vettore x fft(x): calcola la DFT del vettore x fftshift(x): posiziona al centro dello spettro il risultato della fot conv(s1,s2) esegue la convoluzione tra i segnali s1 e s2

10.2.2 Chiamate a funzioni da noi definite in precedenza Chiamate a funzioni esterne utilizzate: gradino: rif. Es 4 – Scheda 2 impulso: rif. Es 3 – Scheda 1


171

10.3 Codice Matlab

10.3.1 Codice Funzione Esercizio % Funzione es25(ampg, ampi, starti, lunghi, xinterv,

ncomp,gradini)

% Argomenti:

% - amp: limite superiore all'ampiezza massima dei gradini

% - starti: inizio dell'impulso

% - lunghi: lunghezza dell'impulso

% - xinterv: intervallo della funzione gradino

% - ncomp: il numero di punti nell'intervallo

% - gradini: il numero di gradini della prima funzione

%


% es25(5,1,60,20,[-1,1] ,1000, 8)

%

% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -


% Data: 20/11/2004

% Versione: 1.0

function y = es25(ampg, ampi, starti, lunghi, xinterv, ncomp,

gradini)

gradini = round(gradini/2);

xmin = xinterv(1,1);

xmax = xinterv(1,2);

% Funzione a Gradino

xgradino = gradino(ampg/gradini, gradini, xinterv, ncomp);

intervallo = interv(xinterv,ncomp);

% Derivatore

ximpulso = impulso(ampi, ncomp, starti, lunghi);

% Viene generato l'intervallo da considerare per la convoluzione

convoluzione = conv(xgradino, ximpulso)/ncomp;

x_conv = interv([xmin,xmax],length(xgradino)+length(ximpulso)-1);

% Generazione Spettro della Funzione Gradino


172

spettro_i = abs(fftshift(fft(xgradino)));

spettro_i = (spettro_i*2)/length(spettro_i);

% Generazione Spettro del Derivatore

spettro_d = abs(fftshift(fft(ximpulso)));

spettro_d = (spettro_d*2)/length(spettro_d);

% Generazione Spettro della Convoluzione

spettro_conv = abs(fftshift(fft(convoluzione)));

spettro_conv = (spettro_conv*2)/length(spettro_conv);




% SEGNALI

subplot(3,2,1)

plot(intervallo, xgradino);

grid on;

title('Funzione a Gradino')

subplot(3,2,3)

plot(intervallo, ximpulso);

grid on;

title('Impulso')

subplot(3,2,5)

plot(x_conv, convoluzione);

grid on;

title('Convoluzione tra i Segnali')

% SPETTRI

subplot(3,2,2)

plot(intervallo, log(spettro_i));

grid on;

title('Spettro della Funzione Gradino')

subplot(3,2,4)

plot(intervallo, spettro_d);


173

grid on;

title('Spettro del Derivatore')

subplot(3,2,6)

plot(x_conv, log(spettro_conv));

grid on;

title('Spettro della Convoluzione')



174

1. Convoluzione tra una funzione con 8 gradini posti a t distanza l’uno dall’altro e un impulso unitario, anche in ampiezza, che inizia all’istante 0.

N° Gradini 8

N. Campioni 1000 Dim. Impulso 1 Inizio Impulso 0

Chiamata es25(5,1,0,1,[-1,1],1000,8)


175

2. Convoluzione tra una funzione con 8 gradini posti a t distanza l’uno dall’altro e un impulso unitario, che inizia all’istante 0 e ha ampiezza pari a due.

N° Gradini 8


Chiamata es25(5,2,0,1,[-1,1],1000,8)


176

3. Convoluzione tra una funzione con 8 gradini posti a t distanza l’uno dall’altro e un impulso di 40 campioni, che inizia all’istante 0. L’ampiezza di tale impulso d’ora in poi, sarà sempre unitaria.

N° Gradini 8


Chiamata es25(5,1,0,40,[-1,1],1000,8)


177

4. Convoluzione tra una funzione con 8 gradini posti a t distanza l’uno dall’altro e un impulso di 80 campioni, che inizia all’istante 0.

N° Gradini 8


Chiamata es25(5,1,0,80,[-1,1],1000,8)


178


N° Gradini 8


Chiamata es25(5,1,0,400,[-1,1],1000,8)


179


N° Gradini 8


Chiamata es25(5,1,0,999,[-1,1],1000,8)


180

7. Convoluzione tra una funzione con 8 gradini posti a t distanza l’uno dall’altro e un impulso unitario, che inizia all’istante 200.

N° Gradini 8


Chiamata es25(5,1,200,1,[-1,1],1000,8)


181


N° Gradini 8


Chiamata es25(5,1,600,1,[-1,1],1000,8)


182


N° Gradini 8


Chiamata es25(5,1,900,1,[-1,1],1000,8)


183

10. Convoluzione tra una funzione con 8 gradini posti a t distanza l’uno dall’altro e un impulso di 500 elementi, che inizia all’istante 250.

N° Gradini 8


Chiamata es25(5,1,250,500,[-1,1],1000,8)


184

10.5 Commenti L’esercizio svolto e gli esperimenti eseguiti ci hanno permesso di individuare quelli che sono gli effetti dell’operazione di convoluzione tra una funzione a gradino e un impulso. In particolare, abbiamo eseguito diversi esperimenti su una funzione a gradino con sette discontinuità (otto gradini quindi) e un impulso di ampiezza unitaria (eccetto che nell’esempio due). Nel primo caso, abbiamo utilizzato un impulso che avesse anche durata unitaria e inizio all’istante 0; abbiamo così potuto verificare come questo lasciasse immutata sostanzialmente la funzione in ingresso; la differenza nello spettro, rispetto a quella iniziale è dovuta alla diversità di intervallo considerato e assume tuttavia una forma che è lecito aspettarsi, anche per quelli che erano gli effetti del teorema di convoluzione (rif. esercizio 3, scheda due). Nel caso seguente, notiamo che ad un aumento di ampiezza dell’impulso, corrisponde un’analogo aumento di ampiezza anche nella convoluzione tra i due segnali (anche qui, differenze nello spettro e nei valori assunti dal segnale ottenuto per mezzo della convoluzione dipendono dalla nuova dimensione del segnale e, sostanzialmente dalla parte destra introdotta, che si presenta sostanzialmente come una sequenza di 0). Nei quattro casi successivi, abbiamo provato ad aumentare progressivamente la dimensione dell’impulso, verificando come ciò portasse ad un restringimento dello spettro dell’impulso stesso (cosa che ci si aspettava per altro, si veda es. 1 scheda 2), ma anche ad un progressivo appiattimento delle discontinuità nel segnale ottenuto per mezzo della convoluzione. A partire dal caso 7, abbiamo invece provato a utilizzare nuovamente un impulso unitario, spostandone però l’inizio ad un istante t, diverso da 0. Abbiamo notato in questi casi, come anche la funzione iniziale venisse traslata di conseguenza. Nell’esempio 10, infine, abbiamo voluto combinare le considerazioni fatte sin qui, utilizzando un impulso di lunghezza non unitaria e il cui inizio non fosse coincidente con l’origine; come ci aspettavamo, abbiamo ottenuto una figura le cui discontinuità erano state smorzate e la cui forma si presentava traslata, in maniera corrispondente a quella che era la traslazione dell’impulso, rispetto all’origine.


185


186

Scheda 3


187


11.1 Traccia Visualizzare le immagini con le primitive Matlab (imshow, imfinfo) anche utilizzando le colormap; visualizzare l’istogramma dei livelli di grigio; calcolare e visualizzare gli spettri di ampiezza e fase tramite la FFT2. Lavorare con immagini sintetiche e reali.

11.2 Svolgimento Per la risoluzione dell’esercizio è stata definita una funzione chiamata es31, che ricevuta un’immagine in ingresso come matrice di valori interi a 8 bit, utilizza le primitive matlab imshow, imhist e mesh per visualizzare rispettivamente l’immagine stessa, l’istogramma dei livelli di grigio di questa, e gli spettri di ampiezza e fase. In particolare, per il calcolo di questi ultimi, è stata utilizzata la primitiva Matlab fft2, che appunto genera la trasformata di Fourier, lavorando in due dimensioni. La scelta di utilizzare per la nostra funzione un parametro in ingresso di tipo immagine (rappresentata da una matrice di interi ad otto bit, ovvero in una scala di grigi, il cui valore minimo è 0 = nero, il massimo è 255 = bianco) è dovuta al fatto che così è possibile decidere di volta in volta quale immagine caricare (o creare). E’ infatti possibile leggere le immagini sia attraverso il comando matlab imread, sia utilizzando la procedura loadimg da noi definita, che permette anche la visualizzazione delle informazioni relative all’immagine (per mezzo di una chiamata alla primitiva Matlab imfinfo). Questo ovviamente per ciò che concerne le immagini reali; per quanto riguarda le immagini sintetiche, abbiamo invece realizzato due procedure distinte, che generano rispettivamente un’immagine secondo i parametri (dimensione e colore della figura da rappresentare) immessi dall’utente, oppure vanno a disegnare due figure di dimensioni e con disposizione pseudo-casuali. Anche in vista delle analisi che si andranno ad effettuare e delle relative considerazioni finali, vale la pena di ricordare che un’immagine, è rappresentabile da una funzione del tipo:

21: DDf a

dove D1 e D2 sono insiemi di osservabili fisiche. Analogamente a quanto visto per i segnali unidimensionali, che rappresentavano i cambiamenti in ampiezza del segnale nel tempo, qui si tiene conto dei cambi di intensità (sempre nella scala di grigi) nello spazio. In particolare, andremo ad operare di seguito su immagini digitali (l’unità di rappresentazione sarà il pixel e le loro dimensioni saranno finite, sia in larghezza, che altezza) unispettrali (vengono


188

considerati i livelli di grigio) e pertanto la descrizione funzionale di un’immagine di questo tipo può essere:

p = p(i,j) dove i e j rappresentano le coordinate del pixel preso in considerazione (con gli 0=<i=<m e 0=<j=<n, di un’immagine, ovvero una matrice di dimensioni m x n) e il valore di p corrisponde alla luminosità (o intensità) del punto considerato. Dell’immagine presa in considerazione, viene quindi eseguita la trasformazione nel dominio degli spettri, per mezzo di una trasformata di Fourier in 2 dimensioni,che consiste semplicemente nell’esecuzione di una trasformata in 1D su ogni riga dell’immagine e poi sull’esecuzione di un’altra trasformazione su ciascuna colonna sull’immagine così ottenuta. Siccome i nostri segnali saranno discreti, necessiteremo quindi di trasformazioni analoghe a quelle utilizzate per i segnali monodimensionali. Tale analogia è soddisfatta dalla seguente coppia di trasformazioni:

Quindi, un’immagine di MxN piexels, avrà anche un set di MxN coefficienti di Fourier complessi. Ognuno di questi punti complessi rappresenta una frequenza presente nell’immagine. Per eseguire tale trasformazione, tuttavia, sarebbe utile disporre di una trasformazione analoga alla FFT, vista nel dominio dei tempi, che ci permetta di computare in modo efficiente i coefficienti della trasformata; è possibile d'altronde separare la DFT 2D in due DFT ad una dimensione, che possono essere così computati con un algoritmo FFT, nella maniera seguente:

I coefficienti della DFT in due dimensioni (in maniera simile a quanto accadeva in 1D), sono disposti tuttavia in maniera non intuitiva, come si vede nello schema seguente:


189

E’ per questo motivo, che utilizziamo nell’esercizio la chiamata fftshift che, analogamente a quanto faceva nel dominio dei tempi, dispone tali coefficienti in maniera sicuramente più intuitiva, addossando i componenti a bassa frequenza attorno all’origine e, disponendo le componenti ad alte frequenza agli esterni dell’immagine. Vengono quindi scambiati rispettivamente il 1° con il 3° e il 2° con il 4° quadrante, in maniera analoga a quanto rappresentato di seguito:

Tale disposizione ci permetterà in seguito di verificare in maniera più intuitive alcune supposizioni e analogie con quanto abbiamo visto in 1D. La trasformata di Fourier produce in output un’immagine formata da valori complessi, che possono essere rappresentati per mezzo di due immagini (sia in termini di parti reali o immaginarie, sia in termini di ampiezza e fase). Nelle nostre visualizzazioni terremo spesso in maggior considerazione l’ampiezza appena citata, considerando comunque il fatto che per ricostruire l’immagine nello spazio reale, avremmo bisogno


190

di entrambe le componenti.

11.2.1 Referenza alle primitive di Matlab utilizzate

imread(x): legge l’immagine x e la carica in una matrice. imfinfo(x): restituisce informazioni sul file immagine x. imshow(x): visualizza l’immagine x. fft2(x): effettua una fft a due dimensioni sulla matrice x. abs(x): restituisce la potenza della trasformata x. angle(x): restituisce la fase della trasformata x. fftshift (x): effettua lo shift di coordinate per quanto riguarda lo spettro x colormap(‘parametro’): imposta a parametro il colormap da utilizzare nelle

rappresentazioni. imhist(x): crea un istogramm dei livelli di grigio per l’immagine x. mesh(x): rappresenta tridimensionalmente una matrice x. L’altezza in particolare è il

corrispondente di un certo valore (intensità) in un punto. Zeros(m,n): costruisce una matrice di dimensione m x n composta da soli 0. Max(a,b): ritorna il valore più grande tra a e b. Round(x): esegue un arrotondamento all’intero più vicino del numero x. Rand(x): genera un vettore di x elementi di valori casuali compresi tra 0 e 1 con

distribuzione gaussiana. im2uint8(X): converte l’immagine x da una rappresentazione floating

point (valori da 0 ad 1) ad una rappresentazione a 8 bit (valori tra 0

e 255).

11.3 Codice Matlab

11.3.1 Codice Programma Principale % Funzione es31(immagine,colmap,scala)

% Argomenti:

% - immagine: immagine (sotto forma di matrice, i cui valori

% rappresentano diverse intensità di grigio)

% - colmap: la colormap da utilizzare per le rappresentazioni

% sono valide colormap quelle riconosciute da Matlab

% di default o quelle definite da utente.

% - scala: specifica la scala per la visualizzazione dello spettro

% scelte valide sono lin (lineare) e log (logaritmica)

%

% Esempi di utilizzo:

% es31(sintetica(40, 40, 30, 10, 0.4, 0.6),’jet’,’log’)


191

% es31(loadimg('prova.tif'),’jet’,’log’)

%

% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -


% Data: 04/12/2004

% Versione: 1.0

function y = es31(immagine,colmap,scala)

format long

% Si ricavano così gli spettri d'ampiezza e fase dell'immagine

ampiezza_fft_img = fftshift(abs(fft2(immagine)));

fase_fft_img = mat2gray(fftshift(angle(fft2(immagine))));




% Numero totale di pixel nell'immagine.

totpixel = size(immagine,1)*size(immagine,2);

% Massimo numero di pixel per ogni intensità

max_livelli = max(imhist(immagine)) + max(imhist(immagine))/20;

subplot(2,2,1)

imshow(immagine,'notruesize');

title('Immagine Originale')

% Visualizzazione Istrogramma

subplot(2,2,2)

imhist(immagine);

axis([-5,260,0,max_livelli]);

grid on

title('Istogramma Livelli di Grigio')

% Visualizzazione spettro d'ampiezza

subplot(2,2,3)

if (scala=='log') funzione = log(ampiezza_fft_img+1);

else funzione = ampiezza_fft_img;

end

mesh(funzione);

title('Spettro di Ampiezza')

% Spettro di Fase

subplot(2,2,4)

imshow(fase_fft_img,'notruesize');

title('Spettro di Fase')


192

colormap(colmap)

11.3.2 Codice Generatore Immagini Sintetiche 1 % Funzione sintetica(larghezza, altezza, dimensioni, bordo, colore1, colore2)

%

% Parametri:

% larghezza: la larghezza dell'immagine generata

% altezza: l'altezza dell'immagine generata

% dimensioni: coefficiente che imposta le dimensioni

% della figura da generare nell'immagine (range 0 - 90)

% Con coefficiente 0 non viene disegnata nessuna

% figura.

% bordo: scostamento in percentuale della figura disegnata dal bordo

% colore1: colore primario della figura

% colore2: colore secondario della figura disegnata

%

% Genera un'immagine di dimensioni larghezza X altezza con una

% rappresentazione sintetica centrata di dimensioni e colori selezionabili

% dall'utente.

%


% sintetica(40, 40, 30, 10, 0.4, 0.6)

%

% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -


% Data: 20/11/2004

% Versione: 1.0

function y = sintetica(larghezza, altezza, dimensioni, bordo, colore1, colore2)

% Creazione di un'immagine vuota (sfondo nero)

immagine = zeros(larghezza,altezza);

% Questo è il caso in cui l'utente vuole disegnare le due figure.

if (dimensioni>0)

% Generazione delle coordinate principali dei due solidi

baseli = max(1,round(larghezza*bordo/100));

baself = (larghezza - baseli +1);

baseai = max(1,round(altezza*bordo/100));

baseaf = (altezza - baseai +1);

% Controllo sul colore

if (colore1=='r') colore1 = max(rand(1),0.1);


193

end

if (colore2=='r') colore2 = max(rand(1),0.1);

end

% Quanto mi stacco dal centro...

aggiuntal = round((larghezza*dimensioni/100)/2);

aggiuntaa = round((altezza*dimensioni/100)/2);

f11 = round(((larghezza/2)-aggiuntal));

f12 = round(baseai);

f13 = round(((larghezza/2)+aggiuntal));

f14 = round(baseaf);

f21 = round(baseli);

f22 = round(((altezza/2)-aggiuntaa));

f23 = round(baself);

f24 = round(((altezza/2)+aggiuntaa));

% Primo Rettangolo

immagine(f11:f13,f12:f14) = colore1;

% Secondo Rettangolo

immagine(f21:f23,f22:f24) = colore2;

% Se l'utente decide di non disegnare alcuna figura, allora

% si usa il primo parametro per il colore immesso e lo si

% utilizza come sfondo.

else

immagine = immagine + colore1;

end

y = im2uint8(immagine');

11.3.3 Codice Generatore Immagini Sintetiche 2 % Funzione generata(larghezza, altezza)

%

% Genera un'immagine di dimensioni larghezza X altezza con due

% rettangoli di dimensione, posizione e colore pseudocasuali.

%

% Il primo di essi è collocato nella metà superiore dell'immagine,

% il secondo in quella inferiore.

%


% generata(larghezza, altezza)

%

% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -


194


% Data: 20/11/2004

% Versione: 1.0

function y = generata(larghezza, altezza)

% Creazione di un'immagine vuota

immagine = zeros(larghezza,altezza);

% La superficie dell'immagine viene divisa in due

% parti [verticalmente]

ba1 = 1;

a1 = round(altezza/2);

ba2 = a1+1;

a2 = altezza;

% Generazione del Primo Rettangolo

baseli = round( 1+ rand(1)*(larghezza/2));

baself = round(min((larghezza/100*95),rand(1)*larghezza+baseli));

baseai = round( ba1 + rand(1)*(a1/2));

baseaf = round(rand(1)*a1+baseai);

immagine(baseli:baself,baseai:baseaf) = max(0.1,rand(1));

% Generazione del Secondo Rettangolo

baseli = round( 1+ rand(1)*(larghezza/2));

baself = round(min((larghezza/100*95),rand(1)*larghezza+baseli));

baseai = round( ba2 + rand(1)*(a1/2));

baseaf = round(min((a2/100*95),rand(1)*a2+baseai));

immagine(baseli:baself,baseai:baseaf) = max(0.1,rand(1));

% Conversione in Integer a 8 Bit

y = im2uint8(immagine');

11.3.4 Codice Loader Immagini Pre-Esistenti % Funzione loadimg(nomefile)

%

% Legge un file di tipo immagine (nei formati e codifiche supportate

% da matlab) e ritorna una matrice con i valori di questa

%

% Argomenti:

% - nomefile: percorso del file immagine che si intende aprire

%


% loadimg('prova.tif')

%

% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -


195


% Data: 04/12/2004

% Versione: 1.0

function y = loadimg(nomefile)

% Visualizzazione delle informazioni relative al file aperto

imfinfo(nomefile)

% Caricamento dei pixel dell'immagine in una matrice

y = imread(nomefile);



196

1. Caso base: utilizzo della funzione sintetica e generazione di un’immagine priva di figure, ovvero di una superficie che presenta un’uniforme intensità di grigio. Nell’istogramma è ovviamente visibile un unico picco (quello relativo al colore 0.4), così come unico è anche il picco visibile nello spettro, ovvero la componente costante, di frequenza 0.

Tipologia Sintetica Dimensione 50 x 50 pixels Immagine

Colori Scala di Grigi (8bit) Colormap ‘gray’

Chiamata es31(sintetica(50, 50, 0, 10, 0.4,0.6),'gray',’log’)


197

2. Viene utilizzata nuovamente la funzione sintetica per la generazione di un’immagine semplice (una figura rettangolare). In questo caso il terzo parametro assume valore diverso da 0 e quindi diventino due i livelli di grigio presenti nell’istogramma (sfondo e colore della figura, ovvero bianco).



Chiamata es31(sintetica(40, 40, 40, 20, 0,1),'gray','log')


198

3. La funzione sintetica genera ancora una volta un’immagine semplice (una figura rettangolare), ma questa volta di dimensioni ridotte rispetto all’esempio precedente.



Chiamata es31(sintetica(40, 40, 8, 48, 0,1),'gray','log')


199

4. Viene utilizzata nuovamente la funzione sintetica per la generazione di un’immagine. In questo caso però il terzo parametro assume valore diverso da 0 e quindi viene creata una figura. Si nota a questo punto come diventino due i livelli di grigio presenti nell’istogramma (sfondo e colore della figura), così come si nota un cambiamento della conformazione spettrale. Per meglio evidenziare tale cambiamento nello spettro, viene impostata la colormap jet. Nonostante lo spettro si presenti comunque più variegato rispetto a prima, si può comunque notare come preponderino le frequenze basse, a causa dei pochi sbalzi di colore tra pixel adiacenti e le diffuse uniformità.


Colori Scala di Grigi (8bit) Colormap ‘jet’

Chiamata es31(sintetica(40, 40, 20, 10, 0.4,0.4),'jet',’log’)


200

5. Utilizziamo ora la funzione generata, così da creare un’immagine sintetica in cui anche il posizionamento delle figure è pseudocasuale. Proviamo anche a variare la colormap utilizzata, richiamando quella denominata “bone”.


Colori Scala di Grigi (8bit) Colormap ‘bone’

Chiamata es31(generata(40, 40),'bone',’log’)


201

6. Utilizziamo ora la funzione generata, come al passo precedente. Ricaviamo quindi nuovamente un’immagine sintetica e vediamo che, per le caratteristiche della funzione generata, questa si presenta diversa da quella del punto precedente. Interessante notare anche le diversità negli spettri. La colormap utilizzata in questo caso è quella denominata “copper”.


Colori Scala di Grigi (8bit) Colormap ‘copper’

Chiamata es31(generata(40, 40),'copper',’log’)


202

7. Utilizziamo nuovamente la funzione generata, per un’ultima volta, così da permetterci di sperimentare colormap diverse. Come ci attendiamo, immagine e spettro differiscono dai due casi precedenti.


Colori Scala di Grigi (8bit) Colormap ‘pink’

Chiamata es31(sintetica(40, 40, 3, 51, 0,0.6),'pink','log')


203

8. Passiamo ora all’analisi della seconda classe di immagini di cui era richiesta l’analisi, ovvero quelle reali. Utilizziamo in questo caso la funzione loadimg, appositamente definita, così da poter riportare in seguito le caratteristiche dell’immagine. Notiamo come ci stiamo trovando in presenza di un’immagine in cui i livelli di grigio dell’istogramma presentano tutti qualche valore (è stato infatti ottimizzato il contrasto). Vi è ovviamente una preponderanza dei valori più elevati, dato che in particolare, il cielo dell’immagine si presenta molto chiaro e ricopre una parte significativa di questa. Anche nello spettro questo è riscontrabile, dato che si nota un’addensamento nelle basse frequenze, dovuto verosimilmente alla vasta area di uniforme intensità nell’immagine, accompagnato comunque da un numero considerevole di componenti alle frequenze più alte, dovuti agli sbalzi piuttosto netti di intensità che si verificano in corrispondenza della linea che idealmente divide il cielo dagli elementi ritratti nella foto.

Tipologia Reale Dimensione 320 x 240 pixels Immagine


Chiamata es31(loadimg('reale1.tif'),'gray',’log’)


204

Informazioni Relative all’Immagine: Filename: [1x57 char] FileModDate: '27-Dec-2004 16:47:02 ' FileSize: 88852 Format: 'tif' FormatVersion: [] Width: 320 Height: 240 BitDepth: 8 ColorType: 'grayscale' FormatSignature: [73 73 42 0] ByteOrder: 'little-endian' NewSubfileType: 0 BitsPerSample: 8 Compression: 'Uncompressed' PhotometricInterpretation: 'BlackIsZero' StripOffsets: 11302 SamplesPerPixel: 1 RowsPerStrip: 240 StripByteCounts: 76800 XResolution: 180 YResolution: 180 ResolutionUnit: 'Inch' Colormap: [] PlanarConfiguration: 'Chunky' TileWidth: [] TileLength: [] TileOffsets: [] TileByteCounts: [] Orientation: 1 FillOrder: 1 GrayResponseUnit: 0.01000000000000 MaxSampleValue: 255 MinSampleValue: 0 Thresholding: 1 NewSubFileType: 0 Make: 'Canon ' Software: 'Adobe Photoshop Elements 2.0 ' DateTime: '2004:27:12 16:47:02 '


205

9. Analizziamo ora una seconda immagine reale. Notiamo come ci troviamo nuovamente di fronte ad un’immagine che presenta pixel ben distribuiti tra tutte le varie intensità di grigio disponibili (vedere istogramma). Diversamente al caso precedente, però, quest’immagine è stata scattata all’interno, in condizioni di luminosità precarie ed è per questo che possiamo notare come la maggior parte dei pixels vadano a distribuirsi nell’intorno più basso dell’istogramma. Per quanto riguarda lo spettro, invece, possiamo notare che non vi è un’area con una distribuzione uniforme e diffusa di una certa intensità (come nel caso precedente) nell’immagine e quindi la preponderanza di basse frequenze, sebbene evidente è meno netta che nell’esempio precedente. La stessa affermazione può essere fatta per le alte frequenze, proprio perché vi è un numero inferiore di bruschi passaggi tra pixel adiacenti, con profonde diversità di intensità.





206

Informazioni Relative all’Immagine: Filename: [1x57 char] FileModDate: '27-Dec-2004 17:01:06' FileSize: 88976 Format: 'tif' FormatVersion: [] Width: 320 Height: 240 BitDepth: 8 ColorType: 'grayscale' FormatSignature: [73 73 42 0] ByteOrder: 'little-endian' NewSubfileType: 0 BitsPerSample: 8 Compression: 'Uncompressed' PhotometricInterpretation: 'BlackIsZero' StripOffsets: 12132 SamplesPerPixel: 1 RowsPerStrip: 240 StripByteCounts: 76800 XResolution: 72 YResolution: 72 ResolutionUnit: 'Inch' Colormap: [] PlanarConfiguration: 'Chunky' TileWidth: [] TileLength: [] TileOffsets: [] TileByteCounts: [] Orientation: 1 FillOrder: 1 GrayResponseUnit: 0.01000000000000 MaxSampleValue: 255 MinSampleValue: 0 Thresholding: 1 NewSubFileType: 0 Make: 'Canon ' Software: 'Adobe Photoshop Elements 2.0 ' DateTime: '2004:27:12 17:01:06 '


207

10. Rimaniamo nell’ambito delle immagini reali e prendiamo ora in considerazione un’immagine buia, leggermente sottoesposta. Notiamo che i valori maggiori nell’istogramma si riscontrano nei livelli inferiori (colori più scuri). Possiamo anche notare che nello spettro vi è una netta preponderanza di basse frequenze , date dalla presenza di aree uniformi; le alte frequenze non appaiono invece in modo evidente, perché assenti sono le differenze sostanziali di livello di intensità tra pixel adiacenti.





208

Informazioni Relative all’Immagine: Filename: [1x57 char] FileModDate: '28-Dec-2004 17:28:36' FileSize: 85436 Format: 'tif' FormatVersion: [] Width: 320 Height: 240 BitDepth: 8 ColorType: 'grayscale' FormatSignature: [73 73 42 0] ByteOrder: 'little-endian' NewSubfileType: 0 BitsPerSample: 8 Compression: 'Uncompressed' PhotometricInterpretation: 'BlackIsZero' StripOffsets: 8590 SamplesPerPixel: 1 RowsPerStrip: 240 StripByteCounts: 76800 XResolution: 72 YResolution: 72 ResolutionUnit: 'Inch' Colormap: [] PlanarConfiguration: 'Chunky' TileWidth: [] TileLength: [] TileOffsets: [] TileByteCounts: [] Orientation: 1 FillOrder: 1 GrayResponseUnit: 0.0100 MaxSampleValue: 255 MinSampleValue: 0 Thresholding: 1 NewSubFileType: 0 Software: [1x29 char] DateTime: '2004:12:28 17:28:36 '


209

11. Trattiamo ora un caso simmetricamente opposto rispetto a quello appena visto, ovvero esaminiamo un’immagine sovraesposta, che appare quindi troppo chiara. Notiamo qui come l’addensamento nell’istogramma compaia nella parte opposta rispetto a quello dell’esercizio precedente (preponderanza di colori chiari) e come invece lo spettro presenti caratteristiche abbastanza simili, con la netta presenza di componenti a basse frequenze e invece componenti alle alte frequenze attenuate, derivanti dalla mancanza di sostanziali differenze di intensità tra pixel adiacenti.





210

Informazioni Relative all’Immagine: Filename: [1x57 char] FileModDate: '28-Dec-2004 18:32:41' FileSize: 86664 Format: 'tif' FormatVersion: [] Width: 320 Height: 240 BitDepth: 8 ColorType: 'grayscale' FormatSignature: [73 73 42 0] ByteOrder: 'little-endian' NewSubfileType: 0 BitsPerSample: 8 Compression: 'Uncompressed' PhotometricInterpretation: 'BlackIsZero' StripOffsets: 9820 SamplesPerPixel: 1 RowsPerStrip: 240 StripByteCounts: 76800 XResolution: 72 YResolution: 72 ResolutionUnit: 'Inch' Colormap: [] PlanarConfiguration: 'Chunky' TileWidth: [] TileLength: [] TileOffsets: [] TileByteCounts: [] Orientation: 1 FillOrder: 1 GrayResponseUnit: 0.0100 MaxSampleValue: 255 MinSampleValue: 0 Thresholding: 1 NewSubFileType: 0 Software: [1x29 char] DateTime: '2004:12:28 18:32:41 '


211

11.5 Commenti L’esercizio svolto ci ha permesso di familiarizzare con alcune delle primitive di Matlab per il trattamento delle immagini e con i relativi spettri, facendo anche emergere alcune analogie con i segnali monodimensionali. Abbiamo anzitutto imparato a caricare immagini da file e lavorare con le codifiche in scala di grigi, intere ad otto bit e in virgola mobile. Ci siamo accorti che il primo dei due formati è molto più parsimonioso in termini di spazio ed è quindi preferibile ove ve ne sia la possibilità di utilizzo (ad esempio per la semplice visualizzazione). Se però si necessita di lavorare sulle immagini e quindi di trattarle come matrici da utilizzare nell’esecuzione di calcoli (ad esempio, se si desidera ricavarne lo spettro, come in questo caso), può essere utile usufruire dell’aumento di precisione offerto dal formato floating point. Per quanto riguarda le operazioni compiute nell’ambito dell’esercizio, come prima cosa abbiamo ricavato l’istogramma di ogni immagine e ciò ci ha permesso di notare come le immagini con una migliore distribuzione dei pixel nei diversi livelli di intensità, risultino all’occhio umano più equilibrate e naturali. Ciò non è determinante nelle immagini sintetiche (soprattutto se ci si trova in presenza di uno scarso numero di colori effettivamente presente), ma risulta invece particolarmente evidente nelle immagini reali, soprattutto alla luce di quanto è possibile notare considerando i risultati dei test 8 e 9 e confrontando gli stessi con i risultati dei test 10 e 11. Nei primi due, le immagini, seppure preponderanti di elementi ai due estremi (pixel ad alta intensità nel primo, a bassa nel secondo, per via del contesto stesso dei soggetti raffigurati) hanno istogrammi ben distribuiti e appaiono quindi naturali e “nitide” a chi le visualizza; nei casi 10 e 11, invece, la distribuzione di colori addensata in una sola parte dell’istogramma giustifica e causa gli effetti “eccessivamente scuro” ed “eccessivamente chiaro” che caratterizzano rispettivamente le immagini. Un’altra operazione applicata in seguito alle immagini è stata le derivazione della rappresentazione spettrale, per mezzo della trasformata di Fourier bidimensionale, introdotta in precedenza. Anche in questo caso, come per i segnali monodimensionali, si è notata la necessità di distinguere tra valori reali e immaginari per la rappresentazione, e di applicare la funzione fftshift, come visto in precedenza, per rappresentare gli spettri in maniera più intuitiva. Abbiamo eseguito dapprima una prova con un’immagine in cui era presente un'unica intensità di grigio equidistribuita (test 1) e ci siamo resi conto di come una componente costante in un’immagine avesse frequenza 0. Ci siamo poi accorti di alcune somiglianze tra gli spettri a due dimensioni e quelli ricavati dai segnali monodimensionali. In particolare, nel test 2, abbiamo notato un’analogia dello spettro di un’immagine sintetica, in cui era rappresentata un’unica figura rettangolare con quanto accadeva in una dimensione con un boxcar. Il risultato nel dominio delle frequenze dell’immagine è stato un segnale di tipo sinc a due dimensioni. Considerando a com’è definita una trasformata di Fourier a 2 dimensioni, era perfettamente lecito supporre un simile comportamento. Nel test 3 abbiamo notato un’ulteriore particolarità della rappresentazione spettrale bidimensionale, utilizzando una figura rettangolare, di dimensione inferiore rispetto a quella dell’esempio


212

precedente. Abbiamo notato come le componenti dello spettro fossero più dilatate, rispetto a quelle ricavate precedentemente. In maniera analoga di quanto avviene nei segnali unidimensionali, quindi, il restringimento di una figura rettangolare nel dominio spaziale, causa un’espansione del segnale sinc nel dominio delle frequenze. Per quanto riguarda gli esempi successivi, invece, abbiamo invece voluto testare le diverse opzioni di colormap offerte da matlab, utilizzando le impostazioni jet (default), bone, copper e pink, rispettivamente per gli esempi 4, 5, 6 e 7; il comportamento spettrale di questi esempi è stato coerente con quello visto in precedenza; soprattutto nell’esempio 7 si è accentuato ulteriormente l’effetto del rimpicciolimento della figura rappresentata sul sinc costituente lo spettro nell’immagine. Gli esempi 5 e 6, invece, così come i casi successivi (8, 9, 10 e 11) ci hanno permesso di notare come fosse tutto sommato abbastanza incisivo l’effetto della doppia trasformazione nel dominio degli spettri. L’esempio 3 ci ha permesso infine di notare come la colormap jet, utilizzata come default da Matlab costituisca probabilmente la migliore alternativa per la visualizzazione degli spettri in 3 dimensioni, mettendo in chiara evidenza (anche cromatica) le forme del segnale spettrale; infatti, anziché essere una scala di intensità per la stessa tonalità (come avviene invece con le colormap pink, gray, copper, ecc), si estende su diverse tonalità, facilitando così la visualizzazione tridimensionale.


213


12.1 Traccia Costruire e visualizzare anche in modalità grafica un’immagine sintetica. Calcolarne e visualizzarne lo spettro anche in presenza di rumore gaussiano additivo a media nulla.

12.2 Svolgimento Per la risoluzione dell’esercizio è stata definita una funzione chiamata es32, che ricevuta un’immagine in ingresso come matrice di valori interi a 8 bit, procede all’aggiunta del rumore gaussiano additivo a media nulla, per mezzo della funzione addnoise, da noi definita. In seguito, per ricavare lo spettro delle due immagini, viene utilizzata la primitiva di Matlab fft2, applicata in modo analogo a quello visto nell’esercizio precedente.

12.2.1 Referenza alle primitive di Matlab utilizzate imshow(x): visualizza l’immagine x. fft2(x): effettua una fft a due dimensioni sulla matrice x. abs(x): restituisce l’ampiezza della trasformata x. fftshift (x): effettua lo shift di coordinate per quanto riguarda lo spettro x colormap(‘parametro’): imposta a parametro il colormap da utilizzare nelle

rappresentazioni. mesh(x): rappresenta tridimensionalmente una matrice x. L’altezza in particolare è il

corrispondente di un certo valore (intensità) in un punto. Max(x): ritorna il valore più grande tra i numeri presenti nel vettore x im2double(X): converte l’immagine X in una matrice di numeri in floating

point (valori da 0 ad 1), aumentando la precisione delle operazioni

eseguibili su di essi. Size(X): Restituisce la dimensione del vettore X (se questo è multidimensionale, viene

restituito un valore per ogni dimensione, in formato array). Randn(x,y) : Restituisce una matrice X x Y di valori casuali normalmente distribuiti, con

media 0 e varianza 1.

12.2.2 Chiamate a funzioni da noi definite in precedenza Chiamate a funzioni utilizzate:

1 sintetica: rif. Es 1 – Scheda 3


214

12.3 Codice Matlab

12.3.1 Codice Programma Principale % Funzione es32(immagine,varianza,colormap,scala)

% Argomenti:

% - immagine: immagine (sotto forma di matrice, i cui valori

% rappresentano diverse intensità di grigio)

% - varianza: varianza del rumore gaussiano aggiunto all'immagine

% - colmap: la colormap da utilizzare per le rappresentazioni

% sono valide colormap quelle riconosciute da Matlab

% di default o quelle definite da utente.

% - scala: specifica la scala da utilizzare per la visualizzazione

% sono valide le scale lin (lineare) e log (logaritmica)

%


% es32(sintetica(25, 25, 60, 20, 0,1),0.04,'jet','log')

%

% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -


% Data: 04/12/2004

% Versione: 1.0

function y = es32(timmagine,varianza,colmap,scala)

format long

immagine = addnoise(im2double(timmagine),varianza);

immagine = abs(immagine ./ max(max(immagine)));

% Si ricavano così gli spettri d'ampiezza delle due immagini

spettro_immagine1 = fftshift(abs(fft2(timmagine)));

spettro_immagine2 = fftshift(abs(fft2(immagine)));

% Predisposizione alla Visualizzazione.



subplot(2,2,1)

imshow(timmagine,'notruesize');


% Visualizzazione Istrogramma.

subplot(2,2,2)

imshow(immagine,'notruesize');

title('Immagine con Rumore')

% Visualizzazione spettro d'ampiezza.


215

subplot(2,2,3)

if (scala=='log') spettro1 = log(spettro_immagine1);

else spettro1 = spettro_immagine1;

end

mesh(spettro1);

title('Spettro Immagine Originale')

% Spettro di Fase

subplot(2,2,4)

if (scala=='log') spettro2 = log(spettro_immagine2);

else spettro2 = spettro_immagine2;

end

mesh(spettro2);

title('Spettro Immagine con Rumore')

colormap(colmap)

12.3.2 Codice Generatore Rumore Gaussiano a Media Nulla % Funzione addnoise(immagine, varianza)

%

% Applica un rumore gaussiano di media nulla e varianza come definita

% da utente attraverso parametro all'immagine fornita come parametro.

%


% addnoise(generata(20, 20),0.4)

%

% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -


% Data: 20/11/2004

% Versione: 1.0

function y = addnoise(immagine,varianza)

immagine = im2double(immagine);

dimensioni = size(immagine);

% Viene generata una matrice con rumore di varianza specificata

rumore = (randn(dimensioni(1),dimensioni(2)) .* sqrt(varianza));

% Il rumore viene aggiunto all'immagine

y = (immagine + rumore);


216



217

1. Rumore con varianza 0.0001, applicato ad un’immagine con una figura sintetica di piccole dimensioni.


Colori Scala di Grigi (8bit) Varianza Rumore 0.0001

Chiamata es32(sintetica(25, 25, 10, 49, 0,1),0.0001,'jet','lin')


218

2. Rumore con varianza 0.01, applicato ad un’immagine con una figura sintetica di piccole dimensioni.



Chiamata es32(sintetica(25, 25, 10, 49, 0,1),0.01,'jet','lin')


219

3. Rumore con varianza 1, applicato ad un’immagine con una figura sintetica di piccole dimensioni.


Colori Scala di Grigi (8bit) Varianza Rumore 1

Chiamata es32(sintetica(25, 25, 10, 49, 0,1),1,'jet','lin')


220

4. Rumore con varianza 0.0001, applicato ad un’immagine con una figura sintetica di piccole dimensioni. Analisi spettrale con scala logaritmica.



Chiamata es32(sintetica(25, 25, 10, 49, 0,1),0.0001,'jet','log')


221

5. Rumore con varianza 0.01, applicato ad un’immagine con una figura sintetica di piccole dimensioni. Analisi spettrale con scala logaritmica.





222

6. Rumore con varianza 1, applicato ad un’immagine con una figura sintetica di piccole dimensioni. Analisi spettrale con scala logaritmica.





223

7. Rumore con varianza 1, applicato ad un’immagine con una figura sintetica di piccole dimensioni. Analisi spettrale con scala logaritmica.



Chiamata es32(sintetica(25, 25, 60, 20, 0,1),1,'jet','lin')


224

8. Rumore con varianza 0.0001, applicato ad un’immagine con una figura sintetica di dimensioni maggiori. Analisi spettrale con scala logaritmica.





225

9. Rumore con varianza 0.01, applicato ad un’immagine con una figura sintetica di dimensioni maggiori. Analisi spettrale con scala logaritmica.





226

10. Rumore con varianza 1, applicato ad un’immagine con una figura sintetica di dimensioni maggiori. Analisi spettrale con scala logaritmica.



Chiamata es32(sintetica(25, 25, 60, 20, 0,1),1,'jet','log')


227

12.5 Commenti L’esercizio ha evidenziato l’effetto del rumore applicato ad un’immagine; le manifestazioni di questo si evidenziano come ci aspettavamo, sia sull’immagine, sia sullo spettro relativo. Abbiamo notato che, in maniera simile a quanto succedeva nella dimensione temporale con i segnali, anche per le immagini, il rumore tende a distribuirsi nelle alte frequenze. E’ tuttavia questo un risultato plausibile, se si tiene in considerazione quella che è l’essenza stessa delle immagini. Abbiamo infatti osservato all’esercizio precedente, come queste siano funzioni, che associano ad una data coordinata spaziale (in 2 dimensioni), un valore di intensità. Tanto più elevata sarà quindi la differenza di intensità tra pixel adiacenti, maggiore sarà la frequenza dell’armonica richiesta per generare tale discontinuità. Notiamo come nell’esempio 1, una varianza bassa (più precisamente uguale a 0.0001), non vada a perturbare in maniera eccessiva né l’immagine, ne lo spettro. Ciò è plausibilmente dovuto al fatto che non si vengono a creare discontinuità evidenti nell’intensità di pixel adiacenti e quindi si nota come le basse frequenze siano ancora preponderanti. Se tuttavia si aumenta la varianza del rumore considerato, come avviene nei casi 2 e 3 (dove questa assume valori rispettivamente 0.01 e 1), notiamo come lo spettro cambi completamente il proprio aspetto, dato che compaiono sempre maggiori punti in cui la discontinuità con i pixel adiacenti è elevata. I casi 4, 5 e 7 presentano gli stessi parametri (immagine dalle medesime dimensioni e raffigurazione; varianza del rumore di 0.0001, 0.01 e 1 rispettivamente), ma lo spettro relativo è raffigurato per mezzo di scala logaritmica, così da evidenziarne ulteriormente il comportamento. Si nota come ad una varianza del rumore di 0.0001 corrisponda uno spettro che poco si discosta da quello originale; ad una varianza di 0.01 lo spettro ha ancora qualche somiglianza con quello originale (così come l’immagine stessa), mentre ad una varianza di 1, come nel caso di scala lineare, è praticamente impossibile individuare analogie sia nello spettro, sia nell’immagine. E’ stato introdotto in questo caso anche un test effettuato con una varianza del rumore di 0.05; la decisione di includere questo esempio nella relazione scaturisce dal fatto che all’occhio umano, in questo caso, l’immagine affetta da rumore ricorda in qualche modo quella originale, mentre non si potrebbe certo dire la stessa cosa confrontando solamente i due spettri. Come ultima cosa, tenendo in considerazione quanto appreso nell’esercizio 3.1, abbiamo provato ad eseguire le stesse operazioni su di un’immagine che presentava una figura generata sinteticamente di dimensioni maggiori. Negli esempi 8,9 e 10 la varianza impiegata è la stessa dei casi 1,2 e 3; nonostante lo spettro dell’immagine iniziale sia diverso (si può già notare una maggiore distribuzione alle alte frequenze), vediamo come l’azione del rumore sia comunque simile a quella che si aveva nei casi 4,5 e 7, con l’immagine di dimensioni inferiori. In particolare, si nota come la presenza esclusiva di elementi orizzontali e verticali nello spettro venga corredata anche da elementi diagonali. Ciò è da attribuirsi al fatto che, mentre nell’immagine originale compaiono solo elementi orizzontali o verticali (in


228

seguito alla conformazione dell’immagine sintetica generata), il rumore va ad inserire anche elementi “obliqui”, proprio a causa della natura “disordinata” che lo caratterizza.


229


13.1 Traccia Filtrare un’immagine attraverso la convoluzione nel dominio spaziale con un filtro passa-basso. Confrontare gli spettri dell’immagine iniziale e di quella filtrata. Eseguire lo stesso filtraggio nel dominio delle frequenze. Lavorare con un’immagine sintetica.

13.2 Svolgimento Per la risoluzione dell’esercizio è stata definita una funzione chiamata es33, che ricevuta in input un’immagine (affetta da rumore a distribuzione gaussiana di media nulla e varianza specificata dall’utente), ne effettua l’operazione di convoluzione e filtraggio nel dominio delle frequenze, per mezzo di filtri costruiti sulla base dei parametri dim_filtro (solo per la convoluzione) e var_filtro, immessi da utente. La convoluzione di cui si è appena fatto accenno, è ovviamente in due dimensioni ed è data dalla formula:

E’ proprio questa operazione che caratterizza il filtraggio nel dominio dei tempi; prima di proseguire con l’analisi del filtraggio nel dominio delle frequenze, vale la pena di confermare quanto avevamo intuito durante l’esecuzione dell’esercizio precedente, ovvero vale la pena di ricordare un’altra analogia riscontrata tra segnali monodimensionali e immagini in scala di grigi. Nei primi, una bassa frequenza caratterizza un segnale che “cambia lentamente” nel tempo, mentre un’alta frequenza caratterizza un segnale con “cambiamenti più repentini”; nelle immagini in scala di grigi tale concetto può essere esteso ad un contesto spaziale; se ci troviamo in presenza di basse frequenze, avremo cambi di intensità tra pixel adiacenti di lieve entità, mentre le alte frequenze caratterizzeranno variazioni più evidenti. Ricordiamo anche che un filtro passa-basso è un filtro lineare, che può essere applicato sia nel dominio spaziale (attraverso la convoluzione), sia nel dominio spettrale (per mezzo della moltiplicazione punto a punto dell’immagine da filtrare con un filtro opportuno); interessante notare che anche qui il procedimento è lo stesso che avevamo impiegato nei segnali ad una dimensione (rif. es. 2.3 – Teorema di convoluzione). Ricordiamo anche che dati due segnali analogici fI(t) ed fO(u), con t e u reali, si definisce filtro (o traduttore) un qualsiasi operatore )(⋅Φ tale che:

)())(( uftf OI =Φ


230

dove gli indici I e O stanno per Input e Output rispettivamente. Definiamo a questo punto due proprietà che tali operatori possono presentare:

• Linearità: Φ è lineare se, dati due segnali di entrata )(1, tf I e )(2, tf I vale

))(())(())()(( 2,21,12,21,1 tfctfctfctfc IIII Φ+Φ=+Φ

per ogni coppia ),( 21 cc di reali, e per ogni valore del reale t.

• Tempo-invarianza: L'operatore Φ si dice tempo-invariante se vale

)())(( ττ +=+Φ tftf OI ,Rt ∈∀ τ reale fissato

Un filtro nel dominio spaziale, come accennavamo, si ottiene per mezzo di una particolare matrice (maschera) applicata attraverso la convoluzione ad un’immagine (attraverso la primitiva matlab conv2); tale matrice viene generata con la primitiva di matlab fspecial; Questa è quadrata, dal numero di elementi dispari (in genere 3, 5, 7 o 9) e la distribuzione dei valori al suo interno (la cui somma, per mantenere la luminosità dell’immagine originale, dev’essere 1) caratterizza l’accentuazione del filtro; maggiore sarà la varianza di questa matrice, più saranno presi in considerazione, nella media pesata, i pixel che circondano quello preso in considerazione; avremo invece un filtro man mano meno incisivo, stringendo la nostra campana, sino ad arrivare a quello che è il “restringimento massimo”, che si verifica quando la matrice presenta un solo 1 al centro e 0 in ogni altra posizione. In questo caso, ovviamente, l’immagine rimane immutata, se sottoposta a convoluzione con tale matrice. Notiamo che nel dominio delle frequenze, il comportamento della gaussiana è invece diametralmente opposto e, ad una campana più stetta, corrisponde un filtraggio più netto (la cosa non dovrebbe sorprendere, sia da un punto di vista teorico – facendo mente locale sul teorema di convoluzione e su quella che è la trasformata di una gaussiana in 2D, ovvero un’altra gaussiana, la cui dilatazione è inversamente proporzionale a quella che si ha nel domino dei tempi – sia dal punto di vista pratico, considerando la moltiplicazione di un segnale, per una campana, con picco al centro e valori decrescenti mano a mano che ci si allontana da questo). I filtri passa-basso così realizzati avranno la caratteristica di lasciare inalterate nell’immagine le componenti che genereranno lievi discrepanze di intensità tra pixel adiacenti, attenuando invece le manifestazioni attribuibili alle alte frequenze; dato che come abbiamo appena visto nell’esercizio precedente è caratteristico del rumore presentare componenti alle alte frequenze, possiamo subito dedurre come uno degli obiettivi dei filtri passa-basso può essere proprio l’attenuazione del rumore; una possibile conseguenza collaterale di un filtro di questo tipo, è ovviamente il peggioramento dei dettagli (soprattutto i contorni) e l’effetto “sfuocato” che si viene a creare sull’immagine.


231

13.2.1 Referenza alle primitive di Matlab utilizzate imshow(x): visualizza l’immagine x. fft2(x): effettua una fft a due dimensioni sulla matrice x. abs(x): restituisce il valore assoluto della variabile x. fftshift (x): effettua lo shift di coordinate per lo spettro x colormap(‘parametro’): imposta a parametro il colormap da utilizzare nelle

rappresentazioni. mesh(x): rappresenta tridimensionalmente una matrice x. L’altezza in particolare è il

corrispondente di un certo valore (intensità) in un punto. max(x): ritorna il valore più grande tra i numeri presenti nel vettore x rand(x): genera un vettore di x elementi di valori casuali compresi tra 0 e 1 con

distribuzione gaussiana. im2uint8(x): converte l’immagine X in una matrice di interi a 8 bit (valori tra 0 e 255). im2double(x): converte l’immagine X in una matrice di numeri in floating point (valori

da 0 ad 1), aumentando la precisione delle operazioni eseguibili su di essi. size(x): Restituisce la dimensione del vettore X (se questo è multidimensionale, viene

restituito un valore per ogni dimensione, in formato array). randn(x,y) : Restituisce una matrice X x Y di valori casuali normalmente distribuiti, con

media 0 e varianza 1. fspecial: costruisce un filtro… conv2(a,b): esegue la convoluzione in 2D tra le immagini a e b mean(x): restituisce il valore medio dell’array x ifft2: esegue la trasformata inversa…


2 sintetica: rif. Es 1 – Scheda 3 3 addnoise: rif. Es 2 – Scheda 3


232

13.3 Codice Matlab

13.3.1 Codice Programma Principale % Funzione es33(immagine,varianza,colormap,scala)

% Funzione es33(immagine, var_rumore, dim_filtro, var_filtro)

% Argomenti:

% - immagine: l'immagine da considerare sotto forma di matrice

% - var_rumore: varianza del rumore che verrà sommato all'immagine

% - dim_filtro: la dimensione del filtro per la convoluzione (es. 3,5,7)

% - var_filtro: la varianza del filtro

%


% es33(es33(generata(5, 5),0.001,5,0.4))

%

% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -


% Data: 24/11/2004

% Versione: 1.1

function y = es33(immagine, var_rumore, dim_filtro, var_filtro)

format long

% IMMAGINE ORIGINALE

% Aggiungo del rumore all'immagine inizialmente fornita

% L'operazione mi ritorna un'immagine in formato double

immagine = abs(addnoise(immagine, var_rumore))

%spettro (reale)

fft_img = fftshift(fft2(immagine));

% creiamo questo per la visualizzazione

fft_img_abs = abs(fft_img);

% DOMINIO SPAZIALE

% Generazione del filtro nel dominio dei tempi per la convoluzione

% Attenzione la varianza non è la stessa che per il filtro nel dominio

% delle frequenze...

filtro = fspecial('gauss',dim_filtro,var_filtro)

sum(sum(filtro));

% Applicazione del filtro passa basso nel dominio spaziale

img_conv = abs(conv2(immagine,filtro,'same'));


233

% Controllo Sulla Luminosità

mean(mean(img_conv))

mean(mean(immagine))

% spettro della convoluzione

fft_conv = fftshift(abs(fft2(img_conv)));

differenzaf1 = abs(double(fft_img_abs) - double(fft_conv));

differenzas1 = abs(double(immagine)-double(img_conv));

% DOMINIO DELLE FREQUENZE

% Genero un filtro passa basso (nel dominio delle frequenze)

dimensione = size(immagine);

dima = dimensione(1);

dimb = dimensione(2);

filtro_freq = spettro_gauss(dima, dimb,10*var_filtro);

% spettro filtrato

filtrata_f = fft_img .* filtro_freq;

% spettro filtrato per rappresentazione

filtrata_frapp = abs(real(fft_img) .* filtro_freq);

% trasformata inversa di fourier per rit. al dom. degli spazi

filtrata = abs(real(ifft2(ifftshift(filtrata_f))));

differenzaf2 = abs(double(fft_img)-double(filtrata_f));

% RAPPRESENTAZIONI

% Prima Figura - Filtraggio nel Dominio dei Tempi

figure;





subplot(3,2,1)

imshow(immagine)


subplot(3,2,2)

colormap('jet')

mesh(log(fft_img_abs))



234

% CONVOLUZIONE

subplot(3,2,3)

imshow(img_conv,[])

title('Filtrata Spazialmente')

subplot(3,2,4)

mesh(log(fft_conv));

title('Spettro Immagine Filtrata Spazialmente');

% DIFFERENZE

subplot(3,2,5)

mesh(differenzas1)

title('Differenza Immagine')

subplot(3,2,6)

mesh(differenzaf1);

title('Differenza Spettro')

% Seconda Figura - Filtraggio nel Dominio delle Frequenze

figure;





subplot(3,2,1)

imshow(immagine)


subplot(3,2,2)

mesh(log(fft_img_abs))


% DOMINIO SPAZIALE

subplot(3,2,3)

imshow(filtrata,[])

title('Filtrata')

subplot(3,2,4)

mesh(log(filtrata_frapp));

title('Spettro Immagine Filtrata')

% DOMINIO DELLE FREQUENZE

subplot(3,2,5)


235

mesh(filtro_freq)

title('Filtro nel Dom. Frequenze')

subplot(3,2,6)

mesh(differenzaf2);

title('Differenza Spettri');

13.3.2 Codice Generatore Filtro Gaussiano 2D % Funzione spettro_gauss(dima, dimb, varianza)

% Argomenti:

% - dima: larghezza desiderata

% - dimb: altezza desiderata

% - varianza: varianza della gaussiana

%


% spettro_gauss (200, 100, 0.01)

%

% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -


% Data: 05/12/2004

% Versione: 1.1

function y = spettro_gauss(dima, dimb, varianza)

fgauss = gaussiana([-dima/2,dima/2],dima,varianza);

fgauss2 = gaussiana([-dimb/2,dimb/2],dimb,varianza)';

y = fgauss * fgauss2;

13.4 Risultati Si riportano di seguito alcuni dei risultati più significativi, selezionati fra i test effettuati in laboratorio. Nota: la varianza assunta dal rumore sarà di 0.2, salvo dove diversamente specificato.


236

1. Utilizziamo in questo caso due filtri passa-basso, uno per la convoluzione nel dominio spaziale, l’altro per il filtraggio (per mezzo di moltiplicazione) nel dominio degli spettri. La gaussiana di entrambi ha varianza 1 (il significato di utilizzare uno stesso parametro per varianze che hanno significati diversi sarà spiegato nelle considerazioni finali dell’esercizio). Per quanto riguarda la dimensione della matrice da utilizzarsi nel dominio spaziale, eseguiamo in questo caso tre test diversi, con matrici da 3, 5 e 9 elementi.

Tipologia Sintetica Immagine

Dimensione 40 x 40 pixels

Varianza del Filtro 1

Dimensioni Filtro nel Dominio Spaziale

3, 5, 9

Chiamata es33(sintetica(40, 40, 40, 20, 0,1), 0.2, 3, 1) es33(sintetica(40, 40, 40, 20, 0,1), 0.2, 5, 1) es33(sintetica(40, 40, 40, 20, 0,1), 0.2, 9, 1)

Fig. 1.1 Filtraggio nel dominio delle frequenze.


237

Fig. 1.2

Filtraggio nel dominio spaziale, attraverso convoluzione dell’immagine con la seguente matrice:


238

Fig. 1.3



239

Fig. 1.4

Filtraggio nel dominio spaziale, attraverso convoluzione dell’immagine con la seguente matrice: filtro =


240

2. Utilizziamo in questo caso un filtro passa-basso, aumentando la varianza a 2; tutti gli altri parametri restano invariati rispetto al punto precedente.





3, 5, 9


Fig. 2.1

Filtraggio nel dominio delle frequenze.


241

Fig. 2.2



242

Fig. 2.3



243

Fig. 2.4



244

3. Aumentiamo ora ulteriormente la varianza delle due campane, lasciando invece invariati gli altri parametri per l’esercizio.





3, 5, 9




245

Fig. 3.2



246

Fig. 3.3



247

Fig. 3.4



248

4. Proviamo ora dunque a diminuire la varianza, ovvero a “restringere” le nostre due campane.





3, 5, 9

Chiamata es33(sintetica(40, 40, 40, 20, 0,1), 0.2, 3, 0.75) es33(sintetica(40, 40, 40, 20, 0,1), 0.2, 5, 0.75) es33(sintetica(40, 40, 40, 20, 0,1), 0.2, 9, 0.75)



249

Fig. 4.2



250

Fig. 4.3



251

Fig. 4.4



252

5. Diminuiamo ulteriormente la varianza delle nostre campane, sino a che, nel domino spaziale, la matrice di convoluzione “degenera” nella matrice sotto raffigurata; l’immagine in questo caso rimane uguale (indipendentemente ovviamente dalle dimensioni della matrice di convoluzione).



Matrice usata per la convoluzione

Chiamata es33(sintetica(40, 40, 40, 20, 0,1), 0.2, 5, 0.1)


253


254

6. Diminuiamo in questo caso ulteriormente la varianza. Il risultato non cambia rispetto all’esempio precedente nel domino spaziale, mentre otteniamo una “degenerazione” in questo caso anche nel dominio delle frequenze, dove viene fatto passare un unico valore e più precisamente quello della componente costante di frequenza 0.



Matrice usata per la convoluzione



255

7. In questo esempio diminuiamo la presenza rumore, abbassando la varianza di questo. Per quanto riguarda la dimensione della matrice da utilizzarsi nel dominio spaziale, date le relativamente piccole dimensioni dell’immagine considerata, utilizzeremo d’ora in poi una matrice da 3 elementi; la varianza del filtro sarà di 1, sia nel dominio dei tempi (la matrice sarà quella del test 1.2), che delle frequenze.


Dimensione 40 x 40 pixels Varianza del Filtro 1 Dimensioni Filtro 3

Varianza Rumore 0.01

Chiamata es33(sintetica(40, 40, 40, 20, 0,1), 0.01, 3, 1)



256


257

8. In questo esempio perturbiamo ulteriormente l’immagine, aumentando la presenza di rumore (var. 0.4); per quanto riguarda la matrice di convoluzione e il filtro nel dominio degli spettri, valgono le considerazioni fatte al punto precedente.


Dimensione 40 x 40 pixels Varianza del Filtro 1 Dimensioni Filtro 3 Varianza Rumore 0.4




258


259

9. Impieghiamo in questo caso un’immagine di dimensioni inferiori, con bassa varianza del rumore. Per quanto riguarda la matrice di convoluzione e il filtro nel dominio degli spettri, valgono le considerazioni fatte al punto 7.






260


261

10. Utilizziamo sempre un’immagine “piccola”, aumentando tuttavia la varianza del rumore. Per quanto riguarda la matrice di convoluzione e il filtro nel dominio degli spettri, valgono le considerazioni fatte al punto 7.






262


263

13.5 Commenti Questo esercizio ci ha permesso di introdurre alcuni concetti legati al filtraggio delle immagini, sia nel dominio spaziale, che in quello delle frequenze. Nello specifico, dopo avere generato sinteticamente un’immagine, abbiamo aggiunto a questa rumore, con diverse varianze e abbiamo infine applicato ad essa azioni di filtraggio, sia nel dominio delle frequenze, sia nel dominio spaziale, per mezzo dell’operazione di convoluzione a due dimensioni, tra l’immagine originale affetta da rumore e una matrice opportunamente preparata. Inizialmente abbiamo fatto variare anche le dimensioni di tale matrice quadrata da 3 a 9 elementi, per vedere qual’era l’influenza di tale parametro sul risultato. Dopo alcuni esempi, abbiamo adottato una matrice 3x3 quale standard per la generazione del filtro nel dominio spaziale; per questo motivo, dato che eseguirne una trasformata di Fourier avrebbe sicuramente introdotto errori dovuti alla scarsa precisione, per quanto riguarda il dominio delle frequenze, abbiamo preferito costruire un filtro da noi, anziché ricavare lo spettro dalla matrice di convoluzione e poi moltiplicarlo per lo spettro dell’immagine, come sarebbe stato possibile per quelle che sono le premesse del teorema di convoluzione, analizzato in 1D. La nostra scelta ci ha anche permesso di verificare il fatto che in due dimensioni, la trasformata di una gaussiana rimane una gaussiana, la cui dilatazione attorno all’origine è inversamente proporzionale alla gaussiana del dominio spaziale. Abbiamo così costruito due “campane” che si stringevano o dilatavano, sia nel dominio spaziale che in quello delle frequenze, in maniera inversamente proporzionale, e ne abbiamo osservato gli effetti. Nel caso 1 abbiamo filtrato un’immagine sia nel dominio delle frequenze, sia nel dominio spaziale; in questo caso, abbiamo utilizzato matrici da 3, 5 e 9 elementi per la convoluzione. Nel caso 2 abbiamo provato a mantenere costante la varianza del rumore che affliggeva l’immagine e ad aumentare la varianza della gaussiana usata da filtro, sia nel dominio spaziale che in quello delle frequenze. Come ci aspettavamo, abbiamo notato che allargando la gaussiana, si aumentava la potenza del filtro nel dominio spaziale, mentre si diminuiva l’efficienza dello stesso, nel dominio delle frequenze, dato che venivano lasciate passare un numero maggiore di componenti alle alte frequenze. Il confronto con l’esempio precedente ci ha anche permesso di capire come nel dominio spaziale la varianza della gaussiana influisca sul contributo dei pixel adiacenti rispetto a quello considerato (maggiore la varianza, maggiore sarà il contributo dei pixel che circondano quello considerato), mentre le dimensioni della matrice influiscono sulla distanza dei pixel che vengono presi in esame (ad esempio, una matrice di 9x9 considererà un numero maggiore di elementi per calcolare la media del pixel in esame, causando così un effetto più evidente). Il caso 3 è stata una conferma di quanto appena detto e anzi, abbiamo visto come a livello spettrale, la distensione ulteriore del filtro faceva si che diminuisse ulteriormente l’effetto di questo (si confrontino i grafici “differenza spettri” dei test 2 e 3 – e soprattutto i relativi assi verticali - per un’immediata constatazione). Nel caso 4 abbiamo invece ristretto la nostra gaussiana bidimensionale, ottenendo così un filtro più efficace nel dominio degli spettri, come immediatamente verificabile. Per quanto riguarda invece il


264

dominio spaziale, notiamo invece come la diminuzione del peso per gli elementi che circondano il pixel considerato faccia si che la dimensione della matrice per la convoluzione sia meno determinante di quanto non lo fosse quando i valori in essa contenuti avevano maggiore distribuzione. Il caso 5 rappresenta un ulteriore restringimento della gaussiana. Nel dominio dei tempi questo si traduce in una matrice che mantiene identica l’immagine originale, mentre nel dominio delle frequenze, ciò fa si che l’immagine subisca un filtro piuttosto accentuato; il caso 6 è un ulteriore accentramento di tale fenomeno (un caso limite), tale che solo la componente costante viene fatta passare dal filtro. Le considerazioni da farsi per il dominio spaziale sono analoghe a quelle fatte nel punto precedente e verranno quindi qui omesse. A partire dal caso 7, abbiamo invece concentrato la nostra attenzione sugli effetti di un filtro passa-basso su di un’immagine cui era stato applicato del rumore con diversa varianza. Date le dimensioni relativamente piccole delle immagini considerate, si è deciso di adottare una matrice di dimensioni 3x3 per la convoluzione. Nel primo esempio (caso 7) abbiamo impiegato un’immagine contenente una figura di dimensioni relativamente grandi e un rumore di bassa entità (varianza del rumore 0.01); ciò ci ha permesso di verificare gli effetti del filtro sui contorni delle figure e quelle che potrebbero essere gli effetti collaterali di un filtro passa-basso; senza dubbio rappresentativi a tale scopo sono sia l’immagine filtrata nel dominio dei tempi, sia il grafico denominato “Differenza Immagine”; appare qui evidente l’effetto della convoluzione sull’immagine. L’esempio 8 si presenta diametralmente opposto al precedente, dato che viene invece presa in considerazione un’immagine fortemente affetta da rumore; interessante notare in questo caso come, anche in presenza di rumore molto forte, il filtro riesca a far percepire all’occhio umano quella che poteva essere la semplice figura originariamente raffigurata. Ovviamente nel caso di figure più complesse non si potrà sperare di ottenere simili risultati, perché ovviamente i maggiori dettagli subiranno lo stesso trattamento del rumore; verificheremo comunque tale comportamento tra poco. Ritorniamo invece con l’esempio 9 ad un’immagine affetta da rumore con varianza limitata (sempre 0.01), ma questa volta di dimensioni notevolmente inferiori. Notiamo in questo caso come possano essere fatte considerazioni analoghe a quelle fatte al punto 7. Spostando la nostra attenzione sull’esempio 10, invece, non possiamo più fare le medesime affermazioni. Notiamo infatti che l’immagine risulta affetta da forte rumore ed è difficilmente riconoscibile, anche dopo il filtraggio. Ciò ci ha fatto capire che, per quanto potente sia quest’operazione, non è sempre impiegabile per liberare immagini dal rumore; per poter raggiungere risultati accettabili ci dev’essere sempre un certo equilibrio tra la qualità e le dimensioni dell’immagine considerata e l’accentuazione del rumore che la perturba.


265


14.1 Traccia Considerare un’immagine sintetica a cui si è aggiunto rumore, ed estrarre i punti di discontinuità spaziale (edge) attraverso i filtri di Sobel, laplaciano a LoG. Confrontare i risultati per diversi valori della varianza stimata del rumore.

14.2 Svolgimento

Per la risoluzione dell’esercizio è stata definita una funzione chiamata es34, che ricevuta in input un’immagine (affetta da rumore a distribuzione gaussiana di media nulla e varianza specificata dall’utente), applica i diversi metodi di estrazione dei punti di discontinuità (edge detection) citati al punto precedente, così da individuare quelli che possono essere considerati come i contorni di eventuali figure presenti nell’immagine. Il primo metodo impiegato è quello che si basa sull’operatore di Sobel, un cosiddetto filtro derivativo di primo ordine, ovvero un filtro direzionale, alla cui base vi è l’idea del gradiente. Tanto più evidente sarà l’edge, maggiore sarà la magnitudine del gradiente, definito come il vettore

le cui componenti sono rappresentate dalle derivate parziali:

⎟⎟⎠

⎞⎜⎜⎝

⎛∂∂

∂∂

=∇ ),(),,(),( 000000 yxyfyx

xfyxf

Le derivate prime assumono che vi sia un Massimo in presenza di un edge, così per un’immagine continua ),( yxf , dove x e y sono rispettivamente le coordinate orizzontali e verticali, vengono

considerate di solito le due derivate direzionali: ),( yxfx∂ e ),( yxfy∂ .

A questo punto, il gradiente può essere rappresentato in termini di ampiezza e fase , rispettivamente per mezzo delle seguenti formule:

22 )),(()),((),( yxfyxfyxf yx ∂+∂=∇ e ⎟⎟⎠

⎞⎜⎜⎝

⎛∂∂

=∠∇),(),(

),(yxfyxf

ArcTanyxfx

y

La descrizione del vettore gradiente avviene quindi per mezzo di due componenti: - La forza dell’edge, ovvero la magnitudine del gradiente - La direzione dell’edge, ovvero l’angolo del gradiente Nella pratica, quindi viene eseguita una convoluzione tra l’immagine e gli operatori orizzontale e verticale (ottenuti per mezzo della chiamata matlab fspecial), prima di combinare i due risultati. Ricordiamo che tale combinazione è resa necessaria dal fatto che tale operatore è direzionale.


266

Supponendo dunque che siano Sx e Sy rispettivamente le nostre due maschere per la convoluzione, definite come segue:

Sx = Sy =

Per ottenere il risultato finale sull’immagine, eseguiamo due convoluzioni distinte, prima di calcolare l’ampiezza e la fase del nostro gradiente:

22

21 SSS += ⎥

⎦

⎤⎢⎣

⎡= −

2

11tanSSα

In particolare noi, negli esempi che seguiranno, chiameremo S l’immagine risultante dal nostro filtraggio. Il secondo filtro utilizzato (Laplaciano) è invece un filtro derivativo di secondo ordine; ciò significa che il filtro ricava gli edge in corrispondenza dei cambiamenti di segno della derivata seconda (zero crossing) e non è direzionale. L’operatore Laplaciano ed è definito come:

{ } { } ),(),(),( 2,2, yxfyxfyxf yx ∂+∂=∆

Viene creata in questo caso una sola maschera di convoluzione, applicata una sola volta; nel nostro specifico caso, sarà utilizzata la seguente matrice 3x3, ovvero quella predefinita di Matlab per la chiamata fspecial:

Come detto, il filtro è derivativo di secondo ordine, non direzionale e quindi il risultato del nostro filtraggio sarà la sola convoluzione tra l’immagine e la matrice 3x3 sopra raffigurata. Il terzo filtro rappresenta un’evoluzione di tale operatore o, per meglio dire, una combinazione di questo, con un filtro passa basso, di tipo gaussiano. L’acronimo LoG significa infatti Laplaciano di Gaussiano e consiste in pratica nella sottoposizione di un’immagine filtrata con un operatore gaussiano ad un filtraggio Laplaciano; tale filtro sarà quindi sempre di secondo ordine e non direzionale. Vi sono in pratica due modi matematicamente equivalenti per utilizzare questo operatore:


267

• eseguire la convoluzione tra un filtro gaussiano e l’immagine e poi eseguire una nuova convoluzione con un operatore Laplaciano.

• eseguire la convoluzione con un filtro lineare composto, che è proprio il Laplaciano di un Gaussiano.

Nello specifico, adotteremo questo secondo approccio, sempre per mezzo della primitiva Matlab fspecial, che ci genererà la seguente matrice di convoluzione nel caso standard 5x5:

Capita l’idea alla base di tale filtro, proveremo anche a modificare i parametri che concorrono alla generazione di tale matrice, per vedere l’effetto che essi hanno su un’immagine, in presenza di rumore con diverse intensità.

14.2.1 Referenza alle primitive di Matlab utilizzate o imshow(x): visualizza l’immagine x. o Size(X): Restituisce la dimensione del vettore X (se questo è multidimensionale, viene

restituito un valore per ogni dimensione, in formato array). o Fspecial: costruisce un filtro… o Conv2(a,b): esegue la convoluzione in 2D tra le immagini a e b


4 sintetica: rif. Es 1 – Scheda 3 5 addnoise: rif. Es 2 – Scheda 3


268

14.3 Codice Matlab

14.3.1 Codice Programma Principale % Funzione es34(immagine, var_rumore, dimensione, var)

% Argomenti:

% - immagine: l'immagine da considerare sotto forma di matrice

% - var_rumore: varianza del rumore che verrà sommato all'immagine

% - dimensione: dimensione del filtro gaussiano per il LoG

% - var: varianza del filtro gaussiano per il LoG

%


% es34(sintetica(40, 40, 40, 20, 0,1), 0.1, 5, 4)

%

% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -


% Data: 03/01/2004

% Versione: 2.0

function y = es34(iimmagine, var_rumore, dimensione, var)

format long


% Aggiungiamo del rumore all'immagine inizialmente fornita

% L'operazione ritorna un'immagine in formato double

immagine = abs(addnoise(iimmagine, var_rumore));

% FILTRAGGIO CON SOBEL

% Generazione del filtro Sobeliano nel dominio dei tempi per la convoluzione

filtro = fspecial('sobel');

% Applicazione del filtro passa basso nel dominio spaziale

img_sobel1 = abs(conv2(immagine,filtro,'same'));

img_sobel2 = abs(conv2(immagine,filtro','same'));

img_sobel = sqrt(img_sobel1.^2 + img_sobel2.^2);

% FILTRAGGIO LAPLACIANO

% Generazione del filtro Sobeliano nel dominio dei tempi per la convoluzione

filtro_lap = fspecial('laplacian');

% Applicazione del filtro tramite convoluzione

img_laplace = abs(conv2(immagine,filtro_lap,'same'));

% CREAZIONE DELL'ANTEPRIMA DEL FILTRAGGIO GAUSSIANO PER IL LOG


269

filtro_ant = fspecial('gaussian',dimensione,var);

anteprima = abs(conv2(immagine,filtro_ant,'same'));

% FILTRAGGIO LOG

% Generazione del filtro LoG

filtro_log = fspecial('log',dimensione,var);

% Applicazione del filtro

img_log = abs(conv2(immagine,filtro_log,'same'));

% RAPPRESENTAZIONI




% RIGA 1

% IMMAGINE DI PARTENZA

subplot(3,2,1)

imshow(iimmagine,[])


% IMMAGINE CON RUMORE

subplot(3,2,2)

imshow(immagine,[])

title('Immagine Con Rumore')

% RIGA 2

% FILTRAGGIO CON SOBEL

subplot(3,2,3)

imshow(img_sobel,[])

title('Filtro Sobeliano')

% FILTRAGGIO GAUSSIANO PER ELIMINARE IL RUMORE (ANTEPRIMA)

subplot(3,2,4)

imshow(anteprima,[])

title('Effetto del solo Passabasso')

% RIGA 3

% FILTRAGGIO LAPLACIANO

subplot(3,2,5)

imshow(img_laplace,[])

title('Laplaciano')

% FILTRAGGIO LOG

subplot(3,2,6)

imshow(img_log,[])


270

title('Laplaciano di Gaussiano')



271

1. Utilizziamo in questo caso un rumore con varianza 0.0001, applicato ad un’immagine sintetica; varianza della gaussiana e dimensione della matrice dell’operatore LoG sono quelli predefiniti di Matlab, ovvero 0.5 e 5 rispettivamente.


Dimensione 40 x 40 pixels Varianza Rumore 0.0001

Dimensioni Matrice LoG 5 Varianza Matrice LoG 0.5



272

2. Utilizziamo in questo caso un rumore con varianza 0.001, applicato ad un’immagine sintetica. Gli altri parametri presentano considerazioni simili a quelle fatte nel caso precedente.






273

3. Utilizziamo in questo caso un rumore con varianza 0.01, applicato ad un’immagine sintetica. Gli altri parametri presentano considerazioni simili a quelle fatte nei casi precedenti.






274







275







276







277

7. Proviamo ad aggiungere nuovamente un rumore con varianza 0.1 all’immagine sintetica. Applichiamo in seguito un filtraggio con operatore LoG la cui matrice ha dimensione 9 e varianza predefinita (0.5).






278

8. Riproviamo ad applicare un rumore con varianza 0.1 all’immagine sintetica, filtrandola in seguito con un filtraggio, il cui operatore LoG ha matrice di dimensione predefinita (5) e varianza 1.


Dimensione 40 x 40 pixels Varianza Rumore 0. 1

Dimensioni Matrice LoG 5 Varianza Matrice LoG 1



279

9. Riproviamo ad applicare un rumore con varianza 0.1 all’immagine sintetica, filtrandola in seguito con un filtraggio, il cui operatore LoG ha matrice di dimensione predefinita (5) e varianza 1.4.






280

10. Aggiungiamo ancora una volta un rumore con varianza 0.1 all’immagine sintetica, filtrandola in seguito con un filtraggio, il cui operatore LoG ha matrice di dimensione predefinita (5) e varianza 2.






281

11. Proviamo ora ad aggiungere nuovamente un rumore con varianza 0.01 all’immagine sintetica, lasciando invariati gli altri parametri, rispetto all’esempio precedente.






282

12. Dopo aver perturbato l’immagine sintetica con rumore dalla varianza di 0.01 all’immagine sintetica, la filtriamo con un operatore LoG, la cui matrice ha dimensione 3 e varianza 1.4.






283

13. Dopo aver perturbato un’immagine sintetica di dimensioni superiori a quelle degli esempi precedenti, con rumore dalla varianza di 0.01, la filtriamo con un operatore LoG, la cui matrice ha dimensione 3 e varianza 1.4, come nel caso precedente.






284

14. Dopo aver perturbato l’immagine sintetica dell’esempio precedente con rumore dalla varianza di 0.1 , la filtriamo con un operatore LoG, la cui matrice ha dimensione 15 e varianza 2.






285

14.5 Commenti In questo esercizio abbiamo confrontato tre diversi operatori di edge detection e ne abbiamo notato il diverso comportamento, in presenza di rumore, su un’immagine sintetica. Abbiamo iniziato a considerare casi in cui aumentavamo progressivamente la varianza del rumore che perturbava l’immagine (casi 1-6). In questa sede abbiamo notato come il comportamento dei tre operatori si manifesti simile nel caso in cui il rumore mantenga soglie contenute (esempi 1 e 2), mentre questo si diversifichi in maniera piuttosto netta nei casi seguenti (3-6); abbiamo notato in particolare come apparentemente (almeno alle condizioni in cui sono stati applicati durante le nostre sperimentazioni), il filtro di Sobel risulti meno condizionato del rumore, rispetto ai filtri che si basano sulla derivata seconda. In particolare, mentre l’operatore di Sobel permette di distinguere (se non altro all’occhio umano) contorni su figure presenti in immagini, in cui il rumore ha varianza anche di 0.1 (esempio 4, anche se, a dire il vero, anche nell’esempio 5, con varianza del rumore di 0.2 si riesce ancora a percepire qualcosa, anche se a malapena), ben diverso è il comportamento per i filtri Laplaciano e LoG, utilizzati con i parametri predefiniti di Matlab. Notiamo in questo caso che già ad una varianza di 0.1 (esempio 4) risulta impossibile distinuguere la minima traccia di contorno di qualche figura. La cosa che però ci ha sorpreso di più in questo caso, è l’apparente uguaglianza tra i risultati ottenuti tra il filtro Laplaciano e il filtro LoG nei primi 6 casi. Indubbiamente l’operatore di Laplace soffre in modo abbastanza evidente il rumore ed è questo il motivo per il quale viene utilizzato il LoG, che cerca, anteponendo un filtro gaussiano, di limitare la presenza di rumore, prima di ricorrere al filtro Laplaciano. Dato che però, almeno nei primi 6 esempi, non abbiamo trovato miglioramenti significativi del LoG rispetto al Laplaciano, abbiamo in seguito cercato di capire come e in che modo si manifestasse questo miglioramento. Abbiamo dapprima provato ad aumentare la dimensione della matrice della convoluzione (esempio 7) ma, come ci aspettavamo, non abbiamo riscontrato miglioramenti significativi. Nel caso 8 siamo allora ritornati ad una matrice di dimensione standard da 5 elementi, aumentando però la varianza di questa a 1; abbiamo notato come l’azione del filtro gaussiano fosse più evidente e come ciò iniziasse ad influire sull’operatore LoG, che iniziava a discostarsi rispetto a quello Laplaciano semplice. Abbiamo dunque provato ad aumentare ancora la varianza (esempio 9), sino ad un valore di 1.4; ci siamo accorti in questo caso che, sebbene il rumore risultasse ora meno evidente, anche i contorni iniziavano ad essere troppo offuscati, per poter essere determinati con precisione; sapendo che un modo per provare a risolvere tale problema poteva essere l’utilizzo di una matrice di convoluzione di dimensioni inferiori, dopo un ulteriore prova di convincimento (esempio 10, con varianza a 2), abbiamo provato a diminuire la dimensione della matrice a 3 (caso 11). Ad ogni modo, si noti che tale risultato si ottiene sempre con una varianza del rumore piuttosto limitata (0.001); pertanto, possiamo affermare che, se non altro quando parliamo di immagini di dimensioni simili a quelle considerate, nonostante il LoG risulti migliore di quello Laplaciano


286

semplice, comunque questo risentirà del rumore in maniera più forte di quanto non accada con il Sobeliano. Abbiamo sin qui accennato alle dimensioni dell’immagine considerata, poiché se le dimensioni cambiano (crescono, in particolare, come negli esempi 13 e 14, possiamo far emergere un altro interessante aspetto. Aumentando infatti le dimensioni dell’immagine, notiamo come, in presenza di certe condizioni, mentre il filtro Sobeliano presenti più o meno le stesse caratteristiche che presentava nelle immagini di dimensioni minori, diverso e interessante è soprattutto il comportamento del filtro LoG. Infatti notiamo qui che se scegliamo opportunamente la dimensione e la varianza della matrice per la convoluzione (come nel caso 14, dove invece il caso 13 presentava parametri non opportuni), possiamo raggiungere risultati decisamente migliori; lo stesso non si può dire del filtro Laplaciano! A questo punto, siamo in grado di mostrare e capire qual è l’utilità della combinazione dei due filtri e quali siano i risultati a cui portano. Pertanto, abbiamo notato come, per immagini di dimensioni relativamente piccole, il filtro Sobeliano sia quello che tollera in maniera più decisa il rumore; al crescere delle dimensioni dell’immagine, invece ciò non è più vero e, scegliendo opportunamente i parametri per la generazione della matrice del filtro LoG, notiamo come questo possa fornire i risultati migliori per quanto riguarda l’edge detection.

Teoria e Tecniche di Elaborazione...

Documents

Transcript of Teoria e Tecniche di Elaborazione...