Appuntidelcorsodi prof.ssa BiancamariadellaVecchia Trasformata di ... · Università“Sapienza”...

Università “Sapienza” Roma, a.a. 2008-2009

Appunti del corso di “Matematica Computazionale”

prof.ssa Biancamaria della Vecchia

Trasformata di Fourier&

Trasformata WaveletDiscreta

Antonio Petrella 1159178 [email protected]

1

Indice

1 Algoritmi di Compressione 41.1 Metodo dell’accetta . . . . . . . . . . . . . . . . . . . . . . . . . . 41.2 Rappresentazione Funzionale . . . . . . . . . . . . . . . . . . . . 4

2 Trasformata di Fourier 62.1 Sintesi e Analisi di Fourier . . . . . . . . . . . . . . . . . . . . . . 7

2.1.1 Risoluzione del sistema . . . . . . . . . . . . . . . . . . . 82.2 Trasformata Veloce di Fourier (FFT) . . . . . . . . . . . . . . . . 9

2.2.1 Aumentare la compattezza dell’analisi di Fourier (2.2) . . 102.2.2 Suddividere la trasformata . . . . . . . . . . . . . . . . . 112.2.3 Velocizzare la trasformata di Fourier (solo se N = 2L) . . 13

2.3 Considerazioni . . . . . . . . . . . . . . . . . . . . . . . . . . . . 152.3.1 Ridurre la dimensione dell’immagine . . . . . . . . . . . . 152.3.2 Trasformata di Fourier con numero di punti diverso da 2L 152.3.3 Trasformata di Fourier nell’algoritmo JPEG . . . . . . . . 152.3.4 Problemi della Trasformata di Fourier . . . . . . . . . . . 17

3 Trasformata Wavelet Discreta 183.1 Esempio Pratico . . . . . . . . . . . . . . . . . . . . . . . . . . . 18

3.1.1 Primo passo della trasformata . . . . . . . . . . . . . . . . 193.1.2 Passi successivi della trasformata . . . . . . . . . . . . . . 21

3.2 Wavelet di Daubechies . . . . . . . . . . . . . . . . . . . . . . . . 243.2.1 Convoluzione con Wavelet da 4 punti . . . . . . . . . . . . 24

3.3 Confronto con la Trasformata di Fourier . . . . . . . . . . . . . . 263.3.1 Localizzazione . . . . . . . . . . . . . . . . . . . . . . . . 263.3.2 Complessità . . . . . . . . . . . . . . . . . . . . . . . . . . 263.3.3 Compressione . . . . . . . . . . . . . . . . . . . . . . . . . 26

4 Trasformata Bidimensionale 27

2

Introduzione

Fotografia Analogica e DigitaleNella macchina fotografica analogica il segnale, considerato continuo, passaattraverso l’otturatore e impressiona la pellicola.

Per quanto riguarda la macchia fotografica digitale, nell’obbiettivo sono di-sposti dei sensori in grado di associare un valore all’intensità di luce che licolpisce.

Di conseguenza nel momento in cui viene scattata una foto verrà scritta unamatrice con i valori rilevati dai sensori (1 byte per sensore), c’è un problema: ladimensione della matrice potrebbe essere molto grande.

Supponiamo di avere una macchina fotografica da 1 megapixel, quindi 1.000.000di pixel, equivalenti a 3.000.000 di sensori nell’obbiettivo, abbiamo 3 sensori perogni pixel (uno per ogni colore della terna RGB1), ciò signifca scrivere unamatrice da 3.000.000 byte (circa 3 Mb) ogni qualvolta si scatta una foto.

Immaginiamo di voler inviare la foto appena scattata ad un amico usandola nostra connessione a 56k2 (trasferimento a cira 6.000 byte/s):

3:000:000

6:000= 500 secondi = 8 minuti

un po’ troppo per mandare una singola foto, anche perché il caso preso inconsiderazione utilizza una fotocamera da appena 1 megapixel (immaginiamodi voler inviare una foto da 10 megapixel).

L’immagine è troppo grande, è necessario ridurre la sua dimensione e perfarlo si utilizzano degli appositi algoritmi di compressione.

1in realtà nelle macchine fotografiche è presente un numero inferiore di sensori.2al momento vengono comunemente usate connessioni nettamente superiori al 56k, ma il

problema si pone comuqnue, ad esempio: nel momento in cui si decide di inviare un album dacentinaia di foto.

3

Capitolo 1

Algoritmi di Compressione

Nella compressione delle immagini vengono usati algoritmi di tipo lossy1, sfrut-tando il fatto che la perdita d’informazione non è sempre visibile.

1.1 Metodo dell’accettaSi tratta di un algoritmo molto semplice e con un rapporto di compressionemolto alto, l’unico difetto è che viene persa troppa informazione, infatti essoopera eliminando tutte le righe e tutte le colonne pari (o dispari) dalla matricedi riferimento, riducendo la risoluzione dell’immagine.

1.2 Rappresentazione Funzionale

Figura 1.1: Grafico della riga

Osserviamo una riga della matrice, sela rappresentiamo su un grafico in cuil’ascissa è la posizione nella riga e l’or-dinata è l’intensità associata a quellaposizione (fig. 1.1), possiamo cercaredi trovare una funzione che passi per ipunti del grafico. Per immagini moltosemplici (un raggio di luce obliquo inun ambiente buio) si può usare la fun-zione della retta (con f(x) si intendel’intensità nel punto x):

f(x) = ax+ b

Questo non è un caso realistico, infatti la maggiorparte delle immagini sonomolto più complesse di un raggio di luce.

1esistono due tipi di compressione: lossy (con perdita d’informazione) e lossless (senzaperdita d’informazione).

4

Potremmo provare un approccio polinomiale con grado n non maggiore delnumero di punti della riga della matrice:

f(x) = a0 + a1x1 + a2x

