Features extraction nel dominio della
frequenza
• Fourier Transform (FT)
• Short-Time Fourier Transform (STFT)
• Continuous Wavelet Transform (CWT)
• Discrete Wavelet Transform (DWT)
• Matching Pursuit (MP)
[email protected]
Fourier Transform
Durante il 19-esimo secolo Fourier dimostrò che ogni funzione periodica può essere espressa
come somma di funzioni esponenziali complesse:

X(f ) 
 j 2ft
x
(
t
)

e
dt



x(t ) 
2 jft
X
(
f
)

e
df


Il segnale x(t) è moltiplicato per un termine esponenziale complesso caratterizzato da una
particolare frequenza f, quindi integrato su tutti gli istanti di tempo.
Se il segnale presenta componenti di frequenza f, tale componente “matcherà” con il termine
sinusoidale ed il prodotto tra di essi avrà valore (relativamente) alto.
[email protected]
Segnale stazionario
FT
Segnale non stazionario
FT
[email protected]
Fourier Transform
FT non è in grado di distinguere i due segnali in maniera efficiente (lo spettro è molto simile,
nonostante i segnali siano profondamente differenti!). Per FT i segnali sono molto simili dal
momento che essi presentano le medesime componenti frequenziali.
FT evidenzia la presenza delle componenti armoniche, ma non permette di avere informazioni
su dove tali frequenze siano posizionate (temporalmente).
Pertanto, FT non è uno strumento valido nello studio di segnali non-stazionari (la maggior
parte dei segnali reali fisici), ovvero segnali le cui componenti cambiano nel tempo.
Se non siamo interessati agli istanti di tempo in cui le componenti frequenziali si presentano,
ovvero siamo interessati a conoscere se tali componenti sono presenti (o meno) nel segnale,
allora FT può tornare utile.
[email protected]
Considerazioni:
La funzione di base e-j2πft è valutata lungo tutto l’asse temporale, per cui ci si limita ad
un’analisi globale del segnale.
D’altra parte, un’analisi locale del segnale permette una localizzazione nel tempo delle
componenti spettrali.
Per i segnali non stazionari occorre inserire nella trasformazione una dipendenza dal tempo. Il
modo più immediato consiste nel rendere locale la FT, non operando su tutto il supporto del
segnale x(t) ma su porzioni di esso, ottenute moltiplicando x(t) per una finestra che trasla nel
tempo.
[email protected]
Short-Time Fourier Transform

STFT ( , f )   x(t ) g * (t   )e j 2ft dt
g
x
•
•
•
•
•

