Dyskretna transformata Fouriera
Dyskretna transformata Fouriera
Discrete Fourier Transform – DFT
Discrete Fourier Transform – DFT
Fourierowskie przedstawienie
ciągów o skończonej długości
(lub ciągów okresowych dla próbek w jednym okresie)
Dyskretna transformata Fouriera
Dyskretna transformata Fouriera
Discrete Fourier Transform – DFT
Discrete Fourier Transform – DFT
Fourierowskie przedstawienie
ciągów o skończonej długości
(lub ciągów okresowych dla próbek w jednym okresie)
Interpretacja
Interpretacja
Ciąg x(n) o skończonej długości N traktowany jako ciąg okresowy
x
N
(n) o okresie N, który w ramach jednego okresu jest identyczny
z ciągiem o skończonej długości (pomijamy na razie wycinanie
fragmentu sygnału – tzw.
okienkowanie
, gdy
x
N
(n) = x(n)w
N
(n)
).
W przeciwieństwie
do szeregu Fouriera
okresowej funkcji ciągłej
x(t) =
f
T
(t)
f
T
(t) =
Σ
F
k
exp(jk
ω
0
t), gdzie
ω
0
= 2
π
/ T
T
∞
∞
∞
∞
(2a)
k = -
∞
∞
∞
∞
2
szereg Fouriera
okresowego ciągu (szeregu)
ma tylko N różnych
zespolonych wyrazów wykładniczych, gdyż
k-ta dyskretna harmoniczna jest okresowa z okresem N
e
j(2
π/Ν)
kn
= e
j(2
π/Ν) (
k+rN)n
F
k
= 1/T
∫
f
T
(t) exp(-jk
ω
0
t
)
d
t
,
T
0
(2b)
k = -
∞
∞
∞
∞
k
ω
0
t
t=nTs
= k
nT
s
2
π
/T
= knT
s
2
π
/ NT
s
= 2
π
kn / N
Rozważmy wyrażenie
exp(-jk
ω
0
t)dt
we wzorze (2b) na
współczynniki Fouriera rozkładu funkcji w przedziale T dla
przypadku czasu dyskretnego, tzn.
t = nT
s
i
dt=T
s
.
Wtedy całka (2b) , przy podstawieniu dt =T i T = NT
3
Wtedy całka (2b) , przy podstawieniu dt =T
s
i T = NT
s
przechodzi w szereg (1)
(transformatę prostą przyjęło się definiować z pominięciem parametru 1/N,
co musimy uwzględnić w transformacie odwrotnej, mnożąc wynik przez 1/N.
F
k
= 1/T
∫
f
T
(t) exp(-jk
ω
0
t)dt
T
0
⇒
⇒
⇒
⇒
Uwagi
Uwagi
• X(k) nazywamy dyskretną transformatą Fouriera
skończonego ciągu x(n) o długości N próbek (wyrazów).
•
X(k) również długość N próbek w dziedzinie dyskretnej
pulsacji (częstotliwości).
Stosując notację DFT można zapisać
N
j
e
W
/
2π
−
=
4
•
Stosując notację DFT można zapisać
w postaci
N
j
N
e
W
/
2π
−
=
(3)
1
0
,
]
[
]
[
1
0
−
≤
≤
∑
=
−
=
N
k
W
n
x
k
X
N
n
n
k
N
Próbek widma jest tyle, ile próbek sygnału. Zespolone, czy to 2 razy więcej?
Właściwości symetrii (właściwie antysymetrii) wynikające z okresowości funkcji exp
(2
π
k/N), to głównie wykorzystywane zależności
X(k) = X*(-k) = X*(N-k).
Fizyczna częstotliwość k-tego
prążka widma to
f = kf
s
/N
, a widmo jest
antysymetryczne względem f=0 i f=fs/2.
Dla sygnału
rzeczywistego
o długości N próbek należy zatem wyznaczać tylko N
niezależnych wartości, mimo że widmo ma N zespolonych prążków, czyli 2N
5
niezależnych wartości, mimo że widmo ma N zespolonych prążków, czyli 2N
wartości (N opisujących części rzeczywiste i N dla części urojonych),
co ilustruje poniższa tabela.
Liczba N próbek
rzeczywistego
sygnału
Zakres
wyznaczanych
prążków k
Liczba
prązków
k
Razem części
Re[X(k)] i
Im[X(k) ]
Łączna liczba wartości
do wyznaczenia
parzysta
0, N/2
N/2+1
N+2
N, ponieważ
Im[X(0)]=Im[X(N/2)]=0
nieparzysta
0, (N-1)/2
(N-1)/2 + 1
N+1
N, ponieważ
Im[X(0)]=0
Rozdzielczość widma (dodawanie zerowych próbek do x(n) )
Rozdzielczość widma DTF, (odległość między kolejnymi prążkami)
wyrażona w Hz, wynosi
∆
f = f
s
/N i jest ona równa częstotliwości
pierwszego prążka (pierwszej harmonicznej).
Rozdzielczość
rośnie
, gdy
∆
f
maleje
, czyli gdy rośnie N.
6
7
Odwrotna transformata DFT
Odwrotna transformata DFT -- IDFT
IDFT
1
0
,
]
[
1
]
[
1
0
−
≤
≤
∑
=
−
=
−
N
n
W
k
X
N
n
x
N
k
kn
N
(4)
Odwrotne przekształcenie do DFT dane jest wzorem
8
Aby wykazać słuszność wzoru (4) (chodzi też o współczynnik N) wykonamy
DFT postaci (4). W tym celu pomnóżmy obustronnie (4) przez
dla i zsumujmy dla .
n
N
W
l
1
0
−
≤
≤
N
n
Otrzymamy
∑
∑
∑
−
=
−
=
−
−
=
=
1
0
1
0
1
0
1
N
n
n
N
N
k
kn
N
N
n
n
N
W
W
k
X
N
W
n
x
l
l
]
[
]
[
∑ ∑
−
−
−
−
=
1
1
1
N
N
n
k
W
k
X
)
(
]
[
l
9
?
∑ ∑
=
=
−
−
=
0
0
1
n
k
n
k
N
W
k
X
N
)
(
]
[
l
∑ ∑
−
=
−
=
−
−
=
1
0
1
0
1
N
k
N
n
n
k
N
W
k
X
N
)
(
]
[
l
Wykorzystując zależność
=
∑
−
=
−
−
1
0
N
n
n
k
N
W
)
(
l
,
,
0
N
,
rN
k
=
− l
for
otherwise
r
całkowite
otrzymujemy, że prawa strona poprzedniej równości wynosi
]
[l
X
10
(
wybierając
k = l
z przedziałów i
)
,
]
[
]
[
l
l
X
W
n
x
N
n
n
N
=
∑
−
=
1
0
co jest poprawną zależnością na DFT ciągu x(n) otrzymanego z (4).
1
0
−
≤
≤
N
k
a zatem
N
= 8
0
1,5
1
3
4
5
6
7
0, 4
2, 6
itd.
itd.
11
N
= 8
k - l =
3
(dodatnie, w lewo co
¾
π
)
k - l = -
6
(ujemne, w prawo 3/2 π)
=
∑
−
=
−
−
1
0
N
n
n
k
N
W
)
(
l
,
,
0
N
,
rN
k
=
−
l
for
otherwise
r
całkowite
N
j
N
e
W
/
2π
−
=
2
5
7
3, 7
Interpretacja sygnału w czasie przez DFT
Interpretacja sygnału w czasie przez DFT
12
Okresowość widma DFT
Okresowość widma DFT
13
Przyczyna okresowości DFT
Przyczyna okresowości DFT
Dla jakich pulsacji w całym przedziale t
∈
(-
∞
,
∞
) zachodzi exp(j
ω
1
t)=exp(j
ω
2
t)?
Tylko dla
ω
2
=
ω
1
, co nie wymusza okresowości widma sygnału ciągłego w czasie.
Dla jakich pulsacji (liczb k) dla całkowitych n w całym przedziale
n
∈
(-
∞
,
∞
) zachodzi exp(j2
π
k
1
n/N)=exp(j2
π
k
2
n/N)?
Dla wszystkich k
2
= k
1
+
rN , gdzie r
∈
C, co powoduje
okresowości widma
dyskretnego
dla sygnału dyskretnego w czasie (szeregu, ciągu próbek czasowych).
dyskretnego
dla sygnału dyskretnego w czasie (szeregu, ciągu próbek czasowych).
Dla jakich pulsacji ciągłych i dla całkowitych n w całym przedziale
n
∈
(-
∞
,
∞
) zachodzi exp(j
ω
1
nT
s
)= exp(j
ω
2
nT
s
) – dotyczy tzw. DTFT ?
Dla wszystkich
ω
2
=
ω
1
+
r
ω
s
=
ω
1
+
r2
π
/Ts , co powoduje
okresowość widma
ciągłego
dla sygnału dyskretnego w czasie
.
Przyczyną okresowości widma jest dyskretyzacja w czasie.
Przykłady
Przykłady DFT
DFT
1
1
0
0
1
−
≤
≤
=
N
n
n
,
,
=
]
[n
x
Dyskretny impuls w zerze –
δ
(n)
w skończonym przedziale czasu n
15
i jego DFT
1
0
0
1
0
=
=
=
∑
−
=
N
N
n
kn
N
W
x
W
n
x
k
X
]
[
]
[
]
[
dla
1
0
−
≤
≤
N
k
Przykłady
Przykłady DFT
DFT cd
cd..
Przesunięty impuls
=
]
[n
y
1
1
1
0
0
1
−
≤
≤
+
−
≤
≤
=
N
n
m
m
n
m
n
,
,
,
i jego DFT
km
km
N
kn
−1
dla
16
km
N
km
N
N
n
kn
N
W
W
m
y
W
n
y
k
Y
=
=
=
∑
−
=
]
[
]
[
]
[
1
0
dla
1
0
−
≤
≤
N
k
Dla stałej c widmo X(k) jest zerowe, za wyjątkiem prążka X(0) = Nc,
gdyż…..
Przykłady
Przykłady DFT
DFT cd
cd..
Ciąg N próbek cosinusa dla o okresie N/r (r-ta harmoniczna)
dla
1
0
−
≤
≤
N
n
1
0
),
/
2
cos(
]
[
−
≤
≤
π
=
N
r
N
rn
n
g
Korzystając ze wzoru Eulera, otrzymamy dla przyjętej notacji
N
j
N
e
W
/
2π
−
=
17
(
)
N
rn
j
N
rn
j
e
e
n
g
/
2
/
2
2
1
]
[
π
−
π
+
=
DFT ciągu g(n) obliczamy zatem jako
∑
−
=
=
1
0
N
n
kn
N
W
n
g
k
G
]
[
]
[
,
2
1
1
0
)
(
1
0
)
(
∑
+
∑
=
−
=
+
−
=
−
−
N
n
n
k
r
N
N
n
n
k
r
N
W
W
a wykorzystując już poprzednio stosowaną tożsamość
−1
N
,
N
,
rN
k
=
− l
for
r
18
=
∑
−
=
−
−
1
0
N
n
n
k
N
W
)
(
l
,
,
0
N
,
rN
k
=
− l
for
otherwise
r
całkowite
otrzymujemy
−
=
=
=
otherwise
0
for
2
for
2
,
,
/
,
/
]
[
r
N
k
N
r
k
N
k
G
Dla N =16 i r =3 wykresy DTFT (
linia ciągła
) i DFT (prążki
o
) są następujące
f=3/16 = 0,1875
f=13/16 = 0,812
brak, gdyż…..
jest natomiast minus 0,1875
Wykresy są przedstawiane w funkcji
znormalizowanej
częstotliwości f (przy czym 1= f
s
),
znormalizowanej
pulsacji
ω
(przy czym T
s
=1s,
ω
s
= 2
π
)
lub numerów harmonicznych
k .
Najczęściej rysowana jest tylko pierwsza połowa wykresu od 0 łącznie z
punktem środkowym
.
f
ω
2
π
k
0 1 2 3 4 ···
N
0,5
π
N
/2
Jeżeli w okresie obserwacji N nie mieści się całkowita liczba okresów
harmonicznej o okresie L próbek, to pulsacja tej harmonicznej wypada
pomiędzy punktami
ω
k
= 2
π
k / N.
20
W sygnale z rys b) nie istnieje składowa 12/8 =
1.5
Obserwujemy tutaj tzw. zjawisko przecieku (czestotliwości)
DFT
DFT w zapisie macierzowym
w zapisie macierzowym
Wyrażenie na wartości transformaty X(k) w postaci definicji (1)
możemy zapisać w postaci macierzowej jako
x
D
X
=
(5)
21
x
D
X
N
=
(5)
gdzie odpowiednie wektory i macierz D są zdefiniowane jako
[
]
T
N
X
X
X
]
[
.....
]
[
]
[
1
1
0
−
=
X
[
]
T
N
x
x
x
]
[
.....
]
[
]
[
1
1
0
−
=
x
Dla zapisu (5), przekształcenie
odwrotne (4) możemy zapisać
jako
X
D
x
1
−
=
N
(6)
gdzie przyjęto oznaczenie
n 0 1 2 . . . N-1
k
22
gdzie przyjęto oznaczenie
*
N
N
N
D
D
1
1
=
−
n 0 1 2 . . . N-1
k
0
1
.
.
.
N-1
Tworzenie macierzy D
N
z wektorów na kole jednostkowym dla N=8
Skok o kąt 2π/N = π/4
• obrót w kierunku
ujemnym (
w prawo
)
dla prostego DFT,
w celu konstrukcji
macierzy
D
N
23
N
• obrót w kierunku
dodatnim (w lewo)
dla odwrotnej DFT,
w celu konstrukcji
macierzy D
N
-1
Tworzenie macierzy D
N
z wektorów na kole jednostkowym cd.
24
próbki widma - prążki
Próbki sygnału
Fast Fourier
Fast Fourier Transform
Transform ((FFT
FFT))
Idea wykorzystania cząstkowych wyników mnożenia wartości próbek przez elementy macierzy
D
N
oraz zależności
W
N
(rN + k)
= W
N
k
i W
N
(N/2 + k)
= -W
N
k
wynikających z okresowości
W
N
k
.
MOTYLEK
X[5] = x
0
W
N
0
+ x
1
W
8
5
+ x
2
W
8
10
+ x
3
W
8
15
+ x
4
W
8
20
+ x
5
W
8
25
+ x
6
W
8
30
+ x
7
W
8
35
=
=
x
0
W
N
0
−−−−
x
1
W
8
1
+
x
2
W
8
2
−−−−
x
3
W
8
3
−−−−
x
4
W
8
0
+
x
5
W
8
1
−−−−
x
6
W
8
2
+
x
7
W
8
3
Próbkowanie transformaty
Próbkowanie transformaty DTFT
DTFT a oryginał
a oryginał
• rozważamy nieograniczony ciąg x(n) i jego DTFT w postaci ,
•
próbkujemy w N równomiernie rozłożonych punktach
dla otrzymując zbiór próbek
,
)
(
ω
j
e
X
)
(
ω
j
e
X
N
k
k
/
2π
=
ω
1
0
−
≤
≤
N
k
)}
(
{
k
j
e
X
ω
26
•
rozważamy otrzymany ciąg jako DFT ciągu czasowego y(n),
•
wyznaczając IDFT porównujemy otrzymany ciąg y(n) z oryginalnym
x(n).
ω
k
∑
=
∞
−∞
=
ω
−
ω
l
l
l
j
j
e
x
e
X
]
[
)
(
(7)
Skoro
a zatem
)
(
)
(
]
[
/
2
N
k
j
j
e
X
e
X
k
Y
k
π
ω
=
=
∑
=
∑
=
∞
∞
π
−
l
l
l
l
k
N
N
k
j
W
x
e
x
]
[
]
[
/
2
27
∑
=
∑
=
−∞
=
−∞
=
π
−
l
l
l
l
l
l
k
N
N
k
j
W
x
e
x
]
[
]
[
/
2
i z odwrotnej dyskretnej transformaty IDFT otrzymujemy
∑
=
−
=
−
1
0
]
[
1
]
[
N
k
n
k
N
W
k
Y
N
n
y
czyli
∑ ∑
−
=
−
∞
−∞
=
=
1
0
1
N
k
n
k
N
k
N
W
W
x
N
n
y
l
l
l
]
[
]
[
=
∑
∑
−
=
−
−
∞
−∞
=
1
0
1
N
k
n
k
N
W
N
x
)
(
]
[
l
l
l
a posługując się tożsamością
28
a posługując się tożsamością
=
∑
−
=
−
−
1
0
1
N
n
r
n
k
N
W
N
)
(
,
,
0
1
mN
n
r
+
=
for
otherwise
otrzymujemy zależność na y(n) w zależności od oryginału x(n) w postaci
1
0
−
≤
≤
+
=
∑
∞
−∞
=
N
n
mN
n
x
n
y
m
],
[
]
[
(8)
Z zależności (8) wynika, że jeżeli ciąg x(n) ma ma skończoną liczbę
wyrazów M mniejszą niż liczba próbek widma N, to wypadkowy ciąg
okresowy
1
0
−
≤
≤
+
=
∑
∞
−∞
=
N
n
mN
n
x
n
y
m
],
[
]
[
(8)
będzie w ramach każdego okresu N powtórzeniem ciągu nieokresowego
29
y
[n] = x[n]
dla i .
N
M ≤
1
0
−
≤
≤
N
n
(9)
W przypadku, gdy
M
> N
zachodzi tzw.
time-domain aliasing.
T
Time
ime--domain aliasing
domain aliasing -- przykład
przykład
Niech
}
{
]}
[
{
5
4
3
2
1
0
=
n
x
↑
Próbkując DTFT w punktach i stosując
czteropunktową IDFT otrzymujemy sekwencję y(n) w postaci
)
(
ω
j
e
X
4
/
2 k
k
π
=
ω
30
czteropunktową IDFT otrzymujemy sekwencję y(n) w postaci
]
[
]
[
]
[
]
[
4
4
−
+
+
+
=
n
x
n
x
n
x
n
y
dla
3
0
≤
≤ n
czyli
}
{
]}
[
{
3
2
6
4
=
n
y
↑
co przedstawiono schematycznie na kolejnym rysunku
n -4 -3 -2 -1 0 1 2 3 4 5 6 7 8 9
x(n)
......
0 1 2 3 4 5
.......
x(n-4)
0 1 2 3 4 5
x(n+4)
0 1 2 3 4 5
____________________________________________________
31
____________________________________________________
y(n) =
Σ
kolumn
4 6 2 3
4 6 2 3
Ciąg x(n) nie może być odtworzony z ciągu y(n).
Inny przykład
y
12
(n)
32
y
7
(n)
Kołowe przesunięcie ciągu
Kołowe przesunięcie ciągu –
– splot kołowy
splot kołowy
Rozpatrzmy ciąg x(n) o długości N zdefiniowany dla .
Próbki dla
n
< 0
i
są równe zeru.
1
0
−
≤
≤
N
n
N
n ≥
Definiujemy kołowe przesunięcie ciągu jako operację modulo N
]
[
]
[
N
o
c
n
n
x
n
x
〉
−
〈
=
(10)
33
N
o
c
Dla
0
>
o
n
(jest to tzw. kołowe przesunięcie ciągu w prawo) otrzymujemy
+
−
−
=
],
[
],
[
]
[
n
n
N
x
n
n
x
n
x
o
o
c
o
o
n
n
N
n
n
<
≤
−
≤
≤
0
for
1
for
Dla przykładu, dla N = 6
x
c1
[n]
x
c4
[n]
34
n
0
= 1
n
0
= 4
x
c
[n]
Inny przykład dla N = 6 i n
0
= -2
(jest to tzw. kołowe przesunięcie ciągu
w lewo
)
35
Splot liniowy
ciągów o długości N-próbek dany jest wzorem
(zakłada się, że poza tym przedziałem ciągi przyjmuję zerowe wartości)
2
2
0
1
0
−
≤
≤
−
=
∑
−
=
N
n
m
n
h
m
g
n
y
N
m
L
],
[
]
[
]
[
(11)
Pytanie: Z czego wynika taka długość splotu liniowego?
36
Splot kołowy
(circular convolution)
definiuje się jako
1
0
],
[
]
[
]
[
1
0
−
≤
≤
∑
〉
−
〈
=
−
=
N
n
m
n
h
m
g
n
y
N
m
N
C
(12)
i oznacza symbolem
y
[n] = g[n] h[n]
N
Przykład splotu kołowego
Niech dane będą czteropunktowe ciągi g(n) i h(n) w postaci
},
{
]}
[
{
1
0
2
1
=
n
g
}
{
]}
[
{
1
1
2
2
=
n
h
↑
↑
2
]
[n
g
2
]
[n
h
37
n
3
2
1
0
1
]
[n
g
n
3
2
1
0
1
Splot kołowy o długości N = 4 zapisujemy dla jako
,
]
[
]
[
]
[
]
[
]
[
3
0
4
∑
〉
−
〈
=
=
=
m
C
m
n
h
m
g
n
h
n
g
n
y
4
3
0
≤
≤ n
a kolejne wartości splotu kołowego wynoszą
∑
〉
〈
−
=
=
3
0
4
]
[
]
[
]
0
[
m
C
m
h
m
g
y
]
1
[
]
3
[
]
2
[
]
2
[
]
3
[
]
1
[
]
0
[
]
0
[
h
g
h
g
h
g
h
g
+
+
+
=
6
2
1
1
0
1
2
2
1
=
×
+
×
+
×
+
×
=
)
(
)
(
)
(
)
(
∑
〉
−
〈
=
3
4
]
1
[
]
[
]
1
[
C
m
h
m
g
y
38
∑
〉
−
〈
=
=0
4
]
1
[
]
[
]
1
[
m
C
m
h
m
g
y
]
[
]
[
]
[
]
[
]
[
]
[
]
[
]
[
2
3
3
2
0
1
1
0
h
g
h
g
h
g
h
g
+
+
+
=
n
3
2
1
0
1
2
]
[n
g
n
3
2
1
0
1
2
]
[n
h
itd.
7
1
1
1
0
2
2
2
1
=
×
+
×
+
×
+
×
=
)
(
)
(
)
(
)
(
Ostatecznie
Inny przykład
39
Animacje
Jeszcze inny przykład
40
Splot liniowy a kołowy
Splot liniowy a kołowy
Rozpatrzmy ciągi o długości L= N / 2 uzupełnione zerami do długości N
41
Obliczanie splotu kołowego metodą
Obliczanie splotu kołowego metodą DFT
DFT
Rozpatrzmy ponownie ciągi g(n) i h(n)
n
3
2
1
0
1
2
]
[n
g
n
3
2
1
0
1
2
]
[n
h
oraz ich splot kołowy y (n) w postaci
42
oraz ich splot kołowy y
C
(n) w postaci
0
1
2
3
g(m)
1
2
0
1
0
h (0-m)
2
1
1
2
6
1
h (1-m)
2
2
1
1
7
2
h (2-m)
1
2
2
1
6
3
h (3-m)
1
1
2
2
5
funkcje
m
y(n)
n
Cykliczny rejestr przesuwający
Czteropunktowe DFT ciągów odpowiednio wynoszą
4
2
1
0
/
]
[
]
[
]
[
k
j
e
g
g
k
G
π
−
+
=
4
6
4
4
3
2
/
/
]
[
]
[
k
j
k
j
e
g
e
g
π
π
−
−
+
+
3
0
2
1
2
3
2
≤
≤
+
+
=
−
−
k
e
e
k
j
k
j
,
/
/
π
π
43
,
]
[
2
1
2
1
2
−
=
−
−
=
G
,
]
[
j
j
j
G
−
=
+
−
=
1
2
1
1
j
j
j
G
+
=
−
+
=
1
2
1
3]
[
,
]
[
4
1
2
1
0
=
+
+
=
G
4
2
1
0
/
]
[
]
[
]
[
k
j
e
h
h
k
H
π
−
+
=
4
6
4
4
3
2
/
/
]
[
]
[
k
j
k
j
e
h
e
h
π
π
−
−
+
+
3
0
2
2
2
3
2
≤
≤
+
+
+
=
−
−
−
k
e
e
e
k
j
k
j
k
j
,
/
/
π
π
π
,
]
[
6
1
1
2
2
0
=
+
+
+
=
H
44
,
]
[
6
1
1
2
2
0
=
+
+
+
=
H
,
]
[
0
1
1
2
2
2
=
−
+
−
=
H
,
]
[
j
j
j
H
−
=
+
−
−
=
1
1
2
2
1
j
j
j
H
+
=
−
−
+
=
1
1
2
2
3]
[
Poprzednie transformaty łatwiej obliczyć używając zapisu macierzowego.
+
−
−
=
−
−
−
−
−
−
=
=
j
j
j
j
j
j
g
g
g
g
G
G
G
G
1
2
1
4
1
0
2
1
1
1
1
1
1
1
1
1
1
1
1
1
3
2
1
0
3
2
1
0
4
]
[
]
[
]
[
]
[
]
[
]
[
]
[
]
[
D
45
+
−
=
−
−
−
−
−
−
=
=
j
j
j
j
j
j
h
h
h
h
H
H
H
H
1
0
1
6
1
1
2
2
1
1
1
1
1
1
1
1
1
1
1
1
3
2
1
0
3
2
1
0
4
]
[
]
[
]
[
]
[
]
[
]
[
]
[
]
[
D
Iloczyn widm (
Jest to tzw. mnożenie tablicowe a nie wektorowe !
Polega na mnożeniu wyrazów na tej samej pozycji.
)
3
0
≤
≤
=
k
k
H
k
G
k
Y
C
],
[
]
[
]
[
wynosi
24
0
0
0
H
G
Y
C
]
[
]
[
]
[
46
−
=
=
2
0
2
24
3
3
2
2
1
1
0
0
3
2
1
0
j
j
H
G
H
G
H
G
H
G
Y
Y
Y
Y
C
C
C
C
]
[
]
[
]
[
]
[
]
[
]
[
]
[
]
[
]
[
]
[
]
[
]
[
a czteropunktowa odwrotna DFT
=
]
[
]
[
]
[
]
[
*
]
[
]
[
]
[
]
[
3
2
1
0
4
1
3
2
1
0
4
C
C
C
C
C
C
C
C
Y
Y
Y
Y
y
y
y
y
D
6
24
1
1
1
1
47
=
−
−
−
−
−
−
−
=
5
6
7
6
2
0
2
24
1
1
1
1
1
1
1
1
1
1
1
1
4
1
j
j
j
j
j
j
Splot liniowy wyznaczany metodą
Splot liniowy wyznaczany metodą DFT
DFT
• Niech g[n] i h[n] będą ciągami o długościach
odpowiednio N i M próbek
• Splot liniowy ma długość
• Zdefiniujmy odpowiednie wydłużone sekwencje
uzupełnione zerami do długości L próbek
1
−
+
=
M
N
L
48
uzupełnione zerami do długości L próbek
wtedy
]
[
]
[
]
[
]
[
]
[
]
[
n
h
n
g
n
y
n
h
n
g
n
y
C
L
=
=
=
*
L
Zero-padding
with
Zero-padding
with
zeros
1)
(
−
M
−
−
+
)
(
1
M
N
point DFT
×
g
[n]
h
[n]
Length-N
]
[n
g
e
]
[n
h
e
−
−
+
)
(
1
M
N
−
−
+
)
(
1
M
N
point IDFT
]
[n
y
L
e
e
49
with
zeros
1)
( −
N
point DFT
e
Length-M
Length
-
)
(
1
−
+M
N
Problem – w przypadku liczenia tą metodą funkcji korelacji
należy odpowiednio przyporządkować próbki
z tablicy IDFT odpowiednim wartościom argumentu m
korelacji – wartości dla m
≥
0 leżą prawidłowo, dla m
<
0 znajdują
się na pozycji (m+L)
(dla splotu wszystkie próbki są prawidłowo ułożone) - przykład.
Splot liniowy metodą
Splot liniowy metodą DFT
DFT –
– przykład
przykład
Dane są ciągi x(n) o długości N =2 i y(n) długości M =3.
Ich splot linowy wyznaczony w dziedzinie czasowej o długości L=N+M-1 = 4
przedstawiono poniżej.
50
-2
-1
0
1
2
3
x (m)
1
1
0
y (0-m)
2
2
1
1
1
y (1-m)
2
2
1
3
2
y (2-m)
2
2
1
4
3
y (3-m)
2
2
1
2
x(n)
∗
∗
∗
∗ y(n)
m
ciągi
n
Wyznaczmy teraz DFT ciągów g(n) i h(n) uzupełnionych zerami do długości
splotu liniowego L=4 i oznaczonych jako x
e
(n) i y
e
(n).
X
e
(k)
51
Y
e
(k)
Splot liniowy (równy kołowemu) wyznaczamy jako
IDFT [X
e
(k)Y
e
(k)]
=
52
Wartości próbek otrzymanego wektora odpowiadają wartościom splotu liniowego x y.
∗∗∗∗
Estymator funkcji korelacji wyznaczany metodą
Estymator funkcji korelacji wyznaczany metodą DFT
DFT
Definicje korelacji wzajemnej (przypomnienie):
korelacja sygnału x(n) z y(n)
R
xy
(m) = E[x(n+m) y(n)] = E(x(n) y(n-m)]
korelacja sygnału y(n) z x(n)
(13a)
53
R
yx
(m) = E[y(n+m) x(n)] = E[y(n) x(n-m)]
(13b)
Zajmiemy się estymatą korelacji w postaci (13b), rzadziej prezentowaną
w podręcznikach, a wykorzystywaną przecież do identyfikacji odpowiedzi
impulsowej układu liniowego przy pobudzeniu szumem białym.
Estymatą odpowiedzi impulsowej
g(n)
jest wtedy
estymata korelacji wzajemnej
sygnału na wyjściu układu y(n) z sygnałem wejściowym u(n).
Estymator funkcji korelacji wyznaczany metodą
Estymator funkcji korelacji wyznaczany metodą DFT
DFT cd
cd..
Estymator obciążony
korelacji wzajemnej lub skrośnej (13b)
ciągów y(n) z x(n) o długościach równych N próbek dany jest wzorem
1
N
__
Σ
n=0
N-m-1
x(n) y(n + m)
dla m
≥
0 ( y przesuwane w lewo)
__
1
Σ
N-1
x(n) y(n+m)
dla m
<
0 ( y przesuwane w prawo)
54
(14)
__
1
N
Σ
n=m
x(n) y(n+m)
dla m
<
0 ( y przesuwane w prawo)
lub najczęściej
1
N
__
Σ
n=0
N-1
x(n) y(n+m)
dla -(N-1)
≤
m
≤
(N-1)
kiedy to zakłada się, że brak próbki y(i + m) dla występującej w sumie (14)
kolejnej próbki x(i) oznacza, że jest ona równa zeru.
←
←
←
← y(n)
m
y(n) →
→
→
→
R
yx
=
∧
∧∧
∧
Bez względu na współczynnik przed sumą w (13) – decydujący o obciążeniu
estymaty, zawsze należy wyznaczyć (2N-1) razy sumę
Σ
N-1
n=0
x(n) y(n+m) =
C
yx
(m) =
(15)
Wyznaczenie wartości (14) w dziedzinie czasu charakteryzuje się
złożonością obliczeniową O(N
2
). Wyznaczenie wartości (14) w dziedzinie
widmowej poprzez zastosowanie DFT (a właściwie jej szybkiego algorytmu
Fast Fourier Transform (FFT) obniża złożoność do postaci Nlog
2
N.
n=0
N-1
Σ
y(n) x(n - m)
55
Porównując (15) ze wzorem na splot liniowy (16)
s
L
(m) = y(m)
∗
x(m) =
Σ
n=0
N-1
y(n) x(m - n)
(16)
C
yx
(m) = y(m)
∗
x(-m)
(17a)
otrzymujemy bardzo wygodną obliczeniowo zależność
C
xy
(m) =
(17b)
oraz analogicznie
x(m)
∗
y(-m)
Obliczenie wartości (16) w dziedzinie widmowej odbywa się według procedury
C
yx
(m) = IDFT [Y
e
(k) X
e
*(k)]
,
(18)
gdzie i oznaczają odpowiednio DFT uzupełnionego zerami
ciągu y(n) i sprzężoną DFT uzupełnionego zerami ciągu x(n), na przykład
jak na rysunku poniżej (uzupełnienie do długości L = M + N – 1 = 4)
X*
e
(k)
Y
e
(k)]
56
M =2
N =3
L = 4
L = 4
Estymator funkcji korelacji metodą
Estymator funkcji korelacji metodą DFT
DFT -- przykład
przykład
Estymaty korelacji (14) wyznaczana w dziedzinie czasu wynosi
-2
-1
0
1
2
3
x (n)
1
1
w prawo
-1
y (n-1)
1
2
2
1
bez
0
y (n)
1
2
2
3
w lewo
1
y (n+1)
1
2
2
4
c
yx
(m)
przesunięcie
m
ciągi
n
57
Natomiast dyskretne transformaty ciągów y
e
(n) i x
e
(n) są równe
X
e
(k)
Y
e
(k)
X*
e
(k)
w lewo
2
y (n+2)
1
2
2
2
Ostatecznie
IDFT [Y
e
(k) X
e
*(k)]
=
Odczyt do m=N/2
PRZYCZYNOWOŚĆ
58
W kolumnie kolejno występują: c
yx
(0) = 3
c
yx
(1) = 4
c
yx
(2) = 2
c
yx
(-1) = IDFT(L-1) = IDFT(3) = 1
Sposób uporządkowania próbek znajdujących się w rejestrze
odwrotnej DFT o wymiarze N próbek, aby odpowiadały one
próbkom wyznaczanej estymaty funkcji korelacji
(często N oznaczane jest w tym przypadku jako L – długość wyciętego okienka z oryginału funkcji
w dziedzinie czasu lub przesunięcia (ang. lag).
Cxy(3)
w rzeczywistości
jest równe zeru (!),
gdyż m
∈
(-1, 2).
Zatem otrzymane
Cxy(3) = Cxy(N-1)=
Cxy(-1),
N
2N
N/2
Cxy(-1),
DFT iloczynu ciągów
DFT iloczynu ciągów –
– splot
splot
kołowy
kołowy
widm DFT
widm DFT
Rozpatrzmy iloczyn ciągów g(n) i h(n) w postaci
n
3
2
1
0
1
2
]
[n
g
n
3
2
1
0
1
2
]
[n
h
i ich iloczyn y(n) = g(n)h(n) , przedstawione w Tab.1
60
i ich iloczyn y(n) = g(n)h(n) , przedstawione w Tab.1
n
0
1
2
3
g(n)
1
2
0
1
h(n)
2
2
1
1
y(n)
2
4
0
1
Tab. 1 Wartości próbek ciągów i ich iloczynu
Wyznaczmy splot kołowy uprzednio wyznaczonych widm G(k) i H(k)
zgodnie ze wzorem definiującym splot kołowy widm w postaci
61
Y(k) = G(m) H[(k-m)
N
] .
1
__
N
Σ
m =0
m =N-1
Dla N = 4 otrzymujemy następującą tabelkę wyznaczania
splotu
kołowego widm
G(k) i H(k) oraz wynik IDFT ze splotu, równy y(n).
zgodnie ze wzorem definiującym splot kołowy widm w postaci
Transformata odwrotna IDFT z wyniku splotu kołowego widm Y(k) wynosi
0
1
2
3
G(m)
4
1 - j
-2
1 + j
0
H (0-m)|4
6
1 + j
0
1 - j
7
1
H (1-m)|4
1 - j
6
1 + j
0
2 - 3j
2
H (2-m)|4
0
1 - j
6
1 + j
-3
3
H (3-m)|4
1 + j
0
1 - j
6
2 + 3j
funkcje
m
Y(k)
k
62
Otrzymany wynik równy jest próbkom iloczynu y(n)=g(n)h(n) z tabeli 1.
Podstawowe własności DFT
Podstawowe własności DFT
63