2 + a3x3 + :::+ an�1x

n�1

Purtroppo i polinomi non vanno bene, un grado molto alto porta a funzionimolto oscillanti e a problemi di stabilità numerica.

Un altra soluzione potrebbe essere la combinazione lineare di funzioni notemoltiplicate per alcuni parametri:

f(xi) = a0'0(xi) + a1'1(xi) + a2'2(xi) + :::+ an�1'n�1(xi)

Una volta scelte le 'i, essendo le f(xi) note, è necessario calcolare i coeffi-cienti, ad esempio, nel primo pixel l’intensità è:

f1 = a0'0(x1) + a1'1(x1) + :::+ an�1'n�1(x1)

nel secondo:

f2 = a0'0(x2) + a1'1(x2) + :::+ an�1'n�1(x2)

...

...

...

nell’n-esimo:

fn = a0'0(xn) + a1'1(xn) + :::+ an�1'n�1(xn)

Si ottiene un sistema di equazioni lineari a n incognite che risolto fornisce lasoluzione.

Un problema, da non sottovalutare, di questa rappresentazione è la scel-ta delle 'i, infatti una scelta sbagliata potrebbe portare ad una complessi-tà dell’algoritmo molto elevata oppure ad una rappresentazione non correttadell’informazione di partenza.

5

Capitolo 2

Trasformata di Fourier

L’idea è di utilizzare le funzioni sin e cos per le 'i, ottenendo:

T (x) = a0 cos(0�x)+a1 cos(2�x)+a2 cos(4�x)+: : :+b0 sin(0�x)+b1 sin(2�x)+b2 sin(4�x)+: : :

Possiamo scrivere la formula di sopra in modo compatto:

T (x) =X

aj cos(2j�x) + bi sin(2j�x) (2.1)

Osservazioni: Proviamo a risolvere la 2.1 per T (0) e per T (1).Ricordando che sin(2�) = 0 e che cos(2�) = 1, è facile verificare che:

� sin(2j�0) = sin(2j�1) = 0;

� cos(2j�0) = cos(2j�1) = 1;

da cui otteniamo:

� T (0) = a0 + a1 + a2 + a3 + : : :+ 0;

� T (1) = a0 + a1 + a2 + a3 + : : :+ 0;

di conseguenza agli estremi dell’intervallo [0; 1] abbiamo gli stessi risultati.Proviamo a risolverla ora per T (x+ 1):

T (x+ 1) =X

aj cos(2j�(x+ 1)) + bj sin(2j�(x+ 1)) =

=X

aj cos(2j�x+ 2j�) + bj sin(2j�x+ 2j�)

sapendo che:

� cos(h+ 2k�) = cos(h);

� sin(h+ 2k�) = sin(h);

otteniamo:

6

Figura 2.1: Periodo della trasformata di Fourier

T (x+ 1) =X

aj cos(2j�x+ 2j�) + bj sin(2j�x+ 2j�) =

=X

aj cos(2j�x) + bj sin(2j�x) = T (x)

Da queste osservazioni notiamo che la formula 2.1 è di periodo 1, ciò significache la trasformata si ripete a destra e a sinistra dell’intervallo [0; 1) (fig. 2.1).

2.1 Sintesi e Analisi di FourierLa formula che esprime l’analisi di Fourier è la seguente:

TN (xk) =

NXj=0

aj cos(2j�xk) + bj sin(2j�xk) (2.2)

dove:

� 2N è il numero di pixel dell’immagine;

� xk è il pixel, essi sono equidistanti nell’intervallo [0; 1], quindi:

xk =k

2N0 � k � 2N

per rappresentare un’immagine in questo modo il nostro obbiettivo è calcolaregli aj e bj .

Osservazioni: Il numero totale di pixel viene rappresentato con 2N , mentrela sommatoria va da j = 0 a N , questa scelta è necessaria per poter calcolarei parametri (aj e bj), infatti abbiamo due parametri per ogni j, quindi 2N + 2parametri.

In base a quanto detto sopra, per poter risolvere il sistema e trovare i 2N+2parametri, dovremmo avere 2N + 2 equazioni, ricordando che il polinomio èperiodico in [0; 1), contiamole:

7

TN (x0) = f0 =

NXj=0

aj cos(2j�0

2N) + bj sin(2j�

0

2N) x0 =

0

2N

f1 =

NXj=0

aj cos(2j�1

2N) + bj sin(2j�

1

2N) x1 =

1

2N

...

...

...

f2N�1 =

NXj=0

aj cos(2j�2N � 1

2N) + bj sin(2j�

2N � 1

2N) x2N�1 =

2N � 1

2N

quindi abbiamo 2N equazioni.C’è un problema: abbiamo 2N+2 parametri da calcolare e solo 2N equazioni,

ma osserviamo meglio il sistema:

� per b0 abbiamo:

b0 sin2 � 0�k2N

= b0 sin 0;

� per bN abbiamo:

bN sin2N�k

2N= bN sin�k;

Ricordando che sin 0 = 0 e sin�k = 0, si nota che nei casi contenenti b0 e bNabbiamo dei risultati indefiniti, dato che vengono moltiplicati per 0, per cuipossiamo escluderli dalla formula.

Ricapitolando, abbiamo 2N + 2 incognite, ma solo 2N equazioni, abbiamoanche visto però che non è necessario calcolare b0 e bN , di conseguenza il numerodi incognite è sceso a 2N , rendendo il sistema risolvibile.

2.1.1 Risoluzione del sistemaIl sistema è ortogonale:

� P cos 2j�xk � cos 2l�xk =