Moltiplicando per g(t-τ) si limita lo spettro del segnale all'interno della finestra temporale
definita dalla funzione g(t-τ) stessa. Si ottiene un'informazione sul suo contenuto armonico
in un intorno dell'istante di tempo τ.
g*(t) è il complesso coniugato della funzione finestra (impulso gaussiano, impulso
rettangolare, etc).
STFT mappa la funzione 1D x(t) in un dominio 2D (tempo/frequenza).
STFT fornisce lo spettro del segnale alterato dalla presenza della finestra (la finestra ha
dimensioni prefissate -> la dimensione della finestra è un parametro critico!).
Vantaggi:
# STFT fornisce una descrizione tempo/frequenza del segnale.
# STFT supera le difficoltà della FT grazie al segnale finestra.
Per ogni τ ed f è calcolato un coefficiente della STFT.
[email protected]
g*(t-τ1)
τ1
g*(t-τ2)
τ2
g*(t-τ3)
τ3
La funzione Gaussiana rossa individua la finestra centrata in t=τ1, la blu quella in t=τ2, la verde
in t=τ3. A ciascuna finestra corrisponde una STFT (in 3 istanti di tempo differenti). In questo
modo otteniamo una vera rappresentazione tempo/frequenza del segnale.
[email protected]
Nel segnale ci sono 4 componenti in
frequenza. Il range 0-250ms presenta una
sinusoide di 300 Hz, i successivi intervalli
(tutti di 250ms) presentano sinusoidi di 200
Hz, 100 Hz e 50 Hz. Apparentemente è un
segnale non-stazionario.
Il grafico è simmetrico rispetto all’asse delle frequenze (la FT
di un segnale reale è sempre simmetrica); dal momento che
la STFT è la versione “finestrata” della FT, non ci sorprende
che la STFT sia simmetrica rispetto all’asse delle frequenze.
Notare la presenza dei 4 picchi corrispondenti alle 4
componenti in frequenza. Diversamente dalla FT, questi 4
picchi sono localizzati a differenti intervalli di tempo
lungo l’asse del tempo (il segnale originale presenta 4
componenti spettrali localizzati in intervalli di tempo differenti).
Contributo in bassa frequenza
Contributo in alta frequenza
[email protected]
STFT :: considerazioni
Grandi finestre: Buona risoluzione in frequenza (separazione delle armoniche), limitata
risoluzione nel tempo.
Finestre strette: Buona risoluzione nel tempo (capacità di localizzazione), limitata risoluzione
in frequenza. Una finestra stretta è abile a rivelare componenti in alta frequenza, ma fornisce
una risposta insoddisfacente per le componenti di bassa frequenza.
Dal principio di indeterminazione di Heisenberg: se la finestra g(t) ha un supporto “largo”, allora
la FT ha un supporto stretto (e viceversa) -> E’ impossibile ottenere contemporaneamente
da una sola finestra una buona risoluzione nel dominio del tempo e nel dominio della
frequenza!
Quindi si possono conoscere gli intervalli di tempo in cui certe bande di frequenze
esistono: Problema della Risoluzione
[email protected]
g(t)=exp(-a(t2)/2)
Finestra “piccola”
(a determina la lunghezza della finestra)
Finestra “grande”
Finestra “media”
[email protected]
Continuous Wavelet Transform
Supera i problemi della STFT usando finestre di lunghezza variabile (finestre piccole per
“catturare” il contenuto di alta frequenza garantendo buona localizzazione temporale; finestre
grandi per “catturare” il contenuto in bassa frequenza garantendo buona risoluzione
frequenziale)
1
* t 
CWT ( , s) 
x(t ) (
)dt

s
|s| t

x
con
x(t) = segnale
 = parametro di traslazione
s = parametro di scaling = 1/f
 = wavelet madre se e solo se sono soddisfatti i seguenti vincoli:

  (t )dt  0


2
|

(
t
)
|
dt  


[email protected]
Continuous Wavelet Transform
Le funzioni kernel usate nella trasformata wavelet sono tutte ottenute da una funzione
prototipo (wavelet madre), mediante operazione di scaling e/o traslazione
1
t b
a ,b (t ) 
(
)
a
a
1, 0 (t )   (t )
Con
a = parametro di scaling
b = parametro di traslazione

Continuous Wavelet Transform:
W ( a, b) 
*
x
(
t
)

a ,b (t ) dt


