Indice - Angewandte Numerische Analysis - M15/Allgemeines ... · i capitoli introduttivi...

Indice

Introduzione alle note iii

1 Serie e trasformate di Fourier 11.1 Spazi di Hilbert . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21.2 Serie e trasformate di Fourier . . . . . . . . . . . . . . . . . . . . 41.3 L’algoritmo della Fast Fourier Transform . . . . . . . . . . . . . . 12

2 Teoria del campionamento 172.1 Relazioni tra le trasformate di Fourier continua e discreta . . . . 172.2 Teoria del campionamento . . . . . . . . . . . . . . . . . . . . . . 202.3 Sistemi ortogonali per L2(Rd): frame di Gabor . . . . . . . . . . 24

3 Analisi tempo-frequenza: frame di Gabor 273.1 Frame per spazi di Hilbert . . . . . . . . . . . . . . . . . . . . . . 283.2 Frame di Gabor . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34

3.2.1 Ridondanza della STFT . . . . . . . . . . . . . . . . . . . 343.2.2 Principio di indeterminazione . . . . . . . . . . . . . . . . 343.2.3 Frame di Gabor e il gruppo di Weyl-Heisenberg . . . . . 353.2.4 Condizioni necessarie per avere una frame di Gabor . . . 363.2.5 Condizioni sufficienti per avere una frame di Gabor . . . . 383.2.6 Frame di Gabor discrete . . . . . . . . . . . . . . . . . . . 41

3.3 Compressione di un segnale . . . . . . . . . . . . . . . . . . . . . 43

4 2D pattern matching: armoniche circolari 474.1 Descrizione del problema . . . . . . . . . . . . . . . . . . . . . . . 474.2 La rappresentazione matematica . . . . . . . . . . . . . . . . . . 494.3 La definizione del matching . . . . . . . . . . . . . . . . . . . . . 554.4 Un’ipotesi di parallelizzazione . . . . . . . . . . . . . . . . . . . . 584.5 Appendice al capitolo . . . . . . . . . . . . . . . . . . . . . . . . 60

5 Elementi della teoria delle ondine (wavelet) 655.1 Trasformate tempo-frequenza e tempo-scala . . . . . . . . . . . . 655.2 Frame di ondine . . . . . . . . . . . . . . . . . . . . . . . . . . . 68

5.2.1 Condizioni necessarie . . . . . . . . . . . . . . . . . . . . . 68

i

ii INDICE

5.2.2 Condizioni sufficienti . . . . . . . . . . . . . . . . . . . . . 695.3 Analisi multirisolutiva e ondine ortonormali . . . . . . . . . . . . 715.4 Ondine a supporto compatto . . . . . . . . . . . . . . . . . . . . 745.5 La trasformata wavelet discreta (DWT) . . . . . . . . . . . . . . 76

6 Problemi inversi con vincoli di sparsita 796.1 Problemi inversi lineari . . . . . . . . . . . . . . . . . . . . . . . . 806.2 Problemi inversi con vincoli di sparsita . . . . . . . . . . . . . . . 826.3 Algoritmi iterativi di soglia . . . . . . . . . . . . . . . . . . . . . 836.4 Algoritmo del gradiente proiettato . . . . . . . . . . . . . . . . . 85

6.4.1 Proiezioni su BR e operatori di soglia . . . . . . . . . . . 886.4.2 Indicazioni di convergenza . . . . . . . . . . . . . . . . . . 90

6.5 Due applicazioni rilevanti . . . . . . . . . . . . . . . . . . . . . . 916.5.1 Ricolorazione di immagini e restauro d’immagini . . . . . 926.5.2 Magnetoencefalografia (MEG) . . . . . . . . . . . . . . . 94

Introduzione alle note

Queste note sono una estensione delle lezioni [47] del corso di Matematica perl’analisi e il trattamento di immagini tenutosi presso l’Universita degli Studi diPadova nel A.A. 2002-2003. La prima edizione si limitava ad illustrare elemen-ti di base dell’analisi armonica numerica con un linguaggio piu moderno, conriferimenti al trattamento numerico di segnali ed immagini. In questa edizionevengono aggiunti argomenti piu avanzati, come nozioni della teoria delle ondine(wavelets) e recenti risultati di ricerca nel trattamento numerico di problemiinversi con vincoli di sparsita. Le applicazioni che vengono aggiunte seguono ariferirsi a problemi di ricostruzione di segnali ed immagini. Rimangono dunquei capitoli introduttivi sull’analisi di Fourier, la teoria del campionamento, e l’a-nalisi tempo-frequenza.

I prerequisiti per una comprensione adeguata delle note sono i corsi di AnalisiMatematica del biennio di ingegneria, matematica o fisica, un corso di introdu-zione alla Analisi Numerica (con particolare riferimento a elementi di algebralineare computazionale) ed elementi di Analisi funzionale (spazi di Hilbert edoperatori lineari, nozioni di serie e integrali di Fourier). Tuttavia, si e cercatodi rendere la lettura consistente ed agevole e rivolta all’utilizzo degli strumentimatematici, limitandosi a volte ad una descrizione puramente intuitiva, sebbe-ne sempre rigorosa ed accompagnata da riferimenti di approfondimento. Perlo svolgimento di alcuni esercizi al calcolatore si presuppone una conoscenza dibase dei linguaggi C/C++ o Matlab e verranno presentati dei tutorial comenotebook di Mathematica, che garantisce la possibilta di operare in modo suf-ficientemente intuitivo e di disporre di una semplice interfaccia grafica per lavisualizzazione dei risultati.

Lo scopo del corso e delle note e quindi di fornire una panoramica deglistrumenti di lavoro nell’ambito dell’analisi armonica numerica, una introduzioneal loro utilizzo e una conoscenza delle loro potenzialita per le applicazioni.

Analisi armonica numerica per immagini

Un’immagine puo venir rappresentata come una funzione matematica, che adogni punto associa un corrispondente valore di colore, spesso espresso come

iii

iv INTRODUZIONE ALLE NOTE

combinazione di colori fondamentali. Tale rappresentazione consente quindiuna interpretazione puramente matematica ed astratta di una immagine. Ildominio di un’immagine, interpretata come funzione matematica, puo esserepensato come continuo, ma al fine di poter eseguire una eleborazione numericaoccorrera definire una rappresentazione anche su un dominio discreto. Tali rap-presentazioni, continua e discreta, devono essere equivalenti, ossia si deve poterpassare dall’una all’altra senza perdita di informazione. Questo garantisce chel’elaborazione che eseguiamo numericamente corrisponda alla fine ad un tratta-mento dell’immagine anche come funzione definita nel continuo. Un processo didigitalizzazione di un segnale f e quindi l’associazione di f ad una successione

di dati→c= (cn(f))n∈N, rappresentanti f . Tali dati

→c potranno essere trattati

numericamente, inviati attraverso un mezzo o trasmessi mediante comunicazionisenza cavo. Tutte queste operazioni richiedono tuttavia due condizioni/principi:

(I) la rappresentazione digitale deve essere stabile: questo significa che piccoleperturbazioni sui dati sono legate a piccole perturbazioni delle funzioneoriginaria. Questa e una naturale richiesta di continuita;

(II) robustezza al rumore e tolleranza dell’errore: possibili errori importantisui dati non dovrebbero modificare troppo la rappresentazione di una fun-zione f . Infatti la trasmissione attraverso un mezzo, ma anche semplicioperazioni numeriche possono produrre errori e la loro propagazione.

Uno degli esempio piu semplici che si possono considerare come processo didigitalizzazione e quello di associare a una funzione continua f la sequenza deisuoi campioni su un reticolo discreto o insieme di nodi cn(f) = f(τn), doveτ > 0 e il passo di campionamento. Se si assume che f appartenga ad unospazio normato di funzioni, si chiami esso (B, ‖ · ‖B), e che la sequenza dei datiappartenga ad uno spazio normato di successioni (Y, ‖ · ‖Y ), una relazione diequivalenza e di continuita e data dalla equivalenza delle norme:

‖(cn(f))n∈N‖Y ≍ ‖f‖B, (1)

ove per “≍” si intende che esistono due costanti A,B > 0, indipendenti da f ,tali che

A‖(cn(f))n∈N‖Y ≤ ‖f‖B ≤ B‖(cn(f))n∈N‖Y .

Assumiamo che la funzione f venga affetta da un errore

f = f + δf.

Allora da (1) si ha che

‖f‖B ≍ ‖(cn(f) + cn(δf))n∈N‖Y ,

e ‖(cn(δf))n∈N‖Y ≤ A−1‖f − f‖B. Quindi con tale equivalenza di norme siassicura la validita del principio/condizione (I).

v

Per esempio se f ∈ Lp(Rd)1 per 1 < p <∞ e una funzione a banda limitata,cioe per cui la sua trasformata di Fourier ha supporto compatto2, allora esisteτ0 > 0 piccolo a sufficienza tale che per ogni 0 < τ ≤ τ0

‖(f(τn))n∈Zd‖ℓp ≍ ‖f‖Lp,

Questa equivalenza di rappresentazione continua e discreta/numerica/digitalesi realizza in particolare per la possibilita di riprodurre f dall’insieme dei suoicampioni. Infatti per un opportuna successione di funzioni a banda limitatafn ∈ L1(Rd), n ∈ Zd, si puo scrivere:

f =∑

n

f(τn)fn.

Cosı il processo di analisi o digitalizzazione (mappa dei coefficienti)

fC→ (cn(f))n∈N

puo essere invertito sulla sua immagine da un processo di sintesi o riproduzione(mappa di recupero)

(cn(f))n∈N

R→ f.

L’equivalenza di norme espressa in (1) e equivalente alla continuita (o allalimitatezza) degli operatori lineari di analisi e sintesi rispettivamente.

Quindi, data una immagine e una funzione f che la rappresenti, si puoconsiderare la decomposizione astratta in pacchetti di informazione numerici ediscreti:

fւ ↓ ց

... ck−1(f) ck(f) ck+1(f) ...

magari indipendenti o correlati tra loro. Ma una decomposizione allo scopo diessere “utile” e fedele (ossia che garantisca che le operazioni numeriche equival-gano a operazioni sulla rappresentazione continua) devono essere associate aduna riproduzione/sintesi:

... ck−1(f) ck(f) ck+1(f) ...ց ↓ ւ

f

E inoltre piuttosto intuitivo che se i pacchetti di informazione sono totalmenteindipendenti allora possibili errori sono individualmente interpretati da ciascunelemento cn(f) di informazione con possibili errori importanti locali/globali nella

1Per p ∈ [1,∞) si ha Lp(X) = f :R

X |f(x)|pdµ(x) < ∞ lo spazio delle funzioni µ-