(0 se j 6= l

6= 0 se j = l;

8

� P sin 2j�xk � sin 2l�xk =

(0 se j 6= l

6= 0 se j = l;

� P sin 2j�xk � cos 2l�xk = 0;

ciò ci consente di risolverlo analiticamente, senza usare costosi metodi numerici:8>><>>:aj = 1

2N �2N�1Pk=0

fk cos 2j�xk

bj = 12N �

2N�1Pk=0

fk sin 2j�xk

(2.3)

dove fk = TN (xk).La 2.3 è l’analisi di Fourier.

Osservazioni: Assumendo che i sin e i cos siano già stati calcolati, analizziamoil costo computazionale della 2.3:

� dobbiamo fare prima 2N moltiplicazioni (fk � cos 2j�xk);� in seguito abbiamo 2N somme;

� 1 moltiplicazione ( 12N );

� tutto ciò per 2N volte (le operazioni vanno ripetute per tutti gli aj e bj);

� per un totale di 2N(2N + 2N + 2) = 8N2 + 4N = O(N2).

2.2 Trasformata Veloce di Fourier (FFT)

Riepilogo sui Numeri ComplessiUn numero complesso è rappresentato da una parte reale e da una parte imma-ginaria:

(a+ {b)

Proprietà importanti

{ è l’unità immaginaria:

{ =p�1 {2 = �1

La somma di due numeri complessi è così definita:

(a+ {b) + (c+ {d) = a+ {b+ c+ {d = (a+ c) + {(b+ d)

La moltiplicazione:

9

(a+ {b) � (c+ {d) = ac+ {ad+ {bc+ {2bd =

= ac+ {ad+ {bc� bd =

= (ac� bd) + {(ad+ bc)

Esponenziale:

e{x = cosx+ { sinx e�{x = cos(�x) + { sin(�x) = cosx� { sinx

Formule d’Eulero:

e{x + e�{x = 2 cosx cosx =e{x + e�{x

2(2.4)

e{x � e�{x = 2{ sinx sinx =e{x � e�{x

2{(2.5)

Infine:

1

{= �{

2.2.1 Aumentare la compattezza dell’analisi di Fourier (2.2)Partendo dalla 2.2:

TN (xk) =

NXj=0

aj cos(2j�xk) + bj sin(2j�xk)

e utilizzando le forumule di Eulero (2.4, 2.5), otteniamo:

TN (xk) =

NXj=0

aj � e{2j�xk + e�{2j�xk

2+ bj � e

{2j�xk � e�{2j�xk

2i=

=

NXj=0

aj � e{2j�xk + e�{2j�xk

2� {bj � e

{2j�xk � e�{2j�xk

2=

=

NXj=0

aje{2j�xk + aje

�{2j�xk � {bje{2j�xk + {bje

�{2j�xk

2

mettendo in evidenza rispetto a e{ e e�{:

TN (xk) =

NXj=0

e{2j�xk(aj � {bj)

2� e�{2j�xk (aj + {bj)

2=

10

=

NXj=�N

cj � e{2j�xk (2.6)

con:

cj =

8>>><>>>:

aj�{bj2 0 < j < N

aj+{bj2 �N < j < 0

aN2 j = N o j = �Na0 j = 0

La sintesi di Fourier è data da:

Cj =1

N�N�1Xk=0

fk � e�{2j�k

N (2.7)

2.2.2 Suddividere la trasformataPartendo dalla 2.7, vogliamo scrivere la somma come coppie di termini:

Cj =1

N�N2 �1Xk=0

fk � e�{2j�k

N + fk+N2� e

�{2j�(k+N2)

N =

=1

N�N2 �1Xk=0

fk � e�{2j�k

N + fk+N2� e�{2j�kN � e

�{2j� N2

N

ponendo w = e�{2�N :

Cj =1

N�N2 �1Xk=0

fk � wjk + fk+N2� wjk � wjN2 (2.8)

analizziamo il termine wjN2 :

wjN2 = (e�{2�N )j

N2 = e

�{2�j N2

N =

e�{2�j

2 = e�{�j

ricordando che sin�k = 0 e che e{x = cosx+ { sinx, abbiamo:

e�{�j = cos(��j) + { sin(��j) = cos�j � { sin�j = cos(�j)

cos�j =

(1 per j pari

�1 per j dispari

11

quindi:cos�j = (�1)j

Possiamo riscrivere la 2.8 così:

Cj =1

N�N2 �1Xk=0

fk � wjk + fk+N2� wjk � (�1)j

A questo punto possiamo suddividere due casi, a seconda che j sia pari odispari: 8>>><

>>>:C2m1

= 1N�N2 �1Pk=0

(fk + fk+N2)w2m1k

C2m1+1 = 1N�N2 �1Pk=0

(fk � fk+N2)w(2m1+1)k

=

=

8>>><>>>:C2m1

= 1N�N2 �1Pk=0

(fk + fk+N2)w2m1k

C2m1+1 = 1N�N2 �1Pk=0

(fk � fk+N2)wk � w2m1k

ponendo ora:

w2 = w1; (fk+fk+N2) = f

(1)k ; (fk�fk+N

2)wk = f

(1)

k+N2

;1

N=

1

2� 1N2

otteniamo: 8>>><>>>:C2m1

= 12 [

1N2

�N2 �1Pk=0

f(1)k � wm1k

1 ]

C2m1+1 = 12 [

1N2

�N2 �1Pk=0

f(1)

k+N2

� wm1k1 ]

(2.9)

Osservazioni: Assumendo di aver già calcolato gli esponenziali, analizziamola complessità della 2.9:

� per i numeri pari:

– occorre fare N2 somme per calcolare gli f (1)k ;

– in seguito vanno fatte N2 moltiplicazioni;

– in più altre N2 somme per risolvere la sommatoria;

– 2 moltiplicazioni;– il tutto ripetuto per N

2 volte;

– per un totale di N2 (

N2 + N

2 + N2 + 2) = 3

4 �N2 +N ;

12

� per i numeri dispari:

– occorre fare N2 + N

2 per calcolare gli f (1)k+N

2

;

– le restanti operazioni sono identiche a quelle per i numeri pari;

– per un totale di N2 (

N2 + N

2 + N2 + N

2 + 2) = N2 +N ;

Il numero totale di operazioni necessarie è: 34 �N2 +N +N2 +N = O(N2), di

conseguenza non abbiamo guadagnato molto rispetto alla 2.3.

2.2.3 Velocizzare la trasformata di Fourier (solo se N = 2L)

Confrontando la 2.7 con la 2.9, si nota che sono molto simili tra loro, ambeduesono delle sommatorie tra valori (fk, f

(1)k ) moltiplicati per degli esponenziali, di

conseguenza si può ulteriormente suddividere la 2.9, ottenendo un’altra formuladella stessa forma a sua volta suddividibile.

Si può continuare questo procedimento fino ad eliminare completamente lasommatoria, infatti all’inizio si avranno N somme da N termini, in seguito N

2

somme da N2 termini e così via, graficamente:

nN somme da N termini )

(N2 somme da N

2 termini )N2 somme da N

2 termini )

8>>><>>>:

N4 somme da N

4 termini � � �N4 somme da N



4 termini � � �

esempio con 8 punti:

0BBBBBBBBBB@

f0f1f2f3f4f5f6f7

1CCCCCCCCCCA)

0BBBBBBBBBBBBBB@

f0 + f4 = f(1)0

f1 + f5 = f(1)1

f2 + f6 = f(1)2

f3 + f7 = f(1)3

� � � � � � � � � � � �f0 � f4 = f

(1)4

(f1 � f5) � w = f(1)5

(f2 � f6) � w2 = f(1)6

(f3 � f7) � w3 = f(1)7

1CCCCCCCCCCCCCCA

)

0BBBBBBBBBBBBBBBBBB@

f(1)0 + f

(1)2 = f

(2)0

f(1)1 + f

(1)3 = f

(2)1

� � � � � � � � � � � �(f

(1)0 � f

(1)2 ) = f

(2)2

(f(1)1 � f

(1)3 ) � w1 = f

(2)3

� � � � � � � � � � � �f(1)4 + f

(1)6 = f

(2)4

f(1)5 + f

(1)7 = f

(2)5

� � � � � � � � � � � �(f

(1)4 � f

(1)6 ) = f

(2)6

(f(1)5 � f

(1)7 ) � w1 = f

(2)7

1CCCCCCCCCCCCCCCCCCA

)

0BBBBBBBBBBBBBBBBBBBBBBBBBB@

f(2)0 + f

(2)1 = C 0

0

� � � � � � � � � � � �f(2)0 � f

(2)1 = C 0

4

� � � � � � � � � � � �f(2)2 + f

(2)3 = C 0

2

� � � � � � � � � � � �f(2)2 � f

(2)3 = C 0

6

� � � � � � � � � � � �f(2)4 + f

(2)5 = C 0

1

� � � � � � � � � � � �f(2)4 � f

(2)5 = C 0

5

� � � � � � � � � � � �f(2)6 + f

(2)7 = C 0

3

� � � � � � � � � � � �f(2)6 � f

(2)7 = C 0

7

1CCCCCCCCCCCCCCCCCCCCCCCCCCA

13

Osservazioni: Analizziamo ora il numero di operazioni necessarie per calco-lare tutti i C 0

k:

� è necessario fare N somme per ogni passaggio;

� per un totale di log2N passaggi;

di conseguenza il costo complessivo dell’algoritmo è O(N � log2N), un enormeguadagno rispetto agli O(N2) precedenti.

Si nota che con questo procedimento viene perso l’ordine dei coefficienti(es: il quinto coefficiente che calcoliamo è C 0

1 invece di C 0

4), per ricostruire taleordine: si prende la rappresentazione in binario fino a k e la si riflette, come sefosse allo specchio.

Con k = 8 (come nell’esempio precedente) sarà così:0BBBBBBBBBB@

0 0001 0012 0103 0114 1005 1016 1107 111

1CCCCCCCCCCA()

0BBBBBBBBBB@

000 0100 4010 2110 6001 1101 5011 3111 7

1CCCCCCCCCCA

Attenzione: Nell’esempio precedente, i C 0

k non sono i coefficienti della tra-sformata di Fourier, questo perché la suddivisione della 2.9 viene fatta solo sullaparte tra parentesi quadre, quindi:

� è necessario moltiplicare 12 per ogni passaggio;

� avendo un totale di log2N passaggi, otteniamo che bisogna moltiplicare iltotale per 1

2log2 N= 1

N;

di conseguenza:

Ck =C 0

k

N

stavolta i Ck sono esattamente i coefficienti della trasformata di Fourier.Da notare che l’aggiunta della moltiplicazione per 1

Nnon cambia la comples-

sità asintotica dell’algoritmo, questo perché:

� sono necessari O(N � log2N) passaggi per calcolare i C 0

k;

� vengono fatti altri N passaggi per i Ck;

� per un totale di O(N � log2N) +N = O(N � log2N).

È necessario tener presente che questo procedimento può essere applicato solose N è una potenza di 2.

14

2.3 ConsiderazioniGrazie alla sua bassissima complessità, la trasformata di Fourier viene imple-mentata a livello hardware, consentendo di fare la compressione di immagini intempo reale.

2.3.1 Ridurre la dimensione dell’immagineUna volta applicata la Trasformata di Fourier abbiamo un elevato numero dicoefficienti (uno per ogni punto), però non tutti i coefficienti sono necessari perottenere una rappresentazione accettabile dell’immagine, infatti ce ne sarannomolti con valori vicini allo 0 (es: 0; 001), la cui eliminazione non porta ad unvisibile degrado dell’immagine.

In altre parole, vengono eliminati i coefficienti “inutili”, riducendo drastica-mente il numero di coefficienti che vengono salvati.

Il motivo per cui possiamo eliminare dei coefficienti essendo sicuri di nonvariare troppo il risultato della trasformata è dato dal teorema di Parseval:

ˆT 2(x) =

Xa2k + b2k (2.10)

che dice che l’energia del segnale è pari alla somma dei quadrati dei coeffi-cienti.

Osservando la 2.10, notiamo anche che i coefficienti con un valore maggioredi 1 contribuiscono all’energia molto di più di quelli compresi tra 0 e 1, adesempio: se abbiamo un coefficiente pari a 1

1:000 il suo quadrato è 11:000:000 ,

quindi 11:000 < 1

1:000:000�= 0.

2.3.2 Trasformata di Fourier con numero di punti diversoda 2

L

In letteratura esistono trasformate di Fourier basate su potenze diverse da 2,però tali algoritmi non possono essere usati a livello hardware, quindi vengonousati solo a livello software.

Esistono diverse soluzioni al problema, implementabili via hardware, cioè:

� Aumentare la dimensione dell’immagine fino a raggiungere la potenza di2 più vicina (ad es: aggiungendo degli 0 a fine immagine).

� Suddividere l’immagine in sezioni di dimensione pari ad una potenza di2 e applicare la trasformata su ogni singola sezione (questo approccio èutilizzato dall’algoritmo JPEG).

2.3.3 Trasformata di Fourier nell’algoritmo JPEGAbbiamo visto che per aumentare la velocità della trasformata di Fourier siutilizzano le formule di Eulero (2.4, 2.5), che si basano sulle proprietà dei numericomplessi, in particolare:

15

Figura 2.3: Segnale Riflesso

eix = cosx+ { sinx

Figura 2.2: Segnale Base

tuttavia, lavorare con la parte immaginaria diun numero complesso può rappresentare dei pro-blemi, questa situazione sconveniente viene risol-ta dall’algoritmo JPEG utilizzando la trasformatacoseno, ovvero una variante della trasformata diFourier che utilizza solo funzioni coseno.

L’eliminazione delle funzioni seno dalla tra-sformata è effettuata basandosi sulle caratteristi-che delle funzioni trigonometriche.

Sappiamo che la funzione coseno è pari1(cosx = cos(�x)), mentre la funzione seno è di-spari2 (sin(�x) = � sinx), quindi il “trucco” pereliminare le funzioni seno è rendere il segnale del-l’immagine simmetrico rispetto all’asse delle ordinate, ciò è ottenuto facilmenteriflettendo i punti del segnale prima di eseguire la trasformata, ricapitolando:

1. ottengo il segnale (fig. 2.2);

2. rifletto a sinistra il segnale (fig. 2.3);

3. applico la trasformata di Fourier.

A questo punto siamo sicuri che la trasformata utilizzerà solo funzioni cosenodato che deve rappresentare un segnale derivante da una funzione pari.

Per quanto riguarda la dimensione dell’immagine, l’algoritmo JPEG suddivi-de l’immagine in blocchi da 8�8 e applica la trasformata su di essi (dopo averliriflessi), di conseguenza la dimensione di un immagine deve essere un multiplodi 8.

1una funzione è pari se f(x) = f(�x), queste funzioni sono simmetriche rispetto all’assedelle ordinate.

2una funzione è dispari se f(�x) = �f(x), queste funzioni sono simmetriche rispettoall’origine.

16

2.3.4 Problemi della Trasformata di FourierLa trasformata è ottima per individuare le frequenze all’interno di un segnale,ma non ci dice nulla riguardo alla loro posizione, ciò la rende utilizzabile soloper segnali stazionari (segnali con frequenze periodiche che non variano con iltempo).

Figura 2.4: Funzione di Gauss

Il motivo è che la funzione coseno è de-finita su tutto l’intervallo del segnale, [0; 1),quindi non può fornire nessuna informazionedi localizzazione (come varia nel tempo) su diesso.

Una possibile soluzione a questo problemaè l’utilizzo della STFT (Trasformata di Fou-rier Short-Time), che consiste nel considerareil segnale localmente e in seguito applicare latrasformata.

La localizzazione viene fatta moltiplican-do il segnale per alcune funzioni, che possono essere (ad esempio):

� la funzione box (vale 1 in un intervallo 0 altrimenti), il problema di questafunzione è che porta ad una rappresentazione con discontinuità (fig. 2.5);

� la funzione “campana” (funzione di Gauss), fig. 2.4, in questo caso ilproblema è la scelta della dimensione della campana, una campana troppostretta porta ad un’individuazione poco precisa delle frequenze, mentreuna troppo larga non risolve il problema della localizzazione del segnale.

Figura 2.5: Applicazione di una funzione box

17

Capitolo 3

Trasformata Wavelet Discreta

La Trasformata di Fourier presenta dei problemi di localizzazione, ciò signifi-ca che per alcune immagini si ottengono dei risultati poco gradevoli, un’altratrasformata applicabile alle immagini è la Trasformata Wavelet.

La Trasformata di Fourier è una combinazione di funzioni lineari note:

f(xi) = a0'0(xi) + a1'1(xi) + a2'2(xi) + :::+ an�1'n�1(xi)

e come funzioni base, 'i, si utilizzano le funzioni sin e cos, ma queste nonsono le uniche funzioni utilizzabili in questo campo, nel caso della TrasformataWavelet si utilizzano le funzioni wavelet1, queste sono delle funzioni molto inte-ressanti che hanno la caratteristica di partire da 0, oscillare e ritornare a 0, essepossono essere sia molto localizzate, sia molto regolari.

Le funzioni wavelet vengono utilizzate tramite l’applicazione di due filtriwavelet:

� uno definito da L coefficienti hk, 0 � k < L, detto filtro passabasso, chesi riferisce a una funzione di scaling;

� un’altro definito da L coefficienti gk, 0 � k < L, detto filtro passaalto, chesi riferisce alla funzione wavelet e i cui elementi sono determinati a partiredagli hk secondo la seguente relazione:

gk = (�1)k � hL�k�1; 0 � k < L;

3.1 Esempio PraticoCome per la Trasformata di Fourier, rappresentiamo il segnale dell’immagine suun grafico, supponiamo di avere le seguenti coppie di punti (x; y):

(0; 1); (1; 3); (2; 3); (3; 2); (4; 4); (5; 5); (6; 4); (7; 3):

1In particolare, nella Trasformata Wavelet Discreta si utilizzano soprattutto le wavelet diHaar e quelle di Daubechies.

18

Fissato L = 2 possiamo scegliere di usare una funzione box, come filtro pas-sabasso, e una wavelet di Haar (fig: 3.1.1), come filtro passaalto, con i seguentifiltri:

(h0 =1

2; h1 =

1

2) (g0 =

1

2; g1 = �1

2)

3.1.1 Primo passo della trasformata

Figura 3.1: Wavelet di Haar

Dobbiamo applicare alle coppie di punti ilfiltro passabasso:

� data una generica coppia di punti[(x0; y0); (x1; y1)], che rappresenta unintervallo [x0; x1+1), l’applicazione con-siste nel generare una nuova coppia dipunti dove:x = x0+x1+1

2 e y = y0 � h0 + y1 � h1;� avremo quindi:

(0; 1) (1; 3) (2; 3) (3; 2) (4; 4) (5; 5) (6; 4) (7; 3) valori del segnale

& . & . & . & .(h0; h1) (h0; h1) (h0; h1) (h0; h1) filtro passabasso

# # # #(1; 2) (3; 52 ) (5; 92 ) (7; 72 ) (x; ci)

Per il filtro passaalto avremo:

� analogamente all’applicazione del filtro passabasso dobbiamo generare unanuova coppia di punti dove: x = x0+x1+1

2 e y = y0 � g0 + y1 � g1;� quindi:

(0; 1) (1; 3) (2; 3) (3; 2) (4; 4) (5; 5) (6; 4) (7; 3) valori del segnale

& . & . & . & .(g0; g1) (g0; g1) (g0; g1) (g0; g1) filtro passaalto

# # # #(1;�1) (3; 12 ) (5;� 1

2 ) (7; 12 ) (x; di)

Osservazioni: I filtri vengono applicati solo alla y delle coppie di punti.Il filtro passabasso genera un nuovo insieme di punti chiamati coefficienti di

approssimazione o di scaling (ci).L’applicazione del filtro passabasso porta ad un’approssimazione del segnale

di partenza, però non consente di fare il procedimento inverso, per questo motivoviene applicato anche il filtro passaalto ottenendo i coefficienti di dettaglio (di),che consentono di ricostruire il segnale di partenza partendo dai coefficienti diapprossimazione.

In pratica, dopo questo primo passo, otteniamo:

19

� il segnale di partenza approssimato (ci);

� il dettaglio perso (di).

Si nota che sono necessarie circa N operazioni per ogni filtro, di conseguenza ilcosto per il primo passo è circa: 2N = O(N).

L’applicazione dei filtri è detta convoluzione.

Rappresentazione Matriciale: È possibile scrivere l’applicazione dei filtriwavelet sotto forma di prodotto tra matrici, ad esempio, se si hanno 8 punti:0

BBBBBBBBBB@

c0d0c1d1c2d2c3d3

1CCCCCCCCCCA

=

0BBBBBBBBBB@

h0 h1 0 0 0 0 0 0g0 g1 0 0 0 0 0 00 0 h0 h1 0 0 0 00 0 g0 g1 0 0 0 00 0 0 0 h0 h1 0 00 0 0 0 g0 g1 0 00 0 0 0 0 0 h0 h10 0 0 0 0 0 g0 g1

1CCCCCCCCCCA�

0BBBBBBBBBB@

y0y1y2y3y4y5y6y7

1CCCCCCCCCCA

chiameremo la matrice che consente la convoluzione W .La rappresentazione matriciale non è utile dal punto di vista computaziona-

le2, ma ci consente di notare che la trasformata inversa può essere fatta sem-plicemente applicando l’inversa della matrice W al vettore dei coefficienti diapprossimazione e dettaglio.

Inoltre si nota che, se la matriceW è ortogonale3, la sua inversa è uguale allasua trasposta4, W�1 = WT , quindi è possibile fare il passo inverso applicandoWT ai coefficienti di scaling e di dettaglio:0

BBBBBBBBBB@

y0y1y2y3y4y5y6y7

1CCCCCCCCCCA

=

0BBBBBBBBBB@

h0 g0 0 0 0 0 0 0h1 g1 0 0 0 0 0 00 0 h0 g0 0 0 0 00 0 h1 g1 0 0 0 00 0 0 0 h0 g0 0 00 0 0 0 h1 g1 0 00 0 0 0 0 0 h0 g00 0 0 0 0 0 h1 g1

1CCCCCCCCCCA�

0BBBBBBBBBB@

c0d0c1d1c2d2c3d3

1CCCCCCCCCCA

Attenzione: Perché la matrice W sia ortogonale i filtri della funzione boxdovrebbero essere:

(h0 =1p2; h1 =

1p2) (g0 =

1p2; g1 = � 1p

2)

2Il costo di un prodotto tra matrici è O(N2), mentre abbiamo visto che è possibile fare ilprimo passo spendendo O(N).

3Una matrice A è ortogonale se A �AT = I, quindi se e solo se le sue righe e le sue colonneformano tra loro una base ortonormale.

4La trasposta di una matrice A, AT , è la matrice stessa, in cui le colonne diventano lerighe e le righe diventano le colonne.

20

per semplificare i calcoli non sono stati usati in questo esempio.

3.1.2 Passi successivi della trasformataRicapitolando, al primo passo della trasformata abbiamo la seguente situazione:0

BBBBBBBBBBBBBBBBBBBBBBBBBB@

c0 = 2

d0 = �1

c1 = 52

d1 = � 12

c2 = 92

d2 = � 12

c3 = 72

d3 = 12


) riordinando )


c0 = 2

c1 = 52

c2 = 92

c3 = 72

� � � � � � � � �d0 = �1

d1 = � 12

d2 = � 12

d3 = 12


Sappiamo inoltre che i coefficienti di approssimazione forniscono il segnaleapprossimato a un numero inferiore di punti, di conseguenza è possibile applicarenuovamente la matrice di convoluzione sui soli ci:0

BBB@c(1)0

d(1)0

c(1)1

d(1)1

1CCCA =

0BB@

h0 h1 0 0g0 g1 0 00 0 h0 h10 0 g0 g1

1CCA �

0BB@

c0c1c2c3

1CCA

avremo quindi:

21

0BBBBBBBBBBBBBBBBBBBBBBBBBBB@

c(1)0 = 9

4

d(1)0 = � 1

4

c(1)1 = 4

d(1)1 = 1

2� � � � � � � � �d0 = �1

d1 = � 12

d2 = � 12

d3 = 12

1CCCCCCCCCCCCCCCCCCCCCCCCCCCA

) riordinando )


c(1)0 = 9

4

c(1)1 = 4� � � � � � � � �d(1)0 = � 1

4

d(1)1 = 1

2� � � � � � � � �d0 = �1

d1 = � 12

d2 = � 12

d3 = 12


Abbiamo approssimato ulteriormente il segnale, quindi possiamo applicaredi nuovo la matrice di convoluzione, stavolta su due soli punti:

c(2)0

d(2)0

!=

�h0 h1g0 g1

��

c(1)0

c(1)1

!

A questo punto abbiamo un unico coefficiente d’approssimazione:0BBBBBBBBBBBBBBBBBBBBBBBBBBB@

c(2)0 = 25

8� � � � � � � � �d(2)0 = � 7

8� � � � � � � � �d(1)0 = � 1

4

d(1)1 = 1

2� � � � � � � � �d0 = �1

d1 = � 12

d2 = � 12

d3 = 12


quindi possiamo fermarci.

22

Osservazioni: Nel complessivo i passaggi della Trasformata Wavelet sonostati:

0BBBBBBBBBBBBBBBBBBBBBBBB@

y0 = 1

y1 = 3

y2 = 3

y3 = 2

y4 = 4

y5 = 5

y6 = 4

y7 = 3

1CCCCCCCCCCCCCCCCCCCCCCCCA

)


c0 = 2

c1 = 52

c2 = 92

c3 = 72

� � � � � � � � �d0 = �1

d1 = � 12

d2 = � 12

d3 = 12


)


c(1)0 = 9

4

c(1)1 = 4� � � � � � � � �d(1)0 = � 1

4

d(1)1 = 1

2� � � � � � � � �d0 = �1

d1 = � 12

d2 = � 12

d3 = 12


)


c(2)0 = 25

8� � � � � � � � �d(2)0 = � 7

8� � � � � � � � �d(1)0 = � 1

4

d(1)1 = 1

2� � � � � � � � �d0 = �1

d1 = � 12

d2 = � 12

d3 = 12


analizziamone la complessità per un numero arbitrario di punti:

� al primo passo i filtri vengono applicati su N punti;

� al secondo passo su N2 punti;

� al terzo su N4 e così via;

� quindi i filtri vengono applicati su N2i punti, dove i indica l’i-esimo passo.

Inoltre, ad ogni passo, il numero di operazioni necessario è il doppio del numerodi punti usati e vengono fatti un totale di logN passi.

Il numero di operazioni effettivo sarà:

logNXi=0

2N

2i= 2N �

logNXi=0

1

2i= 2N � 1�

12logN+1

1� 12

= 2N � 1� ( 12logN

� 12 )12

=

= 2N � (1� 1

2N) � 2 = 2N � (2N � 1

2N) � 2 = 2 � (2N � 1) =

= 4N � 2 = O(N)

quindi il calcolo della Trasformata Wavelet può essere eseguito in tempolineare O(N), che è inferiore a O(N � logN).

23

Figura 3.2: Wavelet di Daubechies da 4 e 12 punti

3.2 Wavelet di DaubechiesNella pratica hanno riscosso un enorme successo, per la Trasformata WaveletDiscreta, le funzioni di Daubechies (fig: 3.2).

Esse sono definite per un numero arbitrario di punti e vengono costruite inmodo che i filtri rispettino le seguenti proprietà:

L�1Pk=0

hk =p2;

L�1Pk=0

h2k = 1;

L�1Pk=0

gk = 0;L�1Pk=0

k � gk = 0:

(3.1)

che implicano:

L�1Xk=0

hk � gk = 0

garantendo un sistema ortogonale, quindi risolvibile analiticamente, senzacostosi metodi numerici.

Tra le proprietà indicate nella 3.1, quelle riguardanti le gk servono a garantireche i filtri approssimino un polinomio di grado p = L

2 e sono dette condizioni diapprossimazione.

3.2.1 Convoluzione con Wavelet da 4 puntiConsideriamo un segnale da 8 punti, indichiamo con gli apici il livello di tra-sformata a cui i coefficienti si riferiscono:

24

� Livello 3:0BBBBBBBBBBBB@

c(2)0

d(2)0

c(2)1

d(2)1

c(2)2

d(2)2

c(2)3

d(2)3

1CCCCCCCCCCCCA

=

0BBBBBBBBBB@

h0 h1 h2 h3 0 0 0 0g0 g1 g2 g3 0 0 0 00 0 h0 h1 h2 h3 0 00 0 g0 g1 g2 g3 0 00 0 0 0 h0 h1 h2 h30 0 0 0 g0 g1 g2 g3h2 h3 0 0 0 0 h0 h1g2 g3 0 0 0 0 g0 g1

1CCCCCCCCCCA�

0BBBBBBBBBB@

y0y1y2y3y4y5y6y7

1CCCCCCCCCCA

(3.2)

�c(2)0 c

(2)1 c

(2)2 c

(2)3

... d(2)0 d

(2)1 d

(2)2 d

(2)3

�� Livello 2:

0BBB@

c(1)0

d(1)0

c(1)1

d(1)1

1CCCA =

0BB@

h0 h1 h2 h3g0 g1 g2 g3h2 h3 h0 h1g2 g3 g0 g1

1CCA �

0BBB@

c(2)0

c(2)1

c(2)2

c(2)3

1CCCA (3.3)

�c(1)0 c

(1)1

... d(1)0 d

(1)1

... d(2)0 d

(2)1 d

(2)2 d

(2)3

�� Livello 1:

c(0)0

d(0)0

!=

�2h0 2h12g0 2g1

��

c(1)0

c(1)1

!(3.4)

�c(0)0

... d(0)0

... d(1)0 d

(1)1

... d(2)0 d

(2)1 d

(2)2 d

(2)3

�

Osservazioni: Nel nostro esempio abbiamo supposto che il segnale sia perio-dico, infatti la 3.2 equivale a:

0BBBBBBBBBBBB@

c(2)0

d(2)0

c(2)1

d(2)1

c(2)2

d(2)2

c(2)3

d(2)3

1CCCCCCCCCCCCA

=

0BBBBBBBBBB@

h0 h1 h2 h3 0 0 0 0 0 0g0 g1 g2 g3 0 0 0 0 0 00 0 h0 h1 h2 h3 0 0 0 00 0 g0 g1 g2 g3 0 0 0 00 0 0 0 h0 h1 h2 h3 0 00 0 0 0 g0 g1 g2 g3 0 00 0 0 0 0 0 h0 h1 h2 h30 0 0 0 0 0 g0 g1 g2 g3

1CCCCCCCCCCA�

0BBBBBBBBBBBBBBBB@

y0y1y2y3y4y5y6y7� � �y0y1

1CCCCCCCCCCCCCCCCA

25

lo stesso ragionamento vale per la 3.3.Per quanto riguarda la 3.4 sono stati utlizzati dei filtri ad-hoc, in quanto il

numero di punti rimanenti è inferiore alla dimensione dei filtri.L’aumento della dimensione dei filtri non cambia la complessità dell’algorit-

mo, che rimane O(N).

3.3 Confronto con la Trasformata di Fourier

3.3.1 LocalizzazioneAbbiamo visto che la Trasformata di Fourier non offre una localizzazione delsegnale, cosa che nella Trasformata Wavelet Discreta non avviene perché:

� il segnale viene considerato localmente (come nella STFT);

� le fuzioni Wavelet hanno la caratteristica di partire da 0 e tornare a 0;

� le funzioni Wavelet sono continue.

Dato che il segnale viene considerato localmente, può esserci un problema dicontinuità ai bordi, ciò può essere risolto in vari modi, ad esempio:

� considerando il segnale periodico (es: 3.2);

� utilizzando dei filtri differenti solo per i bordi (es: 3.4).

3.3.2 ComplessitàCome la Trasformata di Fourier, anche la Trasformata Wavelet Discreta, perla sua semplicità, può essere implementata su un hardware, consentendo dicalcolarla in tempo reale.

Ricordiamo che la Trasformata di Fourier Veloce ha complessità O(N �logN),mentre la trasformata Wavelet ha complessità O(N), quindi è più performante.

Un’altra cosa importante da notare è che nella Trasformata Wavelet è possi-bile fermarsi ad ogni passo ed essere comunque in grado di ricostruire il segnale,nella Trasformata di Fourier Veloce, invece, è necessario completare tutti i passiper avere dei coefficienti validi.

3.3.3 CompressioneAnalogamente alla Trasformata di Fourier, nella Trasformata Wavelet per com-primere basta eliminare i coefficienti di dettaglio “piccoli”.

26

Capitolo 4

Trasformata Bidimensionale

Finora abbiamo visto come applicare una trasformata su una riga della matriceche rapppresenta l’immagine, vediamo ora come applicarla su tutta la matrice:

� si applica la trasformata su tutte le righe della matrice, ottenendo unanuova matrice composta dai coefficienti prodotti;

� a questo punto, sulla nuova matrice, si applica la trasformata su tutte lecolonne.

Questo procedimento si può estendere anche a più di due dimensioni.In modo del tutto analogo si applica la trasformata inversa (prima sulle righe

e poi sulle colonne).

Osservazioni: Supponendo una matrice con N righe e M colonne si scelgonodelle basi per le x e per le y:�

'x0(x0) 'x1(x1) 'x2(x2) � � � 'xM�1(xM�1)

�'y0(y0) '

y1(y1) '

y2(y2) � � � '

yN�1(yN�1)

L’applicazione della trasformata prima sulle righe e poi sulle colonne produce

un prodotto incrociato tra le basi:

'xi (xi) � 'yj (yj)) 'i;j(xi; yj) i = 0; 1; : : : ;M � 1; j = 0; 1; : : : ; N � 1

tale prodotto è separabile, ciò significa che non è importante l’ordine con cuivengono fatte le operazioni (prima sulle righe e poi sulle colonne o viceversa),garantendo di poter fare la trasformata inversa senza preoccuparsi dell’ordinecon cui sono stati eseguiti i calcoli al momento della compressione.

27

Appuntidelcorsodi prof.ssa BiancamariadellaVecchia Trasformata di ... · Università“Sapienza”...

Documents

Transcript of Appuntidelcorsodi prof.ssa BiancamariadellaVecchia Trasformata di ... · Università“Sapienza”...