Valori alti del parametro di scaling “a” permettono un’analisi grossolana (non di dettaglio) del
segnale. Viceversa, bassi valori per “a” permettono di “catturare” l’informazione di dettaglio del
segnale. D’altra parte, nel dominio di Fourier, a basse frequenze (alta scala) corrispondono
informazioni globali, mentre ad alte frequenze (bassa scala) è associata l’informazione di
dettaglio (o il rumore).
[email protected]
Analisi di dettaglio (bassa scala)
In blu wavelet a supporto
compatto con parametro di
scaling unitario (bassa scala,
alta frequenza). In figura la
funzione wavelet è centrata in
to=2, 40, 90 e 140. Per ogni
versione traslata, la wavelet è
moltiplicata per il segnale (in
giallo): il prodotto è non nullo
solo dove il segnale cade
nella regione di supporto
della wavelet; altrimenti è
zero. Shiftando la wavelet il
segnale è localizzato nel
tempo;
cambiando
il
parametro di scaling, il
segnale è localizzato in
frequenza (in scala).
[email protected]
Analisi grossolana (alta scala)
Se il segnale ha componenti spettrali in
corrispondenza del parametro di scaling, il
prodotto della wavelet per il segnale è
relativamente alto. Viceversa, se la
componente spettrale corrispondente al
parametro di scaling non è presente nel
segnale, il valore del prodotto sarà
relativamente basso.
[email protected]
Segnale composto da 4 componenti frequenziali a 30 Hz, 20 Hz, 10 Hz e 5 Hz.
NB: Bassa scala -> alta frequenza (e viceversa)
Informazione
in bassa
frequenza
Informazione
in alta
frequenza
[email protected]
A piccoli valori del parametro di scaling corrispondono alte frequenze (la frequenza decresce
all’aumentare del valore di scaling). La porzione di grafico con parametro di scaling vicino
allo zero corrisponde alla più alta componente frequenziale, mentre ad alti valori di scala
corrispondono contenuti di bassa frequenza. Osserviamo che il segnale presenta prima la
componente di 30 Hz (la più alta frequenza) e ciò si manifesta a valori bassi di scaling e alla
traslazione tra 0 e 30. Successivamente è presente il contenuto a 20 Hz, la seconda
frequenza più alta, e cosi via. La componente a 5 Hz si manifesta verso la parte terminale
dell’asse di traslazione, con alti valori di scala (bassa frequenza).
[email protected]
Inverse Wavelet Transform - Synthesis
E’ possibile ricostruire il segnale originario partendo dai coefficienti wavelet con:
Soggetta alla seguente condizione di ammissibilità della costante:
FT di (t)
[email protected]
Esempi di Wavelets madri (t)
[email protected]
Analisi multirisoluzione - Discrete Wavelet Transform
Il segnale originale x[n] subisce un
filtraggio passa-alto (filtro g[n]) e passabasso (filtro h[n]). L’output dei due
filtraggi è sottocampionato (di 2). Tale
attività completa la decomposizione
wavelet di primo livello.
[email protected]
Analisi multirisoluzione - Discrete Wavelet Transform
L’analisi Wavelet è un’estensione dell’analisi di Fourier. Si trasforma il segnale di input in un
sottospazio rappresentato da versioni scalate e traslate della funzione di base denominata
mother wavelet.
S
cD
S
~ 30 coeff.
cA1
60 Samples
cA
~ 30 coeff.
cA2
cD1
cD2
Haar wavelet
 1 0  t  1/ 2

 (t )   1 1 / 2  t  1
 0 altrimenti

cA3
cD3
0 t  0  t  1
 (t )  
t  [0,1]
1
 j ,k (t )   2 (t  k )  2 j / 2  (2 j t  k )
j
[email protected]
Matching Pursuit e
dizionari overcompleti
[email protected]
Lo scopo dell’attività di rappresentazione dei dati (o segnali) è quello di
individuare un appropriato insieme di coefficienti (features) capaci di
descrivere numericamente (e possibilmente in maniera compatta) i
dati di input
Processo di
CLUSTER
classificazione
Trasformazione
SEGNALE
COEFFICIENTI
(FEATURES)
Processo di
sintesi
SEGNALE
SINTETIZZATO
[email protected]
Normalmente, nell’ambito delle attività legate al Pattern Recognition,
sono usate PCA, Trasformata di Fourier e Trasformata Wavelet al fine
di ottenere una rappresentazione del problema in questione
Ora siamo interessati a studiare sistemi non-ortogonali di funzioni
(over-completi, non basi) in cui la rappresentazione di ciascun
elemento dello spazio di input non sia necessariamente unica. Si
richiede che la rappresentazione possa essere ottenuta in modi diversi
partendo da un insieme di funzioni linearmente dipendenti
[email protected]
Consideriamo un dizionario di poche migliaia di vocaboli per cui un
particolare concetto può essere descritto usando i termini del
vocabolario a fronte di ricorrere a frasi particolarmente lunghe.
Viceversa, usando un dizionario di molti più vocaboli (ad esempio
100.000) lo stesso concetto può essere descritto brevemente, con frasi
molto corte o, addirittura, con un solo vocabolo
In generale, “dizionari” di grandi dimensioni (over-completi) permettono
di descrivere lo stesso concetto in modi differenti per cui, tra le diverse
rappresentazioni, è possibile individuare quella più opportuna per il
particolare contesto applicativo.
Si vuole individuare un insieme di funzioni che non siano
necessariamente tra loro ortonormali o linearmente indipendenti, ma che
possano costituire un operatore di rappresentazione usando espansioni
di funzioni ridondanti
[email protected]
Rappresentazione di un segnale monodimensionale usando un insieme finito
di funzioni linearmente dipendenti (ridondanti) {j}
Sia c il vettore ottenuto applicando
l’operatore di frame F al segnale f
(equazione di analisi)
Il segnale f può essere recuperato
usando l’equazione di sintesi
1
f   c j~ j   c j F T F   j  F †c

