L’algoritmo FFT
Corso di Metodi di Trattamento del Segnale
Se N = 2M
N −1
Fk = ∑ fn e
−
2 π ikn
N
n=0
M −1
= ∑ f2n e
−
2 π ik (2n)
N
n=0
M −1
= ∑ f2n e
n=0
M −1
+ ∑ f2n+1e
−
2 π ik (2n+1)
N
n=0
2 π ikn
−
M
⎛
+⎜e
⎝
trasformata dei campioni
con indice pari
W =e
2π i
−
N
−
2π i
N
⎞
⎟⎠
k
M −1
∑f
e
2n+1
−
2 π ikn
M
n=0
trasformata dei campioni
con indice dispari
Lemma di Danielson-Lanczos
Fk = Fk(e) + W k Fk(o)
trasformata originale
(periodicità N)
trasformata dei campioni
con indice pari (periodicità N/2)
trasformata dei campioni con indice
dispari (periodicità N/2)
Se N è una potenza di 2 si può applicare ripetutamente il lemma:
dopo log2N applicazioni del lemma si devono calcolare N
sottotrasformate, ma ciascuna di queste sottotrasformate si
riferisce ad un sottoinsieme che contiene un solo campione, e
quindi la trasformata coincide con il campione stesso.
Per calcolare la trasformata dobbiamo allora fare le seguenti
operazioni:
1. catalogare i campioni nell'ordine giusto, in modo che i
campioni riordinati ci diano direttamente le trasformate dei
sottoinsiemi di lunghezza 1. Questo passo iniziale lo si può fare
una volta per tutte con O(N) operazioni.
2. partire dal livello più basso e costruire le sottotrasformate di
ordine più elevato.
Ciascuna ricostruzione richiede log2N passi, e poiché la
trasformata ha N componenti, ci vogliono in tutto O(Nlog2N)
operazioni.
Esempio: DFT di un insieme di 8 = 23 campioni
primo livello:
Fk = Fk(e) + W k Fk(o)
3
(e)
k
F
= ∑e
−
2 π ikn
4
−
2 π imn
4
n=0
3
Fk(o) = ∑ e
n=0
W =e
−
2π i
8
campioni in posizione pari
f2n
f2n +1
campioni in posizione dispari
secondo livello
Fk(e) = Fk(ee) + W k Fk(eo)
Fk(o) = Fk(oe) + W k Fk(oo)
Fk(ee) = f0 e
(eo)
k
F
(oe)
k
F
−
= f2 e
= f1e
2 π i·0
·k
2
2 π i·1
−
·k
2
−
2 π i·0
·k
2
−
2 π i·1
·k
2
Fk(oo) = f3e
+ f4 e
+ f6 e
+ f5 e
−
2 π i·2
·k
2
2 π i·3
−
·k
2
−
2 π i·2
·k
2
−
2 π i·3
·k
2
+ f7 e
W =e
−
2π i
4
campioni in posizione pari tra i
campioni in posizione pari
campioni in posizione dispari
tra i campioni in posizione pari
campioni in posizione pari tra i
campioni in posizione dispari
campioni in posizione dispari tra i
campioni in posizione dispari
terzo livello
Fk(ee) = Fk(eee) + W k Fk(eeo)
Fk(eo) = Fk(eoe) + W k Fk(eoo)
Fk(oe) = Fk(oee) + W k Fk(oeo)
(oo)
k
F
=F
W =e
(ooe)
k
−
2π i
2
+W F
= −1
k
(ooo)
k
Fk(eee) = f0
Fk(eeo) = f4
Fk(eoe) = f2
Fk(eoo) = f6
Fk(oee) = f1
Fk(oeo) = f5
Fk(ooe) = f3
Fk(ooo) = f7
assegnazione dei campioni alle trasformate del livello più basso
Fk(eee) = f0
(eeo)
k
F
= f4
Fk(eoe) = f2
Fk(eoo) = f6
Fk(oee) = f1
Fk(oeo) = f5
Fk(ooe) = f3
Fk(ooo) = f7
e→ 0
o→1
000 → 0;
001 → 4;
010 → 2;
011 → 6;
100 → 1;
101 → 5;
110 → 3;
111 → 7;
Esercizi con LabView:
• Realizzazione di un programma che calcola la DFT a
partire dalle operazioni elementari.
• Utilizzo del programma per mostrare che il tempo di
calcolo cresce proporzionalmente a N2, e confronto
con il VI di LabView che calcola la FFT.
Time plot of Wolf Sunspot Number, 1700-2007. This time series is known
to have an irregular cycle with period near 11 years. The long-term mean
is 49.9. Data source: http://sidc.oma.be/sunspot-data/
Densità spe*rale nel caso della DFT nel caso continuo la densità spettrale
1
S(ω ) = lim
T →∞ T
+T /2
∫
2
f (t)e−iω t dt
−T /2
rappresenta la potenza media del segnale ad una certa frequenza (angolare), o
anche la fluttuazione quadratica media a questa stessa frequenza
nel caso della DFT possiamo ragionare in modo analogo
flu*uazione quadra7ca media del segnale campionato 1 N −1 2
W = ∑ fn Δt
T n=0
1
=
N
T = N Δt
N −1
∑
n=0
1
f =
N
2
n
N −1
1
N
N −1
∑ ∑F e
n=0
m=0
m
2 π imn 2
N
2 π imn
1 N −1 N −1 * − 2 πNikn
= 3 ∑ ∑ Fk e
Fm e N
N n = 0 k, m = 0
1 N −1 * N −1 2 π i(mN− k )n
= 3 ∑ Fk Fm ∑ e
N k, m = 0
n=0
1 N −1 *
= 3 ∑ Fk Fm Nδ m, k
N k, m = 0
1
= 2
N
N −1
∑F
k=0
k
Periodogramma
2
2
Fk
Sk = 2
N
Simmetrie di DFT e periodogramma nel caso di segnali reali N −1
2 π ikn
*
N
n
Fk* = ∑ f e
n=0
N −1
= ∑ fn e
n=0
2 π ikn
N
N −1
= ∑ fn e
2 π ikn
2 π iNn
−
N
N
e
n=0
Sk = S N − k
inoltre (se N è pari):
F0 ∈!
FN* 2 = FN 2
⇒ FN 2 ∈!
N −1
= ∑ fn e
n=0
−
2 π in( N − k )
N
= FN − k
Esempio: spe*ro di un segnale sinuisoidale (
A
Fk = N δ k, k0 + δ k, N − k0
2
2
)
(
Fk
A2
Sk = 2 =
δ k, k0 + δ k, N − k0
N
4
)
Definizione di spettro unilatero
2
⎧
F
⎪ S = 0
0
⎪
N2
⎪⎪
1
2
S
=
F
+ FN − k
⎨
k
k
2
N
⎪
2
⎪
F
⎪ SN /2 = N /22
N
⎪⎩
(
2
)
(k ≠ 0
e k ≠ N 2)
Teorema di Wiener-‐Kintchine nel caso della DFT 1
Rm =
N
N −1
∑
n=0
fn* fn + m
N −1
∑R e
m=0
funzione di autocorrelazione nel caso di segnali campiona7 −
m
2π i
mk
N
⎛1
= ∑⎜
m=0 ⎝ N
N −1
N −1
∑
n=0
2π i
−
mk
⎞
*
N
fn fn + m ⎟ e
⎠
2π i
2π i
nk
−
(n + m) k
1
* N
N
= ∑ fn e
fn + m e
N n, m
la DFT della funzione di autocorrelazione è uguale a N volte la densità spe*rale 1
=
N
N −1
∑f
n=0
2
2π i
N −1
nk
* N
n
m=0
e
Fk
=
= NSk
N
∑f
m
e
−
2π i
mk
N
X k = DFT { xn }k
xn = IDFT { X k }n
⎛ 2π ink ⎞
= ∑ xn exp ⎜ −
⎟⎠
⎝
N
n=0
N −1
1
=
N
1
= DFT { X k }− n
N
Implementazione “ovvia” della IDFT: • DFT • inversione della sequenza • divisione per N ⎛ 2π ink ⎞
∑ Xk exp ⎜⎝ N ⎟⎠
k=0
N −1
operazione di scambio xn = an + ibn
ixn* = i ( an − ibn ) = bn + ian
N −1
1
⎛ 2π ink ⎞
*
*
ixn = ∑ iX k exp ⎜ −
=
⎟
⎝
N k=0
N ⎠
1
*
= DFT iX k n
N
( )
{ }
Implementazione efficiente della IDFT: • scambio Re-‐Im • DFT • scambio Re-‐Im • divisione per N ... We tried to assemble the 10 algorithms with the greatest influence on the
development and practice of science and engineering in the 20th century.
Following is our list (here, the list is in chronological order; however, the
articles appear in no particular order):
•
•
•
•
•
•
•
Metropolis Algorithm for Monte Carlo
Simplex Method for Linear Programming
Krylov Subspace Iteration Methods
The Decompositional Approach to Matrix Computations
The Fortran Optimizing Compiler
QR Algorithm for Computing Eigenvalues
Quicksort Algorithm for Sorting
• Fast Fourier Transform
• Integer Relation Detection
• Fast Multipole Method
(from Dongarra and Sullivan: “The Top-Ten Algorithms”, Comp. Sci. Eng.
(2000) n. 1, p. 22)
To(Qp) =digital processing for the purposes of pattern recognition and Wiener
as
REFERENCES
with that of a class of orthogonal
filtering. Its performance is compared
/2M-i= cos (k
TkQp)
the
Karhunenof
that
to
closely
and
found
to
compare
is
transforms
Gx(O)"Z X(m)
[1] S. Winograd, "A new algorithm for inner product," IEEE Trans.
Discrete Cosine Lo'eve
Transform
of
The performances
optimal. July
which isvol.
known
transform,
1968.
C-1 7,topp.be693-694,
(Corresp.),
Comput.
m=o
the Karhunen-Lo'eve and discrete cosine transforms are also found to where Tk(Qp) is the kth Che
2 M-i
chosen
Now, in (2),
criterion. 1974
rate-distortion JANUARY
theCOMPUTERS,
with respect toON
compare
tp isX(m)Co
IEEEclosely
TRANSACTIONS
by Gx(k)=[31
M m0
Index Terms-Discrete cosine transform, discrete Fourier transform,
transform,byrate
an disKarhunen-Loeve
Haar transform,
feature
the filter is represented
Transfonn
Cosinefiltering
required, then it can be
applications,
Wiener
InDiscrete
discreteselection,
kth(2D
where Gx(k) istpthe Cos
filtering.
scalar
and
vector
Wiener
transform,
Walsh-Hadamard
ew algorithm will require (M Xtortion,
M) matrix G. The estimate X of data vector X is given by GZ,
the set of basis vectors {
RAOapproxiK. R.that
ANDimplies
N. and
T. NATARAJAN,
AHMED,
noise vector. This
N is the
where Z = X + N
class of discrete Chebysh
2A
in (2), on
(3) polynomia
Substituting
of
Use
X.
INTRODUCTION
to
are
arithmetic
2M
required
compute
operations
mately
Abstract-A discrete cosine transform (DCT) is defined and an algo- that Chebyshev
(16a)
)
of Ittois
number
substantial
atransform
in anwhich
G fast
orthogonal
developed.
Fourier
using
ityields
compute
rithm
totransforms
respect
interestis with
increasing
been
hasathe
years there
In recent
1
A
can
hence
and
equal
small incosine
magnitude,
relatively
elements
insetthe
usedbearea
of area
digitalof
transform
the be
general
shown
theofdiscrete
in can
=
transforms
orthogonal
class
athat
using are
To (P) =_1
at
is
realized
load
in
computation
reduction
Thus
a
significant
to
zero.
To(Qp)
and Wiener
two
recognition
itself towards
purposes of pattern
addresses
the correspondence
processing forThis
digital
signal processing.
(16b)
error.
estimation
in
the
increase
mean-square
of
a
small
the filtering.
expense
of orthogonal
of a class
recogniwith that
pattem
namely,
compared
processing,
performance
withis image
problems Itsassociated
A
(2
transform
Fourier
discrete
(WHT),
transform
Walsh-Hadamard
The
hm will generally be more
= cos
= cos
]
TkQp)
the
Karhunenof
that
to
closely
and
found
to
compare
2].
is
and
Wiener
filtering
[
tion
transforms
[1
Tk(p)
have of
slant transform
(ST),
and tothetransforms
transform
Haar recognition,
(DFT),
a
noninvertible
The
performances
enable
be
optimal.
which (HT),
isorthogonal
known
Inthe
Lo'eve
transform,
pattern
[4] - [ 9 since these
2],a reduced
[ 1 ], to[transforms
various applications
for from
considered
beentransformation
also found to where Tk(Qp) is the kth C
aredimensionality
cosine
space
pattern
discrete
the
and
Karhunen-Lo'eve
(17a) are the
be computed using fast algorithms.
From
that
th
(4) init follows
transforms that can
orthogonal
is chos
(2),
Now,
criterion.
implemented
rate-distortion
to
be
the
to
scheme
respect
closely
with
This
a
classification
compare
allows
feature
space.
tp
that
with
The performance of these transforms is generally compared in classificaincrease
1
only a issmall
with which
features,(KLT)
with
substantially less
by [31
to be optimal
known
transform
Karhunen-Loeve
of the
transform,
Fourier
discrete
transform,
cosine
Index
Terms-Discrete
tion
error.to the following performance measures: variance distribu(17b) withfeature
respect
transform, rate disselection, Haar transform, Karhunen-Loeve
error criterion [ 2], [4], and
tp Cos
tion [1 ], estimation using the mean-square
filtering.
and scalar
vector
Wiener
transform,
Walsh-Hadamard
tortion,
(2m
is
there
KLT
is
optimal,
functionJanuary
[5]. Although
the rate-distortion
revised June
7, 1973.
29, 1973;the
Manuscript received
= cos
Tk(m)
and
of Electrical
the Departments
N. Ahmed
is with
its fast computation
[ 1] Engineering
that enables
algorithm
ion of a block of sampled no general
Kans.
Manhattan,
University,
Kansas State
Science,
Computer
Substituting (3) in (2),
INTRODUCTION
cosineoftransform
(DCT) is introa discrete
correspondence,
In this
Kansas
Electrical Engineering,
is with the
hod which trades an inT. Natarajan
Department
It is to
Kans.
Manhattan,
its fast interest
University,
that
enables
computation.
withyears
an algorithm
along
with respect
umber of multiplications. ducedState
an increasing
been
hasDepartment
therethe
InK.
recent
1
A
UniEngineering,
of
Electrical
with
R.
Rao
is
that
to
more
DCT
of
the
closely
the
compares
that
performance
shown
is
ime) of a multiplication
general area of digital Comparing (5)Towith
in the
transformsTex.
of atorthogonal
76010.
using
Arlington, Arlington,
Texas
versitya ofclass
(P) =(1) _we
HT.
and
of
the
WHT,
to
the
DFT,
relative
performances
KLT
ithm is always more com- of the
two
towards
signal processing. This correspondence addresses itself
f the correlation, and it is
processing, namely, pattem recogniwith COSINE
image TRANSFORM
problems associated
A
(
DISCRETE
r processing 128 or fewer
= cos
]
tion [1 and Wiener filtering [ 2].
Tk(p)
(M - 1)aisnoninvertible
defmed
0, 1, * *, enable
X(m), m =transforms
sequenceorthogonal
lags" for L < 10 log2 2N.
of a data
TheInDCT
recognition,
pattern
as transformation from the pattern space to a reduced dimensionality
/2M-i This allows a classification scheme to be implemented From (4) it follows that t
feature space.
=
=
Discrete Cosine Transform
2.0
1.5
1.0
0.5
0.0
-2
-1
0
1
2
Estensione periodica indotta dalla DFT
3
Discrete Cosine Transform
2.0
1.5
1.0
0.5
0.0
-3
-2
-1
0
1
2
Estensione periodica che si assume nella DCT
3
2Ck =
N −1
∑
n=− N
=
N −1
∑
n=− N
dove
fn e−2 π i( n+1/2 )k/2 N
N −1
⎛ 2π ( n + 1 2 ) k ⎞
⎛ 2π ( n + 1 2 ) k ⎞
fn cos ⎜
⎟⎠ + i ∑ fn sin ⎜⎝ −
⎟⎠ =
2N
2N
⎝
n=− N
fn = f− n−1
N −1
∑
n=− N
se
N −1
∑
n=− N
⎛ π (n + 1 2) k ⎞
fn cos ⎜
⎟⎠
N
⎝
−N ≤ n < 0
⎛ 2π ( n + 1 2 ) k ⎞
fn cos ⎜
⎟⎠ =
2N
⎝
=
−1
∑
n=− N
−1
∑
n=− N
⎛ π ( n + 1 2 ) k ⎞ N −1
⎛ π (n + 1 2) k ⎞
fn cos ⎜
+
f
cos
∑
ℓ
⎟⎠
⎜⎝
⎟⎠
N
N
⎝
n=0
⎛ π ( n + 1 2 ) k ⎞ N −1
⎛ π (n + 1 2) k ⎞
f− n−1 cos ⎜
+
f
cos
⎟⎠ ∑ n
⎜⎝
⎟⎠
N
N
⎝
n=0
⎛ π ( n − 1 2 ) k ⎞ N −1
⎛ π (n + 1 2) k ⎞
= ∑ fn−1 cos ⎜
+
f
cos
⎟⎠ ∑ n
⎜⎝
⎟⎠
N
N
⎝
N
n=1
n=0
⎛ π ( n + 1 2 ) k ⎞ N −1
⎛ π (n + 1 2) k ⎞
= ∑ fn cos ⎜
+
f
cos
⎟⎠ ∑ n
⎜⎝
⎟⎠
N
N
⎝
N −1
n=0
⎛ π (n + 1 2) k ⎞
= 2∑ fn cos ⎜
⎟⎠
N
⎝
n=0
N −1
n=0
Discrete Cosine Transform
N −1
⎛ π (n + 1 2) k ⎞
1 N −1
−2 π i( n+1/2 )k/2 N
C k = ∑ fn e
= ∑ fn cos ⎜
⎟⎠
2 n=− N
N
⎝
n=0
N −1
C 0 = ∑ fn
n=0
⎛ π ( n + 1 2 ) k ⎞ N −1
⎛ π (n + 1 2) k ⎞
C− k = ∑ fn cos ⎜ −
= ∑ fn cos ⎜
= Ck
⎟
⎟
N
N
⎝
⎠
⎝
⎠
N −1
n=0
N −1
C N = ∑ fn cos ⎡⎣π ( n + 1 2 ) ⎤⎦ = 0
n=0
n=0
Inverse Discrete Cosine Transform
2Ck =
N −1
∑
fn e−2 π i( n+1/2 )k/2 N = e− π ik/2 N Fk(2 N )
Fk(2 N ) = 2eπ ik/2 N Ck
n=− N
1 2 N −1 (2 N ) 2 π ikn/2 N 1 N −1 (2 N ) 2 π ikn/2 N 2 N −1
π ik/2 N 2 π ikn/2 N
fn = ∑ Fk e
=
F
e
=
C
e
e
∑
∑
k
k
N k=0
N k=− N
N k=− N
−1
⎤
2 ⎡ N −1
2 π ik( n+1/2 )/2 N
= ⎢ ∑ Ck e
+ ∑ Ck e2 π ik( n+1/2 )/2 N ⎥
N ⎣ k=0
⎦
k=− N
N −1
2 ⎡ N −1
2 π ik( n+1/2 )/2 N
−2 π ik( n+1/2 )/2 N ⎤
= ⎢ ∑ Ck e
+ ∑ Ck e
⎥
N ⎣ k=0
⎦
k=1
N −1
⎤
2⎡
= ⎢C0 + ∑ Ck e2 π ik( n+1/2 )/2 N + e−2 π ik( n+1/2 )/2 N ⎥
N⎣
⎦
k=1
(
N −1
4 ⎡1
⎛π
⎞⎤
= ⎢ C0 + ∑ Ck cos ⎜ ( n + 1 2 ) k ⎟ ⎥
⎝N
⎠⎦
N ⎣2
k=1
)
Compressione JPEG
La compressione di immagini JPEG utilizza l’algoritmo DCT. I passi essenziali
sono
1. Suddivisione dell’immagine in blocchi di 8 x 8 pixel
2. Trasformazione di ciascun blocco con una DCT di tipo II bidimensionale
⎡π
⎤
⎡π
⎤
Cℓm = ∑ I jk cos ⎢ ( 2 j + 1) ℓ ⎥ cos ⎢ ( 2k + 1) m ⎥
⎣8
⎦
⎣8
⎦
j,k=0
N −1
3. Ciascuna componente della DCT viene divisa per un coefficiente che
corrisponde alla capacità dell’occhio di distinguere le variazioni di luminosità alle
diverse frequenze spaziali: tipicamente il coefficiente è piccolo alle basse
frequenze e grande alle alte frequenze (l’occhio non distingue bene le variazioni
di intensità ad alta frequenza spaziale)
4. Le componenti spaziali pesate, vengono arrotondate. Tipicamente questa
operazione di quantizzazione annulla molti valori della DCT, e quindi il blocco 8 x 8
viene in realtà codificato da un numero di valori molto inferiore a 64.
5. I valori restanti vengono quindi compressi con un algoritmo standard come il
runlength encoding oppure l’Huffman coding
La decodifica dei dati viene realizzata rovesciando la sequenza di codifica:
1. Decompressione dati
2. Moltiplicazione dei valori per i coefficienti di quantizzazione
3. Trasformata coseno inversa
© Copyright 2024 ExpyDoc