misurabili su X e p-sommabili che munito della norma ‖f‖p =`R

X |f(x)|pdµ(x)´

1/pdefinisce

uno spazio normato. Se X e discreto si scrivera ℓp(X) anziche Lp(X).2Illustreremo nel seguito il concetto di trasformata di Fourier e studieremo quindi le funzioni

a banda limitata.

vi INTRODUZIONE ALLE NOTE

consequente sintesi/ricostruzione. D’altra parte, possibili correlazioni tra i pac-chetti di informazione possono distribuire gli errori su tutto l’insieme dei dati,magari riducendo cattive distorsioni locali, magari compensando completamen-te un errore prodotto ad un certo punto, recuperandolo in un punto successivo.Questo si realizza e puo essere intepretato come una intrinseca ridondanza dellarappresentazione numerica/digitale: la ridondanza di informazione e piuttostocomune in natura ed e fortemente raccomandabile ogni volta che errori si posso-no presentare e il grado di ridondanza dovrebbe quindi dipendere dalla frequenzadi tali possibili errori. Basti giusto pensare al DNA umano: e un lungo flussodi informazioni ove solo una minima parte e considerata effettivamente “utile”.Il resto e protezione e parte in qualche senso ridondante. Alcune modificazionipossono essere compensate, e, magari, i loro effetti negativi essere controllati oridotti. Questo corrisponde al principio/condizione (II) e alla seguente formu-lazione matematica:

per una medesima operazione di sintesi/ricostruzione R : Y → B defini-ta da un insieme di funzioni (fn)n∈N esistono molteplici operazioni di anali-si/decomposizione C : B → Y tali che IB = R C.

Questo equivale a dire che possono esistere (cn(f))n∈N e (dn(f))n∈N distintitali che

f =∑

n∈N

cn(f)fn =∑

n∈N

dn(f)fn

Quindi un processo di digitalizzazione ridondante puo offrire diverse possibilirealizzazioni e si puo quindi adattare la scelta di una realizzazione, in dipen-denza del particolare uso che se ne vuol fare. A prima vista, la ridondanzapuo considerarsi in contrasto con il concetto di compressione. Come si vedra inseguito, in realta un sistema ridondante di pattern (pacchetti di informazione)puo aumentare la possibilita di identificazione: e vero che due differenti patternpossono condividere parte del medesimo contenuto di informazione, ma uno puocorrispondere meglio alla data funzione: “piu grande e il dizionario, tante piufrasi si e in grado di costruire e tanto piu corte esse saranno”.

Il corso e le note descrivono quindi i fondamenti matematici dell’analisi ar-monica numerica attraverso il linguaggio moderno delle frame, i.e., espansionicomplete, stabile e ridondanti. Verranno mostrate infine applicazioni per larappresentazione di immagini, la compressione, l’identificazione di immagini ameno di trasformazioni affini, e il restauro di immagini mediante la risoluzio-ne di problemi inversi con vincoli di sparsita rispetto a basi e frame. Le notesono accompagnate da esercizi che si consiglia vivamente di svolgere per unamiglior comprensione e per sviluppare una capacita d’uso degli strumenti. Al-cuni esercizi suggeriti sono di carattere pratico, essenzialmente implementazionial computer.

vii

Gli elementi del corso e approfondimenti biblio-

grafici

Il corso si limita ad una analisi in spazi di Hilbert (essenzialmente L2, ℓ2) ede diviso in sei lezioni, i cui argomenti principali sono brevemente richiamati diseguito, assieme ai principali riferimenti bibliografici di approfondimento.

1. Dopo una breve introduzione sugli spazi di Hilbert, descriveremo il con-cetto di sviluppo in serie di Fourier, trasformata di Fourier e sua imple-mentazione discreta. L’algoritmo FFT (Fast Fourier Transform) [42, 43,64, 65, 85, 14, 24, 25, 55, 56].

2. Formula della somma di Poisson. Teoria del campionamento per funzionia banda-limitata e conversioni digitale-analogico. Stime degli errori dialiasing per funzioni non a banda-limitata. Deduzione di un primo esempiodi sistema tempo-frequenza (frame di Gabor) [64, 70, 45, 38].

3. Frame in spazi di Hilbert. Analisi tempo-frequenza e frame di Gabor,continue e discrete e implementazione dei duali. Applicazioni ai segnaliaudio e alla compressione di immagini [11, 13, 19, 22, 32, 39, 40, 41, 57,59, 61, 63, 71, 63, 72, 73, 75, 83, 84, 86, 2].

4. Generalizzazione delle frame di Gabor su geometria circolare: le armonichecircolari discrete. Rappresentazione di immagini secondo frame di armoni-che circolari, proprieta rispetto all’azione di rotazioni. Metodi ed algoritmidi calcolo per l’indentificazione di immagini. Discussione del caso appli-cativo della ricostruzione assistita dal computer di affreschi minutamenteframmentati [10, 54].

5. Analisi tempo-frequenza e tempo-scala, trasformata wavelet continua, fra-me di wavelet, wavelet ortonormali, e trasformata wavelet discreta [11, 19,62, 68, 69].

6. Problemi inversi lineari, problemi inversi con vincoli di sparsita rispetto aframe, algoritmi iterativi di soglia e di gradiente proiettato per la risolu-zione di problemi inversi con vincoli di sparista. Applicazioni alla ricostru-zione di immagini, in particolare alla ricolorazione di immagini e all’analisidi immagini mediche provenienti da magnetoencefalografia [20, 21, 1, 3, 5,6, 7, 8, 23, 26, 27, 28, 29, 30, 31, 34, 50, 51, 53, 66, 74, 76, 77, 78, 80, 81].

Concludo ringraziando il Dott. Fabio Marcuzzi per l’invito a tenere questelezioni come parte specialistica del suo corso presso l’Universita degli Studi diPadova A.A. 2006/2007, il Program in Applied and Computational Mathema-tics e il Department of Mathematics della Princeton University per l’ospitalitadatami durante la preparazione di queste note.

Massimo Fornasier Princeton, maggio 2007.

viii INTRODUZIONE ALLE NOTE

Capitolo 1

Serie e trasformate diFourier

“... il primo tra gli uomini di scienza europei ...”

scrisse Giuseppe Lodovico Lagrangia (Joseph-Louis Lagrange), descrivendo Joseph

Fourier.

S i vuole iniziare questa lezione cercando di mettere in evidenza l’impor-tanza che l’analisi armonica o di Fourier ha avuto nella storia della matematicamoderna e delle sue molteplici attuali applicazioni.

Si possono sintetizzare una varieta di funzioni complicate da pure sinusoidinello stesso modo con cui si produce un accordo di do maggiore premendo itasti (do, mi, sol) di un pianoforte. Una versione geometrica di questa ideae alla base del modello antico di Ipparco e Tolomeo per il moto planetario(Almagets, II secolo). E pero stato Joseph Fourier (http://www-gap.dcs.st-and.ac.uk/∼history/Mathematicians/Fourier.html) che sviluppo i moderni me-todi dell’uso delle serie e degli integrali trigonometrici, nello studio del flussodel calore nei solidi (Analytical Theory of Heat, 1815). Oggi, l’analisi di Fouriere una branca della matematica piuttosto evoluta e ha sviluppato un importantenumero di applicazioni.

Il concetto di funzione, di integrale di Riemann e di Lebesgue, i concetti diconvergenza puntuale, uniforme e in media, la teoria degli insiemi di Cantor,la teoria delle distribuzioni di Schwartz sono attribuite da alcuni ricercatori(per esempio S. Bochner, Fourier series came first, Amer. Math. Monthly, 86(1979), 197-199), in parte come conseguenze dello sforzo di valutare quando unarappresentazione di Fourier e valida e in che senso. Comunque, come vedremo,l’analisi di Fourier o armonica e un campo di ricerca molto ancora aperto, le sue

1

2 CAPITOLO 1. SERIE E TRASFORMATE DI FOURIER

estensioni in matematica pura costituiscano una branca molto ricca dell’analisie anche le possibilita applicative sono lontane dall’esaurirsi. Ne citiamo alcuneche hanno segnato la storia della scienza e della tecnologia:

• l’articolo matematico piu citatoJ. W. Cooley and J. W. Tukey, An algorithm for the machine computationof complex Fourier series, Math. Comp., 19 (1965), 297-301.

• Circa 3/4 dei premi Nobel in fisica sono stati assegnati a lavori svoltiusando strumenti e concetti di analisi di Fourier.F. N. Magill, ed., The Nobel Prize Winners — Physics, Vol 1-3, SalemPress. Englewood Cliffs, NJ. 1989.

• Francis Crick, James Watson e Maurice Wilkins hanno ottenuto il premioNobel per la medicina e la fisiologia nel 1962 per la scoperta della strutturamolecolare del DNA. Questo fu il primo esempio (a cui Huptmann e Karlesi sono ispirati) dell’uso dell’analisi di Fourier sui dati provenienti dalladiffrazione di raggi X .

• Herbert Hauptmann (un matematico) e Jerome Karle hanno condiviso nel1985 il premio Nobel per la chimica per aver mostrato come usare l’ana-lisi di Fourier per determinare la struttura di grandi molecole dai dati didiffrazione di raggi X.W. A. Hendrickson, the 1985 Nobel Prize in Chemistery, Science 231(1986), 362-364.

• Il numero enorme di prodotti tecnologici (CD-DVD, TV ad alta definizio-ne, telefoni digitali, strumenti per indagini biomediche (medical imaging)...) basati sulla teoria del trattamento di segnali digitali, essenzialmentebasata sull’analisi di Fourier.

Prima di addentrarci nello studio delle serie ed integrali di Fourier sarautile richiamare alcune nozioni relative agli spazi di Hilbert astratti. Per ledimostrazioni dei risultati si veda per esempio [42] Capitolo 5. Indichiamo comeulteriori approfondimenti [42, 43, 64, 65, 85, 14, 24, 25, 55, 56].

1.1 Spazi di Hilbert

Definizione 1.1.1. Sia H uno spazio vettoriale. Un prodotto scalare 〈u, v〉 euna applicazione di H×H a valori in C tale che:

(i) 〈au + bv, z〉 = a〈u, z〉+ b〈v, z〉 per ogni u,v,z ∈ H e a, b ∈ C.

(ii) 〈u, v〉 = 〈v, u〉 per ogni u,v ∈ H.

(iii) 〈u, u〉 ∈ R, 〈u, u〉 ≥ 0 per ogni u ∈ H e 〈u, u〉 6= 0 se u 6= 0.

1.1. SPAZI DI HILBERT 3

Ricordiamo che un prodotto scalare verifica la diseguaglianza di Cauchy-Schwarz

|〈u, v〉| ≤ 〈u, u〉1/2〈v, v〉1/2 ∀u, v ∈ H. (1.1)

e che ‖u‖ = 〈u, u〉1/2 definisce una norma per H.

Definizione 1.1.2. Uno spazio di Hilbert e uno spazio vettoriale H munito diun prodotto scalare 〈u, v〉 e che e completo rispetto alla norma 〈u, u〉1/2.

Salvo indicazioni esplicite, nel corso di questo capitolo con H si indicherasempre uno spazio di Hilbert.

Esempio 1.1.3. Si prongono alcuni esempi:

• Sia Ω ⊂ Rn. Lo spazio vettoriale L2(Ω) = f : Ω → C|∫Ω|f |2dx < ∞ e

uno spazio di Hilbert con il prodotto scalare

〈u, v〉 =

∫

Ω

u(x)v(x)dx.

• Sia Ωd ⊂ Zn. Lo spazio vettoriale ℓ2(Ωd) = f : Ωd → C|∑k∈Ωd|f(k)|2 <

∞ e uno spazio di Hilbert con il prodotto scalare

〈u, v〉 =∑

k∈Ωd

u(k)v(k)dx.

In particolare se |Ωd| = d <∞ allora ℓ2(Ωd) = Cd.

Definizione 1.1.4. Sia E uno spazio vettoriale topologico. Diremo E′ lo spazioduale topologico di E, cioe E′=ϕ di E a valori in C | ϕ lineare e continua.Dati due spazi di Hilbert H,K, l’insieme degli operatori lineari T : H → K,cioe T (au + bv) = aT (v) + bT (v), e che sono continui e uno spazio vettoriale

che munito della norma ‖|T |‖H→K ≡ supu∈H,u6=0‖Tu‖K

‖u‖He uno spazio normato

completo, cioe uno spazio di Banach che indicheremo con L(H,K). In effetti, lacontinuita dell’operatore T si dimostra essere equivalente al fatto che la quantita

supu∈H,u6=0‖Tu‖K

‖u‖Hsia finita.

In particolare se H = K porremo L(H) ≡ L(H,H) e se K = C allora si haL(H,K) = H′.

Teorema 1.1.5 (Teorema di rappresentazione di Riesz-Frechet). Sia H spaziodi Hilbert. Per ogni ϕ ∈ H′ esiste una unica f ∈ H tale che

〈ϕ, v〉 ≡ ϕ(v) = 〈v, f〉 ∀v ∈ H.

Inoltre si ha

‖f‖ = ‖ϕ‖H′ ≡ sup‖v‖=1

|〈ϕ, v〉|.


Proposizione 1.1.6 (Disuguaglianza di Bessel). Sia uαα∈A un insieme or-tonormale in H e cioe 〈uα, uβ〉 = δα,β ove δ·,· e il simbolo di Kronecker. Alloraper x ∈ H: ∑

α∈A|〈x, uα〉|2 ≤ ‖x‖2.

In particolare α|〈x, uα〉 6= 0 e numerabile.

Teorema 1.1.7 (Fourier). Sia uαα∈A un insieme ortonormale numerabile.Allora sono equivalenti:

(i) x =∑α∈A〈x, uα〉uα per ogni x ∈ H.

(ii) (Eguaglianza di Parseval) ‖x‖2 =∑

α∈A |〈x, uα〉|2 per ogni x ∈ H.

(iii) (Completezza) Se x ∈ H e se 〈x, uα〉 = 0 per ogni α, allora x = 0.

Se la successione 〈x, uα〉 e tale che∑α∈A |〈x, uα〉|2 < ∞, allora la serie∑

α∈A〈x, uα〉uα e convergente in H e converge a x. Un insieme ortonormaleche verifichi una delle precedenti condizioni e detto base ortonormale per lospazio di Hilbert H.

Che cosa rappresentano i coefficienti 〈x, uα〉α∈A? In effetti, 〈x, uα〉 calcolala “somiglianza” di x a uα nella distanza definita dalla norma di H. Insomma icoefficienti rappresentano una analisi costitutiva di un elemento x.

Teorema 1.1.8. Ogni spazio di Hilbert ammette una base ortonormale.

Definizione 1.1.9. Uno spazio di Hilbert H e separabile se e soltanto se pos-siede una base ortonormale numerabile. E in tale caso ogni base ortonormale diH e numerabile.

1.2 Serie e trasformate di Fourier

Come si e osservato nel Teorema 1.1.7, l’esistenza di una base ortonormale com-pleta uαα∈A, ossia tale per cui i coefficienti 〈x, uα〉α∈A rappresentano tuttae sola l’informazione relativa all’elemento x ∈ H, garantisce la riproducibilitadell’elemento x come sovrapposizione degli elementi di base uαα∈A. Con illinguaggio adottato nell’introduzione alle note, B = H, Y = ℓ2(A) e si haimmediatamente che la mappa

C : H → ℓ2(A), C(x) = 〈x, uα〉α∈A, (1.2)

e un esempio di digitalizzazione stabile, completa ed invertibile dalla mappa disintesi R = Ruαα∈A

data da

R : ℓ2(A) → H, R(cαα∈A) =∑

α∈Acαuα. (1.3)

1.2. SERIE E TRASFORMATE DI FOURIER 5

Infatti, per 1.1.7 (i) si ha che IH = R C. La stabilita, secondo il principio(I) e l’equivalenza di norme (1) e data dall’eguaglianza di Parseval, Teorema1.1.7 (ii). In tal caso, quindi A = B = 1. In particolare, C ∈ L(H, ℓ2(A))

e R ∈ L(ℓ2(A),H). Proviamo a calcolare ora ‖|C‖|H→ℓ2(A): la norma di→c=

(cα)α∈A ∈ ℓ2(A) e data da ‖ →c ‖ℓ2 =

(∑α∈A |cα|2

)1/2e quindi

‖|C‖|H→ℓ2(A) = supx∈H,x 6=0

‖C(x)‖ℓ2‖x‖ = sup

x∈H,x 6=0

(∑α∈A |〈x, uα〉|2

)1/2

‖x‖ .

Ma per l’eguaglianza di Parseval 1.1.7 (ii) tale quantita e ‖|C‖|H→ℓ22(A) = 1.Per esercizio si provi a calcolare ‖|R‖|ℓ22(A)→H.

-2 -1 1

-2

-1

1

2

Figura 1.1: Orbita di Marte vista dalla Terra secondo il modello tolemaico

Esercizio 1.2.1. Dimostrare che il sistema di digitalizzazione dovuto ad unabase ortonormale completa non e ridondante. Cioe si dimostri che esite unauna unica mappa C tale per cui IH = R C, ove R e fissata come in formula(1.3). Cosa succede ora se alla base ortornormale aggiungo un nuovo ulterioreelemento, chimiamolo u∗ ∈ H? Per la mappa di sintesi R = Ruαα∈A∪u∗esiste ancora una mappa di analisi C tale che IH = R C. E tale mappa unica?E possibile stimare le costanti A,B in questo caso?

Ma allora la ricerca di basi ortonormali permette di costruire dei sistemi perdigitalizzare segnali in maniera stabile e completa, anche se non neccessariamen-te ridondante. Inoltre data una base ortonormale e sempre possibile costruireun sistema che e anche ridondante.Un esempio, classico di base ortonormale e proprio quello suggerito da J. Fouriernel 1815 per la risoluzione del problema del calore.

Diamo ora una semplice interpretazione geometrica di come si costruisce unaserie di Fourier che si ispira al modello tolemaico della rivoluzione dei pianeti. SerP e la distanza media di un pieneta P dal sole in unita astronomiche e τP > 0e il tempo di rivoluzione in anni e se S e un satellite di P che orbita attornoa P con una distanza media rS e con un tempo di rivoluzione τS > 0 allorala traiettoria di S e semplicemente data da zS(t) = rP e

2πit/TQ + rSe2πit/TS .


Analogamente la traiettoria di un secondo pianeta Q visto dal pianeta P saradata da zP→Q(t) = rQe

2πit/TQ − rP e2πit/TP . Se per esempio si considera il

pianeta P come la Terra e il pianeta Q Marte si ha che rP = 1 = τP erQ = 1.55, τQ = 1.88 ≈ 2. La traiettoria di Marte vista dalla Terra si muo-ve nel cielo notturno nella stessa direzione della Luna. Tuttavia ci sono circatre mesi ogni anno in cui il pianeta si muove in direzione opposta (Fig. 1).

Una serie trigonometrica (o una serie di Fourier) in forma complessa diperiodo τ > 0 e allora una serie di funzioni del tipo

T (x) =1√τ

∑

n∈Z

cne2πinx/τ (1.4)

Posto

a0 = 2c0, an = cn + c−n, bn = i(cn − c−n) (1.5)

si ha, per le formule di Eulero (eiw = cos(w) + i sin(w)), che

T (x) =1√τ

(a0

2+

∞∑

k=1

(ak cos(2πkx/τ) + bk sin(2πkx/τ))

)(1.6)

Si osservi ora che 1√τe2πinx/τn∈Z e

√τ2 cos(2πkx/τ)k∈N

⋃√

τ2 sin(2πkx/τ)k∈N,k 6=0

sono due sistemi di funzioni ortonormali nello spazio di Hilbert H = L2(0, τ).Per esempio (vedi 1.1.3)

〈 1√τe2πinx/τ ,

1√τe2πimx/τ 〉 =

1

τ

∫ τ

0

e2πinx/τe2πimx/τdx =1

τ

∫ τ

0

e2πi(n−m)x/τdx =

=

∫ 1

0

e2πi(n−m)xdx = δm,n.

Ci si chiede se costituiscano basi ortonormali complete. Si consideri ora lafunzione data da

f(x) =πx

2− x2

4.

0 1 2 3 4 5 60

0.5

1

1.5

2

2.5

0 1 2 3 4 5 6

0.5

1

1.5

2

2.5

Figura 1.2: f(x) = πx2 − x2

4 = π2

6 −∑∞k=1

cos (kx)k2 .


Si osservi che la funzione e simmetrica e che quindi i coefficienti bk(f) =〈f,√

τ2 sin(2πkx/τ)〉 sono nulli per ogni k. Calcolando (per esercizio! Si integra

due volte per parti) invece i coefficienti ak(f) = 〈f,√

τ2 cos(2πkx/τ)〉 si ottiene

che la serie trigonometrica relativa alla funzione f e data da

Tf(x) =π2

6−

∞∑

k=1

cos (kx)

k2. (1.7)

Come si vede da Fig. 1.2 la serie converge alla funzione f .

Teorema 1.2.2 (Convergenza puntuale). Sia f : R → C periodica di periodoτ > 0 localmente sommabile. Supponiamo che in x0 ∈ R esistano finiti i limitif(x0)

+ e f(x−0 ) e che inoltre i rapporti incrementali destro e sinistro

f(x0 + h) − f(x+0 )

h,

f(x0 − h) + f(x−0 )

h,

siano limitati per h > 0 piccolo. Allora la serie di Fourier Tf di f converge nelpunto x0 alla media dei valori dei limiti destro e sinistro di f in x0:

Tf (x0) =1

τ

∑

n∈Z

〈f, e2πinξ/τ 〉e2πinx/τ =f(x+

0 ) − f(x−0 )

2.

In particolare se f ′(x0) esiste allora Tf (x0) = f(x0).

Dimostrazione. Si veda per esempio [25] VI.6.

Questo fa supporre che 1.1.7 (i) sia verificabile e che quindi il sistema siacompleto. Vediamo come si procede per una traccia della dimostrazione (even-tualmente da percorrere come esercizio).

Un insieme S ⊂ H si dice denso in H se per ogni x ∈ H e per ogni ε > 0esiste p ∈ S tale che ‖x− p‖ ≤ ε. Chiaramente questo equivale a dire che esisteuna successione pnn∈N ⊂ S tale che pn → x per n → ∞. Ma allora se sidimostra che 1.1.7 (ii) e vera per un sottoinsieme denso S allora e vera ancheper H. Infatti, la mappa x→ ‖x‖ e continua in H: se xn → x allora deve essereche ‖xn‖ → ‖x‖. Calcolando 〈xn−x, x+xn〉 ed utilizzando la disuguaglianza diSchwarz si ottiene il limite. Inoltre, per p ∈ H fissato, l’operatore x → 〈x, p〉 elineare (vedi 1.1.1 (i)) e continuo (dalla disuguaglianza di Schwarz e Definizione1.1.4). Ma allora

0 = ‖pn‖2 −∑

α∈A|〈pn, uα〉|2 → ‖x‖2 −

∑

α∈A|〈x, uα〉|2.

Per H = L2(0, τ), lo spazio vettoriale delle funzioni a scalino S e un insiemedenso. E allora sufficiente mostrare 1.1.7 (ii) per funzioni del tipo χ[0,a] per0 < a < τ . Calcolando i coefficienti di Fourier complessi cn(χ[0,a]) si ottiene che

∑

n∈Z

|〈χ[0,a],1√τe2πinx/τ 〉|2 =

a2

τ+

τ

π2

(π2

6−

∞∑

n=1

cos (2πnaτ )

n2

)


Dal Teorema 1.2.2 e da (1.7) sappiamo qual e il valore della serie all’internodelle parentesi, cioe f(2πa/τ). Da cio si conclude facilmente che

a = ‖χ[0,a]‖22 =

∑

n∈Z

|〈χ[0,a],1√τe2πinx/τ 〉|2.

Come volevasi.

Ma allora per f ∈ L2(0, τ) si ha che

f =∑

n∈Z

〈f, 1√τe2πinx/τ 〉 1√

τe2πinx/τ .

Per funzioni f ∈ L2(0, τ) si definisce la trasformata di Fourier F : L2(0, τ) →ℓ2(Z) definita da

Ff(n) = 〈f, 1√τe2πinx/τ 〉 =

1√τ

∫

(0,τ)

f(x)e−2πinx/τdx. (1.8)

0.5 1 1.5 2 2.5 3

0.2

0.4

0.6

0.8

1

Figura 1.3: Serie di Fourier relativa a f(x) = x− [x]−.

Cosa succede se facciamo tendere τ → ∞?

∑

n∈Z


τe2πinx/τ =

1

τ

∑

n∈Z

(∫ −τ/2

τ/2

f(ξ)e−2πinξ/τdξ

)e2πinx/τ .

L’ultima somma e in effetti una somma di Riemann. Questo fa pensare che sef e una funzione sommabile su R allora si possa scrivere qualcosa del tipo

f(x) = limτ→∞

f(x)χ[τ/2,−τ/2] =

= limτ→∞

∑

n∈Z


τe2πinx/τ =

∫

R

(∫

R

f(ξ)e−2πiwξdξ

)e2πiξwdw.

Quando e in che senso tale scrittura sia lecita e pero da provare e da valutare.In verita si ha il seguente


0 50 100 150 200 250-1.5

-1-0.5

00.5

1

0 50 100 150 200 2500

2

4

6

8

Figura 1.4: Segnale sinusoidale affetto da rumore: sono evidenti due componentiin frequenza

Teorema 1.2.3 (Trasformata di Fourier e formula di Plancherel). Si definisceper f ∈ L1(R) la trasformata di Fourier, data da:

Ff(w) =

∫

R

f(x)e−2πiwxdx. (1.9)

Inoltre, se anche Ff ∈ L1(R) allora la trasformata di Fourier e invertibile.

f(x) =

∫

R

Ff(w)e2πiwxdw. (1.10)

Se f ∈ L1⋂L2 allora F : L1

⋂L2 → L2 e un operatore lineare e continuo che

si estende in maniera unica su tutto L2. In particolare, F e un isomorfismoisometrico, ossia vale la formula di Plancherel (eguaglianza di Parseval):

‖f‖2 = ‖Ff‖2. (1.11)

Inoltre, se f, g ∈ L2 allora 〈f, g〉 = 〈Ff,Fg〉.Dimostrazione. Le dimostrazioni si possono trovare in [42] Capitolo 8.

Remark 1.2.4. Questa definizione della trasformata di Fourier e normalmen-te utilizzata nella letteratura relativa all’analisi tempo-frequenza, si veda peresempio [41, 40, 61]. E normalmente utilizzata la definizione

f(ω) =1

2π

∫

R

f(x)e−iωxdx

nella letteratura relativa alle ondine (wavelets), si veda [19, 68]. In queste noteutilizzeremo entrambe per essere consistente con la corrispondente letteratura.In particolare i Capitoli 1, 2 e 3 utilizzeranno la prima, mentre i Capitoli 4 e5 la seconda. In ogni caso, si osservi che Ff(ω) = 2πf(2πω) e che dunque irisultati, salvo opportuna scalatura, possono sempre essere ricondotti all’una oall’altra definizione.

Qual e il significato della trasformata di Fourier? A che cosa serve? Latrasformata di Fourier rappresenta il contenuto in frequenza di una funzio-ne/segnale. Ci dice cioe quali sono i costituenti (importanti) oscillanti di un


segnale e le loro frequenze distintive di oscillazone. La fortuna dell’analisi diFourier risiede essenzialmente nel fatto che e in grado di descrivere uno deicomportamenti piu frequenti in natura: i fenomeni ondulatori, molti dei qualisono regolati da sovrapposizioni di leggi del tipo

yα,s0(t) =

e−αte2πis0t, t > 00, t < 0.

La trasformata di Fourier di yα,s0(t) e detta Lorenziana. Quando, per esempio,delle molecole sono investite da una radiazione elettromagnetica, vengono in-dotte sovrapposizioni di oscillazioni smorzate come descritte dalla legge yα,s0(t).Ciascun costituente della molecola ha proprie oscillazioni distintive ed uniche.Le Lorenziane vengono dette in questo caso lo spettro molecolare. Da qui l’ideache e “costata” il premio Nobel a Richard Earnst per la chimica (1991) perlo sviluppo di uno strumento potente per determinare la struttura di molecoleorganiche complesse. Ma dove si trovano i costituenti? Dove sono localizzatigli atomi che caratterizzano una molecola? La trasformata di Fourier e unostrumento potente ma purtroppo ci dice solamente che alcuni costituenti sonopresenti, ma non come si realizza la loro geometria in una molecola. Sta poial chimico comprendere quale possa essere una forma ammissibile della mole-cola, in base ai possibili legami. Pero, purtroppo non vi sono spesso univocheconfigurazioni. Occorrerebbe quindi uno strumento capace di dare sia l’infor-mazione spaziale sia l’informazione di frequenza. Come vedremo nel seguito unostrumento matematico in grado di svolgere questo servizio e possibile portan-do la nostra capacita di analisi nello spazio delle fasi ( spazio/tempo-frequenza).

Vogliamo ora discutere come sia possibile utilizzare nella pratica gli strumen-ti della analisi di Fourier. Gli esercizi proposti su funzioni semplici, ci mostranoche non e sempre banale poter calcolare espansioni in serie di Fourier. Ne e fa-cile poter calcolare integrali di Fourier, per valutare una trasformata di Fourier.Per funzioni relativamente semplici come f(x) = 1

a2+x2 il calcolo della trasfor-mata di Fourier e tutt’altro che banale (si veda Esercizio 11 Cap. X.7 [25]).In tal caso, si ha Ff(w) = π

a e−|w|a. Una funzione complessa f(z) olomorfa

(ossia derivabile) su un dominio Ω salvo che in un punto a, detto singolarita, siespande in una serie di potenze detta serie di Laurent f(z) =

∑n∈Z

cn(z − a)n.Il coefficiente c−1 e detto residuo e si scrive Res(f, a) = c−1. Si dimostra che sef e olomorfa in un dominio Ω contenente il semipiano Im(z) ≥ 0 eccetto perun numero finito di singolarita non reali e se lim|z|→∞,Im(z)≥0 f(z) = 0 allora siha che per ogni α > 0,

limr→∞

∫ r

−rf(x)eiαxdx = 2πi

∑Res(f(·)eiα·, a), Im(a) > 0.

Analogo risultato per α < 0 e per Im(z) < 0, ma cambiando il segno dell’in-tegrale. Il problema si riduce quindi a verificare tutte le ipotesi e al calcolo deiresidui. Comunque e una procedura piuttosto laboriosa e non ragionevolmenteapplicabile in casi molto complessi. Inoltre, cosa si puo fare con dati provenienti


da sperimentazioni pratiche? Si potrebbero interpolare i dati con funzioni chesono olomorfe nei semipiani superiore ed inferiore, salvo un numero finito disingolarita, calcolare tutti i residui e sperare che la nostra interpolazione siastata sufficientemente buona da poterci assicurare che anche la trasformata diFourier sia ragionevolmente compatibile con il comportamento dei dati. Risul-ta evidente che un simile approccio in casi in cui il numero dei dati sia moltogrande (come in tutti i segnali naturali!) non sia attuabile.

Risulta invece piu ragionevole tentare una approssimazione su un dominiodiscreto e cercare una definizione opportuna di trasformata di Fourier discretache, in qualche senso, possa approssimare bene il comportamento della tra-sformata di Fourier continua. Inoltre, dovremmo assicuraci che il calcolo dellatrasformata di Fourier discreta (DFT) possa essere eseguito in tempi ragionevoli:altrimenti non ne sarebbe valso lo sforzo. Cerchiamo ora di definire quindi unatrasformata di Fourier discreta e rinviamo la discussione sulla approssimazionedella trasformata di Fourier continua alla prossima lezione, ove introdurremo lateoria del campionamento e le opportune stime di errore.

Procederemo, come al solito, con il linguaggio degli spazi di Hilbert. Siconsideri il gruppo ciclico Zn = Z

nZ , n ∈ N\0 delle classi resto modulo n. Inpratica si puo scrivere

Zn = 0, 1, . . . , n− 1,ricordandosi che tali elementi rappresentano i resti dei numeri interi rispetto alladivisione per n. Per esempio, in Z3 l’elemento h = 4 corrisponde all’elemento 1poiche 3 = 3·1+1. Quindi, in particolare, 2+2 = 1 in Z3. Si ha inoltre, ℓ2(Zn) =

(c0, . . . , cn−1)|∑n−1

k=0 |ci|2 <∞ = Cn. Si provi per induzione (esercizio!) che

1 + z + z2 + · · · + zn−1 =

n, z = 1(zn − 1)/(z − 1), altrimenti.

(1.12)

Ma allora non e difficile mostrare che per z = e2πi(k−l)/n si ha

〈 1√n

(e2πikl/n)l∈Zn ,1√n

(e2πilm/n)m∈Zn〉 =

n−1∑

m=0

e2πim(k−l)/n = δk,l.

Quindi, 1√n(e2πikl/n)l∈Znk∈Zn e una base ortonormale completa per lo spazio

di dimensione finita ℓ2(Zn). Ma allora si ha immediatamente che, ogni segnalef di lunghezza n si scrive come:

f =1

n

n−1∑

k=0

〈f , (e2πikl/n)l∈Zn〉(e2πikl/n)l∈Zn . (1.13)

Inoltre anche

Ff(k) =1√n〈f , (e2πikl/n)l∈Zn〉 =

1√n

n−1∑

l=0

f(l)e−2πikl/n, (1.14)


e a sua volta un segnale di lunghezza n. L’operatoreF : ℓ2(Zn) → ℓ2(Zn) e dettotrasformata di Fourier discreta o DFT e realizza tutte le proprieta analoghe aquelle del Teorema 1.2.3. Per esempio, si dimostri per esercizio che

‖f‖ℓ2 = ‖Ff‖ℓ2.

Quanto ci costa calcolare questa trasformata? Assumendo che una operazioneequivalga ad una somma o ad una moltiplicazione allora il costo computazionalee 2n somme e moltiplicazioni per n volte e cioe 2n2. Costando una operazionecomplessa il doppio di una singola reale, il costo totale e C(DFT )(n) = 4n2.Oggi un PC e in grado di eseguire 107 operazioni/secondo e quindi e in gradodi produrre una DFT di un segnale di lunghezza n = 1000 in

4 · 10002 · 1

107= 0.4 sec .

Ma gia per un segnale di lunghezza n = 16384 = 214 il costo e di 107 sec.Un semplice segnale audio di pochi secondi, adeguatamente campionato, puosuperare lunghezza n = 214 senza difficolta. Questo fa comprendere quanto siaimportante poter ridurre ulteriormente questo tempo di calcolo al fine di poterveramente rendere possibile un trattamento di segnali digitali.

1.3 L’algoritmo della Fast Fourier Transform

Allo scopo di comprendere l’algoritmo della Fast Fourier Transform (una manie-ra molto intelligente di svolgere il conto della DFT con un costo computazionalemolto inferiore proposto con grande successo in [14]) occorre approfondire alcu-ne proprieta della DFT. Dato un segnale f di lunghezza n, diremo operatore ditraslazione, operatore definito da

Tmf(k) = f(k −m), m ∈ Zn. (1.15)

e l’operatore di modulazione come

Mmf(k) = e2πimk/nf(k), m ∈ Zn. (1.16)

Inoltre definiremo gli operatori di upsampling e di duplicazione come

U f(h) =

f(h/2), mod(h, 2) = 00, altrimenti,

(1.17)

Df =1

2(f , f), (1.18)

per h ∈ Z2n. L’azione della DTF rispetto a questi operatori e data da

F · Tmf(k) = M−m · Ff(k), F ·Mmf(k) = Tm · Ff(k). (1.19)

F · U f(h) = D · Ff(h). (1.20)

1.3. L’ALGORITMO DELLA FAST FOURIER TRANSFORM 13

Si consideri ora un vettore di lunghezza n = 22 = 4 dato da f = (f0, f1, f2, f3).Vediamo come si possa assemblare f dai singoli fi con le operazioni semplici(1.15) e (1.17)

f0 f2 f1 f3↓ U ↓ U ↓ U ↓ U

(f0, 0) (f2, 0) (f1, 0) (f3, 0)↓ I ↓ T−1 ↓ I ↓ T−1

(f0, 0) (0, f2) (f1, 0) (0, f3)ց ↓ ↓ ւ

(f0, f2) (f1, f3)↓ U ↓ U

(f0, 0, f2, 0) (f1, 0, f3, 0)↓ I ↓ T−1

(f0, 0, f2, 0) (0, f1, 0, f3)ց ւ

f

(1.21)

Osservato che Ffi = fi per ogni i = 0, ...fn−1 si ha che applicando la DFT alprecedente diagramma e sostituendo a U, T−1 rispettivamente D,M1 come datodalle regole di commutazione (1.19) (1.20), si genera un algoritmo ricorsivo dicalcolo della Ff .

f0 f2 f1 f3↓ D ↓ D ↓ D ↓ D(∗, ∗) (∗, ∗) (∗, ∗) (∗, ∗)↓ I ↓M1 ↓ I ↓M1

(∗, ∗) (∗, ∗) (∗, ∗) (∗, ∗)ց ↓ ↓ ւ

(∗, ∗) (∗, ∗)↓ D ↓ D

(∗, ∗, ∗, ∗) (∗, ∗, ∗, ∗)↓ I ↓M1

(∗, ∗, ∗, ∗) (∗, ∗, ∗, ∗)ց ւ

Ff

(1.22)

Abbiamo usato un asterisco per denotare i vettori componenti degli statiintermedi del processo. Cerchiamo di calcolare quanto costa l’algoritmo: innanzitutto assumiamo che il costo effettivo del calcolo di I e di D sia ininfluente.Occorre quindi calcolare solamente il costo di ogni singola M1 che applicata adun vettore di lunghezza l costa infatti l − 1. Si assuma n = 2m. Partendodal basso del diagramma (1.21) viene effettuata una sola M1 su una lunghezzapari al vettore originario e quindi un costo dato da 20( n20 − 1). Tale costova sommato a quello dello stadio superiore, dove occorre eseguire 2(n2 − 1)operazioni corrispondenti a 2 volte M1 calcolata su vettori di meta lunghezza


dell’originale. E cosı via, ottenendo che il costo totale e dato da

C(FFT )(n) =

m−1∑

k=0

2k(n

2k− 1) =

=m−1∑

k=0

(2m − 2k

)= m2m − 2m + 1 = n log2(n) − n+ 1. (1.23)

20 40 60 80 100

0.01

0.02

0.03

0.04

Figura 1.5: Il rapporto C(FFT )(n)C(DFT )(n) per n crescente

Quindi un PC odierno e in grado di produrre una FFT di un segnale dilunghezza n = 16384 in

(21414 − 214 + 1) · 1

107= 0.021 sec,

contro i 107 sec che ci si dovrebbe aspettare dalla DFT!

Lo svolgimento completo dell’esercizio seguente, garantisce la piena com-prensione della lezione e permette di sviluppare una competenza pratica nell’usodegli strumenti teorici sino a qui illustrati.

0 2550751001251500

50

100

150

200

Figura 1.6: Immagine digitale e sua trasformata di Fourier

1.3. L’ALGORITMO DELLA FAST FOURIER TRANSFORM 15

Esercizio 1.3.1. Si implementi in codice C/C++ o Matlab

(DFT) l’algoritmo di calcolo della DFT, eseguita direttamente;

(H-DFT) l’algoritmo di Horner per il calcolo di polinomio; un polinomio p(x) =a0 + a1z + ...+ an−1z

n−1 puo anche essere scritto come

p(x) = a0 + z(...z(an−3 + z(an−2 + zan−1)).

Tale schema e piu vantaggioso del calcolo diretto del polinomio, ma fun-ziona solo per la valutazione puntuale. Si osservi che, posto w = e2πi/n ez = wk si ha che

Ff(k) = f(0) + f(1)z + · · · + f(n− 1)zn−1.

Quindi e possibile applicare lo schema di Horner al calcolo puntuale(!)della DFT,

(FFT) l’algoritmo ricorsivo della FFT, come qui descritto per n = 2m. Imple-mentare inoltre l’antitrasformata di Fourier, cioe quell’operatore F−1 taleper cui

F−1Ff = f .

Si confrontino i tempi di calcolo di queste diverse procedure. Utilizzanto l’al-goritmo di calcolo della FFT unidimensionale, implementare un algoritmo ingrado di calcolare la FFT di un segnale bidimensionale (una immagine).

Si e visto che la trasformata di Fourier continua della funzione f(x) = 1a2+x2

e data da Ff(w) = πa e

−|w|a. Considerare quindi il campionamento di passoτ > 0 di f dato da

f(h) = f(τh), h = −n/2,−(n− 1)/2, . . . , n/2 − 1

Calcolare la FFT di f per diversi valori di τ, n e a. Si confronti il risultato conla trasformata di Fourier continua campionata: si assomigliano? La strategiadiscreta per il calcolo delle trasformate di Fourier e buona?

Abbiamo illustrato l’algoritmo della FFT per n = 2m. Sia w = e2πi/n e sidefinisca la matrice

Fn =

1 1 . . . 11 w . . . wn−1

1 w2 . . . w2·(n−1)

1 w3 . . . w3·(n−1)

· · · . . . . . . . . .1 wn−1 . . . w(n−1)·(n−1)

. (1.24)

La DFT si calcola anche eseguendo

Ff = F ∗ f ,


ove per “∗” si intende il prodotto matrice-vettore. La matrice Fn si fattorizzacome segue:

Fn =

(In/2 Qn/2In/2 Pn/2

)(Fn/2 0

0 Fn/2

)(Hn/2

Ln/2

),

ove In/2 e la matrice identica e

Qn/2 = diag(1, w, w2, . . . , wn−1), Pn/2 = −Qn/2,

ove diag(v) e la matrice diagonale con diagonale principale v. Inoltre

Hn/2 =

1 0 0 0 0 . . . . . . 00 0 1 0 0 . . . . . . 00 0 0 0 1 . . . . . . 00 0 0 0 0 . . . . . . 0· · · . . . . . . . . . . . . . . . . . . . . .0 0 0 0 . . . . . . 1 0

, Ln/2 = rot(Hn/2, 1),

ove rot(M,m) e la matrice generata da M per rotazione destra di m colonne:cioe le ultime m colonne prendono il posto delle prime m colonne che si sonospostate di indice verso destra di m volte. Questo implica che il calcolo della Fnpuo essere risolto ricorsivamente attraverso il calcolo di matrici Fn/2k di ordine

n/2k. Tale interpretazione equivale a (1.22). In effetti tale fattorizzazione dellamatrice (1.24) puo essere generalizzata e dipende dalla fattorizzazione di n intermini di numeri primi n = pm1

1 · pm22 · . . . pml

l . Molte implementazioni dellaFFT realizzano tali fattorizzazioni della matrice (1.23) per primi piccoli, come2, 3, 5 e 7. Uno dei pacchetti piu efficienti e completi per il calcolo della FFTe la Fastest Fourier Transform in the West (FFTW1) sviluppata da M. Frigoe S. G. Johnson [55, 56]. La caratteristica di questa implementazione e datadalla sua adattabilita: infatti la FFTW non usa un algoritmo fissato per ilcalcolo della trasformata, ma adatta l’algoritmo della DFT ai dettagli dellastruttura hardware sottostante allo scopo di raggiungere le migliori prestazioni.Il calcolo ricorsivo cerca di lavorare su piccole strutture di dati che sono in gradodi permanere nella cache, di piu rapido accesso rispetto alla memoria RAM.Sfruttando in modo adattativo le caratteristiche hardware l’algoritmo tenta lamiglior performance. In effetti, il calcolo della trasformata si spezza in due fasi:la FFTW planner costituisce un pre-calcolo in grado di apprendere quale sia ilmodo piu veloce in quel momento per calcolare la trasformata e in seguito il planviene passato alla FFTW executor che, in base alle informazioni del plan calcolala trasformata attuale. Calcolo di FFT mulidimensionali, reali e complessee calcolo parallelo (per MPI) delle stesse sono disponibili nel pacchetto. Permaggiori informazioni per l’uso di tale libreria si veda [56]. La libreria e infattiliberamente scaricabile da internet da http://www.fftw.org/.

1La FFTW e una libreria C/C++, ma puo essere richiamata anche da altri linguaggi diprogrammazione, come il Fortran. Inoltre, la FFT attualmente in uso nelle piu recenti versionidi Matlab, dalla 6.1 in poi, monta la FFTW come libreria di calcolo.

Capitolo 2

Teoria del campionamento

N ella lezione precedente abbiamo discusso la costruzione di sistemi di digi-talizzazione stabili per basi ortonormali. In particolare, si sono definite le basidi Fourier e si e visto come sia possible rappresentare rispetto a tali basi unavarieta di segnali continui e discreti. In particolare, si e mostrato che i segnalidiscreti possono essere trattati in modo estremamente efficiente usando l’algorit-mo della FFT. In questa lezione si vuole discutere la qualita di approssimazioneche la rappresentazione discreta offre rispetto alla rappresentazione continua ele relazioni tra le differenti versioni della trasformata di Fourier (1.8,9,14). In-dichiamo come approfondimenti [64, 70, 45, 38]. Avendo utilizzato il medesimosimbolo F per le differenti definizioni per domini continui limitati, illimitati ediscreti, si invita il lettore a prestare attenzione all’interpretazione del simbolo aseconda della funzione a cui e applicato, se una funzione periodica o a supportolimitato, oppure definita su tutto Rd o discreta.

Nel seguito assumeremo di lavorare in dimensione d ≥ 1. La teoria cheabbiamo sviluppato nella precedente lezioni si estende infatti al caso multidi-mensionale in maniera piuttosto naturale. Per τ = (τ1, ..., τd) ∈ Rd+ si definiscedet(τ) = τ1 · . . . τd e ω = 1

τ . Si ponga T d = (0, 1)d. Si denoti con (v, w) ilprodotto scalare in Rd.

2.1 Relazioni tra le trasformate di Fourier con-

tinua e discreta

Teorema 2.1.1. Se f ∈ L1(Rd) allora la serie∑

k∈Zd f(x − ωk) converge in

L1(ωT d) a una funzione periodica Pf tale che ‖Pf‖1 ≤ ‖f‖1. Inoltre per k ∈ωZd F(Pf)(k) (trasformata di Fourier su ωT d) e uguale a (det(ω))−1/2Ff(k/ω)(trasformata di Fourier su Rd).

Dimostrazione. Sia Q =∏di=1 ωi[− 1

2 ,12 ). Si ha

∫Q

∑k∈Zd |f(x − ωk)|dx =∑

k∈Zd

∫Q+ωk

|f(x)|dx =∫

Rd |f(x)|dx < ∞. Ma allora per la convergenza do-

17

18 CAPITOLO 2. TEORIA DEL CAMPIONAMENTO

minata applicata alle serie (cfr. [42] Folland, pag. 53) si ottiene che la serieconverge in L1 ad una funzione Pf tale che ‖Pf‖1 ≤ ‖f‖1. Per la secondaparte dell’asserto si noti che per (1.8):

F(Pf)(h) =1

(det(ω))1/2

∫

Q

∑

k∈Zd

f(x− ωk)

e−2πi(h,x)/ωdx =

=1

(det(ω))1/2

∑

k∈Zd

∫

Q+ωk

f(x)e−2πi(h,x+ωk)/ωdx =

=1

(det(ω))1/2

∫

Rd

f(x)e−2πi(h,x)/ωdx = (1.9)

= (det(ω))−1/2Ff(h/ω).

Teorema 2.1.2 (Formula della somma di Poisson). Sia f ∈ C(Rd), |f(x)| ≤C(1 + |x|)−d−ε, e |Ff(x)| ≤ C(1 + |x|)−d−ε per un qualche C, ε > 0. Siaω = (ω1, ..., ωd) ∈ (0,∞)d. Allora vale la formula:

(det(ω))1/2∑

k∈Zd

Ff(x− ωk) = (det(1

ω))1/2

∑

n∈Zd

f(n/ω)e−2πi(x,n/ω).

ove entrambe le parti dell’equazione convergono uniformemente in ωT d. Inparticolare per x = 0 si ha

(det(ω))1/2∑

k∈Zd

Ff(ωk) = (det(1

ω))1/2

∑

n∈Zd

f(n/ω).

Dimostrazione. La convergenza assoluta e uniforme delle serie segue dall’ap-plicazione del criterio dell’integrale per la serie

∑k∈Zd(1 + |k|)−n−ε. Allora

Pf(x) =∑

k∈Zd f(x− ωk) e in C(ωT d) e quindi in L2(ωT d). Per il teorema diFourier 1.1.7 e per il precedente teorema 2.1.1 si ha:

Pf(x) =1

(det(ω))1/2

∑

h∈Zd

F(Pf)(h)e2πi(h,x)/ω =

=1

det(ω)

∑

h∈Zd

Ff(h/ω)e2πi(h,x)/ω.

Scambiando ora i ruoli di f e Ff e osservando che FFf = f(−x) si ottienequanto voluto.

Definizione 2.1.3. Sia K ⊂⊂ Rd compatto. Porremo

C(K) = f ∈ C(int(K))|esiste f ∈ C(Rd), f|K = f.

2.1. RELAZIONI TRA LE TRASFORMATE DI FOURIER CONTINUA E DISCRETA19

Si definisca inoltre:

C0(K) = f ∈ C(K)|f|∂K = 0,sottospazio di C(K). C(K) munito della norma ‖ · ‖∞ definisce uno spazio diBanach.

Remark 2.1.4. Si noti che, nelle notazioni della precedente definizione, l’ope-ratore J da C0(K) su Cc(R

n), definito da J(f) = f su K, J(f) = 0 altrove, eun isometria di C0(K) sull’immagine J(C0(K)).

Definizione 2.1.5. Siano f ∈ C(Rn) e τ = (τ1, ..., τd) ∈ (0,∞)d. Diremo cam-pionamento di passo τ l’operatore che a f associa f c ≡ f|τZd. Se f ∈ Cc(R

d),

ovviamente f c risulta nulla fuori da un compatto di Zd. Inoltre, salvo trasla-zioni, si puo pensare f c come definita su un toro discreto d-esimo τZd/τmZd

isomorfo in effetti a Zdm Il campionamento induce quindi un operatore di Cc(Rd)

in ℓ2(Zdm), per un qualche m = (m1, ...,md) ∈ Nd\0.Cerchiamo di studiare alcune proprieta dell’operatore ·c.

Remark 2.1.6. Ci si potrebbe chiedere che relazione intercorra tra Ffc (intesacome DFT) e Ff (intesa come in (1.9)). Una risposta e data dal seguente:

Corollario 2.1.7 (Teorema di Shannon). Sia f ∈ Cc(Rd), |f(x)| ≤ C(1 +

|x|)−d−ε, e |Ff(x)| ≤ C(1+ |x|)−d−ε per un qualche C, ε > 0. Se, per esempio,

supp(f) ⊂∏di=1[0, ai], allora vale

Ffc(l) =(det(m))1/2

(det(a))

∑

k∈Zd

Ff(x− ωk)

c

(l/a), (2.1)

ove il campionamento a sinistra e eseguito con passo τ = 1/ω e quello a destrae eseguito con passo 1/a, in modo tale che ωa = m ∈ Nd.

Dimostrazione. Si assuma supp(f) ⊂∏di=1[0, ai] e si ponga a = (a1, . . . , ad). Si

assuma di aver scelto a > 0 in modo tale che ωa = m ∈ Nd. Ma allora

(det(1

w))1/2

∑

n∈Zd

f(n/ω)e−2πi(x,n/ω) =(det(a))1/2

(det(m))1/2

∑

n∈Zdm

f(n/ω)e−2πi(xa,n)/m.

Se assumiamo che xa = l ∈ Zd allora, per 2.1.2 si ha

(det(m))1/2

(det(a))1/2

∑

k∈Zd

Ff(l/a− ωk) =(det(a))1/2

(det(m))1/2

∑

n∈Zdm

f(n/ω)e−2πi(l,n)/m.

Posto τ = 1/ω, si ha da (1.14) quindi che

Ffc(l) =(det(m))1/2

(det(a))

∑

k∈Zd

Ff(l/a−ωk) =(det(m))1/2

(det(a))

∑

k∈Zd

Ff(x− ωk)

c

(l/a),

ove l’ultimo campionamento e effettuato con passo 1/a.


2.2 Teoria del campionamento

Lemma 2.2.1. Poniamo sinc(x) ≡ sin(πx)πx , x ∈ R e Q =

∏di=1 ωi[−1/2, 1/2).

Allora vale

1

det(ω)1/2

∫

Q

e2πi(k,ξ)/ωe2πi(ξ,x)dx = det(ω)1/2d∏

i=1

sinc(wixi − ki), x ∈ Rd.

Dimostrazione. Usare Fubini-Tonelli.

Teorema 2.2.2 (Del campionamento perturbato in L2). Siano Q =∏di=1 ωi[−1/2, 1/2)

e f ∈ C(Rn)∩L2(Rn) tale che f|τZd ∈ ℓ2. Si scriva f = η+ ǫ, ove Fη = Ff suQ. Allora vale

f(x) =∑

k∈Zn

(f c(τk) − ǫc(τk))

d∏

i=1

sinc(τ−1i xi − ki) + ǫ(x), in L2(Rd). (2.2)

Dimostrazione. Si scriva f = η + ǫ, ove Fη = Ff su Q =∏di=1 ωi[−1/2, 1/2).

Consideriamo P [Fη](ξ) =∑

k∈Zd Fη(ξ − ωk). Sappiamo che P [Fη] ∈ L2(ωT d)e quindi:

P [Fη](ξ) =1

det(ω)1/2

∑

k∈Zd

F(P [Fη])(k)e2πi(k,ξ)/ω ,

con

F(P [Fη])(k) =1

det(ω)1/2

∫

ωTd

P [Fη](x)e−2πi(k,x)/ωdx.

Ovviamente si ha che P [Fη]|Q = Fη. Ma allora si vede subito che F(P [Fη])(k) =1

det(ω)1/2 (f c(−k/ω) − ǫc(−k/ω)), in quanto Fη e non nulla solo su un compatto

e quindi Fη ∈ L1(Rd). Ma allora sostituto P [Fη] a Fη su Q ed estendendo a 0fuori da Q:

f(x) =1

det(ω)1/2F−1

∑

k∈Zd

F(P [Fη])(k)e2πi(k,ξ)/ω

(x) + ǫ(x) =

=∑

k∈Zn

F(P [Fη])(k)F−1

(1

det(ω)1/2e2πi(k,ξ)/ω

)(x) + ǫ(x) =

(2.2.1)

=∑

k∈Zn

(f c(τk) − ǫc(τk))

d∏

i=1

sinc(τ−1i xi − ki) + ǫ(x),

ove tali uguaglianze valgono in L2.

Definizione 2.2.3. Per Q =∏di=1 ωi[−1/2, 1/2), porremo

L2Q(Rd) = f ∈ L2(Rd) : supp(Ff) ⊂ Q.

Diremo inoltre che f ∈ L2Q(Rd) e una funzione di banda ω-limitata.

2.2. TEORIA DEL CAMPIONAMENTO 21

Corollario 2.2.4 (Whittacker-Shannon). Se f ∈ L2(Rd) e una funzione dibanda ω-limitata esiste τ0 > 0 tale che per ogni 0 < τ ≤ τ0

f(x) =∑

k∈Zd

f(τk)

d∏

i=1

sinc(τ−1i xi − ki). (2.3)

Dimostrazione. Si usi il Teorema 2.2.2 osservando che ε = 0 e τ0 = ω−1.

Proposizione 2.2.5. Se f ∈ Cc(Rd) allora la sua trasformata Ff su Rd non

e a supporto compatto.

Remark 2.2.6. Il teorema di Whittacker-Shannon ci garantisce la completariproducibilita di una funzione a banda ω-limitata a partire dai suoi campioni.Tuttavia dalla precedente proposizione si deduce che per le funzioni di Cc(R

d)non c’e speranza di poter applicare il teorema del campionamento e di poter rico-struire la funzione continua direttamente dai suoi campioni! Queste considera-zioni assieme al Corollario 2.1.7 dicono inoltre che l’applicazione dell’operatoredi campionamento seguito dalla applicazione della DFT

·c FCc → ℓ2(Zdm) → ℓ2(Zdm)f → fc → Ffc

equivale a eseguire una trasformata di Fourier continua, una periodicizzazionee un campionamento

F P ·cCc → L2(Rd) → L2(ωT d) → ℓ2(Zd)f → Ff → P [Ff ] → (P [Ff ])c

Inoltre SE Ff e una funzione che decade in modo sufficientemente rapido si hache

Ff(w) ≈ P [Ff ](w), ∀w ∈ Q

e che quindi, per 2.2.4 e pensabile che

Ff(w) ≈∑

l∈Zd,l/a∈Q(P [Ff ])c(l/a)

d∏

i=1

sinc(aiwi − li)

=det(a))

(det(m))1/2

∑

l∈Zd,l/a∈QFfc(l)

d∏

i=1

sinc(aiwi − li).

QUINDI se f e una funzione a supporto compatto la cui trasformata di Fou-rier decade piuttosto rapidamente (e questo e garantito non appena la funzionef e sufficientemente regolare, ossia derivabile a sufficienza!) allora la sua tra-sformata di Fourier puo essere calcolata con buona approssimazione applicandouna DFT ai campioni di f e la formula del Teorema 2.2.4. Quindi per funzionicomplicate, ove sia difficile applicare tecniche di calcolo come quelle dei residui,si puo tentare una approssimazione mediante le trasformate di Fourier discrete.


Assumiamo di seguito che ℓ2(τZd) sia munito del prodotto scalare 〈f, g〉ℓ2 :=det(τ)

∑k∈Zd f(τk)g(τk).

Teorema 2.2.7 (Prodotti scalari). Si assumano f, g ∈ C(Rn)∩L2(Rn) tali chef|τZd, g|τZd ∈ ℓ2. Allora si approssima il prodotto scalare di f e g come :

|〈f, g〉L2 − 〈f|τZd , g|τZd〉ℓ2 | ≤

≤ ‖εf‖ℓ2‖g‖ℓ2 + ‖εg‖ℓ2‖f‖ℓ2 + ‖εf‖ℓ2‖εg‖ℓ2 + ‖εf‖L2‖εf‖L2. (2.4)

Dimostrazione. Per il Teorema 2.2.2 si ha

〈f, g〉L2 =

∫

Rd

f(x)g(x)dx =

∫

Rd

∑

k∈Zd

(f(τk) − εf (τk))

d∏

i=1

sinc(τ−1i xi − ki) + εf (x)

×

×

∑

h∈Zd

(g(τh) − εg(τh))

d∏

i=1

sinc(τ−1i xi − hi) + εg(x)

dx. (2.5)

Dal momento che le funzioni del tipo η∗ appartengono a sottospazi ortogonalidi L2 con rispetto alle funzioni del tipo ε∗, si ha:

〈f, g〉L2 =∑

h,k∈Zd

(f(τk) − εf (τk))(g(τh) − εg(τh))×

×(∫

Rd

d∏

i=1

sinc(τ−1i xi − ki) sinc(τ−1

i xi − hi)dx

)+ 〈εf , εg〉L2 =

dalla ortogonalita dei sinc

= det(τ) ·∑

k∈Zd

(f(τk) − εf (τk))(g(τk) − εg(τk)) + 〈εf , εg〉L2 =

= 〈f|τZd , g|τZd〉ℓ2+〈f|τZd , (εg)|τZd〉ℓ2+〈(εf )|τZd , g|τZd〉ℓ2+〈(εf )|τZd , (εg)|τZd〉ℓ2+〈εf , εg〉L2

Quindi usando la disuguaglianza di Cauchy-Schwarz, si ha la seguente disugua-glianza.

|〈f, g〉L2 − 〈f|τZd , g|τZd〉ℓ2 | ≤

≤ ‖εf‖ℓ2‖g‖ℓ2 + ‖εg‖ℓ2‖f‖ℓ2 + ‖εf‖ℓ2‖εg‖ℓ2 + ‖εf‖L2‖εf‖L2 .

Remark 2.2.8. Si noti che per funzioni f, g a banda ω-limitata si ha che

∫

Rn

fg dx = det(τ)∑

k∈Zn

f c(τk)gc(τk).

2.2. TEORIA DEL CAMPIONAMENTO 23

Remark 2.2.9. Risulta anche evidente ora che se f, g ∈ Cc(Rn) allora in

generale si avra: ∫

Rn

fg dx 6=∑

k∈Zn

f c(k)gc(k).

L’unica speranza di limitare gli errori dovuti al campionamento sui prodotti sca-lari e di scegliere f, g in modo che le code delle trasformate abbiano “energia” (=integrate quadrato) minimo. Per una trattazione piu dettagliata ed approfonditadi veda [?, 70].

In questa sezione si e studiato in particolare il sistema di digitalizzazione percampionamento. Si e in particolare dimostrato che per ogni f ∈ L2

Q la mappa

di analisi C : L2Q → ℓ2 data da C(f) = f(τk)k∈Zd puo essere invertita dalla

mappa di sintesi R : ℓ2 → L2Q definita da R(

→c ) =

∑k∈Zd ck

∏di=1 sinc(τ−1

i xi −ki). Infatti, per 2.2.4 si ha I = R C. Inoltre 2.2.8 assicura che il sistema estabile. Il sistema e ridondante?

Proposizione 2.2.10. L’operatore di traslazione e dato da Txg(t) = g(t − x)per ogni g su Rd. Data una funzione a banda limitata g ∈ L1(Rd) tale cheFg 6= 0 on Ω ⊂ Rd compatto, allora esiste τ > 0 tale che per ogni f ∈ L2

Ω(Rd)

f =∑

k∈Zd

ck(f)Tτkg, (2.6)

per opportuni coefficienti (ck(f))k∈Zd ∈ ℓ2.

Dimostrazione. E noto che τd/2e2πiτkxk∈Zd e ONB per L2([−1/(2τ), 1/(2τ)]d).Per τ > 0 piccolo a sufficienza si ha Ω ⊂ [−1/(2τ), 1/(2τ)]d. Esiste g1 ∈ L1

funzione a banda limitata tale che Fg1 · Fg ≡ 1 su Ω. Quindi per f ∈ L2Ω(Rd)

e per τ > 0 piccolo a sufficienza si ha

Ff(w) = [(FfFg1)Fg](w) =∑

k∈Zd

τd〈FfFg1, e2πiτkx〉e2πiτkwFg(w). (2.7)

Applicando la trasformata di Fourier inversa si ha

f =∑

k∈Zd

τd(f ⋆ g1)(τk)Tτkg, (2.8)

ove ⋆ e il simbolo della convoluzione:

f ⋆ g(x) =

∫

Rd

f(x− y)g(y)dy.

Si puo dimostrare che

∑

k∈Zd

τd|(f ⋆ g1)(τk)|2 = ‖f ⋆ g1‖22 ≤ ‖f‖2

2‖g1‖21. (2.9)


Quindi posto ck(f) = τd(f ⋆ g1)(τk) si ha:

∑

k∈Z

|ck(f)|2 ≤(τd‖g1‖2

1

)‖f‖2

2. (2.10)

In particolare,∑

k∈Z|〈f, Tτkg〉|2 =

∑k∈Z

|〈Ff, e2πiτkwFg〉|2 ≤ τ−d‖f‖22‖Fg‖2

∞.

Il seguente esercizio e difficile: richiede una comprensione approfondita ditutta la lezione e uno sforzo personale per realizzare l’uso concreto degli stru-menti matematici illustrati.

Esercizio 2.2.11. Si implementi in C/C++ o Matlab il sistema di digitalizza-zione descritto dalla precedente proposizione. In particolare si trovi

• una funzione g come richiesta dalle ipotesi;

• una funzione g1 e la si calcoli numericamente;

• un modo per calcolare efficientemente i coefficienti (ck(f))k∈Zd ∈ ℓ2 peruna funzione f a banda ω-limitata;

• fare alcuni esempi di digitalizzazione per funzioni f a banda ω-limitata emostrare la loro sintesi mediante le formule (2.3) e (2.6). Che differenzeci sono tra le due?

2.3 Sistemi ortogonali per L2(Rd): frame di Ga-

bor

Nella precedente lezione abbiamo costruito sistemi e basi ortogonali per funzio-ni definite su insiemi limitati di Rd. Abbiamo visto come tale costruzione ciha permesso di dedurre l’estensione della trasformata di Fourier su tutto Rd,come descritto nel Teorema 1.2.3. Non abbiamo pero mai costruito un sistemaortogonale per tutto L2(Rd). Nel seguito costruiamo quindi sistemi ortogonaliper tutto L2(Rd). Si osservi che anche per funzioni di L2(Rd) valgono le pro-prieta descritte in (1.15-16). Infatti, si definiscono gli operatori di traslazione emodulazione come:

Txf(t) = f(t− x), Mwf(t) = e2πiwtf(t). (2.11)

e vale

FTx = M−xF , FMw = TwF . (2.12)

Inoltre se la convoluzione di f, g ∈ L2(Rd) e la funzione di L2(Rd) definita da

f ⋆ g(x) =

∫

Rd

f(x− y)g(y)dy. (2.13)

2.3. SISTEMI ORTOGONALI PER L2(RD): FRAME DI GABOR 25

si ha cheF(f ⋆ g) = F(f) · F(g). (2.14)

Si dimostri per esercizio che esiste una analoga forma di convoluzione per segnalidiscreti e che vale (2.14) anche in tal caso, utilizzando la DFT.

Si consideri ora Q =∏di=1 ωi[−1/2, 1/2). Certamente ωj+Qj∈Zd definisce

un ricoprimento disgiunto di Rd, e cioe Rd =⋃j∈Zd(ωj + Q). Ma allora ogni

f ∈ L2(Rd) si scrive come

f = F−1

∑

j∈Zd

Ff · χωj+Q

=

∑

j∈Zd

F−1 (Ff · χωj+Q) . (2.15)

Ciascuna funzione F−1 (Ff · χωj+Q) e a banda limitata e supp((Ff · χωj+Q)) ⊂ωj +Q. Ma allora supp(Tωj (Ff · χωj+Q)) ⊂ Q. Da (2.12) si ha quindi che

M−ωjF−1 (Ff · χωj+Q) = M−ωjfj,

e una funzione a banda ω-limitata, ove si e posto fj = F−1 (Ff · χωj+Q). DalTeorema 2.2.4 si ha subito quindi che

M−ωjfj =∑

k∈Zd

fj(τk)d∏

i=1

sinc(τ−1i xi − ki).

Ma allora si ha

fj =∑

k∈Zd

fj(τk)Mωj

d∏

i=1

sinc(τ−1i xi − ki). (2.16)

Combinando (2.15) e (2.16) si ottiene subito che

f =∑

j∈Zd

∑

k∈Zd

fj(τk)Mωj

d∏

i=1

sinc(τ−1i xi − ki). (2.17)

Inoltre, posto g(x) =∏di=1

sin(πxi/τ)πxi/τi

=∏di=1 sinc(xi/τi) la formula (2.17) si

riscrive:f =

∑

j∈Zd

∑

k∈Zd

fj(τk)MωjTτkg. (2.18)

Esercizio 2.3.1. Si dimostri che il sistema 1(det(τ))1/2MωjTτkgj,k∈Zd e una

base ortonormale per L2(Rd).

Si noti che Proposizione 2.2.10 ci dice anche che esistono molte g per cuii termini fj possono essere riprodotti da formule come (2.16). Ma per talig il sistema non e piu in generale ortonormale e, soprattutto, e in generaleridondante. Un sistema ridondante del tipo MωjTτkgj,k∈Zd e detto frame diGabor.


Esercizio 2.3.2. Si ripercorra la costruzione della formula (2.18) per unag come in Proposizione 2.2.10. Si usi quindi l’Esercizio 2.2.11 per calcolarenumericamente i coefficienti fj(τk)j,k.

Capitolo 3

Analisi tempo-frequenza:frame di Gabor

Come si e mostrato nella lezione precedente, e possibile costruire dei sistemiortonormali per L2(Rd) generati per modulazione e traslazione della funzione

g(x) =∏di=1 sinc(xi/τi).

1 2 3 4 5 6

x 104

−1

−0.8

−0.6

−0.4

−0.2

0

0.2

0.4

0.6

0.8

1

Time

Am

plitu

de

Figura 3.1: Un segnale test

I coefficienti di una espansione relativa ad una tale base ortonormale possonoessere visti come il campionamento della seguente trasformata (detta trasfor-

27

28 CAPITOLO 3. ANALISI TEMPO-FREQUENZA: FRAME DI GABOR

mata di Gabor o SFTF - Short Time Fourier Transform) su un reticolo d × ddimensionale:

Vg(f)(b, w) = 〈f,MwTbg〉 =

∫

Rd

f(t)e−2πiw·tg(t− b)dt (3.1)

Infatti, per τ > 0 fissato e ω = 1/τ si ha

f =1

det(τ)

∑

j,k∈Zd

Vg(f)(τk, ωj)MωjTτkg (3.2)

Un segnale f e mappato sul piano tempo-frequenza attraverso tale trasfor-mazione, realizzata valutando il prodotto scalare in L2 di f rispetto ad unafunzione analizzatrice g modificata per applicazione di una traslazione e unamodulazione. La figura 3.1 mostra un segnale test e figura 3.2 mostra la corri-spondente STFT, definita sul piano tempo-frequenza. Il risultante spettrogram-ma assomiglia ad uno spartito musicale, dove ad un dato tempo si puo riferireun insieme di frequenze simultanee (le note se ci si riferisce alla musica). Pertale interpretazione espansioni di segnali per mezzo di sistemi del tipo MωjTτkgcostituiscono uno strumento importante per l’analisi (tempo-frequenza) di suo-ni. Come si e discusso nella sezione 2.3 possono esistere numerose g capaci digenerare espansioni come in (2.18). Normalmente tali espansioni, non avendounici coefficienti (mappe di analisi) sono ridondanti. A questo punto diventa ne-cessaria una formalizzazione matematica per quei sistemi gnn∈N che generanodigitalizzazioni stabili e ridondanti.

Indichiamo come approfondimenti [11, 13, 44, 19, 22, 32, 39, 40, 41, 57, 59,61, 63, 71, 63, 72, 73, 75, 83, 84, 86, 2].

3.1 Frame per spazi di Hilbert

Non e ancora stato coniato un termine italiano per il concetto matematico di fra-me. Sollevandoci dalla responsabilita di essere i primi a coniarlo, ci accontente-remo di riferirci a frame come ad un termine femminile, in spontanea assonanzacon quello di “base” (ortonormale ad esempio).

Le frame sono state introdotte per la prima volta presumibilmente nel 1952da Duffin e Shaeffer in [32] nel loro studio di serie non armoniche di Fourier.Successivamente Daubechies, Grossmann, Meyer, Coifman, Janssen e Morletsono stati gli attori principali dello sviluppo della teoria delle frame.Possiamo pensare ad una frame come ad una generalizzazione di una base peruno spazio. Infatti, dato uno spazio di Banach (per esempio Lp(R)), conoscereuna base per lo spazio, ovvero un fissato insieme di vettori gnn∈Z, ci permettedi scrivere ogni vettore f dello spazio come f =

∑n∈Z

cngn per un’unica sceltapossibili di scalari cn. Per molti spazi che si utilizzano comunemente sappiamoche esistono tali basi, ma avremmo bisogno di poter trovare facilmente gli ele-menti della base, o che questi elementi soddisfino a certe proprieta, e di calcolarefacilmente i cn. Ma tutte queste condizioni sono difficili da ottenere contempo-raneamente.

3.1. FRAME PER SPAZI DI HILBERT 29

Time

Fre

quen

cy

50 100 150 200 250

200

400

600

800

1000

1200

1400

1600

1800

2000

Figura 3.2: Spettrogramma o |STFT |2 del segnale test

In spazi di Hilbert, (per esempio L2(R)), sappiamo che esistono basi ortonor-mali, ma e un problema tuttora studiato quello di cercare basi ortonormali chesoddisfino anche ulteriori proprieta.Nasce cosı il concetto di frame. Infatti una frame permette di scrivere f comef =

∑n∈Z

cngn dove gli scalari cn sono noti. Non richiediamo che i gn sianoortogonali o che i cn siano unici, ma vogliamo poter avere un buon controllodei cn e del comportamento della somma. Quindi richiedere che un insieme siauna frame non e cosı restrittivo come richiedere che sia una base ortonormale,ma ci permette di richiedere che altre condizioni vengano soddisfatte dai gn.Noi tratteremo solamente di frame in spazi di Hilbert, ma bisogna ricordare chequesto concetto e stato esteso da Grochenig [60] ad una larga classe di spazi diBanach e che questa estensione non e banale.Certi tipi di frame dette coerenti (per esempio wavelet e frame di Gabor) so-no, come vedremo, piu facili da costruire rispetto a basi ortogonali, perche sipossono creare partendo da un’unica funzione, detta ondina madre, che vieneopportunamente trasformata mediante operatori unitari, per esempio dilata-zioni, traslazioni e modulazioni. A questo minor sforzo per la definizione deign, si aggiunge il vantaggio che i coefficienti cn che ci permettono di ricostrui-re la f non sono unici. Grazie a cio vengono soddisfatte delle richieste moltoimportanti: la stabilita, la robustezza al rumore e la tolleranza all’errore. Ov-vero si richiede che possibili errori importanti sui cn non modifichino “troppo”la corretta riproducibilita di f . Infatti la trasmissione attraverso un mezzo, osemplici operazioni numeriche come l’arrotondamento, possono produrre errori


importanti sui cn.Notiamo che se i cn sono totalmente indipendenti (come nel caso di una baseortonormale), allora possibili errori sulla f sono individualmente interpretatida ciascun elemento cn di informazione, con possibili errori importanti, globa-li, quando vogliamo ricostruire la f . D’altra parte, possibili correlazioni trai pacchetti di informazione cn possono distribuire gli errori su tutto l’insiemedei dati, magari riducendo cattive distorsioni locali, magari compensando com-pletamente un errore prodotto ad un certo punto, recuperandolo in un puntosuccessivo. Questo si realizza effettivamente usando le frame e puo essere in-terpretato come una intrinseca ridondanza della rappresentazione numerica. Laridondanza di informazione e piuttosto comune in natura ed e fortemente rac-comandabile ogni volta che si possono presentare errori e il grado di ridondanzadovrebbe quindi dipendere dalla frequenza di tali possibili errori. Basti pensa-re ad esempio al DNA umano: e un lungo filamento di informazioni ove solouna minima parte e considerata effettivamente utile. Il resto e protezione ede una parte in qualche senso ridondante. Alcune modificazioni possono esserecompensate, e, magari, i loro effetti negativi essere controllati o ridotti. Questoequivale a dire che possono esistere successioni (cn)n∈Z e(dn)n∈Z distinte taliche f =

∑n∈Z

cnfn =∑

n∈Zdnfn: il processo di analisi ridondante puo offrire

diverse possibili scelte di cn e si puo quindi adattare la scelta di una famigliadi coefficienti rispetto ad un’altra, in dipendenza al particolare uso che se nevuol fare. Tale liberta lascia spazio alla creazione di metodi numerici adat-tativi per la scelta della miglior rappresentazione possibile della soluzione diinteresse, tipicamente quella che richiede minor sforzo computazionale e quellapiu compressa. A prima vista, la ridondanza puo considerarsi in contrasto conil concetto di minima informazione utile. In realta un sistema ridondante dipacchetti di informazione gn puo aumentare la possibilita di identificazione: evero che due differenti gn possono condividere parte del medesimo contenuto diinformazione, ma uno puo corrispondere meglio alla data funzione: come direche “tanto piu grande e il dizionario, tante piu frasi si e in grado di costruire” oanche “tanto piu grande e il dizionario, tanto piu corte sono le frasi che scrivero”.

Data la loro “semplicia d’uso”, flessibilita e robustezza, le frame sono utiliz-zate in diversi campi applicativi, ad esempio, nella simulazione numerica adat-tativa (i.e., soluzione di equazioni integrali e differenzali) [17, 15, 18, 12, 82],nella codifica, compressione e trasmissione di segnali digitali [41, 40, 61, 68] enell’analisi di immagini e suoni [2, 54]. Possiamo raccogliere i metodi associatiall’uso di frame in questi ambiti applicativi sotto il nome di “analisi armonicanumerica”.Da un punto di vista teorico, le frame compaiono in modo naturale come meto-di di discretizzazione di trasformazioni generate da rappresentazioni irriducibili,unitarie e quadrato sommabili di gruppi localmente compatti, come la trasforma-ta di Gabor o la trasformata wavelet (analisi armonica). Lo studio delle funzioniriproducibili da una determinata frame (teoria della approssimazione) conducealla cosiddetta teoria degli spazi coorbita [37, 60, 16] e alla caratterizzazione dispazi funzionali di Banach mediante espansioni frame [48, 46, 35, 49, 52, 16].


Operatori definiti su spazi funzionali vengono anch’essi caratterizzati dalle pro-perieta di struttura e spettrali delle corrispondenti matrici associate a determi-nate frame (analisi funzionale) [61, 69]. Lo studio di frame in spazi di dimensionefinita per il loro uso in teoria dei codici richiede strumenti di geometria di spazidi Banach in dimensione finita. Entrando poi nello specifico delle strutture al-gebriche associate a certe frame, citiamo, infine, di grande interesse teorico, larecente osservazione che la dualita di Wexler-Raz nelle frame di Gabor non ealtro che una espressione della equivalenza (forte) di Morita dovuta a Rieffel [75]della C∗-algebra generata da π(λ) := MwTx (modulazione e traslazione), conλ = (x, ω) su un reticolo Λ del piano tempo-frequenza R × R con la C∗-algebragenerata da π(µ) con µ = (x, ω) sul reticolo Λ, detto reticolo aggiunto [67](geometria non commutativa).

Data la capacita della teoria delle frame di intersecare ambiti molto diversidelle matematica e delle sue applicazioni, e presumibile un loro ulteriore grandesviluppo sia teorico sia applicativo. Basti pensare ai recenti risultati di TerenceTao (medaglia Fields 2006, http://www.math.ucla.edu/∼tao/) [5, 8, 7, 6, 9] inquesto settore.

Sia H uno spazio di Hilbert separabile.

Definizione 3.1.1. Una sequenza gnn∈N ⊂ H e una frame per H se esistonodue costanti A,B > 0 tali che

A · ‖f‖2 ≤∑

n∈N|〈f, gn〉|2 ≤ B · ‖f‖2, ∀f ∈ H. (3.3)

Si osservi che tale condizione (3.3) rappresenta una forma piu debole dellaeguaglianza di Parseval 1.1.7 ii). Inoltre garantisce che il sistema di digitalizza-zione Cg : f → 〈f, gn〉 e stabile. Ci si chiede se tale sistema e fedele, ossia seesiste una mappa di sintesi R che inverta la mappa C.

In effetti questa inversione di attua attraverso un operatore, detto operatoredi frame definito da S : H → H

Sf =∑

n∈N〈f, gn〉gn. (3.4)

In particolare, (3.3) implica che S e positivo, autoaggiunto e invertibile. Cosı siha che

f = SS−1f =∑

n∈N〈f, S−1gn〉gn = S−1Sf =

∑

n∈N〈f, gn〉S−1gn. (3.5)

Il sistema gn = S−1gnn∈N e ancora una frame detta duale canonica dignn∈N con corrispondente operatore di frame S−1. Dal momento che unaframe e tipicamente ridondante, nel senso che non esiste un unica mappa deicoefficienti cnn∈N tale che

f =∑

n∈Ncn(f)gn,


esistono quindi molti possibili duali gnn∈N ⊂ H per cui

f =∑

n∈N〈f, gn〉gn.

In particolare la mappa Rg :→c→ ∑

n∈N cn(f)gn e tale che I = Rg Cg, cioe euna inversa a sinistra della mappa Cg. Si noti inoltre che S = Rg Cg.

Quando gn = S−1gn allora Rg e un inverso sinistro di Cg e vale che

‖|Cg Rg − I‖|ℓ2→ℓ2

e minima. In tal caso Rg e detto l’operatore pseudo-inverso di Cg secondoMoore-Penrose. Quindi per calcolare il duale canonico si puo

• invertire l’operatore di frame S;

• calcolare l’operatore pseudo-inverso di Cg.

Nel caso particolare in cui A = B la frame e detta stretta e si ha che gn =S−1gn = A−1gn. Ovviamente una base ortonormale e una frame stretta conA = B = 1.

Esercizio 3.1.2. Se gnn∈N e una frame stretta con A = B = 1 e se ‖gn‖ = 1per ogni n allora gnn∈N e una base ortonormale.

Si puo dimostrare che un sistema finito Gn = g0, ..., gn−1 e sempre unaframe per H = spanGn = ∑n

i=0 cigi : ci ∈ C.Esempio 3.1.3. Si consideri H = R2, f = (−1, 3) e g0 = (1,−1), g1 = (0, 1),g2 = (1, 1). I coefficienti sono dati da

cnn=0,1,2 = −4, 3, 2

e la frame duale canonica g = (12 ,− 1

3 ), (0, 13 ), (1

2 ,13 ) che da la ricostruzione di

f come:

f =2∑

i=0

cngn = (−2,4

3) + (0, 1) + (1,

2

3) = (−1, 3).

Per esercizio si trovi un altro possibile duale g, ma non canonico.

Come si calcola una frame duale canonica di una frame di dimensione finita?Si consideri una base ortonormale eim−1

i=0 , m ≤ n, di H e si consideri la matriceGm×n = (〈gk, ei〉)i=0,...,m−1,k=0,...n−1. In particolare si ha che per ogni f ∈ H e

considerato→f= (〈f, ei〉)i=0,...,m−1

Cgf = GTm×n∗→f . (3.6)

La procedura cosiddetta SVD (Singular Value Decomposition) applicata aduna matrice M ∈ Cn×m e un elemento importante dell’algebra lineare numerica.


L’idea di base e di poter scrivere la matrice M nella forma M = U∗ ∗D ∗ V oveD e una matrice diagonale e U , V sono matrici le cui righe sono ortonormali eU∗ e la trasposta Hermetiana di U . Gli elementi non nulli di D sono noti comei valori singolari di M . Il rapporto tra il il piu grande valore singolare e il piupiccolo rappresenta il numero di condizionamento di M . Ma allora una inversaa sinistra di M si puo scrivere come

M † = V T ∗D−1 ∗ U.

Infatti M † ∗M = (U∗ ∗D ∗ V ) ∗(V T ∗D−1 ∗ U

)= I. In effetti si dimostra che

la somma delle entrate al quadrato della matrice (M ∗M †−I) e minima e quindisi ha che M † e proprio la matrice dell’operatore pseudo-inverso dell’operatorerappresentato da M . Si noti inoltre che

M † ∗M = MT ∗ (M †)T = MT ∗ (MT )†.

Quindi per una matrice N ∈ Cm×n con n ≥ m si definisce N † = (NT )† e si haN ∗N † = I.

Esercizio 3.1.4. Si provi che le colonne della matrice(G†m×n

)Trappresenta-

no in coefficienti 〈gk, ei〉i=0,...,m−1 del duale canonico di Gn. Se prendo un

insieme di coefficienti→c= (ck)k=0,...,n−1 si ha che

(G†m×n

)T∗ →c restituisce un

vettore di lunghezza m. In particolare si ha che

Rg→c=

((G†m×n

)T∗ →c

)T· eim−1

i=0 . (3.7)

Quindi, combinando (3.6) e (3.7) si ha che

f = Rg Cgf =

((G†m×n

)T∗ GTm×n∗

→f

)T· eim−1

i=0 =

=

(→fT

∗Gm×n ∗ G†m×n

)· eim−1

i=0 (3.8)

Inoltre, si ha che

Rg→c=

(Gm×n∗

→c)T

· eim−1i=0 .

Quindi posto Sg = Gm×n ∗ GTm×n si ha

Sf = Rg Cgf =

(Gm×n ∗ GTm×n∗

→f

)T· eim−1

i=0 =

(Sg∗

→f

)T· eim−1

i=0 .

Ma allora si puo calcolare il duale canonico anche come

gk =(S−1g ∗ →

g k

)T· eim−1

i=0 .


3.2 Frame di Gabor

Per semplicita si assuma d = 1, anche se tutto cio che segue vale in dimensionearbitraria.

3.2.1 Ridondanza della STFT

Con la STFT rappresentiamo il segnale 1-dimensionale f nel piano tempo-frequenza, 2-dimensionale. E quindi presumibile che non sia necessaria conosceretutta la STFT per ricostruire f , ma che sia possibile considerare un campio-namento della STFT su un reticolo nel piano in modo da avere un numerosufficiente di valori della STFT da poter ricostruire la f . Scriviamo la Vg comesegue

Vg(f)(x, ω) := 〈f,MωTxg〉,e invece di considerare tutta la STFT ne consideriamo solo un campionamentosu un reticolo di passi a, b > 0,

(〈f,MbmTang〉)n,m∈Z.

E possibile ricostruire f dai campioni (〈f,MbmTang〉)n,m∈Z?

Il Corollario ?? suggerisce una risposta elegante a questa domanda. Se fosseMbmTangn,m∈Z una frame per L2(R), allora potremmo ricostruire f me-diante una frame duale (canonica). La questione e quindi quando e possibileche MbmTangn,m∈Z sia una frame. In particolare, quanto densamente (cioequanto piccoli devono essere a, b > 0) devo campionare per avere sufficienteinformazione?

3.2.2 Principio di indeterminazione

La trasformata di Gabor distribuisce l’energia della funzione (cioe ‖f‖2) sututto il piano tempo-frequenza R × R. Sebbene Vg non sia suriettiva vogliamoora dimostrare che l’energia viene distribuita in modo uniforme e che non epossibile trovare “concentrazioni” di energia.

Proposizione 3.2.1 (Principio di indeterminazione debole per la STFT). Siano∥∥f∥∥2

2=∥∥g∥∥2

2= 1, U ⊂ R × R e C > 0 costante, tali che

∫ ∫

U

| Vgf(a, b) |2 dadb ≥ C.

Allora |U | ≥ C.

Dimostrazione. . La disuguaglianza di Cauchy-Schwarz implica come in (??)che | Vgf(a, b) |≤ 1. Quindi

C ≤∫ ∫

U

| Vgf(a, b) |2 dadb ≤ ‖Vgf‖2∞|U | ≤ |U |. .

3.2. FRAME DI GABOR 35

In altre parole, se si vuole raccogliere una parte di informazione sulla fdal piano tempo-frequenza, e necessario considerare un sottoinsieme di R × R

di misura non nulla e limitata dal basso, e questa e una naturale alternativainterpretazione del principio di indeterminazione di Heisenberg:

Teorema 3.2.2 (Principio di indeterminazione di Heisenberg). Siano f ∈L2(R) e siano a, b ∈ R due punti arbitrari. Allora

(∫(x− a)2|f(x)|2dx

)1/2(∫(ω − b)2|Ff(ω)|2dω

)1/2

≥ ‖f‖22

4π.

Se f e molto localizzata attorno a a (cioe la gran parte dell’energia e localiz-zata attorno ad a) allora necessariamente Ff deve essere molto poco localizzataattorno a qualsivoglia punto b e viceversa.

A questo punto dovrebbe essere intuitivamente chiaro che le costanti direticolo a, b > 0 di una frame di Gabor MbmTangm,n∈Z devono essere suf-ficientemente piccole per poter campionare sufficiente informazione sul pianotempo-frequenza.

3.2.3 Frame di Gabor e il gruppo di Weyl-Heisenberg

Feichtinger e Grochenig hanno dimostrato in [36, 37, 60] che la trasformatadi Gabor e profondamente collegata alla teoria di rappresentazione dei gruppi.In particolare, la teoria sopra esposta sulla trasformata di Gabor continua puoessere generalizzata a funzioni di L2(G) oveG e un gruppo localmente compatto,i.e., uno spazio topologico localmente compatto con una operazione del gruppo· : G × G → G che mappa (x, y) 7→ x · y e tale che l’operazione da G in Gx 7→ x−1 siano entrambe continue.Se intendiamo con H = T × R × R il gruppo di Weyl-Heisenberg, possiamodefinire una rappresentazione W di H su L2(R) con

W (t, a, b)f(x) = te2πibxf(x− a) = tMbTaf(x)

ove (t, a, b) ∈ H e f ∈ L2(R).Allora si ha

W (t1, a1, b1)W (t2, a2, b2)f(x)

= W (t1, a1, b1)(t2e2πib2xf(x− a2))

= t1e2πib1xt2e

2πib2(x−a1)f(x− a1 − a2)

= t1t2e−2πib2a1e2πi(b1+b2)xf(x− a1 − a2)

= W (t1t2e−2πib2a1 , a1 + a2, b1 + b2)f(x)

E quindi, dato che W e una rappresentazione delle operazioni del gruppo su H,deve essere:

(t1, a1, b1) · (t2, a2, b2) = (t1t2e−2πib2a1 , a1 + a2, b1 + b2).


E facile verificare che questa e infatti una operazione del gruppo ovvero e as-sociativa, con identita (l’elemento (1,0,0)) e inverso (l’elemento (t, a, b)−1 =(t−1e2πiab,−a,−b)). Poiche ogni W (t, a, b) e un operatore unitario su L2(R), siha che W e una rappresentazione unitaria di H su L2(R).Quindi possiamo generalizzare la trasformata di Gabor ponendoVgf(a, b) := 〈f,W (1, a, b)g〉 e osservando che il Teorema ?? vale anche nellaforma piu generale:

Teorema 3.2.3. Se f, g ∈ L2(R) allora

∫

R

∫

R

∫

T

| 〈f,W (t, a, b)g〉 |2 dtdadb =∥∥f∥∥2

2

∥∥g∥∥2

2

E dal teorema precedente deduciamo che∫

R

∫

R

∫

T

| 〈f,W (t, a, b)g〉 |2 dtdadb =

∫

R

∫

R

| 〈f,W (1, a, b)g〉 |2 dadb

e si possono concludere le stesse affermazioni della Sezione ??. Infatti nonsi perde nulla se si ignora la componente torica della rappresentazione W delgruppo, e quindi ci si puo ricondurre alla definizione sopra data di Vg. La teoriadi Feichtinger e Grochenig offre un contesto molto generale che garantisce cheun opportuno campionamento MbmTangm,n∈Z generi una frame per L2(R).Non entriamo nei dettagli di questa teoria e studiamo invece delle condizioninecessarie e sufficienti e direttamente verificabili perche MbmTangm,n∈Z generiuna frame.

3.2.4 Condizioni necessarie per avere una frame di Gabor

I seguenti risultati sono dovuti essenzialmente al lavoro di Daubechies, Gross-man [19] e Meyer [69].

Definizione 3.2.4. Data g ∈ L2(R) e siano a, b > 0, diciamo che (g,a,b) generauna frame di Gabor per L2(R) se MmbTangm,n∈Z e una frame per L2(R).La funzione g viene chiamata ondina madre o wavelet di analisi.I numeri a, b si dicono parametri di frame ed a e il parametro di traslazione,mentre b e il parametro di modulazione.

E chiaro che MmbTangm,n∈Z e una frame per L2(R) sse TanMbngm,n∈Z

e una frame per L2(R). E successivamente useremo una o l’altra definizione aseconda della convenienza.

Proposizione 3.2.5. Sia g ∈ L2(R) e a, b > 0 tali che:

1. esistono costanti A,B t.c. 0 < A ≤∑n∈Z| g(x− na) |2≤ B <∞ q.o.

2. g ha supporto compatto, con supp(g) ⊂ I ⊂ R, con I intervallo di lun-ghezza 1/b

Allora (g, a, b) genera una frame di Gabor per L2(R) con frame bound b−1A, b−1B


Dimostrazione. . Fissato n, osserviamo che la funzione fTnag ha supporto inIn = I − na = x − na : x ∈ I, un intervallo di lunghezza 1/b. Segue dallacondizione 1. che g e limitata, quindi fTnag ∈ L2(In). Ma l’insieme di funzionib1/2e2πimbχIn m∈Z = b1/2MmbχInm∈Z e una base ortonormale per L2(In),quindi∑m∈Z

| 〈fTnag,MmbχIn〉 |2= b−1∫| f(x) |2| g(x− na) |2 dx. Quindi :

∑

m,n∈Z

| 〈f,MmbT−nag〉 |2=∑

m,n∈Z

| 〈fTnag,MmbχIn〉 |2

= b−1∑

n∈Z

∫

R

| f(x) |2| g(x− na) |2 dx

= b−1

∫

R

| f(x) |2∑

n∈Z

| g(x− na) |2 dx

E il risultato ora segue usando la condizione 1. e notando che si somma suZ.

Proposizione 3.2.6. Sia g una funzione continua con supporto contenuto inun intervallo I di lunghezza L > 0 che non assume il valore 0 nei punti interni diI. Allora (g, a, b) genera una frame per L2(R) per ogni 0 < a < L e 0 < b ≤ 1/L.

Dimostrazione. . Ovviamente g ∈ L2(R). Poiche 1/b ≥ L vediamo che ilsupporto di g e contenuto in un intervallo di lunghezza 1/b. Definiamo G(x) =∑n∈Z

| g(x − na) |2. Dal teorema precedente segue che (g, a, b) genera unaframe per L2(R), se dimostriamo che G e limitata dal basso e dall’alto. Poicheg ha supporto compatto, la somma che definisce G e una somma finita, che haal piu 1/ab termini perch, dato che l’intervallo e contenuto in uno di lunghezza1/b gli addendi diversi da zero sono quelli t.c. R < x− na < R + 1/b dove conR indichiamo l’estremo inferiore di I. Quindi G e limitata dall’alto dato che glo e, essendo continua con supporto in I. Sia J un sottointervallo di I con lostesso centro, ma di lunghezza a. Dato x ∈ R, esiste un n ∈ Z t.c. x− na ∈ J ,cosı infx∈R G(x) ≥ infx∈J | g(x) |2> 0.

Proposizione 3.2.7. Sia che g abbia supporto compatto o meno, e necessaria lacondizione 1. del Teorema 3.2.5 perche (g, a, b) generi una frame e in particolareg deve essere limitata.

Dimostrazione. . Sia G(x) =∑n∈Z

| g(x − na) |2 . Dobbiamo dimostrare chenon e limitato dall’alto e dal basso. Assumiamo che ess infx∈RG(x) = 0. Datoδ > 0 possiamo trovare un insieme ∆ ⊂ I ⊂ R, dove I e in un intervallo dilunghezza 1/b t.c. | ∆ |> 0 e G(x) ≤ bδ su ∆. Se poniamo f = χ∆, allora∥∥∥f∥∥∥

2

2=| ∆ | e, come nella dimostrazione del Teorema 3.2.5,

∑

n,m∈Z

| 〈f,MmbTnag〉 |2= b−1∑

n∈Z

∫

R

| f(x) |2| g(x− na) |2 dx


= b−1

∫

∆

∑

n∈Z

| g(x− na) |2 dx ≤ b−1 | ∆ | bδ = δ∥∥f∥∥2

2

Poiche δ e arbitrario, (g, a, b) non genera una frame. Una dimostrazione similemostra che G deve essere limitato dall’alto.

Illustriamo di seguito un risultato generale e non semplice da dimostrare [61,Capitolo 7], legato a risultati profondi di Rieffel [75].

Teorema 3.2.8. Per g ∈ L2(R), a, b > 0, se (g, a, b) genera una frame di Gaborper L2(R), allora ab ≤ 1.

3.2.5 Condizioni sufficienti per avere una frame di Gabor

Poiche F(MmbTnag) = TmbM−naFg e la trasformata di Fourier e una mappaunitaria di L2(R) in L2(R), vediamo che (g, a, b) genera una frame di Gaborper L2(R) se e solo se (Fg, a, b) genera una frame di Gabor per L2(R). E facilevedere che se g soddisfa la condizione 2. del Teorema 3.2.5 e se ab > 1, alloraess infx∈R

∑n∈Z

| g(x− na) |2= 0, e cosı g non puo generare una frame. Infattil’insieme MmbTnagn,m∈Z non e nemmeno completo, cioe la f ∈ L2(R) ortogo-nale a tutte le (g,a,b) non e la funzione nulla. Infatti

⋃n∈Z

supp(Tnag) non copreR e quindi ogni funzione con supporto nel complementare di

⋃n∈Z

supp(Tnag)risulta essere ortogonale a ogniMmbTnag. Si puo dimostrare che se ab > 1 alloraMmbTnag non puo mai essere completa in L2(R) per ogni g ∈ L2(R). Accan-to al principio di indeterminazione sopra dimostrato, quest’ultima osservazionespiega in modo piu rigoroso la necessita di campionare in modo sufficientementedenso la STFT per poter avere una completa ricostruzione di una funzione f .Il seguente teorema sull’esistenza delle frame di Gabor in L2(R) con g a supportonon compatto e stato formulato da Daubechies [19], usando nella dimostrazionela formula della somma di Poisson. La dimostrazione riportata, invece non lautilizza e generalizza le condizioni del teorema della Daubechies.

Teorema 3.2.9. Sia g ∈ L2(R) e sia a > 0 t.c.

1. esistono costanti A,B t.c. 0 < A ≤∑n∈Z| g(x− na) |2≤ B <∞ q.o.

2. limb→0

∑k 6=0,k∈Z

β(k/b) = 0 ove

β(s) = ess supx∈R

|∑

n∈Z

g(x− na)g(x− s− na) |=∥∥∑

n∈Z

Tnag · Tna+sg∥∥∞

allora esistono b0 > 0 t.c. (g, a, b) genera una frame di Gabor per L2(R) perogni 0 < b < b0.

Dimostrazione. . Prima assumiamo che f sia continua e con supporto compatto.Questo garantisce che possiamo scambiare le seguenti sommatorie e integrali.Per un fissato n consideriamo un funzione 1/b-periodica data da

Fn(t) =∑

k∈Z

f(t− k/b)g(t− na− k/b)


Ora, Fn ∈ L1[0, 1/b] dato che f e g sono limitate, vale

∫

R

f(t)g(t− na)e−2πimbtdt =

∫ 1/b

0

Fn(t)e−2πimbtdt.

Dato che b1/2Mmbχ[0,1/b]dm∈Z e una base ortonormale per L1[0, 1/b], dallaformula di Plancherel si ha

∑

m∈Z

|∫ 1/b

0

Fn(t)e−2πimbtdt |2= b−1

∫ 1/b

0

| Fn(t) |2 dt

Quindi:

∑

n∈Z

∑

m∈Z

| 〈f,MmbTnag〉 |2=∑

n∈Z

∑

m∈Z

|∫

R

f(t)g(t− na)e−2πimbtdt |2

= b−1∑

n∈Z

∫ 1/b

0

|∑

k∈Z

f(t− k/b)g(t− na− k/b) |2 dt

= b−1∑

n∈Z

∫ 1/b

0

∑

l∈Z

f(t− l/b)g(t− na− l/b) ·∑

k∈Z

f(t− k/b)g(t− na− k/b)dt

= b−1∑

n∈Z

∑

l∈Z

∫ 1/b

0

f(t− l/b)g(t− na− l/b) ·∑

k∈Z


= b−1∑

n∈Z

∫

R

f(t)g(t− na) ·∑

k∈Z


= b−1∑

k∈Z

∫

R

f(t)f(t− k/b) ·∑

n∈Z

g(t− na)g(t− na− k/b)dt

= b−1

∫

R

| f(t) |2 ·∑

n∈R

| g(t− na) |2 dt+

b−1∑

k∈Z;k 6=0

∫

R

f(t)f(t− k/b) ·∑

n∈Z

g(t− na)g(t− na− k/b)dt = (∗)

Usando la disuguaglianza di Cauchy-Schwarz otteniamo

(∗) ≤ b−1B∥∥f∥∥2

2+ b−1

∑

k∈Z;k 6=0

β(k/b)

∫

R

f(t)f(t− k/b)dt ≤ B0(b)∥∥f∥∥2

2

e

(∗) ≥ b−1A∥∥f∥∥2

2− b−1

∑

k∈Z;k 6=0

β(k/b)

∫

R

f(t)f(t− k/b)dt ≥ A0(b)∥∥f∥∥2

2

dove

A0(b) = b−1A− b−1∑

k∈Z;k 6=0

β(k/b) e B0(b) = b−1B + b−1∑

k∈Z;k 6=0

β(k/b)


Dalla condizione 2. esiste, per definizione di limite e dato che A > 0 e B < ∞,un b0 > 0 t.c. A0(b) > 0 e B0(b) <∞ per ogni 0 < b < b0.Supponiamo ora f ∈ L2(R) arbitraria. Allora possimo trovare una successionedi funzioni continue a supporto compatto fj, t.c. fj → f in L2(R) per j → ∞.

Dal risultato appena ottenuto abbiamo che vale Ao(b)∥∥fj∥∥2

2≤ ∑

n∈Z

∑m∈Z

|〈fj ,MmbTnag〉 |2≤ Bo(b)

∥∥fj∥∥2

2. Non e difficile mostrare che queste disugua-

glianze valgono anche nel limite per j → ∞, da cui (g, a, b) genera una frame diGabor con frame bound A0(b), B0(b) per ogni 0 < b < b0.

Se a, b > 0 e se un sistema MajTbkgj,k∈Zd e frame di Gabor per tuttoL2(Rd) allora deve essere

A‖f‖22 ≤

∑

j,k

|〈f,MajTbkg〉|2 ≤ B‖f‖22, (3.9)

e si puo dimostrare che necessariamente a · b ≤ 1. In particolare si ha chel’operatore di frame

Sgf =∑

j,k

〈f,MajTbkg〉MajTbkg,

commuta con gli operatori di modulazione e traslazione nel senso che

SgMajTbk = MajTbkSg.

Questo implica che il duale canonico S−1g MajTbkg = MajTbkS

−1g g. Si ponga

g = S−1g g e si ha

f =∑

j,k

〈f,MajTbkg〉MajTbkg. (3.10)

Si osservi in particolare che 〈f,MajTbkg〉 = Vg(f)(bk, aj).

Esercizio 3.2.10. Si provi che se G(g, a, b) = MajTbkgj,k e una frame diGabor allora G(Fg, b, a) = MbjTakFgj,k e frame di Gabor con duale canonicogenerato da F g.

Teorema 3.2.11 (Esistenza delle frame di Gabor). Se g ∈ L2(Rd) e una fun-zione “ben localizzata” nello spazio tempo-frequenza (per esempio una funzio-ne molto regolare con decadimento rapido all’infinito), allora esistono a, b >0 sufficientemente piccoli tali che G(g, a, b) = MajTbkgj,k e G(Fg, b, a) =

MbjTakFgj,k sono frame di Gabor. In particolare se g(x) = e−π|x|2

e laGaussiana si ha che G(g, a, b) e frame di Gabor per ogni a, b > 0 tali che a·b < 1.

Come al solito il calcolo dei coefficienti di una espansione (infinita) come in(3.10) costituisce un problema: soprattutto non e chiaro come si possa calcolareil duale, ossia l’inverso di Sg. Ancora una volta conviene spostarsi su un dominiodiscreto e in dimensione finita ed applicare le considerazioni relative alle framein dimensione finita della precedente sezione.


-10 -5 5 10

-1

-0.5

0.5

1

0 50 100 150 200020406080100

Figura 3.3: Spettrogramma o |STFT |2 del segnale test a sinistra

3.2.6 Frame di Gabor discrete

Sia d = 1 e si consideri a, b, L ∈ N tali che a|L e b|L e a · b ≤ L. Si ponganoN = L/a e M = L/b. Allora si definisce un sistema di Gabor discreto come

gm,n = MmbLTang, m = 0, ...,M − 1, n = 0, ..., N − 1, (3.11)

ove g ∈ ZL. Si noti che N ·M ≥ L. Ma allora si cercano frame per ℓ2(ZL)del tipo G(g, a, b) = gm,nm=0,...,M−1,n=0,...,N−1. Considerata eiL−1

i=0 la basecanonica di H = ℓ2(ZL) si puo definire la matrice

GL×N ·M =(gmod(k,M), k−mod(k,M)

M

)

k=0,...,N ·M−1, (3.12)

le cui colonne sono rappresentate dai valori di G(g, a, b) = gm,nm=0,...,M−1,n=0,...,N−1.Come si e mostrato nella precedente sezione la matrice relativa all’operatore diframe e data da

Sg = GL×N ·M ∗ GTL×N ·M . (3.13)

Infatti si ha che Sg : ℓ2(ZL) → ℓ2(ZL). Allora un sistema come in (3.11) e unaframe di Gabor discreta se la matrice Sg e invertibile, ossia quando det(Sg) 6= 0.

Esercizio 3.2.12. Posto che il prodotto riga-colonna matriciale equivale al cal-colo di un prodotto scalare e posto che il prodotto scalare di una funzione perun’altra traslata equivale ad una convoluzione:

〈f , Takg〉 = f ⋆ g(ak),

si osservi che, per (2.13-14) si ha

〈f , Takg〉 = F−1(Ff · Fg)(ak).

Si usi la FFT per calcolare rapidamente l’operatore di frame (3.13) e si calcolila complessita dell’algoritmo in dipendenza di L, a, b.


0 50 100 150 200-0.2

00.20.40.60.8

1

0 50 100 150 200-0.04-0.02

00.020.040.060.08

0 20 40 60 80 1001200

0.2

0.4

0.6

0.8

1

0 20 40 60 80 100120-0.04-0.02

00.020.040.060.080.1

Figura 3.4: Atomi di Gabor e corrispondenti duali. In alto e rappresentata asinistra la Gaussiana e a destra il corrispondente duale per L = 132, a = b = 11,con ridondanza L/(ab) = 1.09. In basso e rappresentato il seno cardinale e ilrispettivo duale per L = 240, a = b = 15, con ridondanza L/(ab) = 1.06.

Una volta calcolato l’operatore di frame, l’applicazione di una inversione nu-merica della matrice Sg permette il calcolo del duale. Come l’operatore anchela sua matrice Sg e simmetrica, positiva ed invertibile. Inoltre la matrice e spar-sa, nulla fuori da alcune specifiche diagonali ed e circolante. La struttura dellamatrice permette un calcolo estremamente rapido dei prodotti matrice-vettore,ma si rimanda a [?, ?] per maggiori dettagli. In casi reali, come nel trattamentodi un segnale audio di alcuni secondi, sicuramente le dimensioni di tali matricipuo essere anche del tipo 10000× 10000.

Risulta quindi necessario l’uso di metodi iterativi per l’inversione della ma-trice come, per esempio quello del gradiente coniugato (CG).

0 50 1001502000

0.00050.001

0.00150.002

0 50 100 150 200-0.02

00.020.040.06

0 50 100 150 200-0.04-0.02

00.020.040.060.08

0 50 100 150 200-0.04-0.02

00.020.040.060.08

Figura 3.5: Calcolo iterativo di un duale mediante CG

Si noti che tanto piu ci si avvicina alla soglia critica a · b ≈ L tanto piu peg-giora il numero di condizionamento della matrice Sg. In molti casi l’applicazionediretta del semplice CG fallisce non riuscendo a risolvere il sistema. L’uso diopportuni precondizionatori, e quindi del PCG, diventa quindi necessario. Adesempio il calcolo del duale in figura 3.5 viene effettuato con 139 iterazioni edun errore residuo inferiore a 10−7, essendo il numero di condizionamento della

3.3. COMPRESSIONE DI UN SEGNALE 43

matrice circa 9. In seguito al precondizionamento il numero di condizionamentodella matrice scende a 2 comportando che il medesimo calcolo viene eseguito in75 iterazioni.

L’analisi dei segnali mediante frame di Gabor e quindi un esempio importantedi applicazione dove strumenti tra i piu importanti dell’analisi numerica vengonoimpiegati.

3.3 Compressione di un segnale

Una volta costruita una frame di Gabor G(g, a, b) ed un suo corrispondenteduale G(g, a, b) ogni segnale f di lunghezza L puo essere (de)composto comesegue:

f =∑

m,n

〈f , gm,n〉gm,n. (3.14)

0 200 400 600 8000

5000

10000

15000

20000

25000

30000

35000

Figura 3.6: Coefficienti |〈f , gm,n〉|2.

Il numero di coefficienti e pari aN ·M ≥ L e quindi rappresenta una digitaliz-zazione ridondante. In altre parole il numero di elementi di informazione che sistanno utilizzando e piu grande rispetto al numero di campioni che costituisco-no il segnale stesso. Questo a priori sembrerebbe aver complicato solamente larappresentazione persino peggiorando la quantita di dati da trattare. In realta,sebbene due elementi della rappresentazione condividano una certa parte di in-formazione (da cui la ridondanza), puo succedere che uno dei due stia in effettirappresentando meglio (rispetto all’altro) il segnale. Per cui la ridondanza au-menta la possibilita di identificazione: piu il dizionario e grande piu corte sonole frasi che si possono comporre.

Inoltre come si puo vedere da figura 3.6 i coefficienti della espansione (3.14)non sono tutti importanti allo stesso modo. Ve ne sono di molto rilevanti eve ne sono di molto poco rilevanti. Ci si puo quindi aspettare che eliminandoquegli elementi della espansione i cui coefficienti sono al di sotto di una soglia diaccettabilita si possa ridurre di molto l’informazione necessaria a rappresentarecomunque bene il segnale. Figura 3.7 mostra infatti solo i coefficienti che sonosuperiori ad una data soglia.


0 200 400 600 8000

5000

10000

15000

20000

25000

30000

35000

Figura 3.7: Coefficienti |〈f , gm,n〉|2 superiori ad una soglia fissata.

Figura 3.8 descrive un segnale audio compresso al 58% e i rispettivi spet-trogrammi. Come si puo osservare dallo spettrogramma, la riduzione, sebbenesostanziale a livello di informazione utilizzata, non cambia molto la effettivarappresentazione del segnale stesso.

All’indirizzo http://www.math.unipd.it/∼mfornasi/SignalProcessing.zip e pos-sibile scaricare un pacchetto contenente un tutorial come notebook di Mathe-matica (Wolfram Research) che implementa quanto descritto in questa lezione.Si suggerisce di utilizzarlo e in seguito di realizzare il seguente esercizio.

Esercizio 3.3.1. Si mostri che se g1nn∈N e g2

nn∈N sono due frame per L2(R)(o per ℓ2(ZN ) risp.) con rispettivi duali g1

nn∈N e g2nn∈N allora gn(x1, x2) =

g1n(x1) · g2

n(x2)n∈N e una frame (prodotto tensore) per L2(R2) (o per ℓ2(Z2N )

risp.) con duale gn(x1, x2) = g1n(x1) · g2

n(x2)n∈N.

Si usi questo risultato per costruire frame di Gabor per segnali bidimensio-nali in ℓ2(Z2

N ) ed implementare (in C/C++ oppure Matlab) un algoritmo dicompressione di immagini digitali.

3.3. COMPRESSIONE DI UN SEGNALE 45

Conpression rate: 58.3333%

0 200 400 600 8000

100

200

300

400

8218.11 Second, Null<

0 200 400 600 8000

100

200

300

400

837.84 Second, Null<

Figura 3.8: Coefficienti |〈f , gm,n〉|2 superiori ad una soglia fissata.


Figura 3.9: Confronto di un algoritmo di compressione mediante frame di Gabore JPEG

Capitolo 4

2D pattern matching:armoniche circolari

La presente lezione si propone di illustrare un metodo per rappresentare econfrontare immagini digitali identificate matematicamente con funzioni dellospazio ℓ2(Ω), Ω cerchio di Z2, a meno di rotazioni.I coefficienti di un particolare sviluppo di Fourier delle immagini, inizialmenteconsiderate su dominio continuo, secondo una base costituita da autofunzio-ni di ogni operatore di rotazione, definiscono, infatti, un insieme di “pseudo-invarianti” che caratterizzano l’immagine indipendentemente dal suo orienta-mento. Tale sviluppo, identificato da un’ equazione che coinvolge i cosiddetticoefficienti di Bessel di prima specie Jm, si mantiene, per una sua buona parte,rappresentazione efficace anche restringendo il dominio al reticolo Z2, determi-nando un carattere di invarianza per rotazioni anche di un’immagine campio-nata. Mediante l’opportuna definizione di un coefficiente di matching, vienesuggerita una particolare applicazione di tale metodo al problema della rico-struzione virtuale di affreschi minutamente frammentati, dei quali si dispongadi un congruo numero di parti e di un’immagine antecedente alla distruzione.

Indichiamo come approfondimenti [10, 54].

4.1 Descrizione del problema

Lo sviluppo tecnologico, caratterizzato negli ultimi anni in modo particolare dalcrescente utilizzo dei calcolatori elettronici e dal potenziamento delle loro possi-bilita, ha fornito nuovi strumenti per la risoluzione di questioni che, per la lorocomplessita, erano considerate insormontabili, aprendo, cosı, nuove frontiere eproblematiche. Uno dei campi di applicazione e stato quello connesso alla rap-presentazione, al trattamento e alla trasmissione delle immagini. E noto, infatti,il notevole impatto che l’informazione visiva riesce ad ottenere, costituendo mo-tivo di forte interesse. Nuovo impulso e stato dato anche alla matematica al fine

47

48 CAPITOLO 4. 2D PATTERN MATCHING: ARMONICHE CIRCOLARI

di occuparsi di problemi concretamente legati alle immagini con le piu moderneteorie dei segnali [?, ?].

Uno degli strumenti che, sempre di piu, e affiancato al calcolatore e lo scan-ner, dispositivo usato al fine di digitalizzare immagini, per poterle rappresentaree trattare con facilita con i moderni programmi di grafica. Il campionanamentoeffettuato dagli scanner e, generalmente, su reticoli a geometria quadrata: in ta-le contesto, un’immagine campionata e ruotata non puo, in generale, coincidere,salvo interpolazioni e quindi errori, all’immagine stessa ruotata e campionata,tranne per le rotazioni di angoli multipli di π

2 , che rispettano la geometria delreticolo di campionamento. Di conseguenza, si evidenzia il problema di un con-fronto di immagini digitali a meno di rotazioni.

Richiedere un tale confronto presuppone una ricerca di quei contenuti diinformazione che siano invarianti rispetto all’orientamento. In particolare, siesige che solo una parte dell’informazione possa venire utilizzata allo scopo delconfronto. Infatti, il fatto di voler lavorare con una geometria circolare, diffe-rente da quella intrinseca del supporto, imporra necessariamente di operare conapprossimazioni e con uno spazio di informazione ridotto e, soprattutto, sullarestrizione ad un dominio isotropo. Nel caso di figure geometriche piane sonostati studiati innumerevoli invarianti, quali, ad esempio, le aree, i perimetri, imomenti angolari, etc; volendo mantenere, invece, una certa generalita, la classedi modelli per il confronto si dovra supporre non facilmente descrivibile comequella delle figure piane e si dovra, quindi, richiedere un metodo di riconosci-mento piu generale.Un metodo capace di un confronto di questo tipo ha innumerevoli possibilitadi applicazione. Una particolare di esse e la ricostruzione virtuale di affreschiframmentati in piccole parti, delle quali sia stato conservato un congruo nume-ro, assieme ad un’immagine dell’opera antecedente alla distruzione [54]. Infatti,non e possibile stabilire quale debba essere a priori l’orientamento di un fram-mento rispetto ad uno prefissato sull’immagine dell’affresco.

Il problema della ricostruzione virtuale e similare al problema di un “puzzle”.Intervengono pero delle ulteriori complicazioni dovute al fatto che i frammenti ingenere non costituiscono la totalita, le loro dimensioni sono piccole relativamentealla superficie intera dell’affresco e i bordi non rappresentano un’informazionesufficientemente precisa ed univoca. I dati che permettono di individuare unframmento dovranno allora essere basati esclusivamente sul contenuto pittoricoe la valutazione dovra essere indipendente dalla eventuale disposizione di al-tri frammenti, che potra costituire fattore discriminante solo successivamente,escludendo posizioni che potrebbero dare luogo a sovrapposizioni. Con tali pre-messe, questo problema offre un luogo adatto all’applicazione di un metodo diconfronto come sopra descritto.

Facciamo riferimento, ad esempio, agli affreschi contenuti nella Chiesa degliEremitani in Padova che furono gravemente danneggiati l’ 11 marzo 1944 a

4.2. LA RAPPRESENTAZIONE MATEMATICA 49

causa di un bombardamento. Parte di questa meravigliosa opera del Mantegna,Guariento, Pizolo e di altri illustri artisti, e stata frantumata in migliaia dipezzi, producendo uno dei danni culturalmente piu gravi subiti dall’Italia nellaseconda guerra mondiale. Si dispone di alcune foto in bianco e nero di questiaffreschi che risalgono alla campagna fotografica degli anni 1920 condotta daifratelli Alinari. Recentemente, con lo sviluppo tecnologico e la percezione dinuove possibilita, sono accresciuti il desiderio e la speranza di vedere ricostruitiquesti inestimabili affreschi [10, 54].

Figura 4.1: Frammenti degli affreschi.

Cercheremo ora di illustrare i dettagli di un modello matematico basatosu frame locali di cosiddette armoniche circolari (una variante delle frame diGabor, ma a geometria circolare) che rappresenti le immagini e ne consenta ilconfronto a meno di rotazioni e concluderemo con la definizione di un opportunocoefficiente di matching (corrispondenza) che ne permetta l’eventuale direttaapplicazione numerica al caso specifico degli affreschi.

4.2 La rappresentazione matematica

L’idea del metodo consiste nello sviluppare le immagini, considerate come fun-zioni di L2(Ω) ove Ω risulta essere un cerchio di R2, secondo una successioneortonormale costituita da invarianti per operatori di rotazione. In questo modol’immagine risulta identificata dai coefficienti dello sviluppo che saranno inva-rianti per rotazioni, a meno di prodotto per un autovalore di modulo unitario


del tipo eimθ. Un opportuno confronto dei coefficienti potrebbe permettere ladefinizione della classe di immagini che siano congrue (sovrapponibili) a menodi una rotazione ad una data.

Dopo aver stabilito l’esistenza di tale sviluppo su di un dominio continuo, lemaggiori difficolta alla applicazione del metodo alle funzioni campionate sonodovute per lo piu al fenomeno dell’ aliasing: piu funzioni possono ammettere unmedesimo campionamento (restrizione ad un reticolo discreto) e, quindi, rap-presentare degli “alias” di uno stesso oggetto. Nel nostro caso si verifica unaperdita di ortonormalita della base, che potrebbe, nel discreto, non avere piuil significato che nel continuo ci e garantito. Infatti, l’incertezza che si ha sullefunzioni necessariamente si propaga anche sui loro prodotti scalari in ℓ2 come sie descritto nel Teorema 2.2.7. In particolare, si rischia di non poter piu fruttarela proprieta di essere autofunzioni di operatori di rotazione. Allo scopo di ope-rare su immagini digitalizzate e, quindi, definite su reticoli discreti, occorreracautelarsi e verificare in quali condizioni sia da considerare lecito il passaggiodal continuo al discreto.

Data, infatti, una funzione continua f e considerata la sua restrizione (cam-pionamento) sul reticolo Z2, diciamola f c, l’incertezza sulla possiblita di poterricostruire in modo diretto ed univoco f a partire da f c e strettamente legata allasua trasformata di Fourier f . Se supponiamo unitario il passo di campionamen-to, allora l’incertezza sulla rappresentazione di f da parte del campionamentof c, dipende dalla “coda” di f , ossia la funzione ǫf che coincide con f fuoridi [−π, π)2 ed e nulla su [−π, π)2. Le funzioni che ammettono trasformata diFourier con supporto (banda) limitato in [−π, π)2 sono invece ben descritte, inquanto direttemente ricostruibili, dai propri campioni. Come conseguenza, an-che i prodotti scalari in ℓ2 saranno affetti da un’incertezza che dipendera dallecode delle trasformate di Fourier delle funzioni in gioco. In particolare, abbiamodimostrato nel Teorema 2.2.7 che la radice quadrata dell’ “energia” (integralequadrato) delle code risulta essere il parametro di giudizio caratteristico di taleincertezza.

Ricordando che per funzioni definite su compatti la trasformata di Fouriere restrizione di una funzione intera ad R2, e, quindi, non a supporto com-patto, si deduce che qualunque sviluppo venisse scelto, non potrebbe, in ge-nerale, mantenere la propria ortonormalita nel discreto. L’unica possibilitaper il campionamento di rispettare il prodotto scalare e di scegliere una ba-se che limiti le incertezze, almeno per una buona parte delle funzioni che lacompongono. La nostra ricerca si dova limitare a quelle basi che possano es-sere costituite da autofunzioni per operatori di rotazione. Dal Lemma 4.6.1riportato in Appendice si deduce che una generica base debba essere del tipodm,n(r, θ) = eimθRm,n(r),m ∈ Z, n ∈ N, ove (Rm,n)n∈N e una successione orto-normale nello spazio L2

r1/2 delle funzioni f tali che r1/2f(r) ∈ L2. Per questotipo di funzioni si riescono a dare delle stime (Appendice, Teoremi 4.6.5 e 4.6.7)dell’andamento delle loro trasformate di Fourier, che non sono abbastanza fi-


ni per garantire, in generale, che, restringendo il dominio al discreto, vengamantenuta l’ortonormalita e l’invarianza per rotazioni, compatibilmente al teo-rema del campionamento. Infatti se Tα e un operatore di rotazione di angolo α(Appendice) allora per ogni f si ha

〈Tαf, dm,n(r, θ)〉 = 〈f, T−αdm,n(r, θ)〉 = eimα〈f, dm,n(r, θ)〉.

Quindi una buona rappresentazione dei prodotti scalari anche nel discreto (Teo-rema 2.2.7) assicura una estensione della invarianza per rotazione anche perfunzioni campionate.

Il Corollario 4.5.8 suggerisce una particolare classe di basi, per cui la tra-sformata di Fourier si riesce a scrivere in forma esplicita e conseguentementefornisce uno strumento di indagine diretto.

Proseguiamo il nostro studio, allora, sulla particolare classe di queste basi,identificata da opportuni insiemi numerabili di reali positivi (ρm,n)m∈Z,n∈N, chesoddisfino all’equazione:

(∗) ρm,n1J′m(ρm,n1)Jm(ρm,n2) = ρm,n2J

′m(ρm,n2)Jm(ρm,n1),

ove Jm, e il coefficinte di Bessel di ordine m, mediante la relazione:

dm,n(r, θ) = gm,neimθJm(ρm,nr),

gm,n costante di normalizzazione. In particolare, la condizione (*) e strettamen-te connessa al fatto che il Laplaciano si comporti come un operatore autoag-giunto (hermetiano), rispetto al prodotto scalare in L2, sullo spazio vettorialegenerato dalle soluzioni di un problema agli autovalori analogo a quello descrittoin Appendice, con la condizione che le derivate di una soluzione nella direzionenormale al bordo siano proporzionali alla stessa. Come ulteriore giustificazione,quindi, per la scelta di tali basi, si ricordi che il Laplaciano e essenzialmentel’unico operatore differenziale a coefficienti costanti che commuti con gli opera-tori di rotazione e che il suo spettro, quando e autoaggiunto, caratterizza unasuccessione di autospazi ortogonali che risultano, cosı, invarianti per rotazioni.Una condizione necessaria e sufficiente affinche si verifichi (*) e che valga:

(∗∗) ρm,nJ′m(ρm,n) + σJm(ρm,n) = 0,

per ogni m ed n e σ costante fissata. Questa seconda equazione mista ammetteun’infinita numerabile di zeri reali positivi (Appendice, proposizione 6.7) laquale rende dm,n(r, θ) = gm,ne

imθJm(ρm,nr) una successione ortonormale che,al tendere di σ → ∞, diviene anche base completa di L2.

Ci si deve chiedere ora come queste particolari funzioni si comportino neiconfronti dell’aliasing. Il fatto che si considerino come definite su un cerchioe, quindi, essenzialmente a supporto compatto, ci porterebbe a concludere che,in generale, non v’e speranza di mantenere l’ortonormalita e l’invarianza perrotazioni. Allo scopo di salvaguardare lo sviluppo si potrebbero usare le cosid-dette “funzioni finestra” che convolute con una funzione a supporto compatto


0 1 2 3 4 5 6

w

-0.02

0

0.02

0.04

0.06

Figura 4.2: Parte radiale della TdF di d1,5,7 per una particolare scelta di ρ1,5

la approssimano e la rendono piu sicura rispetto all’aliasing, limitando le codedelle trasformate: tale soluzione pero risulta non necessaria e addirittura scon-veniente, una volta osservata una peculiarita delle trasformate di Fourier delle(dm,n)m∈Z,n∈N. Esse non sono a supporto compatto, ma sono pressoche nullefuori da un compatto, come mostriamo di seguito.Sia (dm,n)m∈Z,n∈N una generica base come descritta da (*) (o anche da (**)).Allo scopo di considerarla su domini di raggio arbitrario, si ponga per a > 0:

dm,n,a(r, θ) ≡1

adm,n(

r

a, θ).

Dai risultati riportati in Appendice, si ha che |dm,n,a(ρ, ϕ)| ammette un mas-

simo assoluto per ρ = ρm,n/a. Inoltre, |dm,n,a(ρ, ϕ)| decade all’infinito almenocome ρ−3/2 (Appendice) e, dunque, piuttosto rapidamente (Fig. 4.1).In base a queste osservazioni, allo scopo di limitare l’intervento di aliasing e laperdita di ortonormalita, si potrebbe pensare di scegliere tra le (dm,n,a)m∈Z,n∈N,quelle che ammettono il picco caratteristico della trasformata tutto contenutoin [0, π) e cioe con ρm,n+1/a < π. Alcune stime e delle dirette integrazioni nu-meriche confermano questa idea. Infatti, in questo modo, si limita l’ “energia”delle code e, conseguentemente, le incertezze sui prodotti scalari.Ci si potrebbe chiedere, ora, se, tra questi possibili sviluppi, ve ne sia uno chemeglio si adatta a mantenere una legalita nel passaggio dal continuo al discre-to. A questo proposito, si osservi che la distribuzione dei valori (ρm,n)m∈Z,n∈N,non dipende dalla particolare base scelta, in quanto ρm,n ∈ [jm,n, jm,n+1], ove(jm,n)n sono gli zeri positivi di Jm, come si prova studiando il comportamen-

to della funzione Jm(x)xJ′

m(x) (Appendice). Semplici osservazioni sulle distribuzioni

possibili mostrano, tuttavia, che il numero massimo di valori (ρm,n)m∈Z,n∈N checadono nell’intervallo [0, aπ) e raggiunto scegliendo come base le armoniche cir-colari (em,n,a)m∈Z,n∈N, cioe quella per cui ρm,n = jm,n. In particolare il fattoche |em,n,a(ρ, ϕ)| decada piu rapidamente, e cioe come ρ−5/2 (Appendice), sug-gerisce che il numero di elementi della base per cui la perdita di ortonormalita elimitata sotto una certa incertezza sia maggiore per le (em,n,a)m∈Z,n∈N. Questa


potrebbe, quindi, costituire una buona scelta per ottenere una massima rappre-sentazione dell’immagine digitale, senza pero dover subire la gravosa perdita diortonormalita e di invarianza per rotazioni.Una diretta integrazione numerica sullo spettro delle (em,n,a)m∈Z,n∈N mostracome l’errore sui prodotti scalari si evidenzi non appena jm,n+1/a oltrepassi lafrequenza limite di π.

2 4 6 8 10

n

00.10.20.30.40.50.60.7

1-En

a = 7

Figura 4.3: Energia della coda della Tdf di e1,n,7 al variare di n

Un’interessante analisi visuale mostra come, in questo caso, la perdita di or-tonormalita coincida con la perdita della simmetria circolare della trasformata diFourier. La peculiare forma a torre della trasformata di Fourier, caratterizzatada un picco e da un rapido decadimento, si inviluppa a causa della tassellaturaperiodica che la trasformata di Fourier discreta determina sul dominio, compor-tando la sovrapposizione delle code, come descritto dal teorema di Shannon edalla formula della somma di Poisson.

0 50100150200250050100150200250

0 50100150200250050100150200250

0 50100150200250050

100150200250

0 50100150200250050

100150200250

Figura 4.4: Re(e5,n,7) per valori di n crescenti

Risulta allora chiaro che limitare l’incertezza sui prodotti scalari dovuta alcampionamento equivale a mantenere la frequenza di picco caratteristica entroil limite di π e a conservare la simmetria circolare.

Tale equivalenza, puo essere ulteriormente illustrata con un esempio: tra learmoniche circolari quelle per cui

‖ǫem,n,a‖2 ≤ 0.3,


5 10 15 20 25m

246810

n

Coppie M-N a = 11

Figura 4.5: Spazio delle coppie (m,n) per cui ‖ǫem,n,11‖2 e limitato da 0.3

si ottiene una rappresentazione dell’insieme delle armoniche circolari per cuil’incertezza sui prodotti scalari e limitata da 0.3 per a = 11 come un insieme dicoppie (m,n), m > 0, come in Fig. 4.4.Calcolando, invece, l’insieme delle coppie (m,n), m > 0, tali per cui jm,n+1/a <π si ottiene una rappresentazione come in Fig. 4.5 e tali rappresentazioni sonopressoche equivalenti. Studiando i risultati numerici, si deduce che lo spazio sucui effettuare i confronti tra le immagini puo essere approssimativamente benrappresentato dall’insieme:

SpMNa = (m,n) ∈ Z × N : 0 ≤ |m| ≤ 2a, 0 < n ≤ −m/2 + a.

Il fatto di essersi ristretti precedentemente al caso m > 0 e, ora, di averintrodotto il modulo nella definizione di SpMNa e giustificato da:

e−m,n,a = (−1)mem,n,a,

conseguenza della nota equazione per i coefficienti di Bessel:

J−m(r) = (−1)mJm(r).

Le rappresentazione delle Figg. 4.4 e 4.5 sono quindi sufficienti a descriveretutto SpMNa.

5 10 15 20 25m

246810

n

Coppie M-N a = 11

Figura 4.6: Spazio delle coppie (m,n) per cui jm,n/a < π

Si noti che tale spazio ha dimensione circa pari a due volte il quadrato delraggio a del frammento. Di conseguenza il numero di elementi utili della basee direttamente proporzionale all’area del cerchio, per cui l’informazione rappre-sentata cresce ovviamente con il crescere del raggio.

4.3. LA DEFINIZIONE DEL MATCHING 55

Al fine di valutare la capacita delle armoniche circolari di rappresentare le im-magini digitali, proponiamo la definizione di un opportuno coefficiente, comeillustrato nel seguente paragrafo.

Figura 4.7: Spazio delle coppie (m,n) e le corrispondenti armoniche circolari

4.3 La definizione del matching

Allo scopo di poter localizzare le immagini di frammenti di un affresco su diuna sua rappresentazione fotografica, applichiamo i risultati presentati nei pre-cedenti paragrafi, cercando di trovare un criterio per stabilire se un’immagine econgrua ad un’altra a meno di rotazione, ricordando come sia impossibile risa-lire a priori all’orientamento che in origine doveva avere avuto un frammento,rispetto all’affresco, prima della distruzione.Ossevato che le immagini digitali vengono ben rappresentate da matrici bidi-mensionali, si considerino campionate tutte le funzioni che seguono e le si iden-tifichi con la matrice del loro campionamento. In base alle considerazioni fatte,un’immagine supportata su di un cerchio discreto di raggio a > 0 puo essereproiettata su uno spazio di dimensione finita, identificato dall’insieme SpMNa.Tale spazio ha la proprieta di avere come base (quasi) ortonormale alcuni in-varianti rispetto ad operatori di rotazione. In particolare se α ∈ [0, 2π) e Tα el’operatore di rotazione associato, si ha:

Tαem,n,a = eimαem,n,a, ∀(m,n) ∈ SpMNa.

La scelta di tale spazio ci garantisce che il passaggio dal continuo al discreto siaquasi esente da errori. Se FSpMNa rappresenta la proiezione su SpMNa di un

cerchio, diciamolo F , estratto dall’immagine di un frammento e FSpMNa la pro-iezione del suo ruotato di −α, diciamolo F , in modo da rispettare l’orientamento


originale, si ha:

FSpMNa ≈∑

(m,n)∈SpMNa

eimαfm,nem,n,a,

ove fm,n = (FSpMNa , em,n,a) e (·, ·) e il prodotto scalare in ℓ2(Z2). Ma, riscri-vendo la precedente:

FSpMNa ≈∑

0≤|m|≤2a

fm,

ove fm =∑

0<n≤−m/2+a eimαfm,nem,n,a. Se gi,j rappresentasse la proiezio-

ne su SpMNa di un cerchio estratto attorno ad una qualche posizione (i, j)dell’immagine dell’affresco G potremmo analogamente scrivere:

gi,j ≈∑

0≤|m|≤2a

gi,jm ,

ove gi,jm =∑

0<n≤−m/2+a gi,jm,nem,n,a, g

i,jm,n = (gi,j , em,n,a).

Ci si propone di poter confrontare i due sviluppi e di definire un coefficientedi matching invariante per rotazioni e per fattori scalari. Premettiamo, pero,alcuni risultati noti e delle importanti osservazioni.

Di seguito supporremo che le operazioni eseguite sulle matrici agiscano pun-tualmente, salvo differenti indicazioni esplicite. Si consideri, ora, estesa a zero lafm fuori del suo dominio in modo tale che essa sia rappresentata da una matricedelle stesse dimensioni di G e di traslarne il supporto in modo opportuno. Posto⋆ l’operatore di convoluzione per quanto osservato sulla ortonormalita della basesi hanno le seguenti uguaglianze e approssimazioni:

(fm ⋆ G)(i, j) = (fm, gi,j) ≈ (fm, g

i,jm ) ≈

∑

0<n≤−m/2+aeimαfm,n(g

i,jm,n).

Per una proprieta notevole della convoluzione si ha anche:

F−1[(F [fm]F [G])](i, j) = c(fm ⋆ G)(i, j),

ove F [·] e la trasformata di Fourier discreta e c > 0 una qualche costante chedipende dalla definizione di F [·] (in genere potremmo supporre c = 1). Questatecnica e molto usata per trasportare un’immagine lungo la superficie di un’al-tra ed eseguire il confronto con la parte sottostante.

Allo scopo di specificare il metodo al caso concreto, si osservi che le im-magini dei frammenti e dell’affresco, divise per canale (RGB), possono essereconsiderate come funzioni a valori reali. Si deve, quindi, notare che i terminicon m < 0 dello spazio SpMNa sono, per immagini reali, ridondanti ed essen-zialmente ricavabili dalla proiezione sulla parte di SpMNa per cui m > 0. Vale,infatti:

(h, e−m,n,a) = (−1)m(h, em,n,a),

4.3. LA DEFINIZIONE DEL MATCHING 57

per ogni funzione h a valori reali definita sul cerchio. Eviteremo, quindi, diconsiderare le e−m,n,a, ottenendo di ridurre a meta lo spazio su cui fare i con-fronti e conseguentemente, come vedremo, i tempi di calcolo. Alcuni risultatinumerici rivelano, inoltre, che, mediamente, le proiezioni sullo spazio generatoda e0,n,a sono le piu importanti in termini di peso, rispetto a tutte le altre. Purcostituendo e0,n,a l’unica vera e propria componente invariante per rotazionidello sviluppo, definiremo un coefficiente di matching che non terra conto deicontributi dati dalle proiezioni su e0,n,a. Infatti, rischieremmo di far pesare trop-po una componente rispetto a tutte le altre, comunque decisive ed importanti.Osservato questo, diamo ora la seguente

Definizione 4.3.1. Sia k ∈ N\0 e z, w ∈ C\0.Si consideri l’operatore binario non commutativo e non associativo ⊕k definitoda:

z ⊕k w = z + w(z

|z| )1−k.

Diremo somma circolare di ordine m ∈ N\0 di z1, ..., zm ∈ C\0 l’operatoredefinito da: ⊕m

k=1 zk = z1, m = 1⊕mk=1 zk = (

⊕m−1k=1 zk) ⊕m zm m > 1

.

Figura 4.8: Comportamento della somma circolare in presenza di rumore

Siamo ora pronti a definire il matching:

Definizione 4.3.2. Nelle notazioni sin qui portate, si definisce la matrice dimatching tra F e G

MFG =

1∑0<m≤2a ‖fm‖2

2

2a⊕

m=1

(‖fm‖2

(‖gi,jm ‖2)i,jF−1[(F [fm]F [G])],

ove ‖fm‖2 e la norma di fm in ℓ2(Z2) e (‖gi,jm ‖2)i,j e la matrice delle norme inℓ2(Z2) degli gi,jm .

Per ogni (i, j) punto dell’immagine G, la matrice di matching restituisce ilcoefficiente di matching dato da 0 < |MF

G (i, j)| ≤ 1, indice della corrispondenzacon il frammento a meno di rotazioni. Supponendo, infatti, che esista (i, j)posizione in G, tale che:

gi,j = FSpMNa ,


si ottiene:F−1[(F [fm]F [G])](i, j) ≈ eimα‖fm‖2

2,

per ogni 0 < m ≤ 2a.. Ma allora si ha:

MFG (i, j) ≈ 1∑

0<m≤2a ‖fm‖22

2a⊕

m=1

eimα‖fm‖22 = eiα.

Si deduce che condizione necessaria affinche una posizione (i, j) sia quella origi-naria del frammento e che:

|MF

G (i, j)| ≈ 1Arg(MF

G (i, j)) ≈ α

Questo fatto, che non costituisce condizione sufficiente, non e comunquefacile a verificarsi, in quanto si dovrebbe avere

|(fm, gi,jm )| ≈ c‖fm‖22,

per ogni 0 < m ≤ 2a e per un qualche c > 0 costante. Se a e abbastanza grande,lo spazio generato, ad esempio, dalle funzioni (e1,n,a)(1,n)∈SpMNa

ha dimensionea − 1 e due vettori presi casualmente in un tale spazio hanno una probabilitamolto bassa di essere proporzionali. La condizione diventa, dunque, anche mol-to stringente, non appena a >> 0.

Il metodo proposto e stato implementato con linguaggi C e C++, su variepiattaforme ed e stato applicato nel caso degli affreschi della Chiesa degli Eremi-tani in Padova, offrendo dei risultati notevoli in termini di capacita di riconosci-mento. Si rimanda per maggiori dettagli ad una piu completa argomentazione[?].

4.4 Un’ipotesi di parallelizzazione

Le immagini sono uno degli insiemi di dati piu gravosi da un punto di vistadell’eleborazione, richiedendo enormi risorse sia in termini di velocita sia in ter-mini di memoria. I tempi di calcolo sono uno dei fattori importanti di cui tenereconto nei problemi che coinvolgono le immagini.

Allo scopo di parallelizzare il calcolo e conseguentemente di ridurre i tempidi elaborazione, e importante evidenziare quali siano le operazioni effettuabili inmodo indipendente e quali invece no. Il fatto che ⊕k non sia, in generale, com-mutativo e associativo impone un ordine nel calcolo della matrice di matching,almeno per la definizione che si e data in 4.3.2. La somma circolare di ordinek e, quindi, subordinata alla somma circolare di ordine k− 1 dei termini prece-denti. Invece il calcolo di ogni singolo termine non e vincolato da quello deglialtri. Si deve inoltre osservare che i tempi di esecuzione di una convoluzione deltipo F−1[(F [fm]F [G])] e il calcolo della matrice delle norme, sono in generale

4.4. UN’IPOTESI DI PARALLELIZZAZIONE 59

Figura 4.9: Alcuni frammenti ricomposti

preponderanti rispetto al calcolo delle somme e della trasmissione dei dati, percui una parallelizzazione in questo senso permetterebbe di diminuire l’attesadi un fattore che possiamo supporre prossimo a 2a, nel caso si disponga di unpari numero di esecutori. In particolare si potrebbe incaricare il processore checalcola il termine per m = 1 di raccogliere i risultati provenienti dagli altri edeseguire le somme. Questa idea si puo generalizzare ad un numero N < 2a diprocessori che sia divisore intero di 2a. Ciascuno di essi sara incaricato del cal-colo di 2a

N termini della somma, ed uno, in particolare, di raccogliere i risultatiin ordine ed eseguire le somme. In questo caso, pero, le attese per lo scambio deirisultati tra esecutori rischiano di aumentare. Allo scopo di diminuire i tempi sipotrebbero distribuire i termini da calcolare in modo che, nel momento in cuitutti gli esecutori hanno finito di operare un termine, si rendano disponibili iprossimi 2a

N termini della somma nell’ordine opportuno, eseguendo, per cosı dire,un calcolo per righe. In questo caso verranno effettuate (2a− 2a

N ) trasmissioni edunque supporremo ancora che il grado di parallelizzazione si approssimi a N .Si deduce che nel caso limite di N = 2a la parallelizzazione sia massima, macon il maggiore dispendio di tempo in trasmissione. Nel caso si disponga di unamacchina parallela, come ad esempio il Cray T3E del CINECA, si potrebberovalutare i tempi medi effettivi e calcolare quale possa essere la soluzione ottimalein termini del numero processori. Precedenti esperimenti fatti sul Cray, hannomostrato che i tempi di trasmissione crescono mediamente in modo lineare conla quantita di informazione scambiata. Questo potrebbe far supporre che lascelta per N = 2a, sia, in effetti, comunque la migliore.


Figura 4.10: Struttura ad albero per la parallelizzazione

4.5 Appendice al capitolo

La presente Appendice raccoglie brevemente i risultati citati nel corso dellalezione.

Lemma 4.5.1. Sia α ∈ [0, 2π[ e sia Tα operatore di rotazione associato ad αdalla relazione

Tα[u(r, θ)] ≡ u(r, θ + α),

ove u e definita su R2 in coordinate polari. Allora u e autofunzione per Tα, perogni α ∈ [0, 2π[ se e soltanto se u e del tipo:

u(r, θ) ≡ ei(θm)R(r),

per un qualche m ∈ Z, e una R funzione di [0,∞[ a valori in C.

Si consideri la nota descrizione del Laplaciano in coordinate polari:

∆ = ∂2r +

1

r∂r +

1

r2∂2ϕ.

Proposizione 4.5.2. Ci proponiamo di risolvere il problema agli autovaloridato da:

(∗)

(∂2r + 1

r∂r + 1r2 ∂

2ϕ)fm(r, ϕ) = −λmfm(r, ϕ),

fm|∂Ω = 0.

4.5. APPENDICE AL CAPITOLO 61

ove Ω e un cerchio di R2. A tale scopo consideriamo fm(r, ϕ) = Rm(r)Φm(ϕ)a variabili separate e sostituendo fm quale proposta soluzione di (*) si ottiene:

r2R′′m(r)

Rm(r)+ r

R′m(r)

Rm(r)− λmr

2 = c = −Φ′′m(ϕ)

Φm(ϕ).

Si prova che la costante di separazione c e del tipo c = m2 con m ∈ Z. Si trovache

Φm(ϕ) = (2π)−1/2eimϕ,

e che Rm e soluzione dell’equazione:

(∗∗)(∂2r +

1

r∂r + (λm − m2

r2))Rm = 0.

E facile accorgersi della relazione che intercorre tra la precedente equazione (**)e l’equazione di Bessel e possiamo trovare la soluzione generica per (**) datada:

Rm(r) = AmJm(λ1/2m r) +BmYm(λ1/2

m r),

ove Jm e Ym sono le funzioni di Bessel di ordine m, di prima e seconda specie

rispettivamente. Osserviamo ora che Ym(λ1/2m r) ha una singolarita in 0 e che

rY 2m(λ

1/2m r)non e sommabile e, quindi, affinche Rm soddisfi alla richiesta fm ∈

L2(Ω) deve essere Bm = 0. Le condizioni al bordo di annullamento impongono

inoltre che Jm(λ1/2m ) = 0 e conseguentemente che λ

1/2m sia zero della funzione

di Bessel Jm. Come si puo dimostrare Jm ammette una infinita numerabile dizeri semplici reali.

-1-0.5

00.5

1 -1

-0.5

0

0.51

-0.500.51

-1-0.5

00.5

1

Figura 4.11: Parte reale dell’armonica circolare di ordine (1,2)

L’effetto di tale condizione e:

λm,n = (jm,n)2,

ove m ∈ Z, n ∈ N e jm,n e l’n-esimo zero della funzione di Bessel Jm. Questacondizione permette di ottenere gli autovalori e le autofunzioni del laplacianosul cerchio in condizioni di Dirichlet come:

em,n = cm,nJm(jm,nr)eimϕ, m ∈ Z, n ∈ N.


La costante cm,n viene introdotta come normalizzazione rispetto al prodottoscalare in L2 e si prova che:

cm,n = π− 12 [dJm(s)

ds|s=jm,n ]−1.

La successione (em,n)m∈Z,n∈N definisce una base completa per L2(Ω) che sap-piamo essere costituita da invarianti rispetto agli operatori di rotazione. Chia-meremo tali funzioni armoniche circolari. Inoltre, non e difficile dimostrare cheil sistema (ekm,n)m,k∈Z,n∈N definito da

ekm,n(x) = em,n(x− k)

e una frame per L2(R2). Tale sistema generalizza le frame di Gabor bidimen-sionali, ma su una geometria circolare.

0 10 20 30 400

10

20

30

40

0 10 20 30 400

10

20

30

40

Figura 4.12: L’occhio della Gioconda ricostruito con la frame (ekm,n)m,k∈Z,n∈N.

L’intervento degli zeri jm,n della funzione di Bessel Jm e dovuto alla par-ticolare scelta della condizione di annullamento al bordo. Questa condizionepuo essere sostituita dalla condizione piu generale che le derivate parziali nel-la direzione della normale positiva siano proporzionali alla funzione stessa. Inparticolare, nel nostro caso tale richiesta e soddisfatta se e soltanto se

λ1/2m J ′

m(λ1/2m ) + σJm(λ1/2

m ) = 0,

per ogni m. Con tale condizione il Laplaciano si mantiene autoaggiunto e comesoluzione del problema agli autovalori si ottiene un insieme ortonormale deltipo dm,n(r, θ) = gm,ne

imθJm(ρm,nr) dove ρm,n sara una radice dell’equazionerJ ′m(r)+σJm(r) = 0 e gm,n una opportuna costante di normalizzazione. Questa

classe di basi e detta di Dini di cui fanno parte quelle che abbiamo chiamatoarmoniche circolari.

Proposizione 4.5.3. Per ogni σ costante fissata e n ∈ N\0, esiste uno edun solo zero dell’equazione rJ ′

m(r) + σJm(r) = 0 nell’intervallo [jm,n, jm,n+1].

Definizione 4.5.4. Sia f una funzione di [0,+∞) in C. Diremo trasformatadi Hankel di ordine m della funzione f, ogni qual volta sia definita, la funzione:

Hm(f)(ρ) =

∫ +∞

0

rf(r)Jm(ρr)dr.

4.5. APPENDICE AL CAPITOLO 63

Teorema 4.5.5. (Riduzione della Trasformata di Fourier) Siano m ∈ Z Fm ∈L1(R2), Fm(r, θ) = f(r)eimθ con f funzione di [0,+∞) a valori in C. Alloravale l’uguaglianza:

Fm(ρ, ϕ) ≡ (−i)meimϕHm(f)(ρ).

Sia (dm,n)m∈Z,n∈N una generica base di Dini. Allo scopo di considerarla sudomini di raggio arbitrario, si ponga per a > 0:

dm,n,a(r, θ) ≡1

adm,n(

r

a, θ).

Proposizione 4.5.6. Siano a > 0 e (dm,n,a)m∈Z,n∈N una base di Dini definitasu B(0,a). Allora

dm,n,a(ρ, ϕ) = eimϕ(−i)magm,nρm,nJm−1(ρm,n)Jm(aρ) − (ρa)Jm−1(ρa)Jm(ρm,n)

(aρ)2 − ρ2m,n

.

Teorema 4.5.7. Sia f una funzione definita su [0,+∞) a valori in C, tale che

r2−k+1

2 f(r) ∈ Lk([0,+∞)), k = 1, 2, e tale che f(r) = 0 per r > α > 0. Sia Fla funzione di R2 in C definita da:

F (r, θ) = f(r).

Vale allora la stima:

|Hm(f)(ρ)| ≤ ‖F‖22

2π[α2

2(J2m(αρ) − Jm−1(αρ)Jm+1(αρ))].

Corollario 4.5.8. Nelle ipotesi e premesse del precendente teorema si ha chela stima e raggiunta ad un dato ρ = ρ0 se e soltanto se

f(r) =

AJm(rρ0), r ≤ α0, r > α

e

Hm(f)(ρ) = Aαρ0Jm(ρα)Jm−1(ρ0α) − ρJm(ρ0α)Jm−1(ρα)

ρ2 − ρ20

,

ove A e una costante non nulla. In particolare, calcolando il limite per ρ 7→ ρ0

si ha

Hm(f)(ρ0) =‖F‖2

2α2

4π[J2m(αρ0) + J2

m+1(αρ0)] −‖F‖2

2mα

2πρ0Jm(αρ0)Jm+1(αρ0).

Capitolo 5

Elementi della teoria delleondine (wavelet)

I n questa lezione vogliamo presentare elementi di base della teoria delle ondine(wavelets). Le ondine sono funzioni a forma di onda che decadono rapidamenteall’infinito (da cui il nome) le cui dilatazioni e traslazioni generano basi e frameper L2(Rd). La prima parte della lezione e dedicata a chiarire la diffrenza tral’analisi eseguita da questo tipo di espansioni e le frame di Gabor precendente-mente introdotte. In particolare descriveremo la trasformata wavelet continua ediscuteremo l’esistenza di frame di wavelet che la discretizzino. Infine introdur-remo elementi della teoria della analisi multirisolutiva e la costruzione di basiortonormali di wavelet. Indichiamo come approfondimenti [11, 19, 62, 68, 69].

5.1 Trasformate tempo-frequenza e tempo-scala

Per semplicita nel seguito assumiamo d = 1. Spiegeremo in seguito come gene-ralizzare i risultati a dimensioni superiori.

Iniziamo questa sezione richiamando il concetto di trasformata di Gabor:

Vg(f)(b, w) = 〈f,MwTbg〉 =

∫

R

f(t)e−2πiw·tg(t− b)dt. (5.1)

In analogia con la trasformata di Fourier per cui vale 〈Ff1,Ff2〉 = 〈f1, f2〉, perf1, f2 ∈ L2(R), si ha il seguente importante risultato.

Proposizione 5.1.1 (Relazioni di ortogonalita per la STFT). Siano f1, f2, g ∈L2(R). Si ha Vg(fj) ∈ L2(R2), j = 1, 2, e

〈Vg(f1), Vg(f2)〉 = ‖g‖22〈f1, f2〉.

65

66CAPITOLO 5. ELEMENTI DELLA TEORIA DELLE ONDINE (WAVELET)

Dimostrazione. Si assuma fj , g ∈ L1 ∩ L∞. Il risultato generale segue perdensita. Si ha fjTxg ∈ L2 per ogni x ∈ R. Quindi, per la formula di Plancherelis ha

∫

R

∫

R

Vgf1(x, ω)Vgf2(x, ω)dωdx

=

∫

R

(∫

R

F(f1Txg)(ω)F(f2Txg)(ω)dω

)dx

=

∫

R

(∫

R

f1(t)f2(t)|g(t− x)|2dt)dx

Per le nostre assunzioni possiamo applicare Fubini-Tonelli e si ha

〈Vg(f1), Vg(f2)〉 =

∫

R

f1(t)f2(t)

(∫

R

|g(t− x)|2dx)dt

= ‖g‖22〈f1, f2〉.

Corollario 5.1.2. Siano f, g ∈ L2(R). Allora

‖Vg(f)‖2 = ‖g‖2‖f‖2.

In particolare se ‖g‖2 = 1 allora la STFT e un isomorfismo isometrico da L2(R)sulla propria immagine.

Queste proprieta di isomorfismo delle STFT garantiscono che essa rappresen-ti in modo unico e fedele una funzione. Tuttavia, come abbiamo precentementeosservato mediante un’analisi del principio di indeterminazione, la STFT pro-duce un’analisi fortemente ridondante. Inoltre, sempre a causa del principio diindeterminazione, se si utilizza una funzione finestra g con supporto largo, allorala STFT avra una buona risoluzione delle alte frequenze. Infatti si avra che Fge fortemente localizzata. Se viceversa avremo g molto localizzata, la STFT avrauna buona risoluzione del tempo e delle basse frequenze, ma sara “trasparente”alle alte frequenze.

Allo scopo di migliorare la risoluzione dello spazio tempo-frequenza per tuttele frequenze contemporaneamente, si potrebbe volere una finestra il cui supportosi adatta alla frequenza che vogliamo analizzare, cioe vorremo che essa venga“ristretta” propozionalmente alla frequenza. Seguendo questo ragionamentopotremmo modificare la STFT come segue:

V δg (f)(b, w) =

∫

R

f(t)e−2πiw·tg(w(t− b))wdt. (5.2)

A questo punto, introduciamo l’operatore di dilatazione

Daf(t) :=1

|a|1/2 f(t/a), a ∈ R,

5.1. TRASFORMATE TEMPO-FREQUENZA E TEMPO-SCALA 67

facciamo alcuni semplici passaggi e otteniamo

V δg (f)(b, w) =

∫

R

f(t)e−2πiw·tg(w(t− b))dt

= e−2πiw·b∫

R

f(t)TbD1/w(e−2πitg(t))dt.

Se ci dimentichiamo del fattore complesso unitario davanti all’integrale e chia-miamo ψ(t) := e−2πitg(t) si ottiene che f viene analizzata dalla precedentetrasformazione mediante traslazioni e dilatazioni di ψ. Si osservi inoltre cheψ(t) = cos(2πt)g(t) + i sin(2πt)g(t), e dunque ha forma di onda. Le sue dilata-zioni hanno come effetto quello di aumentare la frequenza dell’onda e di variarneil supporto. Sulla base di queste osservazioni si definisce la trasformata waveletcontinua (CWT) data da

Wψ(f)(b, a) :=

∫

R

f(t)TbDaψ(t)dt. (5.3)

In analogia con la trasformata di Fourier e la STFT si ha una corrispondenterelazione di ortogonalita.

Proposizione 5.1.3 (Relazioni di ortogonalita per la CWT). Siano f, g ∈L2(R) and ψ ∈ L2(R) tale che

cψ = 2π

∫

R

|ξ|−1|ψ(ξ)|2dξ <∞. (5.4)

Allora si ha Wψ(g),Wψ(f) ∈ L2a(R

2), e

〈Wψ(f),Wψ(g)〉 = cψ〈f, g〉.

Dimostrazione. Per la formula di Plancherel, si ha

∫

R

∫

R

Wψf(x, a)Wψf(x, a)dadx

a2

=

∫

R

∫

R

(∫f(ξ)|a|1/2e−ixξψ(aξ)dξ

)(∫g(ξ)|a|1/2e−ixξψ(aξ)dξ

)dadx

a2

Se definiamo

Fa(ξ) := |a|1/2f(ξ)ψ(aξ), e Ga(ξ) := |a|1/2g(ξ)ψ(aξ),


si ha, per Plancherel nuovamente e Fubini-Tonelli,∫

R

∫

R

Wψf(x, a)Wψf(x, a)dadx

a2

= 2π

∫

R

∫

R

Fa(x)Ga(x)dadx

a2

= 2π

∫

R

∫

R

Fa(ξ)Ga(ξ)dadξ

a2

= 2π

∫

R

∫

R

f(ξ)g(ξ)|ψ(aξ)|2 dadξ|a|

= 2π

∫

R

f(ξ)g(ξ)

(∫

R

|ψ(aξ)|2 da|a|

)dξ

= cψ〈f, g〉.

Corollario 5.1.4. Siano f ∈ L2(R) and ψ ∈ L2(R) tale che

cψ = 2π

∫

R

|ξ|−1|ψ(ξ)|2 <∞.

Allora si ha Wψ(f) ∈ L2a(R

2), e

‖Wψ(f)‖2,a =√cψ‖f‖2.

5.2 Frame di ondine

Come nel caso della STFT, anche la trasformata wavelet continuaWψ(f) puo es-sere campionata su un insieme discreto e i suoi campioni utilizzati per rappresen-tare univocamente la funzione f . In modo simile, prodecedermo dimostrando l’e-sistenza di frame generate per traslazione e dilatazione della ψ. Vogliamo studia-

re in particolare insiemi di funzioni del tipo (ψ, a0, b0) := a−m/20 ψ(a−m0 x−nb0) :m,n ∈ Z ove ψ ∈ L2(R) e a0, b0 > 0. Come notazione, scriveremo ψm,n(x) :=

a−m/20 ψ(a−m0 x − nb0) = Dam

0Tb0nψ. Valuteremo per quali condizioni sulla fun-

zione ψ, detta ancora ondina madre, e sui parametri a0, b0 l’insieme (ψ, a0, b0)definisce una frame ovvero quando i valori Wψ(f)(a−m0 nb0, a

−m0 ) = 〈f, ψm,n〉

completamente caratterizzano f e rappresentano una digitalizzazione stabileper f .

5.2.1 Condizioni necessarie

Il seguente risultato e dovuto a Daubechies, [19, Theorem 3.3.1].

Proposizione 5.2.1. Se (ψ, a0, b0) definisce una frame per L2(R) con costantidi frame A,B > 0 allora

b0 ln a0

2πA ≤

∫ ∞

0

|ξ|−1|ψ(ξ)|2dξ ≤ b0 ln a0

2πB, (5.5)

5.2. FRAME DI ONDINE 69

e

b0 ln a0

2πA ≤

∫ 0

−∞|ξ|−1|ψ(ξ)|2dξ ≤ b0 ln a0

2πB, (5.6)

5.2.2 Condizioni sufficienti

Con una dimostrazione simile a quella prodotta per il Teorema 3.2.9 si provanole seguenti condizioni sufficienti.

Proposizione 5.2.2. Se ψ e a0 sono tali che

inf1≤|ξ|≤a0

∞∑

m=−∞|ψ(am0 ξ)|2 > 0, (5.7)

sup1≤|ξ|≤a0

∞∑

m=−∞|ψ(am0 ξ)|2 <∞, (5.8)

e se β(s) = supξ∑∞

m=−∞ |ψ(am0 ξ)||ψ(am0 ξ+s)| decade almeno come (1+|s|)1+ε,con ε > 0, allora esite un b0 > 0 tale che (ψ, a0, b0) e un frame per L2(R) perogni 0 < b0 ≤ b0. Per 0 < b0 ≤ b0 si hanno le seguenti espressioni per le costantiA,B > 0 di frame:

A =2π

b0

inf1≤|ξ|≤a0

∞∑

m=−∞|ψ(am0 ξ)|2 −

∑

k 6=0

[β

(2π

b0k

)β

(−2π

b0k

)]1/2

,

B =2π

b0

sup1≤|ξ|≤a0

∞∑

m=−∞|ψ(am0 ξ)|2 +

∑

k 6=0

[β

(2π

b0k

)β

(−2π

b0k

)]1/2

.

Le condizioni sono soddisfatte non appena |ψ(ξ)| ≤ C|ξ|α(1 + |ξ|)−γ con α > 0,γ > α+ 1.

Dimostrazione. Dobbiamo dimostrare che esistono A,B > 0 tali che

A‖f‖22 ≤

∑

m,n

|〈f, ψm,n〉|2 ≤ B‖f‖22,

per ogni f ∈ L2. Dal momento che 1√τe2πinx/τn∈Z e base ortonormale per


H = L2(0, τ), scriviamo

∑

m,n

|〈f, ψm,n〉|2 =∑

m,n

∣∣∣∣∫f(ξ)a

m/20 ψ(am0 ξ)e

ib0am0 nξ

∣∣∣∣2

dξ

=∑

m,n

am0

∣∣∣∣∣

∫ 2πb−10 a−m

0

0

eib0am0 nξ

∑

ℓ∈Z

f(ξ + 2πℓb−10 a−m0 )a

m/20 ψ(am0 ξ + 2πℓb−1

0 )

∣∣∣∣∣

2

dξ

=2π

b0

∑

m

∫ 2πb−10 a−m

0

0

∣∣∣∣∣∑

ℓ∈Z

f(ξ + 2πℓb−10 a−m0 )a

m/20 ψ(am0 ξ + 2πℓb−1

0 )

∣∣∣∣∣

2

dξ

=2π

b0

∑

m,k

∫f(ξ)f (ξ + 2πka−m0 b−1

0 )ψ(am0 )ψ(am0 ξ + 2πkb−10 )dξ

=2π

b0

∫|f(ξ)|2

∑

m

|ψ(am0 ξ)|2 + R(f).

Vogliamo ora limitare dall’alto e dal basso il resto R(f). Per Cauchy-Schwarz eun cambio di variabile si ha

|R(f)| =

∣∣∣∣∣∣2π

b0

∑

m,k 6=0

∫f(ξ)f(ξ + 2πka−m0 b−1

0 )ψ(am0 )ψ(am0 ξ + 2πkb−10 )dξ

∣∣∣∣∣∣

≤ 2π

b0

∑

m,k 6=0

[∫|f(ξ)|2|ψ(am0 )||ψ(am0 ξ + 2πkb−1

0 )|dξ]1/2

×[∫

|f(ξ)|2|ψ(am0 )||ψ(am0 ξ − 2πkb−10 )|dξ

]1/2

Applicando ora Cauchy-Schwarz alla somma su m is ottiene

|R(f)| ≤ 2π

b0

∑

k 6=0

[∫|f(ξ)|2

∑

m

|ψ(am0 )||ψ(am0 ξ + 2πkb−10 )|dξ

]1/2

×[∫

|f(ξ)|2∑

m

|ψ(am0 )||ψ(am0 ξ − 2πkb−10 )|dξ

]1/2

≤ 2π

b0‖f‖2

2

∑

k 6=0

[β

(2π

b0k

)β

(−2π

b0k

)]1/2.

Con queste stime si ha subito che se 0 < A,B < ∞ allora (ψ, a0, b0) e unaframe. Si lascia per esercizio dimostrare che le condizioni su ψ e a0 assicuranoche esiste b0 > 0 tale che per ogni 0 < b0 ≤ b0 le costanti A,B sono davverocome volute.

5.3. ANALISI MULTIRISOLUTIVA E ONDINE ORTONORMALI 71

5.3 Analisi multirisolutiva e ondine ortonormali

Abbiamo visto che per le frame di Gabor esiste una teoria discreta corrispon-dente che permette di realizzare un’analisi numerica efficiente. Per le frame diondine non esiste invece in generale una versione discreta che approssimi benequella continua. Il problema sorge per il fatto the le dilatazioni trasformanofunzioni definite su reticoli in funzioni definite su reticoli diversi e, a meno dinon lavorare su insiemi discreti con topologie “esotiche” (topologie a0-adiche),l’analisi puo risultare assai complessa. Per poter realizzare un’analisi numericaper le wavelet occorre completamente cambiare punto di vista ed introdurre unostrumento in grado di gestire discretizzazioni che vivono su insiemi a diversa ri-soluzione, possibilmente tra loro innestati. Questo strumento innovativo e statointrodotto da Meyer e Mallat [19, 68, 69] e si chiama analisi multirisolutiva(MRA, multiresolution analysis).

Una analisi multirisolutiva consiste di una successione di spazi approssimantiVj . Piu precisamente, i sottospazi chiusi Vj ⊂ L2 soddisfano

. . . V2 ⊂ V1 ⊂ V0 ⊂ V−1 ⊂ V−2 ⊂ . . .

con ⋃

j∈Z

Vj = L2(R),

⋂

j∈Z

Vj = 0.

Si denotano con Pj the proiezioni ortogonali su Vj e dunque si ha f = limj→∞ Pjfper ogni f ∈ L2(R). Ovviamente esistono moltissimi possibili spazi con que-ste proprieta che non hanno nulla a che fare con proprieta di multirisoluzione.Queste ultime sono rese possibili dalle ulteriori condizioni

f ∈ Vj ⇔ f(2j·) ∈ V0, ∀j ∈ Z,

ef ∈ V0 ⇔ f(· − n) ∈ V0, ∀n ∈ Z.

Dal momento che V0 e uno spazio invariante per traslazioni risulta naturale as-sumere che esso abbia una base costituita da traslazioni di una singola funzioneϕ ∈ V0, cioe assumeremo che ϕ(· − n) : n ∈ Z e una base ortonormale per V0.E anche possibile rilassare la condizione di ortonormalita chiedendo “meno”.Di recente si sono studiate analisi multirisolutive che generano frame (e nonbasi), ma questo argomento va oltre gli scopi di questa lezione. Se denotiamoϕj,n(t) = 2−j/2ϕ(2−jt− n), si dimostra facilmente dalle proprieta di multiriso-luzione che ϕj,n : n ∈ Z e una base per Vj .

Per definizione Vj ⊂ Vj−1 e dunque esiste uno spazio complementare (cioeortogonale a Vj in Vj−1) Wj tale che

Vj−1 = Vj ⊕Wj .


Come vedremo, potremmo definire una base per Wj definita da funzioni deltipo ψj,n(t) = 2−j/2ψ(2−jt − n), ove ψ e un’ondina madre. Per le proprieta dimultirisoluzione si mostra facilmente che

L2(R) =⊕

j∈Z

Wj ,

e dunque che ψj,n : j, n ∈ Z e una base per L2(R).

Vogliamo chiarire ora come, a partire da ϕ, si possa costruire una ψ chegeneri per dilatazioni diadiche e traslazioni intere una base per L2(R). Dalmomento che V0 ⊂ V−1 la funzione ϕ deve essere necessariamente rappresentatada una combinazione lineare di φ(2 · −n) e cioe deve essere

ϕ(t) =√

2∑

n

hnϕ(2t− n). (5.9)

La precedente equazione e detta relazione di scala e la funzione ϕ e dunquechiamata funzione scalante. Osserviamo la struttura di ϕ, computandone latrasformata di Fourier. Si ha

ϕ(ξ) = m0(ξ/2)ϕ(ξ/2), (5.10)

ove

m0(ξ) =1√2

∑

n

hne−inξ, (5.11)

che e una funzione 2π-periodica. Vediamone alcune proprietaa.

Lemma 5.3.1. Per quasi ogni ζ ∈ R si ha

|m0(ζ)|2 + |m0(ζ + π)|2 = 1. (5.12)

Dimostrazione. Per l’ortonormalita delle ϕ(· − k) si ha

δk,0 =

∫ϕ(x)ϕ(x − k)dx =

∫|ϕ(ξ)|2eikξdξ

=

∫ 2π

0

eikξ∑

ℓ∈Z

|ϕ(ξ + 2πℓ)|2dξ

Dal momento che 1√τe2πinx/τn∈Z e base ortonormale per L2(0, τ), si ha l’e-

quazione fondamentale ∑

ℓ∈Z

|ϕ(ξ + 2πℓ)|2 =1

2π, (5.13)

con validita quasi ovunque. Se sostituiamo (5.10) nella precedente otteniamo

∑

ℓ∈Z

|m0(ξ + πℓ)|2|ϕ(ξ + πℓ)|2 =1

2π.

5.3. ANALISI MULTIRISOLUTIVA E ONDINE ORTONORMALI 73

Spezzando la somma sugli ℓ pari e sugli ℓ dispari, usando la periodicita di m0 eapplicando (5.13) ancora una volta, si ha

|m0(ζ)|2 + |m0(ζ + π)|2 = 1.

Si lascia per esercizio verificare questi passaggi.

Con questo strumento, siamo ora in grado di caratterizzare gli elementi diW0 = V ⊥

0 . Si ha il seguente

Lemma 5.3.2. Se f ∈W0 allora

f(ξ) = eiξ/2m0(ξ/2 + π)ν(ξ)ϕ(ξ/2), (5.14)

per una qualche ν funzione 2π-periodica.

Dimostrazione. Se f ∈W0 allora f =∑n fnϕ−1,n e

f(ξ) =1√2

∑

n

fne−inξ/2ϕ(ξ/2) = mf (ξ/2)ϕ(ξ/2), (5.15)

ove

mf (ξ) =1√2

∑

n

fne−inξ/2.

Il vincolo f ∈ V ⊥0 implica

0 =

∫f(ξ)ϕ(ξ)eikξdξ

=

∫ 2π

0

eikξ∑

ℓ∈Z

f(ξ + 2πℓ)ϕ(ξ + 2πℓ)dξ.

Sostituendo alla precedente (5.10) e (5.15), raggruppando le somme per ℓ paree dispari, e usando (5.13) si ottiene

mf (ξ)m0(ξ) +mf (ξ + π)m0(ξ + π) = 0. (5.16)

Dal momento che m0(ξ) e m0(ξ+π) non si annullano contemporaneamente per(5.12), esiste una funzione λ(ξ) 2π-periodica tale che

mf (ξ) = λ(ξ)m0(ξ + π), (5.17)

eλ(ξ) + λ(ξ + π) = 0.

Tale equazione ha come uniche soluzioni funzioni del tipo

λ(ξ) = eiξν(2ξ), (5.18)

per ν funzione 2π-periodica. Sostituendo (5.18) e (5.17) in (5.15) si ha

f(ξ) = eiξ/2m0(ξ/2 + π)ν(ξ)ϕ(ξ/2).


In base a quest’ultimo lemma, risulta naturale supporre che la ψ che cer-chiamo sia del tipo

ψ(ξ) = eiξ/2m0(ξ/2 + π)ϕ(ξ/2).

Infatti, non preoccupandoci in dettaglio delle questioni di convergenza, dalprecedente lemma si deduce anche che

f(ξ) = ν(ξ)ψ(ξ) =∑

k

νke−ikξψ(ξ),

da cuif(t) =

∑

k

νkψ(t− k).

Cioe le traslazioni di ψ riproducono ogni funzione di W0 e dunque le funzioniψj,k riproducono ogni funzione di Wj .

Per esercizio si lascia da verificare che ψ(t− k) : kZ sono davvero ortonor-mali.

Sulla base dei precedenti risultati si ha dunque il seguente fondamentalerisultato della teoria delle ondine.

Teorema 5.3.3. Sia (Vj)j∈Z una analisi multirisolutiva con le proprieta sopradescritte. Allora esiste una base ortonormale associata ψj,k : j, k ∈ Z perL2(R), tale che

ψ(ξ) = eiξ/2m0(ξ/2 + π)ϕ(ξ/2),

o equivalentemente

ψ =∑

n

(−1)n−1h−n−1︸︷︷︸:=gn

ϕ−1,n. (5.19)

5.4 Ondine a supporto compatto

Il problema di costruire ondine si riduce dunque alla identificazione di opportunefunzioni scalanti ϕ che risolvono (5.9) e che generano una analisi multirisolutiva.In realta si puo dimostrare che se le ϕ(·−k) sono ortonormali, ϕ risolve (5.9), e ϕha media non nulla, e.g.,

∫ϕ(t)dt 6= 0, allora essa certamente genera una analisi

multirisolutiva. Purtroppo, anche la sola risoluzione della equazione (5.9), che euna equazione di punto fisso, non e finitamente computabile se il supporto dellasuccessione (hn)n non e finito e se la soluzione ϕ non e a supporto compatto.Viceversa, si avrebbe

ϕ(k) =∑

ℓ

h2k−ℓϕ(ℓ), k ∈ Z.

Posto il vettore Φ = (ϕ(k))Tk e la matrice H = (h2k−ℓ)k,ℓ si ha

Φ = HΦ,

5.4. ONDINE A SUPPORTO COMPATTO 75

200 400 600 800 1000

-0.02

0.02

0.04

0.06

0.08

Figura 5.1: Esempio di funzione scalante a supporto compatto

cioe Φ e un semplice autovettore della matrice H . Dunque risolvendo questoproblema agli autovalori si ottengono i valori di ϕ sugli interi. Dalla relazionescalante

ϕ(ℓ

2) =

∑

n

hnϕ(ℓ− n),

si ottengono ricorsivamente i valori di ϕ su un qualsiasi reticolo diadico di nodi2−jℓ. Inoltre la stessa computabilita dell’ondina

ψ =∑

n

(−1)n−1h−n−1ϕ−1,n,

dipende dal fatto che (hn)n abbia supporto finito.

Una classificazione molto ampia di funzioni scalanti a supporto compattoe stata scoperta (costruita?) da Daubechies in un suo celebratissimo articolodel 1988, stabilendo la nascita dell’analisi numerica mediante wavelet, con ri-cadute in moltissimi ambiti applicativi, dal trattamento numerico di segnali edimmagini, alla simulazione numerica, e recentemente alla risoluzione di proble-mi inversi e di discontinuita libera. Non entriamo nei dettagli della costruzionedelle wavelet a supporto compatto e indichiamo come approfondimento il [19,Capiolo 6].


200 400 600 800 1000

-0.05

0.05

0.1

Figura 5.2: Esempio di ondina a supporto compatto

5.5 La trasformata wavelet discreta (DWT)

Data una funzione f j−1 ∈ Vj−1 si ha

f j−1 =∑

k

cj−1k ϕj−1,k

=∑

k

cjkϕj,k

︸︷︷︸Vj

+∑

k

djkψj,k

︸︷︷︸Wj

=∑

k

cj+1k ϕj+1,k

︸︷︷︸Vj+1

+∑

k

dj+1k ψj+1,k

︸︷︷︸Wj+1

+∑

k

djkψj,k

︸︷︷︸Wj

= . . .

Vogliamo mostrare ora come si trasformano i coefficienti cj−1k nella coppia coeffi-

cienti (cjk, djk) e viceversa. Queste trasformazioni sono dette trasformata wavelet

discreta e trasformata wavelet discreta inversa. Permettono dunque di associaread una funzione una digitalizzazione stabile (anche se non ridondante in questocaso). La trasformata diretta e dunque data da

cjk = 〈f j−1, ϕj,k〉 =∑

n

hn−2k〈f j−1, ϕj−1,k〉 =∑

n

hn−2kckj−1

djk = 〈f j−1, ϕj,k〉 =∑

n

gn−2k〈f j−1, ϕj−1,k〉 =∑

n

gn−2kckj−1,

5.5. LA TRASFORMATA WAVELET DISCRETA (DWT) 77

e quella inversa da

cj−1k = 〈f j−1, ϕj−1,k〉

=∑

k

cjk〈ϕj,k, ϕj−1,k〉 +∑

k

djk〈ψj,k, φj−1,k〉

=∑

k

[hn−2kc

jk + gn−2kd

jk

].

Si puo dimostrare che la trasformata wavelet discreta ha complessita O(n) oven = 2−j+1, j < 0. Quindi tale trasformata e persino piu rapida della trasfor-mata di Fourier discreta. Nella pratica numerica, specialmente nel trattamentonumerico dei segnali, dato un segnale f j−1 si assume che cj−1

k ≈ f j−1(2−jk)siano una buona approssimazione del campionamento del segnale su un reticolodiadico, e si procede dunque alla trasformata wavelet veloce per calcolare tuttii coefficienti a scale inferiori. In modo simile a quanto fatto per le frame di Ga-bor discrete, si possono utilizare le espansioni wavelet discrete per comprimeresegnali ed analizzarne proprieta di interesse, come per esempio la regolarita.

Capitolo 6

Problemi inversi con vincolidi sparsita

La lezione che segue presenta alcuni metodi recenti per la risoluzione di pro-blemi inversi con vincoli di sparsita. Un problema inverso e determinato daun modello, matematicamente dato da un operatore lineare o non lineare, cheapplicato ad una soluzione di interesse fornisce una misura. Tipicamente ta-le operatore e non invertibile e dunque la misura data e incompleta per poterfornire una ricostruzione della soluzione cercata. In realta le soluzioni posso-no essere infinite in questo caso. Per circoscrivere e determinare la soluzionedi interesse del problema inverso occorre aggiungere ulteriori assunzioni sullasoluzione che non sono date dalle misure. In molti problemi provenienti dallafisica, e naturale assumere la soluzione come la soluzione di energia minima, cheviene misurata come la norma di un opportuno spazio di Hilbert. Dato che talevincolo aggiuntivo e strettamente convesso, esso tipicamente impone l’unicitadella soluzione. Associati a tale formulazione del problema esistono moltissimialgoritmi che consentono di calcolare la soluzione voluta in modo efficiente. Puressendo tale formulazione di grande interesse applicativo, in questa lezione cioccuperemo di vincoli che promuovono la sparsita della soluzione rispetto aduna frame, piuttosto che l’energia minima. Per soluzione sparsa si intende unasoluzione che, rappresentata rispetto ad una determinata frame, si scriva comecombinazione lineare di pochissimi elementi della frame. Tale assunzione, rispet-to a quella di minimia energia, risulta piu realistica per moltissime applicazioniper cui la soluzione di interesse non ha energia minima. Per esempio, comediscuteremo, le immagini naturali non hanno normalmente minima norma ℓ2.Presenteremo dunque metodi ed algoritmi per il calcolo numerico di soluzioni diproblemi inversi con vincoli di sparsita ed alcune applicazioni nel trattamentonumerico di immagini digitali. Indichiamo come approfondimenti [20, 21, 1, 3,5, 6, 7, 8, 23, 26, 27, 28, 29, 30, 31, 34, 50, 51, 53, 66, 74, 76, 77, 78, 80, 81].

79

80 CAPITOLO 6. PROBLEMI INVERSI CON VINCOLI DI SPARSITA

6.1 Problemi inversi lineari

In molti problemi pratici la quantita di interesse, essa sia una proprieta fisica oun oggetto con caratteristiche geometriche e morfologiche come una immagine,e misurata mediante quantita direttamente osservabili, che sono normalmenterelazionate alla quantita di interesse per una dipendenza lineare. In effettimoltissimi problemi sono in realta non lineari, ma e sorprendente scoprire che laapprossimazione lineare funziona molto bene in pratica. Percio, sebbene ancheil caso non lineare sia trattabile, nella nostra lezione ci limiteremo a considerareproblemi inversi lineari, che possono esssere matematicamente descritti comesegue. Siano K e H due spazi di Hilbert separabili, e A : K → H un operatorelineare e continuo (limitato). Si assume di avere un osservabile g ∈ H,

g = Af.

Il nostro scopo consiste nella ricostruzione dell’incognita f ∈ K. Come descrittonell’introduzione a questa lezione, siamo interessati alla situazione in cui la cor-rispondente mappa dal vettore f al vettore g non e invertibile o la cui inversasia fortemente mal condizionata. Questo e il caso, per esempio, per operatoriA compatti con kernel non nullo, come gli operatori di convoluzione in L2(Rd).In tal caso, soprattutto quando il dato g e determinato in modo “incerto”, acausa di misure imprecise corrotte da cosiddetto rumore, la soluzione f non edeterminabile a meno di aggiungere ulteriori assunzioni sulla soluzione.

Iniziamo la nostra discussione con un esempio concreto. Si consideri unaframe G = gn : n ∈ N per lo spazio di Hilbert H. Allora si consideri la mappadi sintesi (o di ricostruzione) Rg : ℓ2(N) → H data da ~c 7→ ∑

n cngn. Gia sap-piamo che tale mappa ha kernel non nullo non appena la frame e ridondante.Esistono dunque coefficienti ~d e ~c tali che Rg(~d) = Rg(~c).

E dunque lecito chiedersi come si possano calcolare dei “buoni” coefficientiper una data g ∈ H. Come abbiamo detto una possibile scelta e data daicoefficienti di norma minima, cioe:

~c(g) := argmin‖~d‖ℓ2 : g =∑

n

dngn. (6.1)

In tal caso la soluzione del problema (6.1) e unica. Dunque, tra tutte le infinitesoluzioni possibili, ci siamo ristretti a quella di norma minima. Vediamo comesi possano caratterizzare tali coefficienti.

Proposizione 6.1.1. Data una frame gnn∈N in H e f ∈ H, sia an =〈f, S−1gn〉, e f =

∑n∈N

angn. Se e possibile trovare altri scalari cn t.c. f =∑n∈N

cngn allora vale che∑n∈N

| cn |2=∑n∈N| an |2 +

∑n∈N

| an − cn |2.

Dimostrazione. Notiamo che 〈gn, S−1f〉 = 〈S−1gn, f〉 = an. Sostituendo f =∑n∈N

angn e f =∑

n∈Ncngn nei primi termini del prodotto interno ottenia-

mo∑

n∈N| an |2= ∑

n∈Nanan =

∑n∈N

an〈f, S−1gn〉 = 〈f,∑n∈NanS

−1gn〉 =

6.1. PROBLEMI INVERSI LINEARI 81

〈f,∑n∈N〈f, S−1gn〉S−1gn〉 = 〈f, S−1

∑n∈N

〈f, S−1gn〉gn〉 = 〈f, S−1∑

n∈Nangn〉 =

〈f, S−1f〉 = 〈∑n∈Ncngn, S

−1f〉 =∑

n∈Ncn〈gn, S−1f〉 =

∑n∈N

cnan. Quindi,

∑

n∈N

| an |2 +∑

n∈N

| an − cn |2 =∑

n∈N

| an |2 +∑

n∈N

(| an |2 −ancn − ancn+ | cn |2)

=∑

n∈N

| cn |2 .

dato che la somma dei moduli degli an al quadrato e uguale al suo coniugato.

In altre parole, i coefficienti duali canonici hanno norma minima tra tutti ipossibili coefficienti.

Inoltre mostriamo concretamente nel seguente esempio che dati i coefficienticn = 〈f ,gn〉, i vettori gn che ricostruiscono f non sono univocamente determi-nati.

Esempio 6.1.2. Sia H = C2, f = (−1, 3) e gT1 = (1,−1), gT2 = (0, 1), gT3 =(1, 1). I coefficienti (cn)n∈1,2,3 sono

(cn)n∈1,2,3 =(〈f ,g1〉, 〈f ,g2〉, 〈f ,g3〉

)= (−4, 3, 2).

La frame duale canonica gnn∈1,2,3 e data da

gnn∈1,2,3 = S−1gnn∈1,2,3 =

(1/2−1/3

),

(0

1/3

),

(1/21/3

).

Ottengo cosı una ricostruzione di f :

f =

3∑

n=1

cngn =(( −2

43

)+

(01

)+

(123

))=

(−13

)

Prendendo un’altra duale definita da

hnn∈1,2,3 =

(140

),

(01

),

(00

)

abbiamo ancora una ricostruzione di f.

f =

3∑

n=1

cnhn =(( −1

0

)+

(03

)+

(00

))=

(−13

)

Osserviamo che nel caso ricostruissimo il segnale f con la duale hnn∈1,2,3abbiamo bisogno di un coefficiente in meno. Quindi, pur producendo coefficientidi norma minima, i.e.,

min∑

n

|cn|2, soggetto a f =∑

n

cngn,


il duale canonico non ha tutte le virtu possibili. In particolare, il duale canoniconon produce rappresentazioni sparse (cioe quella con il numero minimo di coef-ficienti non nulli) di un vettore. E invece possibile trovare la rappresentazionepiu sparsa risolvendo un altro problema di minimio, dato da:

min∑

n

|cn|, soggetto a f =∑

n

cngn. (6.2)

Tale problema si puo facilmente riformulare come un problema di programma-zione lineare [?]. Tanto piu e ridondate una frame tanto piu sparsa e la soluzionedi quest’ultimo problema.

Per comprendere in modo semplice ed intuitivo il perche la minimizzazionedella norma ℓ1 comporti la promozione di una soluzione sparsa proponiamo lasituazione come raffigurata in Figura 6.10. In tal caso si mostra la soluzione diun problema lineare del tipo Af = g in R2 la cui norma ℓ1 e minima. Come sivede da immediate considerazioni geometriche, a meno che il kernel di A noncontenga vettori del tipo (1, 1) o (1,−1), tale soluzione e su un vertice dellapalla unitaria in norma di ℓ1, che ha necessariamente almeno una componentenulla. Dunque e una soluzione piu sparsa. Tale non sarebbe minimizzando lanorma di ℓ2.

-3 -2 -1 0 1 2 3-1

-0.5

0

0.5

1

1.5

2

2.5

Figura 6.1: La soluzione piu spars del problema Af = g.

6.2 Problemi inversi con vincoli di sparsita

In molti problemi concreti le misure g date da Af = g sono imprecise, cioeg = g0 + ν ove g0 e la misura ideale non affetta da alcuna distorsione e ν e ladistorsione (rumore). Quindi la ricostruzione esatta di una f tale che Af = gnon rappresenta probabilmente la misura cercata e non possiamo voler essere piuprecisi del rumore σ2 = ‖ν‖2. Ci accontenteremo dunque di voler ricostruire unaf tale che ‖Af − g‖2

H = σ2. Come sopra citato assumeremo che la soluzione siasparsa rispetto ad una data frame. Come abbiamo visto nei capitoli precedenti,esistono numerosi tipi di frame, per esempio le frame di Gabor e le frame diondine, che producono diverse analisi. Allora modifichiamo il problema iniziale

6.3. ALGORITMI ITERATIVI DI SOGLIA 83

componendo il modello A con la mappa di ricostruzione Rg ottenendo un nuovooperatore T = A Rg tale che

σ2 = ‖Af − g‖2H = ‖T~u− g‖2

H,

ove f =∑

n ungn. Con tale modifica, trasportiamo il problema in ℓ2(N), chee lo spazio dei coefficienti frame. Si puo dimostrare mediante l’uso di unageneralizzazione del teorema dei moltiplicatori di Lagrange che la soluzione delproblema

argmin∑

n

|un|, soggetto a ‖T~u− g‖2H = σ2, (6.3)

equivale alla risoluzione del seguente problema

arg min~u∈ℓ2(N)

‖T~u− g‖2H + 2τ

∑

n

|un| (6.4)

per una opportuno moltiplicatore τ > 0 che dipende da σ2. Nelle seguenti sezionici proponiamo di studiare algoritmi efficienti per la soluzione di problemi deltipo (6.4).

6.3 Algoritmi iterativi di soglia

La minimizzazione di un funzionale (regolare) avviene mediante lo studio delleequazioni di Eulero-Lagrange associate. Dato un funzionale convesso J : K →R che si desidera minimizzare, si calcolano i suoi minimi come soluzioni delleequazioni

∇J (u) = 0, u ∈ K.L’operatore ∇ rappresenta una derivata generalizzata, che, se esiste, e definitada

〈∇J (u), v〉 = limh→0

J (u+ hv) − J (u)

h, ∀v ∈ K.

Nel nostro caso J (~u) := ‖T~u− g‖2H + 2τ

∑n |un| e si ottiene

〈∇J (~u), en〉 = limh→0

J (u+ hen) − J (u)

h= [2T ∗(T~u− g)]n + sign(un),

per ogni n ∈ N tale che un 6= 0 e en = δn e un elemento della base canonica diℓ2(N). Il problema e che tale derivata fallisce di essere definita per quegli n taliche un = 0. Inoltre anche la diretta soluzione delle equazioni

[2T ∗(T~u− g)]n + sign(un) = 0,

risulta complicata dal fatto che l’operatore T ∗T accoppia tutte le entrate di~u rendendo le equazioni un sistema non lineare accoppiato. Per tali ragioni enecessario escogitare un approccio diverso.


Un importante strumento per l’analisi di funzionali non derivabili e dei lorominimi e il concetto di sottodifferenzale. Si ricorda che per un funzionale con-vesso F su uno spazio di Banach V il suo sottodifferenziale ∂F (x) in un puntox ∈ V con F (x) <∞ e definito come l’insieme

∂F (x) = x∗ ∈ V ′, x∗(z − x) + F (x) ≤ F (z) per ogni z ∈ V ,ove V ′ e lo spazio duale topologico di V . Da questa definizione e ovvio che0 ∈ ∂F (x) se e solo se x e un punto di minimo per F .

Esempio 6.3.1. Sia V = ℓ1(N) e F (x) := ‖x‖1 e la la norma in ℓ1. Si ha

∂‖ · ‖1(x) = ξ ∈ ℓ∞(N) : ξλ ∈ ∂| · |(xλ), λ ∈ N (6.5)

ove ∂| · |(z) = sign(z) se z 6= 0 e ∂| · |(0) = [−1, 1]. La dimostrazione di talecaratterizzazione del suttodifferenziale e lasciata per esercizio.

Per semplicita di notazione nel seguito dimenticheremo ~· nella scrittura divettori in ℓ2, cioe porremo con abuso u = ~u. Ci risultera utile introdurre ilseguente funzionale, detto il surrogato di J ,

J S(u, a) := ‖g − Tu‖2H + 2τ‖u‖1︸︷︷︸

:=J (u)

+‖u− a‖22 − ‖Tu− Ta‖2

H. (6.6)

Un calcolo diretto (per esercizio!) da

J S(u, a) = ‖(a+T ∗(g−Ta))−u‖22+‖u‖1−‖a+T ∗(g−Ta)‖2

2+‖g‖2H−‖Ta‖2

H+‖a‖22.

Nel seguito si assumera ‖T ‖ < 1. Questa condizione e sempre ottenibile me-diante un opportuna scalatura di T e g. Si osservi che

‖u− a‖22 − ‖Tu− Ta‖2

H ≥ C‖u− a‖22, (6.7)

per C = (1 − ‖T ‖2) > 0. Quindi, si ha

J (u) = J S(u, u) ≤ J S(u, a), (6.8)

eJ S(u, a) − J S(u, u) ≥ C‖u− a‖2

2. (6.9)

In particolare, J S e un funzionale strettamente convesso rispetto ad u e dunqueha un unico minimizzatore, non appena sia a fissato. Osservando che ∂(‖T ·−g‖2

H)(u) = 2T ∗(Tu − g) (see [53, Lemma 3.2]) e per applicazione di [?,Proposition 5.2] combinati con l’esempio appena sopra riportato, si ottiene lasequente caratterizzazione dei sottodifferenziali di J e di J S .

Lemma 6.3.2. i) Il sottodifferenziale di J at u e dato da

∂J (u) = 2T ∗(Tu− g) + 2τ∂‖ · ‖1(u)

= ξ ∈ ℓ∞(N) : ξn ∈ [2T ∗(Tu− g)]n + 2τ∂| · |(un).ii) Il sottodifferenziale di J S rispetto alla sola variabile u e dato da

∂uJ S(u, a) = −2(a+ T ∗(g − Ta)) + 2u+ 2τ∂‖ · ‖1(u)

= ξ ∈ ℓ∞(N) : ξn ∈ [−2(a+ T ∗(g − Ta))]n + 2un + 2τ∂| · |(un).

6.4. ALGORITMO DEL GRADIENTE PROIETTATO 85

Posto

Sτ (x) =

x− sign(x)τ, |x| > τ0, altrimenti.

si ha Sτ (z) := (Sτ (zn))n∈N che e l’unica soluzione della inclusione differenziale0 ∈ 2(u− z)+ 2τ∂‖ · ‖1(u), si veda per esempio [53] e [20, Proposition 2.1]. DalLemma 6.3.2 ii) si ottiene immediatamente

arg minu∈ℓ2(N)

J S(u, a) = Sτ (a+ T ∗(g − Ta)).

L’operatore non lineare S e detto operatore di soglia, poiche azzera tutte leentrate di un vettore che non superano in valore assoluto una determinata soglia.Con questi ingredienti si formula il seguente algoritmo iterativo. Per un inizialeu(0) ∈ ℓ2(N) (u(0) = 0 va bene) si itera

u(n+1) = Sτ (u(n) + T ∗(g − Tu(n))), n ≥ 0. (6.10)

Per (6.8) si ha

J (u(n)) = J S(u(n), u(n))

≥ J S(u(n+1), u(n))

≥ J S(u(n+1), u(n+1)) = J (u(n+1)).

Quindi si ha che (J (u(n)))n e una successione decrescente e, in quanto limitatada 0 dal basso, essa e anche necessariamente convergente. Per (6.9) si ha anche

J (u(n)) − J (u(n+1)) ≥ ‖u(n+1) − u(n)‖22,

e dunque si ha ‖u(n+1)−u(n)‖22 → 0. Senza entrare nei dettagli di una dimostra-

zione rigorosa, i conti che abbiamo appena svolto suggeriscono che l’algoritmotende a minimizzare il funzionale J e che la successione (u(n))n in effetti con-verge ad un minimizzatore u∗ per J . La dimostrazione completa di questofatto si puo trovare in []. L’algoritmo (6.10) ha una duplice azione. Il passou(n) 7→ u(n) +T ∗(g−Tu(n)) tende a minimizzare ‖Tu− g‖2

H nella direzione delsuo gradiente (questo passo e detta iterazione di Landweber), mentre il passou(n) + T ∗(g − Tu(n)) 7→ Sτ (u

(n) + T ∗(g − Tu(n))) tende a promuovere solo lecomponenti importanti dell’attuale iterazione, cioe a minimizzare ‖u(n+1)‖1 e amassimizzare la sparsita.

6.4 Algoritmo del gradiente proiettato

Il problema combinatorio di trovare la soluzione piu sparsa per ‖Tu − g‖2 =σ2 e dunque risolvibile mediante la risoluzione di un problema convesso cheimplementa la minimizzazione della norma in ℓ1. Questa e la ragione per cui inalgoritmo del tipo (6.10) si dice veloce. E veloce rispetto alla ricerca esausitivadi forza-bruta per la soluzione piu sparsa, che e un problema cosiddetto NP (dicomplessita non polinomiale).


0 5 10 15 20 25

104

105

log(χ2) vs. ||w||1, 2000 iterations

Path (one choice of τ)Trade−off curve (pre−computed, 300000 it.)

Figura 6.2:

Comunque, per essere piu onesti, occorrerebbe confrontare la velocita di con-vergenza di (6.10) con quella degli usuali metodi iterativi per la risoluzione disitemi lineare, come i metodi di gradiente e gradiente coniugato. In tal casosi trova che in pratica, sebbene esso sia semplicissimo da implementare e dun-que molto attraente per un utilizzatore, l’algoritmo non e affatto competitivo.Allo scopo di presentare una prima possibile accelerazione di tale algoritmo,vogliamo ora descriverne il comportamento, vedi Fig. 6.2: infatti, mentre l’al-goritmo inizialmente converge in modo relativamente veloce, in seguito eccedeil valore limite per ‖u(n)‖1 e richiede molto tempo ed iterazioni per correggerel’approssimazione. Dobbiamo immaginare che, partendo da u(0) = 0, l’algo-ritmo generi un cammino u(n)n∈N che e inizialmente confinato nella pallaBR := u ∈ ℓ2(N) : ‖u‖1 ≤ R, with R := ‖u∗(τ)‖1. In seguito esce dalla pallaper poi ritornarvi (tipicamente su un vertice o su uno spigolo, corrispondenti avettori con poche componenti non nulle) solo al limite. Quindi il primo modointuitivo per evitare questa lunga e non necessaria deviazione esterna e di man-tenere tutte le iterazione sulla palla BR. Un metodo per forzare le iterazionisulla palla e di sostituire l’operatore di soglia con la proiezione ortogonale in ℓ2

sulla palla, data da PR := PBR , che e definita come l’unico punto sulla palla condistanza minima. Si ricordi che tale punto esiste sempre ed e unico perche BRe un insieme chiuso e convesso. Otteniamo dunque il seguente algoritmo:

u(n+1) = PR

[u(n) + T ∗g − T ∗Tu(n)

]. (6.11)

Chiameremo l’iterazione (6.11) l’iterazione di Landweber proiettata.In Fig. 6.3 mostriamo la tipica dinamica dell’algoritmo. Notiamo che infatti

esso non eccede piu R e allo stesso tempo raggiunge una buona approssimazione


0 5 10 15 20 25

104

105

log(χ2) vs. ||w||1, 200 iterations, LW

final χ2= 2679.5832

Path (one choice of target L1 norm)Trade−off curve (pre−computed, 300000 it.)Trade−off curve (pre−computed, 2000 it.)Trade−off curve (pre−computed, 1000 it.)

Figura 6.3:

di u∗(τ) con un minor numero di iterazioni rispetto a (6.10). Comunque, anchequesta modificazione non e ancora ottimale e il sentiero u(n)n∈N puo correrelungo facce e spigoli della palla con passi relativamente piccoli. Allo scopo di farepassi piu lunghi, e potenzialmente saltare da una faccia all’altra attraverso unascorciatoia interna alla palla, potremmo modificare le iterazioni introducendoun parametro di discesa β(n) > 0 come segue:

u(n+1) = PR

[u(n) + β(n)T ∗(g − Tu(n))

]. (6.12)

Ci riferiamo a tale algoritmo come alla iterazione di gradiente proiettata pro-jected. In Fig. 6.4 mostriamo la corrispondente dinamica dove possiamo notareuna significativa accelerazione. La scelta di β(n) e una questione cruciale. Unrisultato recente [] stabilisce la convergenza forte delle (u(n))n∈N prodotte da(6.12) a un minimizzatore di J (per un qualche τ che dipende implicitamenteda R), quando β(n) e determinato come segue:si definisca ΓR(β;x) := PR(x+ βT ∗(g − Tx)),

β(n) := supβ > 0 : ∀ γ ∈ [0, β), γ‖T (ΓR(γ, u(n))−u(n))‖2H ≤ r‖ΓR(γ;u(n))−u(n)‖2,

(6.13)dove r > 0 e un limite superiore per ‖T ∗T ‖ℓ2→ℓ2 , r < 1. Segue che β(n) ≥ 1 perogni n e

β(n)‖T (u(n+1) − u(n))‖2H ≤ r‖u(n+1) − u(n)‖2. (6.14)

Per ragioni tecniche si definisce infine

β(n) := minβ, β(n), (6.15)


0 5 10 15 20 25

104

105

log(χ2) vs. ||w||1, 200 iterations, SD

final χ2= 2066.9603

Path (one choice of target L1 norm)Trade−off curve (pre−computed, 300000 it.)Trade−off curve (pre−computed, 2000 it.)Trade−off curve (pre−computed, 1000 it.)

Figura 6.4:

per un valore β > 1 sufficientemente largo.

Pur non volendo entrare nei dettagli della dimostrazione delle convergenzadel precedente algoritmo, vogliamo dare una descrizione piu precisa di comesi implementino le proiezioni su palle BR e alcune indicazioni intuitive sullaconvergenza dell’algoritmo.

6.4.1 Proiezioni su BR e operatori di soglia

In questa sezione vogliamo discutere alcune proprieta delle proiezioni su palleBR. In particolare, studiamo le loro relazioni con gli operatori di soglia e il lorocalcolo esplicito. Diamo inoltre una stima della complessita computazionaledi tali proiezioni in spazi di dimensione finita, essendo un aspetto cruciale pervalutare l’efficienza dell’algoritmo (6.12).

Innanzi tutto osserviamo una utile proprieta dell’operatore di soglia.

Lemma 6.4.1. Per ogni a ∈ ℓ2(N) e τ > 0, ‖Sτ (a)‖1 e una funzione di τ linearea tratti, continua e decrescente. Inoltre ‖Sτ=0(a)‖1 = ‖a‖1 e ‖Sτ≥maxi |ai|(x)‖1 =0.

Dimostrazione. ‖Sτ (a)‖1 =∑k |Sτ (ak)| =

∑k Sτ (|ak|) =

∑|ak|>τ (|ak|−τ), ove

la somma a destra e finita per τ > 0.

Una illustrazione schematica di questa proprieta di legge in Fig. 6.5.

Il seguente lemma mostra come determinare la proiezione di un vettore asulla palla BR come applicazione di un operatore si soglia.


τ

‖Sτ (a)‖1

‖a‖1

τ

R

maxi |xi|

Figura 6.5: Per a ∈ ℓ2, la funzione τ 7→ ‖Sτ (a)‖1

Lemma 6.4.2. La proiezione di a su BR e data da PR(a) = Sτ (a) ove τ(dipendente da a e R) e scelto tale che ‖Sτ (a)‖1 = R se ‖a‖1 > R. Se ‖a‖1 ≤ Rallora PR(a) = S0(a) = a.

Dimostrazione. Si assuma ‖a‖1 > R. Si scelga τ > 0 tale che ‖Sτ (a)‖ = R(questo e possibile grazie al Lemma 6.4.1). Abbiamo gia menzionato nell’in-troduzione che b = Sτ (a) e l’unico minimizzatore di ‖u − a‖2

2 + 2τ‖u‖1 equindi:

‖b− a‖22 + 2τ‖b‖1 < ‖u− a‖2

2 + 2τ‖u‖1

per ogni u 6= b; in particolare dal momento che ‖b‖1 = R e se ‖u‖1 ≤ R, segueche

‖b− a‖2 < ‖u− a‖2 ∀u ∈ BR

Quindi PR(a) = b = Sτ (a)

In pratica e piu esplicitamente, la proiezione su BR si calcola in dimensionefinita m come segue: si ordinino le entrate di a per valore assoluto |ai1 | ≥ |ai2 | ≥. . . ≥ |aim |. Sia n ∈ 1, . . . ,m il piu grande naturale tale che

|ain | ≥1

n− 1

(n−1∑

k=1

|aik | −R

).

Allora si ha

(PR(a))ij = aij −sign(aij )

n

(n∑

k=1

|aik | − R

), j = 1, . . . , n,

(PR(a))ij = 0, j = n+ 1, . . . ,m.


Queste formule sono ricavate anche in [53, Lemma 4.1 and Lemma 4.2], osser-vando che PR(a) = a− S∞

R (a), dove

S∞R (a) = arg min

x∈Rm(‖a− x‖2

2 + 2R‖x‖∞), x ∈ Rm. (6.16)

Quest’ultimo e ancora un operatore di soglia, ma e riferito alla norma ℓ∞. E im-mediato verificare che la complessita di calcolo di una proiezione e O(m logm).

6.4.2 Indicazioni di convergenza

In modo simile a quanto fatto nel caso di (6.10) e del funzionale surrogato, siha il seguente risultato.

Proposizione 6.4.3. Si assuma ‖T ‖ℓ2→H < 1 e β ≥ 1. Il vettore ΓR(β;u) =PR(u+ βT ∗(y − Tu)) e l’unico minimizzatore del funzionale

Fβ(w;u) = ‖Tw − y‖2H − ‖T (w − u)‖2

H +1

β‖w − u‖2, on BR. (6.17)

Inoltre la successione (‖Tu(n) − g‖2)n∈N e decrescente, e

limn→∞

‖u(n+1) − u(n)‖ = 0. (6.18)

Dimostrazione. Prima di tutto si osservi che da ‖T ‖ < 1 e β ≥ 1, il funzionaleFβ e strettamente convesso, e dunque ha un unico minimizzatore in BR. Sia uil suo minimizzatore, allora per ogni w ∈ BR e t ∈ [0, 1] si ha

Fβ(u;u) ≤ Fβ((1 − t)u+ tw;u)

⇒ 2t

[〈T u− g, T (w − u)〉 − 〈T u− Tu, T (w− u)〉 +

1

β〈u− u,w − u〉

]+t2

β‖w − u‖2 ≥ 0

⇒ [β〈Tu− g, T (w − u)〉 + 〈u − u,w − u〉] + t

2‖w − u‖2 ≥ 0

⇒ 〈w − u, u− u+ βT ∗(Tu− g)〉 ≥ 0

⇒ 〈u+ βT ∗(g − Tu)− u, w − u〉 ≤ 0.

L’ultima implicazione equivale a dire che u = ΓR(β;u).

Per definizione u(n+1) in (6.12) e per γ = 1r − 1 si ha

‖Tu(n+1) − g‖2H ≤ ‖Tu(n+1) − g‖2

H + γ‖T (u(n+1) − u(n))‖2H

≤ ‖Tu(n+1) − g‖2H + (1 + γ)‖T (u(n+1) − u(n))‖2

H − ‖T (u(n+1) − u(n))‖2H

≤ ‖Tu(n+1) − g‖2H − ‖T (u(n+1) − u(n))‖2

H +1

β(n)‖u(n+1) − u(n)‖2

= Fβ(n)(u(n+1);u(n)) ≤ Fβ(n)(u(n);u(n)) = ‖Tu(n) − g‖2H.

6.5. DUE APPLICAZIONI RILEVANTI 91

Inoltre si ha

−Fβ(n+1)(u(n+1);u(n+1)) + Fβ(n)(u(n+1);u(n)) =1

β(n)‖u(n+1) − u(n)‖2 − ‖T (u(n+1) − u(n))‖2

H

≥ 1 − r

β(n)‖u(n+1) − u(n)‖2 ≥ 1 − r

β‖u(n+1) − u(n)‖2.

Cio implica

N∑

n=0

‖u(n+1) − u(n)‖2 ≤ β

1 − r

N∑

n=0

(Fβ(n)(u(n+1);u(n)) − Fβ(n+1)(u(n+1);u(n+1))

)

≤ β

1 − r

N∑

n=0

(Fβ(n)(u(n);u(n)) − Fβ(n+1)(u(n+1);u(n+1))

)

=β

1 − r

(Fβ(0)(u(0);u(0)) − Fβ(N+1)(u(N+1);u(N+1))

)

≤ β

1 − rFβ(0)(u(0);u(0)).

Quindi la serie∑∞

n=0 ‖u(n+1) − u(n)‖2 converge e limn→∞ ‖u(n+1) − u(n)‖ =0.

Il precente risultato suggerisce che l’algoritmo tende a minimizzare ‖Tu−g‖2H

in BR e allo stesso tempo che esso converge. La dimostrazione completa cheesso converga infatti ad un minimizzatore u∗ di ‖Tu− g‖2 in BR si puo trovarein [21]. Si osservi che, per quanto precentendemente osservato,

u∗ = arg minu∈BR

‖Tu− g‖2H = arg min

u∈ℓ2(N)‖Tu− g‖2

H + 2τ‖u‖1,

per un opportuno τ > 0 e dunque che u∗ minimizza anche J .

6.5 Due applicazioni rilevanti

Concludiamo la lezione e questo corso, con due applicazioni. La prima si ri-ferisce al problema delle ricolorazione automatica degli affreschi di Mantegnaframmentati, di cui ci siamo occupati nel Capitolo 4. La seconda e relativa allaricostruzione delle densita di corrente all’interno del cervello durante la sua at-tivita a partire dal dato del campo magnetico generato che viene esternamentemisurato, tale tecnica e detta magnetoencefalografia. Per entrambe queste ap-plicazioni e cruciale poter esprimere la soluzione dei problemi inversi associaticome espansione sparsa con rispetto ad una base/frame di ondine. Ecco dunquecome tutta la teoria delle espansioni frame e dei problemi inversi si incontrano.


6.5.1 Ricolorazione di immagini e restauro d’immagini

Nel Capitolo 4 abbiamo proposto un algoritmo per la localizzazione dei fram-menti degli affreschi del Mantegna frammentati da un bombardamento. L’algo-ritmo ha consentito il restauro parziale di questi affreschi. Il risultato, sebbenedi importanza e storica e culturale, non permette una lettura completa dell’ope-ra. In questa sezione vogliamo mostrare come gli strumenti sin qui sviluppati,i.e., frame di ondine e metodi ed algoritmi per la soluzione di problemi inversi,possano essere usati per stimare in modo accurato i colori degli affreschi nellelacune, a partire dal dato vero dei colori dei frammenti localizzati e i livelli digrigio delle lacune, dati dalle foto Alinari degli affreschi del 1920.

50 100 150 200 250

50

100

150

200

250

Figura 6.6: Stima della curva non lineare L da una distribuzione di punti concoordinate date da combinazioni lineari αr + βg + γb dei colori (r, g, b) deiframmenti e dai corrispondenti livelli di grigio sottostanti forniti dalle foto del1920.

Una immagine digitale a colori puo essere modellata come una funzioneu : Ω ⊂ R2 → R3

+, tale che ad ogni “punto” x dell’immagine si associa il vettoreu(x) = (r(x), g(x), b(x)) ∈ R3

+ dei colori rappresentati per differenti canalidi rosso, verde e blu. In particolare, una digitalizzazione di una immagine ucorrisponde ad un suo campionamento su un reticolo τZ2, τ > 0. Scriveremoancora u : N → R3

+, u(x) = (r(x), g(x), b(x)), per x ∈ N := Ω ∩ τZ2.

I livelli di grigio di una immagine possono essere descritti come una sotto-varieta M ⊂ R3 da

M := Mσ = σ(x) : x = L(r, g, b) := L(αr + βg + γb), (r, g, b) ∈ R3+,

ove α, β, γ > 0, α + β + γ = 1, L : R → R e una funzione non decrescente eσ : R+ → R3

+ e una opportuna sezione L σ = idR+ . La funzione L si assumeregolare, non lineare, e normalmente non covessa e non concava. Per esempioFig. 4.9 illustra una tipica situazione dove questo modello si applica e Fig. 6.5.1descrive la forma tipica di una funzione L. Qui L viene stimata per fitting diuna distribuzione di dati dai colori reali dei frammenti localizzati. Dato che lafunzione L e invertibile, e sempre possibile, salvo una “ri-equalizzazione” deilivelli di grigio, poter fare in modo che L(r, g, b) = 1

3 (r + g + b). In seguito


supporremo dunque di essere di questa situazione lineare.

In questo caso, c’e un area Ω\D del dominio Ω ⊂ R2 dell’immagine doveframmenti con colori sono piazzati e un’informazione completa e disponibile, eun’altra area D (detta lacuna) dove solo i livelli di grigio sono noti.

Figura 6.7: Ricolorazione degli affreschi di Mantegna mediante la risoluzione diun problema inverso con vincoli di sparsita rispetto a wavelet.

Vogliamo ricolorare gli affreschi mediante la risoluzione del seguente proble-ma inverso:

arg min~u∈(ℓ2(Λ))3

J (~u) = µ1

∫

Ω\D|R(~u)(x) − u(x)|2dx +µ2

∫

D

|L(R(~u)(x)) − v(x)|2dx

+2τ

3∑

i=1

∑

λ∈Λ

|ui,λ|,

ove R : (ℓ2(Λ))3 → L2(Ω; R3) data da

(u1,λ, u2,λ, u3,λ)λ∈Λ 7→ (∑

λ

u1,λψλ,∑

λ

u2,λψλ,∑

λ

u3,λψλ)λ∈Λ,

e una mappa di ricostruzione (per canale) mediante una base wavelet Ψ = ψλ :λ ∈ Λ. Qui gli indici λ = (j, k) includono in forma compatta sia la scala siala posizione della wavelet. Le costanti µ1, µ2 > 0 regolano la fedelta ai dati, inaccordo con il livello di rumore.


Non e troppo difficile mostrare che il problema (6.19) e della forma (6.4).Dunque disponiamo di tutti gli ingredienti, metodi ed algoritmi, per la sua riso-luzione. Infatti, e stato dimostrato che le immagini possono essere rappresentatein modo sparso rispetto a decomposizioni wavelet [] (si pensi alla compressioneJPEG2000, basata sulle wavelets) e dunque il nostro modello risulta ragionevole.In Fig. 6.5.1 mostriamo il risulato della ricolorazione dovuta ad una applicazionedell’algoritmo di minimizzazione (6.10).

6.5.2 Magnetoencefalografia (MEG)

Il difficile problema inverso della MEG consiste nel derivare la densita di correnteJ, che fluisce nel cervello durante la sua attivita, a partire dalla misura esternadi alcune componenti del campo neuromagnetico.

Figura 6.8: La strumentazione per la MEG.

Per semplicita si assume che le sorgenti siano confinate in una regione sfericaV0, centrata nell’origine. Quindi il campo neuromagnetico B e dato semplica-mente dalla legge di Biot-Savart

B(x) := B(x,J) =µ0

4π

∫

V0

J(x) × (x− y)

‖x− y‖3R3

d(y) (6.19)

dove J sono le correnti intracellulare che fluiscono nei neuroni. Si assume inoltreche Ω ⊃ V0 e una sfera concentrica a V0 con δ := dist(V0, ∂Ω) > 0. Il camponeuromagnetico e misurato da sensori disposti sulla superficie esterna ∂Ω e nemisurano solo la componente normale alla superficie data da

Br(x) := Br(x,J) =µ0

4π

∫

V0

J(y) × (x− y) · er(x)‖x− y‖3

R3

d(y). (6.20)

Ne consegue che la componente radiale del campo presso il sensore xk ∈ ∂Ω e


data da

Br(xk,J) =µ0

4π

∫

V0

(er(xk) × (xk − y)

‖xk − y‖3R3

)· J(y) d(y)

=3∑

ℓ=1

µ0

4π

∫

V0

(er(xk) × (xk − y)

‖xk − y‖3R3

)

ℓ

Jℓ(y) d(y)

=3∑

ℓ=1

Aℓ,k Jℓ, k = 1, . . . , N , (6.21)

dove Aℓ,k : L2(V0; R) → R e l’operatore

Aℓ,k F :=µ0

4π

∫

V0

(er(xk) × (xk − y)

‖xk − y‖3R3

)

ℓ

F (y) d(y) , (6.22)

ove F ∈ L2(V0; R). Come sopra, assumiamo ora che lo spazio L2(V0; R) disponga

Figura 6.9: Identificazione di aree attive del cervello a seguito di un’analisi dellaMEG.

di una base di wavelet ψλ : λ ∈ Λ il cui operatore di ricostruzione e dato da

R~u :=∑

λ∈Λ

uλ ψλ , (6.23)

dove ~u = (uλ)λ∈Λ ∈ ℓ2(Λ; R). Se definisce Tℓ,k : ℓ2(Λ; R) → R come

Tℓ,k := Aℓ,k R , (6.24)

cosı che

Tℓ,kjℓ =

∑

λ∈Λ

jℓλ

(µ0

4π

∫

V0

(er(xk) × (xk − y)

‖xk − y‖3R3

)

ℓ

ψλ(y) d(y)

), (6.25)


e(T~j)

k:=

3∑

ℓ=1

Tℓ,kjℓ , k = 1, . . . , N, (6.26)

ove ~j e la successione dei coefficienti wavelet associati alle componenti vettorialidi J. Ancora una volta possiamo modellare la ricostruzione delle correnti comeun problema inverso con vincoli di sparsita dato da

argmin~j

‖T~j − g‖22 + 2τ

M∑

ℓ=1

∑

λ∈Λ

|jℓλ|,

ove g = (g1, . . . , gM ) e un insieme di misure del campo.

-1.5 -1 -0.5 0 0.5 1 1.5-1.5

-1

-0.5

0

0.5

1

1.5

Current

Sensor

Figura 6.10: Sinistra: La matrice compressa T ∗T associata all’operatore di

Biot-Savart nel 2D, i.e., B(x, j) :=∫V0

(x−y)⊥‖(x−y)⊥‖2

R2j(y)dy, rispetto alle wavelet

di Daubechies. Abbiamo mantenuto solo le entrate della matrice che eccedonjo10−6. Solo il 5.97% della matrice non e nullo. Destra: Un modello 2D dellaMEG. I sensori sono distribuiti su un semicerchio attorno all’area dove le correntisono distribuite.


0 10 20 30 40 50 60 700

10

20

30

40

50

60

70

0 10 20 30 40 50 60 700

10

20

30

40

50

60

70

0 10 20 30 40 50 60 700

10

20

30

40

50

60

70

0 10 20 30 40 50 60 700

10

20

30

40

50

60

70

0 10 20 30 40 50 60 700

10

20

30

40

50

60

70

0 10 20 30 40 50 60 700

10

20

30

40

50

60

70

0 10 20 30 40 50 60 700

10

20

30

40

50

60

70

0 10 20 30 40 50 60 700

10

20

30

40

50

60

70

0 10 20 30 40 50 60 700

10

20

30

40

50

60

70

0 10 20 30 40 50 60 700

10

20

30

40

50

60

70

0 10 20 30 40 50 60 700

10

20

30

40

50

60

70

0 10 20 30 40 50 60 700

10

20

30

40

50

60

70

Figura 6.11: Nella prima colonna illustramo alcune correnti fittizie come nelmodello 2D in Fig. 6.10. Nella seconda colonna mostriamo il risultato dellaricostruzione dovuto all’algoritmo (6.10). La terza colonna mostra la corrispon-dente ricostruzione mediante la risoluzione del problema inverso la cui soluzioneha norma L2 minima. Si vede come l’algoritmo (6.10) sia decisamente piuaccurato nella localizzazione spaziale delle correnti.

Bibliografia

[1] S. Anthoine, Different Wavelet-based Approaches for the Separation ofNoisy and Blurred Mixtures of Components. Application to AstrophysicalData., Ph.D. thesis, Princeton University, 2005.

[2] G. Assayag, H. G. Feichtinger, J. F. Rodrigues, Mathematics and Music,Springer Berlin, 2002

[3] T. Blumensath, M. Yaghoobi, and M. Davies, Iterative hard thresholdingand ℓ0 regularisation, Proceedings of the IEEE International Conferenceon Acoustics, Speech and Signal Processing, 2007.

[4] E. J. Candes, D. L. Donoho, New Tight Frames of Curvelets and OptimalRepresentations of Objects with Piecewise C2 Singularities, Comm. PureAppl. Math. 57 (2004), no. 2, 219–266

[5] E. J. Candes, J. Romberg, T. Tao. Robust Uncertainty Principles: Exact Si-gnal Reconstruction from Highly Incomplete Frequency Information, IEEETrans. Inform. Theory 52 (2006), no. 2, 489–509. 94A13

[6] E. J. Candes, J. Romberg, T. Tao, Stable Signal Recovery from Incompleteand Inaccurate Measurements, Comm. Pure Appl. Math. 59 (2006), no. 8,1207–1223. 94A12

[7] E. J. Candes, T. Tao, Decoding by Linear Programming, IEEE Trans.Inform. Theory 51 (2005), no. 12, 4203–4215

[8] E. J. Candes, T. Tao, Near Optimal Signal Recovery From RandomProjections: Universal Encoding Strategies?, preprint 2004.

[9] E. J. Candes, T. Tao The Dantzig Selector: Statistical Estimation When pis Much Larger than n, preprint 2005.

[10] R. Cazzato, G. Costa, A. Dal Farra, M. Fornasier, D. Toniolo, D. Tosato,C. Zanuso, Il Progetto Mantegna: storia e risultati (Italian), , in “AndreaMantegna e i Maestri della cappella Ovetari: La ricomposizione virtua-le e il restauro” (Eds. Anna Maria Spiazzi, Alberta De Nicolo’ Salmazo,Domenico Toniolo), Skira, 2006.

99

100 BIBLIOGRAFIA

[11] O. Christensen, An Introduction to Frames and Riesz Bases, Birkhauser,2003.

[12] A. Cohen, Numerical Analysis of Wavelet Methods, Studies in Mathema-tics and its Applications, 32. North-Holland Publishing Co., Amsterdam,2003

[13] L. Cohen, Time-frequency analysis, PTR. PH.,28, Chap 3, Chap 7, 1995

[14] J. W. Cooley and J. W. Tukey, An algorithm for the machine computationof complex Fourier series, Math. Comp., 19 (1965), 297-301.

[15] S. Dahlke, M. Fornasier, T. Raasch, Adaptive Frame Methods for EllipticOperator Equations, apparira in Adv. Comput. Mat.

[16] S. Dahlke, M. Fornasier, H. Rauhut, G. Steid, G. Teschke, Generali-zed Coorbit Theory, Banach frames and the Relation to Alpha-modulationSpaces, preprint 2005.

[17] S. Dahlke, M. Fornasier, K. Grochenig, Optimal Adaptive Computationsin the Jaffard Algebra and Localized Frames, preprint 2006.

[18] S. Dahlke, M. Fornasier, T. Raasch, R. Stevenson, M. Werner, Adapti-ve Frame Methods for Elliptic Operator Equations: The Steepest DescentApproach, apparira in IMA J. Numer. Anal.

[19] I. Daubechies, Ten Lectures on Wavlets, CBMS-NSF Reg. Conf. Series inApplied Math. SIAM,1992.

[20] I. Daubechies, M. Defrise, and C. DeMol, An iterative thresholding algo-rithm for linear inverse problems, Comm. Pure Appl. Math. 57 (2004),no. 11, 1413–1457.

[21] I. Daubechies, M. Fornasier, and I. Loris Accelerated projected gradientmethod for linear inverse problems with sparsity constraints, preprint, 2007.

[22] I. Daubechies, H. Landau, and Z. Landau, Gabor Time-frequency Latticesand the Wexel-Raz Identity, J. Four. Anal. Appl.,1(4):437-478, 1995.

[23] I. Daubechies and G. Teschke, Variational image restoration by means ofwavelets: Simultaneous decomposition, deblurring, and denoising., Appl.Comput. Harmon. Anal. 19 (2005), no. 1, 1–16.

[24] G. De Marco, Analisi Due/1, Decibel-Zanichelli, 1994 e successive ristampe.

[25] G. De Marco, Analisi Due/2, Decibel-Zanichelli, 1994 e successive ristampe.

[26] D. L. Donoho, Superresolution via sparsity constraints., SIAM J. Math.Anal. 23 (1992), no. 5, 1309–1331.

[27] D. L. Donoho, De-noising by soft-thresholding., IEEE Trans. Inf. Theory41 (1995), no. 3, 613–627.

BIBLIOGRAFIA 101

[28] D. L. Donoho, Nonlinear solution of linear inverse problems by wavelet-vaguelette decomposition., Appl. Comput. Harmon. Anal. 2 (1995), no. 2,101–126.

[29] D. L. Donoho, Compressed sensing, IEEE Trans. Inf. Theory 52 (2006),no. 4, 1289–1306.

[30] D. L. Donoho and I. M. Johnstone, Ideal spatial adaptation by waveletshrinkage, Biometrika 81 (1994), no. 3, 425–455.

[31] D. L. Donoho, I. M. Johnstone, G. Kerkyacharian, and D. Picard, Waveletshrinkage: asymptopia?, J. Roy. Statist. Soc. Ser. B 57 (1995), no. 2, 301–369.

[32] R.J. Duffin, A.C. Schaeffer, A Class of Nonharmonic Fourier Series,Trans.Amer.Math.Soc.,72(1952) pp.341-366

[33] B. Efron, T. Hastie, I. Johnstone, and R. Tibshirani, Least angle regression,Ann. Statist. 32 (2004), no. 2, 407–499.

[34] M. Elad, J.-L. Starck, P. Querre, and D.L. Donoho, Simultaneous car-toon and texture image inpainting using morphological component analysis(MCA), Appl. Comput. Harmon. Anal. 19 (2005), 340–358.

[35] H. G. Feichtinger, M. Fornasier, Flexible Gabor-wavelet Atomic Decompo-sitions for L2-Sobolev Spaces, Ann. Mat. Pura Appl. (4) 185 (2006), no. 1,105–131

[36] H. G. Feichtinger, P. Grobner, Banach spaces of distributions defined bydecomposition methods, I, Math.nachr., 123 (1985), pp.97-120

[37] H. G. Feichtinger, K. Grochenig, Banach Spaces Related to Integrable GroupRepresentations and Their Atomic Decomposition, I, J. Funct.Anal., toappear

[38] H. G. Feichtinger, K. Grochenig, T. Strohmer, Frame Theory and SamplingProblems in Time-Frequency Analysis and Wavelet Theory (Part II), Inter-national Journal of Wavelets, Multiresolution and Information Processing,Vol. 3, N. 3, Sept. 2005

[39] H. G. Feichtinger, F. Luef, T. Werther, A Guided Tour from Linear Algebrato the Foundations of Gabor Analysis, NuHAG, Univ. of Vienna, August2005

[40] H. G. Feichtinger, T. Strohmer, Advances in Gabor analysis, Applied andNumerical Harmonic Analysis. Birkhuser Boston, Inc., Boston, MA, 2003

[41] H. G. Feichtinger, T. Strohmer, Gabor Analysis and Algorithms. Theoryand applications, Applied and Numerical Harmonic Analysis. BirkhuserBoston, Inc., Boston, MA, 1998 T

102 BIBLIOGRAFIA

[42] G. B. Folland, Real Analysis. Modern Techniques and their Applications.2nd ed., Pure and Applied Math. A Wiley-Interscience Series of Texts,Monog., and Tract., New York, 1999.

[43] G. B. Folland, A Course in Abstract Harmonic Analysis, CRC Press, 1995.

[44] G. B. Folland, Harmonic Analysis in Phase Space, Princeton Univ. Press,no. 122 in Annals Math. Studies 1989.

[45] M. Fornasier, Function spaces inclusions and rate of convergence ofRiemann-type sums in numerical integration, Numer. Funct. Anal. Opt.,Vol. 24, Nos. 1 & 2, 2003, pp. 45–57.

[46] M. Fornasier, Banach frames for α-modulation spaces, Appl. Comp.Harmon. Anal., Vol. 22, No. 2, 2007, pp. 157–175.

[47] M. Fornasier, Matematica per l’analisi e il trattamento di immagini, noteper il Corso di Master presso l’Universita di Padova A.A. 2002/2003.

[48] M. Fornasier, Quasi-orthogonal Decompositions of Structured Frames, J.Math. Anal. Appl. 289 (2004), no. 1, 180–199.

[49] M. Fornasier, K. Grochenig, Intrinsic Localization of frames, Constr.Approx. 22 (2005), no. 3, 395–415.

[50] M. A. T. Figueiredo and R. D. Nowak, An EM algorithm for wavelet-basedimage restoration., IEEE Trans. Image Proc. 12 (2003), no. 8, 906–916.

[51] M. Fornasier and F. Pitolli, Adatpive iterative thresholding algorithms formagnetoenceophalography (MEG), (2007), preprint.

[52] M. Fornasier, H. Rauhut, Continuous Frames, Function Spaces and theDiscretization Problem, J. Fourier Anal. Appl. 11 (2005), no. 3, 245–287

[53] M. Fornasier and H. Rauhut, Recovery algorithms for vector valued datawith joint sparsity constraints, to appear in SIAM J. Numer. Anal.

[54] M. Fornasier, D. Toniolo, Fast, Robust and Efficient 2D Pattern Recogni-tion for Re-assembling Fragmented Digital Images, Pattern Recognition,Vol. 38, No. 11, 2074-2087 2005.

[55] M. Frigo, S. Johnson, The Fastest Fourier Transform in the West, MIT-LCS-TR-728, September 1997.

[56] M. Frigo, S. Johnson, FFTW user’s manual - 2.1.3, November 1999,http://www.fftw.org/.

[57] D. Gabor, Theory of Communication, J.IEE London,93(III):429-457,November 1946

BIBLIOGRAFIA 103

[58] V. K. Goyal, M. Vetterli, N.T. Tao, Quantized Overcomplete Expansionin RN : Analysis, Synthesis and Algoritms, IEEE Trans. Inform. Theory,44(1998), pp. 1631

[59] K., Grochenig, Aspects of Gabor analysis on locally compact abelian groups.Gabor analysis and algorithms, 211–231, Appl. Numer. Harmon. Anal.,Birkhuser Boston, Boston, MA, 1998.

[60] K. Grochenig, Describing Functions: Atomic Decompositions vs. Frames,Monatshefe fur matematik, 112 (1991), pp 1-41

[61] K. Grochenig, Foundations of Time-frequency Analysis, Applied andNumerical Harmonic Analysis. Birkhuser Boston, Inc., Boston, MA, 2001

[62] C. E. Heil, D. F. Walnut, Continuous and Discrete Wavelet Transform,SIAM Rewiew,31 (1989), 628-666 http://www.math.gatech.edu/

[63] A. J. E. M. Janssen, Duality and Biorthogonality for Weyl-HeisenbergFrames, J. Four. Anal. Appl.,1(4):403-437, 1995.

[64] Y. Katznelson, An Introduction to Harmonic Analysis, J.Wiley, NewYork,1968

[65] D. W. Kammler, A First Course in Fourier Analysis, Prentice Hall, UpperSaddle River, New Jersey 07458, 2000.

[66] I. Loris, G. Nolet, I. Daubechies, and F. A. Dahlen, Tomographic inversionusing ℓ1-norm regularization of wavelet coefficients, (2006), preprint.

[67] F. Luef, Gabor Analysis, Noncommutative Tori and Feichtinger’a Algebra,Lect. Notes Ser.Inst. Math.Sci. Natl.Univ. Singap.,World Sci.Publishing,River Edge, to appear.http://homepage.univie.ac.at/franz.luef/LuSing.pdf

[68] S. Mallat, A Wavelet Tour of Signal Processing, (English summary)Academic Press, Inc., San Diego, CA, 1998

[69] Y. Meyer, Ondelettes et Operateurs, I et II. (French) ActualitsMathmatiques. Hermann, Paris, 1990

[70] A. Papoulis, Signal Analysis, McGraw-Hill, New York, 1977.

[71] P. Prinz, Theory and Algorithms for Discrete 1-dimensional Gabor Frames,Tesi di dottorato, Vienna 1996

[72] S. Qian, Introduction to Time-frequency and Wavelet Transforms, PTR.PH., 2002

[73] S. Qiu, Block Circulant Gabor-matrix Structure and Discrete GaborTransforms, Opt. Engrg., pp 2872-2878,1995.

104 BIBLIOGRAFIA

[74] R. Ramlau and G. Teschke, Tikhonov replacement functionals for iterati-vely solving nonlinear operator equations., Inverse Probl. 21 (2005), no. 5,1571–1592.

[75] M. A. Rieffel, Strong Morita Equivalence of Certain Transformation GroupC∗-algebras, Math. Ann. 222 (1976), no. 1, 7–22.

[76] J.-L. Starck, M. K. Nguyen, and F. Murtagh, Wavelets and curvelets forimage deconvolution: a combined approach, Signal Proc. 83 (2003), 2279–2283.

[77] G. Teschke, Multi-frames in thresholding iterations for nonlinear operatorequations with mixed sparsity constraints, preprint, 2005.

[78] G. Teschke, Multi-frame representations in linear inverse problems withmixed multi-constraints, preprint, 2006, to appear in Appl. and Comput.Harm. Anal.

[79] R. Tibshirani, Regression shrinkage and selection via the lasso, J. Roy.Statist. Soc. Ser. B 58 (1996), no. 1, 267–288.

[80] J. Tropp, Algorithms for simultaneous sparse approximation. Part II:Convex relaxation, Signal Processing 86 (2006), 589–602.

[81] J. Tropp, K. Herrity, and A. C. Gilbert, Sparse approximation via iterativethresholding, Proceedings of the 2006 IEEE International Conference onAcoustics, Speech, and Signal Processing, 2006.

[82] R. Stevenson, Adaptive Solution of Operator Equations Using WaveletFrames, SIAM J. Numer. Anal. 41 (2003), no. 3, 1074–1100

[83] T. Strohmer, A Unified Approach to Numerical Algorithms for DiscreteGabor Expansion, in Proc. SampTA-Sampling Theory and Applications.Aveiro/Portugal, 1997 pp 297-302

[84] T. Strohmer, Numerical algorithms for discrete Gabor expansions. Gaboranalysis and algorithms, 267–294, Appl. Numer. Harmon. Anal., BirkhuserBoston, Boston, MA, 1998

[85] C. F. van Loan, Computational Frameworks for Fast Fourier Transform,SIAM, 1992.

[86] T. Werther, Y. Eldar, N. Subbanna, Dual Gabor Frames:Theory andComputational Aspects, IEEE Trans.Signal Proc., Vol 53, no.11, 2005.

Indice - Angewandte Numerische Analysis - M15/Allgemeines ... · i capitoli introduttivi...

Documents

Transcript of Indice - Angewandte Numerische Analysis - M15/Allgemeines ... · i capitoli introduttivi...