j 1
j 1
L
c j  Ff  j  f ,  j j  1...L
L
~ j
L’ultima relazione mette in evidenza la differenza tra l’operatore di frame F e
l’operatore F†. F associa un vettore di coefficienti c (le features) al segnale
decomponendolo (proiettandolo) in termini di atomi del dizionario. Viceversa,
F† (operatore di sintesi o ricostruzione) permette il recupero del segnale come
sovrapposizione degli atomi del frame duale opportunamente pesati a
coefficienti c. Le colonne di F† sono gli elementi del frame duale.
[email protected]
Funzioni di Gabor monodimensionali
La funzione di Gabor 1D è rappresentata
da una funzione Gaussiana modulata da
una funzione esponenziale complessa
Esponenziale
complesso
 1  x   2 
1
g  , , 0 ( x) 
exp  
  exp  j 0 x
2

2 
 
 
Funzione Gaussiana
La funzione di Gabor permette un’ottima
localizzazione sia nel dominio del tempo
che delle frequenze
La funzione di Gabor può essere usata per
rappresentare segnali e per estrarre
features in un processo di classificazione
[email protected]
Dizionario overcompleto 1D di Gabor
 1  x   2 
1
g  , , 0 ( x) 
exp  
  exp  j 0 x
2 
 2    
[email protected]
Dizionario overcompleto 1D di Gabor
 1  x   2 
1
g  , , 0 ( x) 
exp  
  exp  j 0 x
2 
 2    
[email protected]
Dizionario overcompleto 1D di Gabor
 1  x   2 
1
g  , , 0 ( x) 
exp  
  exp  j 0 x
2 
 2    
[email protected]
Dizionario overcompleto 1D di Haar
[email protected]
L’algoritmo di Matching Pursuit per l’estrazione di features
MPMP
decompone
un segnale
f in una
espansione
è un algoritmo
sub-ottimo
che
raffina
linearel’approssimazione
di funzioni j che del
appartengono
segnale al
dizionario overcompleto
seleziona
gli atomi
progressivamente
mediante euna
procedura
iterativa.
cheottimizza
meglio sil’energia
adattanoacquisita
alla struttura
globale
di f
MP
ad ogni
iterazione
Sia u0 = f ed L la cardinalità del dizionario F = {j}
For k = 0 to n = numero di funzioni di Gabor desiderate
Calcola per ciascun elemento di F:
cj= <uk,j> , j = 1, ..., L
Seleziona l’atomo jk che massimizza jk =argmax cj j =1,…,L
Segnale uk
(k = 0 alla prima iterazione)
Atomi del
dizionario
overcompleto
1
c1 = <uk, 1 >
2
c2 = <uk, 2 >
3
c3 = <uk, 3 >
4
c4 = <uk, 4 >
…
…
…
…
…
…
…
…
…
…
L
cL = <uk, L>
Inserisci jk nel sotto-dizionario F*
Aggiorna il segnale residuo: uk+1 = uk - cjjk
End loop
jk = argmax cj
uk+1 = uk - cj jk
Selezione
dell’atomo migliore
Aggiornamento
del segnale
residuo
[email protected]
L’algoritmo di Matching Pursuit per l’estrazione di features
Lo schema identifica solo un sotto-dizionario di pochi atomi
in grado di rappresentare qualsiasi segnale monodimensionale
Sotto-dizionario F*
SCHEMA MP
(Insieme di 256 funzioni del
dizionario overcompleto)
Rappresentazione di qualsiasi segnale mediante proiezione sul sotto-dizionario F*
Si usa l’equazione di analisi della Teoria dei frame
c j  F * f  j  f ,  j j  1...L
[email protected]
L’algoritmo di Matching Pursuit con l’uso del dizionario di Gabor
Errore di ricostruzione per il segnale usato nel processo di generazione del
sotto-dizionario F* usando lo schema di MP (equazioni di analisi e sintesi)
I primi atomi di F* catturano
le informazioni di bassa
frequenza del segnale
I successivi atomi
descrivono l’informazione
di dettaglio (contenuto in
alta frequenza)
L’errore residuo converge
a zero in un numero finito
di iterazioni
[email protected]
L’algoritmo di Matching Pursuit con l’uso del dizionario di Gabor
Errore di ricostruzione per un segnale non usato dallo schema
MP per la generazione di F* (equazioni di analisi e sintesi)
Si garantisce ancora una
buona approssimazione
del segnale
Il sotto-dizionario F*
spanna bene lo spazio dei
segnali 1D
Alta capacità di sintesi di F*
[email protected]
L’algoritmo di Matching Pursuit con l’uso del dizionario di Haar
Errore di ricostruzione per il segnale usato nel processo di generazione del
sotto-dizionario F* usando lo schema di MP (equazioni di analisi e sintesi)
La convergenza dell’errore
di ricostruzione è più
veloce rispetto al caso di
dizionari di Gabor
[email protected]
L’algoritmo di Matching Pursuit – Versione generalizzata
Let R0p = fp with 1 p m and m=number of training signals
Loop: for i = 0 to n=number of desired basis elements
Compute Cg for each atom { i } in the dictionary
Select the i-th basis element i that maximizes Cg
Update coefficients for each block pi = < Rifp, i >
Update the residue blocks Ri+1fp = Rifp - pi i
Increment the iteration counter i = i + 1
End loop
At iteration i let (p0, …, pi-1) the coefficients computed
through iteration i-1 to represent block fp with 1 p m
Compute the centroid  of (p0, …, pi-1, < Rifp, i >)
Compute the mean distance U between  and
(p0, …, pi-1, < Rifp, i >)

