Dyskretna transformata Fouriera
Dyskretna transformata Fouriera
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
Fourierowskie przedstawienie
ciągów o skończonej długości
(lub ciągów okresowych dla próbek w jednym okresie)
2
Discrete Fourier Transform
Discrete Fourier Transform
-
-
DFT
DFT
Definicja
DFT dla skończonego ciągu x(n) o długości N
próbek dla i jego DTFT , jest to
równomiernie spróbkowana transformata wzdłuż
osi ω, w przedziale w punktach
dla .
1
0
−
≤
≤
N
n
)
(
ω
j
e
X
)
(
ω
j
e
X
π
≤
ω
≤
2
0
N
k
k
/
2π
=
ω
1
0
−
≤
≤
N
k
dla .
1
0
−
≤
≤
N
k
(1)
Definicja DFT
3
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)
szereg Fouriera
ciągu (szeregu) okresowego
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
T
(t) =
Σ
F
k
exp(jk
ω
0
t), gdzie
ω
0
= 2
π
/ T
F
k
= 1/T
∫
f
T
(t) exp(-jk
ω
0
t)dt ,
T
0
∞
∞
∞
∞
k = -
∞
∞
∞
∞
(2a)
(2b)
4
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
s
i T = NT
s
przechodzi w szereg (1).
F
k
= 1/T
∫
f
T
(t) exp(-jk
ω
0
t)dt
T
0
⇒
⇒
⇒
⇒
5
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ć
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?
6
DTFT
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
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
1
0
−
≤
≤
N
n
8
?
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
]
[
]
[
∑ ∑
−
=
−
=
−
−
=
1
0
1
0
1
N
n
N
k
n
k
N
W
k
X
N
)
(
]
[
∑ ∑
−
=
−
=
−
−
=
1
0
1
0
1
N
k
N
n
n
k
N
W
k
X
N
)
(
]
[
9
Wykorzystując zależność
=
∑
−
=
−
−
1
0
N
n
n
k
N
W
)
(
,
,
0
N
,
rN
k
=
−
for
otherwise
r
całkowite
otrzymujemy, że prawa strona poprzedniej równości wynosi
(
wybierając
k = l
z przedziałów i
)
,
]
[
X
]
[
]
[
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
10
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
)
(
,
,
0
N
,
rN
k
=
−
for
otherwise
r
całkowite
N
j
N
e
W
/
2π
−
=
0
1,5
1
2
3
4
5
6
7
0, 4
2, 6
3, 7
itd.
itd.
Interpretacja sygna
Interpretacja sygna
ł
ł
u w czasie przez DFT
u w czasie przez DFT
11
Okresowo
Okresowo
ść
ść
widma DFT
widma DFT
12
Przyczyna okresowo
Przyczyna okresowo
ś
ś
ci DFT
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).
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.
(
)
m
Granice sumowania w DFT
-1
- komentarz
1
0
,
]
[
1
]
[
1
0
−
≤
≤
∑
=
−
=
−
N
n
W
k
X
N
n
x
N
k
kn
N
Dlaczego?
→
x
a
(t)
1
0
,
]
[
1
]
[
1
0
−
≤
≤
∑
=
−
=
−
N
n
W
k
X
N
n
x
N
k
kn
N
m
→x(n)
15
Przyk
Przyk
ł
ł
ady DFT
ady DFT
1
1
0
0
1
−
≤
≤
=
N
n
n
,
,
=
]
[n
x
Dyskretny impuls w zerze –
δ
(n)
w skończonym przedziale czasu n
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
16
Przyk
Przyk
ł
ł
ady DFT
ady 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
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ż…..
17
Przyk
Przyk
ł
ł
ady DFT
ady 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
rn
j
N
rn
j
e
e
n
g
/
2
/
2
2
1
]
[
π
−
π
+
=
N
j
N
e
W
/
2π
−
=
18
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
0
N
n
n
k
N
W
)
(
,
,
0
N
,
rN
k
=
−
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
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=3/16 = 0,1875
f
ω
2
π
k
0 1 2 3 4 ···
N
0,5
π
N
/2
f=13/16 = 0,812
brak, gdyż…..
jest natomiast minus 0,1875
20
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.
W sygnale z rys b) istnieje składowa 12/8 =
1.5
21
DFT w zapisie macierzowym
DFT 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
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
22
Dla zapisu (5), przekształcenie
odwrotne (4) możemy zapisać
jako
X
D
x
1
−
=
N
(6)
gdzie przyjęto oznaczenie
*
N
N
N
D
D
1
1
=
−
n 0 1 2 . . . N-1
k
0
1
.
.
.
N-1
23
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
• obrót w kierunku
dodatnim (w lewo)
dla odwrotnej DFT,
w celu konstrukcji
macierzy D
N
-1
24
Tworzenie macierzy D
N
z wektorów na kole jednostkowym cd.
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
.
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
MOTYLEK
26
Pr
Pr
ó
ó
bkowanie transformaty DTFT a orygina
bkowanie transformaty DTFT 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
,
•
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).
)
(
ω
j
e
X
)
(
ω
j
e
X
N
k
k
/
2π
=
ω
1
0
−
≤
≤
N
k
)}
(
{
k
j
e
X
ω
ω
k
27
∑
=
∞
−∞
=
ω
−
ω
j
j
e
x
e
X
]
[
)
(
(7)
Skoro
a zatem
)
(
)
(
]
[
/
2
N
k
j
j
e
X
e
X
k
Y
k
π
ω
=
=
∑
=
∑
=
∞
−∞
=
∞
−∞
=
π
−
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
28
czyli
∑ ∑
−
=
−
∞
−∞
=
=
1
0
1
N
k
n
k
N
k
N
W
W
x
N
n
y
]
[
]
[
=
∑
∑
−
=
−
−
∞
−∞
=
1
0
1
N
k
n
k
N
W
N
x
)
(
]
[
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)
29
Na podstawie (8) wnosimy, że w wyniku dyskretyzacji
ciągłego widma DTFT ciągu x(n)
otrzymujemy w wyniku odwrotnej dyskretnej
transformaty DFT ciąg y(n) będący
nieskończoną sumą poprzesuwanych o okres N
ciągów pierwotnych x(n).
Jest to analog zjawiska próbkowania sygnału x(t),
którego widmo jest sumą poprzesuwanych widm
sygnału analogowego (przypomnienie).
)
(
ω
j
e
X
30
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
y
[n] = x[n]
dla i .
N
M ≤
1
0
−
≤
≤
N
n
(9)
W przypadku, gdy
M
> N
zachodzi tzw.
time-domain aliasing.
31
T
T
ime
ime
-
-
domain aliasing
domain aliasing
-
-
przyk
przyk
ł
ł
ad
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
π
=
ω
]
[
]
[
]
[
]
[
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
32
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
____________________________________________________
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).
33
Inny przykład
y
12
(n)
y
7
(n)
34
Ko
Ko
ł
ł
owe przesuni
owe przesuni
ę
ę
cie ci
cie ci
ą
ą
gu
gu
–
–
splot ko
splot ko
ł
ł
owy
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)
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
35
Dla przykładu, dla N = 6
]
[n
x
]
1
[
6
〉
−
〈n
x
n
0
= 1
]
4
[
6
〉
−
〈n
x
n
0
= 4
36
Inny przykład dla N = 6 i n
0
= -2
(jest to tzw. kołowe przesunięcie ciągu
w lewo
)
37
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?
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
38
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
↑
↑
n
3
2
1
0
1
2
]
[n
g
n
3
2
1
0
1
2
]
[n
h
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
39
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
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
=
×
+
×
+
×
+
×
=
)
(
)
(
)
(
)
(
40
Ostatecznie
]
[n
y
C
Inny przykład
[ ]
[ ]
(
)
(
)
[
]
∑
−
=
−
=
1
N
0
m
N
2
1
3
m
n
x
m
x
n
x
41
Jeszcze inny przykład
42
Splot liniowy a ko
Splot liniowy a ko
ł
ł
owy
owy
Rozpatrzmy ciągi o długości L= N / 2 uzupełnione zerami do długości N
43
Obliczanie splotu ko
Obliczanie splotu ko
ł
ł
owego metod
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
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
44
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
,
/
/
π
π
,
]
[
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
45
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
,
]
[
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]
[
46
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
+
−
=
−
−
−
−
−
−
=
=
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
47
Iloczyn widm
3
0
≤
≤
=
k
k
H
k
G
k
Y
C
],
[
]
[
]
[
wynosi
−
=
=
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
48
=
]
[
]
[
]
[
]
[
*
]
[
]
[
]
[
]
[
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
=
−
−
−
−
−
−
−
=
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
49
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
−
≤
≤
−
≤
≤
=
1
0
1
0
L
n
N
N
n
n
g
n
g
e
,
],
[
]
[
−
≤
≤
−
≤
≤
=
1
0
1
0
L
n
M
M
n
n
h
n
h
e
,
],
[
]
[
50
wtedy
]
[
]
[
]
[
]
[
]
[
]
[
n
h
n
g
n
y
n
h
n
g
n
y
C
L
=
=
=
*
L
Zero-padding
with
zeros
1)
( −
N
point DFT
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
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.
e
e
51
Splot liniowy metod
Splot liniowy metod
ą
ą
DFT
DFT
–
–
przyk
przyk
ł
ł
ad
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.
-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
52
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)
Y
e
(k)
53
Splot liniowy (równy kołowemu) wyznaczamy jako
IDFT [X
e
(k)Y
e
(k)]
Wartości próbek otrzymanego wektora odpowiadają wartościom splotu liniowego x y.
∗∗∗∗
=
54
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)]
R
yx
(m) = E[y(n+m) x(n)] = E[y(n) x(n-m)]
korelacja sygnału y(n) z x(n)
(13a)
(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).
55
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
(14)
1
N
__
Σ
n=0
N-m-1
x(n) y(n + m)
dla m
≥
0 ( y przesuwane w lewo)
__
1
N
Σ
N-1
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
=
∧
∧∧
∧
56
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.
Porównując (15) ze wzorem na splot liniowy (16)
s
L
(m) = y(m)
∗
x(m) =
Σ
n=0
N-1
y(m) x(m - n)
(16)
C
yx
(m) = y(m)
∗
x(-m)
(17a)
n=0
N-1
Σ
y(n) x(n - m)
otrzymujemy bardzo wygodną obliczeniowo zależność
C
xy
(m) =
(17b)
oraz analogicznie
x(m)
∗
y(-m)
57
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)]
M =2
N =3
L = 4
L = 4
58
Estymator funkcji korelacji metod
Estymator funkcji korelacji metod
ą
ą
DFT
DFT
-
-
przyk
przyk
ł
ł
ad
ad
Natomiast dyskretne transformaty ciągów y
e
(n) i x
e
(n) są równe
Estymaty korelacji (14) wyznaczana w dziedzinie czasu wynosi
X
e
(k)
Y
e
(k)
X*
e
(k)
-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
w lewo
2
y (n+2)
1
2
2
2
c
yx
(m)
przesunięcie
m
ciągi
n
59
Ostatecznie
IDFT [Y
e
(k) X
e
*(k)]
=
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
60
DFT iloczynu ci
DFT iloczynu ci
ą
ą
g
g
ó
ó
w
w
–
–
splot
splot
ko
ko
ł
ł
owy
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
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
61
Wyznaczmy splot kołowy uprzednio wyznaczonych widm G(k) i H(k)
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
62
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
Otrzymany wynik równy jest próbkom iloczynu y(n)=g(n)h(n) z tabeli 1.
63
Podstawowe w
Podstawowe w
ł
ł
asno
asno
ś
ś
ci DFT
ci DFT