U  1/m  μ   0p ,..., ip-1 , R i f p , i

Return U as value of Cg
[email protected]
[email protected]
[email protected]
Esempio: An automated approach for vascular oral network detection
RGB
RED
Prevailing channel
Few information
GREEN
Best vessel/background
contrast
BLUE
Low contrast
[email protected]
GREEN
9000
RESTRICTED
DYNAMICS
8000
7000
6000
5000
4000
3000
2000
1000
0
HE
0
50
100
150
200
250
0
50
100
150
200
250
9000
8000
7000
6000
5000
4000
3000
POOR
QUALITY
CLAHE
2000
1000
0
STRETCHED
DYNAMICS
3000
2500
2000
1500
1000
500
0
0
50
100
150
200
250
[email protected]
2D Gaussian
function
Sinusoidal
modulator
2D Gabor
wavelet
[email protected]
2D Gaussian
function
Sinusoidal
modulator
2D Gabor
wavelet
[email protected]
2D Gaussian
function
Sinusoidal
modulator
2D Gabor
wavelet
[email protected]
CWT
can be values
meantof as
The intensity
thecorrelation
equalized
measure
between
vesselsopixels
green channel
are the
reversed
that
and
a particular
Gabor wavelet
blood
vessels present
values higher
than
the
background
pixels
(otherwise CWT detects background!)
CWT presents stronger responses
when Gabor wavelet is tuned to
vessels pixels (maximum correlation)
A grey-levels image is defined by
evaluating the scaling parameters for
which
the
maximum
modulus
response is realized by the CWT
In white low-frequency information (background) is
shown. Darker regions are related to different size vessels
(high frequency)
[email protected]
Image binarization is realized according to the Otsu scheme[3]
The automated threshold guarantees maximum distance of the
clusters (background, vessel) and minimum intra-cluster variance
18000
16000
14000
12000
VESSEL
10000
8000
BK
THRESHOLD
6000
4000
2000
0
0
50
100
150
200
250
[3] N. Otsu, "A Threshold Selection Method from Gray-Level Histograms," IEEE Transactions on
Systems, Man and Cybernetics, vol.9, no.1, 1979, 62-66
[email protected]
Operatori morfologici
Thinning
Segmentazione finale
[email protected]
Esempio: Classificazione da immagini telerilevate mediante Gabor function
[email protected]