N´ ahodn´ e sign´ aly ˇ Jan Cernock´ y ´ UPGM FIT VUT Brno,
[email protected]
1
N´ ahodn´ e sign´ aly • deterministick´e sign´aly (m˚ uˇzeme je zapsat rovnic´ı) maj´ı jednu z´asadn´ı nev´yhodu – nesou velmi m´alo informace (napˇr. kos´ınusovka: amplituda, frekvence, poˇc´ateˇcn´ı f´aze). • sign´aly re´aln´eho svˇeta se daj´ı popsat deterministicky velmi tˇeˇzce nebo v˚ ubec (napˇr. fyzik´aln´ı model pro ˇreˇcov´y sign´al je velmi sloˇzit´y a stejnˇe se jedn´a o zjednoduˇsen´ı). ⇒ na tyto uˇ ziteˇ cn´ e sign´aly se budeme z hlediska teorie d´ıvat jako na n´ ahodn´ e sign´ aly ˇ e poˇsty, kurs Kˇc/EUR . . . ). (procesy) (napˇr. ˇreˇc, cirkulace z´asilek poboˇckami Cesk´ Podle charakteru ˇcasov´e osy dˇel´ıme na n´ahodn´e sign´aly se spojit´ ym ˇ casem (definov´any pro vˇsechna t) a s diskr´ etn´ım ˇ casem (jen pro diskr´etn´ı n). Sign´aly nem˚ uˇzeme popsat ve vˇsech ˇcasech (to by byly deterministick´e), budeme sp´ıˇse hledat charakteristick´e vlastnosti n´ahodn´ych sign´al˚ u, jako stˇredn´ı hodnota, funkce hustoty rozloˇzen´ı pravdˇepodobnosti, atd. 2
Definice n´ ahodn´ eho procesu • spojit´y ˇcas: syst´em {ξt } n´ahodn´ych veliˇcin definovan´ych pro vˇsechna t ∈ < se naz´yv´a n´ahodn´y proces, oznaˇcujeme ξ(t). • diskr´etn´ı ˇcas: syst´em {ξn } n´ahodn´ych veliˇcin definovan´ych pro vˇsechna n ∈ N se naz´yv´a n´ahodn´y proces, oznaˇcujeme ξ[n]. Mnoˇ zina realizac´ı n´ ahodn´ eho procesu moˇznou reprezentac´ı n´ahodn´eho procesu je nekoneˇcnˇe mnoho jeho r˚ uzn´ych pr˚ ubˇeh˚ u– realizac´ı. Omez´ıme se na koneˇcn´y poˇcet Ω a kaˇzdou realizaci oznaˇc´ıme ξω (t), pˇr´ıpadnˇe ξω [n]. Pokud budeme na souboru realizac´ı n´ahodn´eho procesu odhadovat nˇejak´e jeho parametry, bude se jednat o souborov´ e odhady. Pˇr´ıklad: n´ahodn´y proces je zvukov´y sign´al vody tekouc´ı vodovodn´ı trubkou doma u ˇ Cernock´ ych. Bylo nahr´ano 1068 realizac´ı po 20 ms. Pro v´yklad spojit´ych n´ahodn´ych proces˚ u si je budeme pˇredstavovat jako spojit´e sign´aly ξω (t), pro v´yklad diskr´etn´ıch n´ahodn´ych proces˚ u jako ξω [n]. 3
ξω (t) pro ω = 1, 200, 500, 1000 0.2 0 −0.2 0
0.002
0.004
0.006
0.008
0.01
0.012
0.014
0.016
0.018
0
0.002
0.004
0.006
0.008
0.01
0.012
0.014
0.016
0.018
0
0.002
0.004
0.006
0.008
0.01
0.012
0.014
0.016
0.018
0
0.002
0.004
0.006
0.008
0.01
0.012
0.014
0.016
0.018
0.2 0 −0.2
0.2 0 −0.2
0.2 0 −0.2
4
ξω [n] pro ω = 1, 200, 500, 1000 0.2 0 −0.2 0
50
100
150
200
250
300
0
50
100
150
200
250
300
0
50
100
150
200
250
300
0
50
100
150
200
250
300
0.2 0 −0.2
0.2 0 −0.2
0.2 0 −0.2
5
Distribuˇ cn´ı funkce je definov´ana pro jednu n´ahodnou veliˇcinu: n´ahodn´y proces pro urˇcit´y ˇcas t nebo n je takovou n´ahodnou veliˇcinou. Definice: F (x, t) = P{ξ(t) < x}, F (x, n) = P{ξ[n] < x}, kde P{ξ(t) < x} nebo P{ξ[n] < x} je pravdˇepodobnost toho, ˇze n´ahodn´a promˇenn´a zde nabude hodnoty menˇs´ı neˇz x. Uvˇedomme si pros´ım, ˇze x nen´ı nic n´ahodn´eho, je to pomocn´a promˇenn´a, kterou nasad´ıme na nˇejakou hodnotu a pro tuto hodnotu sledujeme pravdˇepodobnost. Souborov´ y odhad distribuˇ cn´ı funkce: posad´ıme se do urˇcit´eho ˇcasu t nebo n, vezmeme Ω realizac´ı, kter´e m´ame k disposici a pro x odhadujeme: PΩ ω=1 1 pokud ξω (t) < x, 0 jinak ˆ F (x, t) = Ω PΩ ω=1 1 pokud ξω [n] < x, 0 jinak ˆ F (x, n) = Ω 6
1 F(x,0.1ms) F(x,3.1ms) F(x,6.3ms) F(x,9.4ms)
0.9
0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
−0.3
−0.2
−0.1
0
0.1 x
7
0.2
0.3
0.4
0.5
1 F(x,1) F(x,50) F(x,100) F(x,150)
0.9
0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
−0.3
−0.2
−0.1
0
0.1 x
8
0.2
0.3
0.4
0.5
Funkce hustoty rozdˇ elen´ı pravdˇ epodobnosti je opˇet definov´ana pro jednu n´ahodnou veliˇcinu (n´ahodn´y proces pro urˇcit´y ˇcas t nebo n je takovou n´ahodnou veliˇcinou). Definice: p(x, t) =
δF (x, t) δx
p(x, n) =
δF (x, n) δx
Souborov´ y odhad funkce hustoty rozdˇ elen´ı pravdˇ epodobnosti: Funkci m˚ uˇzeme z´ıskat numerick´ym derivov´an´ım z odhadnut´e Fˆ (x, t) nebo Fˆ (x, n) nebo ji odhadnout pomoc´ı histogramu: • posad´ıme se do urˇcit´eho t nebo n • Mus´ıme si zvolit L hodnot x od xmin do xmax , nejl´epe s pravideln´ym krokem −xmin : ∆ = xmaxL−1 x1 = xmin , x2 = xmin + ∆, x3 = xmin + 2∆ . . . . . . xL−1 = xmin + (L − 2)∆, xL = xmin + (L − 1)∆ = xmax 9
∆ dostaneme tak L “chl´ıvk˚ u” o ˇs´ıˇrce ∆, pro xi je dan´y chl´ıvek chi od xi − ∆ 2 do xi + 2 . Lev´y okraj spodn´ıho chl´ıvku (1) nat´ahneme aˇz do −∞, prav´y okraj horn´ıho (L) do +∞.
• odhad pro xi poˇc´ıt´ame jako pˆ(xi , t) =
PΩ
1 pokud ξω (t) ∈ chi , 0 jinak Ω∆
pˆ(xi , n) =
PΩ
1 pokud ξω [n] ∈ chi , 0 jinak Ω∆
ω=1
ω=1
10
3.5
3 2.5
2.5
p(x,3.1ms)
p(x,0.1ms)
3
2 1.5
1.5 1
1
0.5
0.5 0
2
−0.2
0
0.2
0.4
0
0.6
−0.2
0
x
0.2
0.4
0.6
0.2
0.4
0.6
x
3.5 3
3
p(x,9.4ms)
p(x,6.3ms)
2.5 2 1.5
2.5 2 1.5
1
1
0.5
0.5
0
−0.2
0
0.2
0.4
0
0.6
x
−0.2
0 x
11
3.5
3
3
2.5
p(x,50)
p(x,1)
2.5 2 1.5
1.5 1
1
0.5
0.5 0
2
−0.2
0
0.2
0.4
0
0.6
−0.2
0
x
0.2
0.4
0.6
0.2
0.4
0.6
x
3.5 3
3 2.5
p(x,150)
p(x,100)
2.5 2 1.5
2 1.5
1
1
0.5
0.5
0
−0.2
0
0.2
0.4
0
0.6
x
−0.2
0 x
12
F (x, t), p(x, t) a pravdˇ epodobnosti pokud m´ame za u ´kol vypoˇc´ıtat pravdˇepodobnost toho, ˇze se hodnota procesu v ˇcase t nebo n vyskytne v intervalu [a, b], m´ame tyto moˇznosti (budeme ukazovat jen na spojit´em ˇcase, vzoreˇcky pro diskr´etn´ı jsou naprosto stejn´e): • vypoˇc´ıtat tuto pravdˇepodobnost z distribuˇcn´ı funkce, pokud si uvˇedom´ıme jej´ı definici: F (x, t) = P{ξ(t) < x}, pak P{a < ξ(t) < b} = F (b, t) − F (a, t) • vypoˇc´ıtat ji z funkce hustoty rozdˇelen´ı pravdˇepodobnosti. Uˇz jsme si zvykli, ˇze pokud je nˇekde hustota, bude se integrovat: Z b p(x, t)dx P{a < ξ(t) < b} = a
13
Z tˇechto vzoreˇck˚ u vypl´yvaj´ı d˚ uleˇzit´e vlastnosti distribuˇcn´ı funkce a funkce hustoty rozdˇelen´ı pravdˇepodobnosti: • hodnoty n´ahodn´eho procesu budou tˇeˇzko menˇs´ı neˇz −∞, proto F (−∞, t) = P{ξ(t) < −∞} = 0. • hodnoty n´ahodn´eho procesu zˇrejmˇe vˇsechny menˇs´ı neˇz ∞, proto F (+∞, t) = P{ξ(t) < +∞} = 1. • funkce hustoty rozdˇelen´ı pravdˇepodobnosti je d´ana jako derivace distribuˇcn´ı funkce, naopak plat´ı integr´al: Z x
F (x, t) =
14
p(g, t)dg
−∞
jelikoˇz F (+∞, t) = 1, mus´ı b´yt Z
+∞
p(x, t)dx = 1 −∞
• hodnota funkce hustoty rozdˇelen´ı pravdˇepodobnosti pro urˇcit´e x nen´ı pravdˇ epodobnost !!! (ˇcast´y chyt´ak statistik˚ u).
15
Ilustrace F (x, t), p(x, t) – beˇ cka piva. . . na kolej´ıch se od x = 18.00 do 22.00 pije beˇcka piva. Definujeme funkci p(x) jako okamˇzitou spotˇrebu piva (pic´ı funkce) a F (x) jako funkci vypit´eho piva (x je v tomto pˇr´ıkladu v´yjimeˇcnˇe ˇcas):
Chov´an´ı je podobn´e funkci hustoty rozdˇelen´ı pravdˇepodobnosti a distribuˇcn´ı funkci, pouze za 1 dosad’te “1 beˇcka”: • F (x) je nulov´a v ˇcase −∞ (je nulov´a dokonce aˇz do 18.00), protoˇze pivo nebylo. 16
• F (x) je 1 beˇcka v ˇcase +∞ (dokonce uˇz ve 22.00), protoˇze pivo je vypit´e a v´ıc ho nebude. Rx • mnoˇzstv´ı vypit´eho piva v ˇcase x: F (x) = −∞ p(g)dg. R +∞ • celkov´e mnoˇzstv´ı vypit´eho piva: F (+∞) = −∞ p(x)dx = 1 beˇcka.
u • mnoˇzstv´ı vypit´eho pivo od ˇcasu x1 do ˇcasu x2 je moˇzn´e spoˇc´ıtat jako rozd´ıl dvou bod˚ F (x) nebo integrac´ı p(x). • hodnotˇe p(x) (napˇr. p(19.00))nem˚ uˇzeme ˇr´ıkat mnoˇzstv´ı piva – za nekoneˇcnˇe kr´atk´y ˇcasov´y interval ho do nikoho nevtekla ani kapka.
17
Momenty na rozd´ıl od funkc´ı jsou momenty ˇ c´ısla, kter´a charakterizuj´ı n´ahodn´y proces v dan´em ˇcase t nebo n: stˇredn´ı hodnota nebo t´eˇz “Expectation” (oˇcek´av´an´ı), prvn´ı moment: Z +∞ xp(x, t)dx a(t) = E{ξ(t)} = −∞
a[n] = E{ξ[n]} =
Z
+∞
xp(x, n)dx −∞
souborov´ y odhad stˇredn´ı hodnoty je pro kaˇzd´y ˇcas t nebo n d´an jako pr˚ umˇer vzork˚ u pˇres vˇsechny realizace: Ω 1 X ξω (t) a ˆ(t) = Ω ω=1 Ω 1 X ξω [n] a ˆ[n] = Ω ω=1 18
Pro naˇse sign´aly: −3
x 10 5 0 −5 −10 −15 0
0.002
0.004
0.006
0.008
0.01 t
0.012
0.014
0.016
0.018
−3
x 10 5 0 −5 −10 −15 0
50
100
150
200 n
19
250
300
rozptyl (disperze), smˇ erodatn´ a odchylka D(t) = E{[ξ(t) − a(t)]2 } = 2
D[n] = E{[ξ[n] − a[n]] } =
Z Z
+∞
[x − a(t)]2 p(x, t)dx −∞ +∞
[x − a[n]]2 p(x, n)dx −∞
smˇerodatn´a odchylka (std - standard deviation) je odmocninou rozptylu: p p σ[n] = D[n] σ(t) = D(t) souborov´ y odhad rozptylu a smˇ erodatn´ e odchylky je pro kaˇzd´y ˇcas t nebo n d´an: Ω X 1 ˆ [ξω (t) − a ˆ(t)]2 , D(t) = Ω ω=1 Ω X 1 ˆ [ξω [n] − a ˆ[n]]2 , D[n] = Ω ω=1 20
q ˆ σ ˆ (t) = D(t) σ ˆ [n] =
q
ˆ D[n]
Pro naˇse sign´aly: 0.14
0.135
0.13
0.125 0
0.002
0.004
0.006
0.008
0.01 t
0.012
0.014
0.016
0.018
0.14
0.135
0.13
0.125 0
50
100
150
200 n
21
250
300
Korelaˇ cn´ı funkce Ud´av´a podobnost mezi hodnotami n´ahodn´eho procesu v ˇcasech t1 (nebo n1 ) a t2 (nebo n2 ): Z +∞ Z +∞ x1 x2 p(x1 , x2 , t1 , t2 )dx1 dx2 , R(t1 , t2 ) = −∞
R(n1 , n2 ) =
Z
+∞ −∞
−∞
Z
+∞
x1 x2 p(x1 , x2 , n1 , n2 )dx1 dx2 , −∞
kde p(x1 , x2 , t1 , t2 ), resp. p(x1 , x2 , n1 , n2 ) je dvourozmˇern´a funkce hustoty rozdˇelen´ı pravdˇepodobnosti mezi ˇcasy t1 a t2 , resp. n1 a n2 . Teoreticky ji vypoˇc´ıt´ame z dvourozmˇern´e distribuˇcn´ı funkce: F (x1 , x2 , t1 , t2 ) = P{ξ(t1 ) < x1 a ξ(t2 ) < x2 }, F (x1 , x2 , n1 , n2 ) = P{ξ[n1 ] < x1 a ξ[n2 ] < x2 }
22
pomoc´ı derivace podle x1 a x2 : δ 2 F (x1 , x2 , t1 , t2 ) p(x1 , x2 , t1 , t2 ) = δx1 δx2 δ 2 F (x1 , x2 , n1 , n2 ) p(x1 , x2 , n1 , n2 ) = δx1 δx2 N´as bude ale sp´ıˇse zaj´ımat souborov´y odhad, kter´y zaˇr´ıd´ıme pomoc´ı dvourozmˇern´eho histogramu: • podobnˇe jako u standardn´ıho histogramu vytvoˇr´ıme “chl´ıvky”, tentokr´at ale budou ∆ do x + y dvourozmˇern´e (ˇctvereˇcky): chij je pro prvn´ı rozmˇer od xi − ∆ i 2 2 a pro druh´ ∆ rozmˇer od xj − ∆ 2 do xj + 2 .
23
24
• hodnota 2D histogramu pro chl´ıvek chij urˇcen´y hodnotami xi a xj bude: PΩ ω=1 1 pokud ξω (t1 ), ξω (t2 ) ∈ chij , 0 jinak pˆ(xi , xj , t1 , t2 ) = Ω∆2 PΩ 1 pokud ξω [n1 ], ξω [n2 ] ∈ chij , 0 jinak pˆ(xi , xj , n1 , n2 ) = ω=1 Ω∆2
25
Pro naˇse sign´aly (uv´ad´ıme jen diskr´etn´ı ˇcasy, stejnˇe jste uˇz pˇriˇsli na to, ˇze ty se spojit´ym ˇcasem jsou u ´plnˇe stejn´e. . . ): n1 = 0, n2 = 0, 1, 5, 11. 0.4
0.4
100
20
0.2
0.2 80
60
0
x1
x1
15
0 10
40 −0.2
−0.2 5 20
−0.4 −0.4
−0.2
0 x2
0.2
0.4
−0.4 −0.4
0
0.4
−0.2
0 x2
0.2
0.4
0
0.4 14
12
20
0.2
0.2 10 15
x1
x1
8 0
0
6
10
4
−0.2
−0.2 5
2
−0.4 −0.4
−0.2
0 x2
0.2
−0.4 −0.4
0.4
26
−0.2
0 x2
0.2
0.4
Pro n1 = 0, n2 = 0, 1, 5, 11 n´am po numerick´e integraci vyˇsly n´asleduj´ıc´ı autokorelaˇcn´ı koeficienty • R(0, 0) = 0.0188: ve stejn´em bodˇe je s´am sobˇe proces vˇzdy nejv´ıce podobn´y. . . • R(0, 1) = 0.0151: posuneme-li se do vedlejˇs´ıho vzorku, je st´ale jeˇstˇe dost podobn´y ˇcasu n1 = 0. • R(0, 5) = 0.0030: ˇcasy n1 = 0 a n2 = 5 si nejsou v˚ ubec podobn´e. • R(0, 11) = −0.0133: v ˇcasech n1 = 0 a n2 = 11 si je proces podobn´y, ale s opaˇcn´ym znam´enkem ! Je pravdˇepodobn´e, ˇze kdyˇz bude hodnota ξ[n1 ] kladn´a, bude ξ[n2 ] z´aporn´a, a naopak.
27
Korelaˇcn´ı funkce pro n1 = 0, n2 = n1 +k pro k = 0 . . . 40, srovn´an´ı s n1 = 100, n2 = n1 +k pro k = 0 . . . 40: 0.02 0.015
R(0,k)
0.01 0.005 0 −0.005 −0.01 −0.015
0
5
10
15
20 k
25
30
35
40
0
5
10
15
20 k
25
30
35
40
0.02
R(100,100+k)
0.015 0.01 0.005 0 −0.005 −0.01 −0.015
28
´ ´ STACIONARITA NAHODN EHO PROCESU lidovˇe ˇreˇceno, chov´an´ı stacion´arn´ıho n´ahodn´eho procesu se nemˇen´ı v ˇcase. Statistick´e veliˇciny nejsou z´avisl´e na aktu´aln´ım t nebo n. Korelaˇcn´ı funkce nen´ı z´avisl´a na pˇresn´e poloze t1 , t2 nebo n1 , n2 , ale pouze na jejich rozd´ılu: τ = t2 − t1 , k = n2 − n1 . Pro stacion´arn´ı syst´emy se spojit´ym ˇcasem: F (x, t) → F (x)
p(x, t) → p(x)
a(t) → a D(t) → D
σ(t) → σ
p(x1 , x2 , t1 , t2 ) → p(x1 , x2 , τ ) R(t1 , t2 ) → R(τ ) Podobnˇe pro diskr´etn´ı ˇcas: F (x, n) → F (x)
p(x, n) → p(x)
a[n] → a D[n] → D p(x1 , x2 , n1 , n2 ) → p(x1 , x2 , k) 29
σ[n] → σ R(n1 , n2 ) → R(k)
V pˇr´ıkladu s tekouc´ı vodou jsme zˇrejmˇe mˇeli stacion´arn´ı sign´al, protoˇze: • stˇredn´ı hodnota byla pro vˇsechny ˇcasy podobn´a (pokud bychom mˇeli k disposici v´ıce realizac´ı, byla by jeˇstˇe “stejnˇejˇs´ı”). • smˇerodatn´a odchylka tak´e (dtto). • korelaˇcn´ı funkce pro n1 = 0, n2 = n1 + k a pro n1 = 100, n2 = n1 + k vypadala podobnˇe.
30
Stacion´arn´ı vs. nestacion´arn´ı sign´al: 2 1.5 1
k
x (t)
0.5 0
−0.5 −1 −1.5 −2 0
0.5
1
1.5 t
2
2.5
3
0.5
1
1.5 t
2
2.5
3
1 0.8 0.6 0.4
k
x (t)
0.2 0
−0.2 −0.4 −0.6 −0.8 −1 0
31
´ ´ ERGODICITA NAHODN EHO PROCESU nutnost m´ıt k disposici mnoho realizac´ı k jak´emukoliv odhadu je trochu svazuj´ıc´ı. U ergodick´ych n´ahodn´ych proces˚ u m˚ uˇzeme parametry odhadnout z jedin´e realizace.
32
Pˇr´ıklad stacion´arn´ıho a ergodick´eho a stacion´arn´ıho, ale neergodick´eho procesu: 2 1.5 1
k
x (t)
0.5 0
−0.5 −1 −1.5 −2 0
0.5
1
1.5 t
2
2.5
3
0.5
1
1.5 t
2
2.5
3
2 1.5 1
k
x (t)
0.5 0
−0.5 −1 −1.5 −2 0
33
Vˇsechny souborov´ e odhady m˚ uˇzeme nahradit ˇ casov´ ymi odhady, m´ame k disposici interval o d´elce T (pro spojit´e procesy) pˇr´ıpadnˇe o poˇctu vzork˚ u N (pro diskr´etn´ı). Jedinou realizaci, kterou m´ame k disposici, nazveme klasicky x(t), resp. x[n]: • pomoc´ı histogram˚ u m˚ uˇzeme stejn´ym zp˚ usobem jako u souborov´ych odhad˚ u odhadnout distribuˇcn´ı funkci a funkci hustoty rozdˇelen´ı pravdˇepodobnosti. • stˇredn´ı hodnota, rozptyl, smˇerodatn´a odchylka: Z T Z T 1 ˆ = 1 x(t)dt D [x(t) − a ˆ]2 dt a ˆ= T 0 T 0 N −1 1 X x[n] a ˆ= N n=0
N −1 X 1 ˆ = D [x[n] − a ˆ ]2 N n=0
• korelaˇcn´ı funkce 1 ˆ R(τ ) = T
Z
T
x(t)x(t + τ )dt 0
N −1 X 1 ˆ x[n]x[n + k] R(k) = N n=0 34
σ ˆ=
σ ˆ=
p
p
ˆ D
ˆ D
N´ ahodn´ e sign´ aly II. ˇ Jan Cernock´ y ´ UPGM FIT VUT Brno,
[email protected]
1
ˇ Casov´ y odhad autokorelaˇ cn´ıch koeficient˚ u pro ergodick´y n´ahodn´y proces s diskr´etn´ım ˇcasem. •
N −1 X 1 ˆ = x[n]x[n + k], R[k] N n=0
kde N je poˇcet vzork˚ u, kter´e m´ame k disposici, se naz´yv´a vych´ ylen´ y odhad (biased estimation). Kdyˇz totiˇz odsouv´ame sign´aly od sebe, odhadujeme R(k) pouze z N − k vzork˚ u. Dˇel´ıme vˇsak st´ale N , takˇze se hodnoty ke kraj˚ um budou sniˇzovat. •
N −1 X 1 ˆ = R[k] x[n]x[n + k], N − |k| n=0
je nevych´ ylen´ y odhad, kdy se dˇel´ı skuteˇcnˇe pouˇzit´ym poˇctem vzork˚ u. Krajn´ı koeficienty k → N − 1, k → −N + 1 jsou ale zat´ıˇzeny znaˇcnou chybou (protoˇze se k jejich odhadu pouˇz´ıv´a m´alo vzork˚ u!), proto d´av´ame pˇrednost vych´ylen´emu odhadu. 2
−3
10
x 10
8
6
4
2
0
−2
−4
−6 −300
−200
−100
0
100
200
300
−200
−100
0
100
200
300
−3
x 10 14 12 10 8 6 4 2 0 −2 −4 −6
−300
3
Spektr´ aln´ı hustota v´ ykonu – spojit´ yˇ cas I u n´ahodn´ych proces˚ u n´as bude zaj´ımat chov´an´ı ve frekvenˇcn´ı oblasti, ale: ˇ protoˇze n´ahodn´e sign´aly nejsou periodick´e. • nem˚ uˇzeme pouˇz´ıt FR,
• nem˚ uˇzeme pouˇz´ıt ani FT, protoˇze n´ahodn´e sign´aly maj´ı nekoneˇcnou energii (FT na tyto sign´aly ˇsla aplikovat, ale pouze ve speci´aln´ıch pˇr´ıpadech) Budeme uvaˇzovat pouze ergodick´e n´ahodn´e sign´aly, realizaci x(t). Odvozen´ı spektr´ aln´ı hustoty v´ ykonu (power spectral density – PSD) • definujeme interval o d´elce T , vezmeme pouze u ´sek od −T /2 do T /2: x(t) pro |t| < T /2 xT (t) = 0 jinde 4
x
-0.5T
0
+0.5T
• definujeme Fourier˚ uv obraz: XT (jω) =
Z
+∞
xT (t)e−jωt dt
−∞
• definujeme spektr´aln´ı hustotu energie (viz pˇredn´aˇska o FT): Z ∞ Z ∞ Z ∞ 1 LT (jω)dω XT (jω)XT (−jω)dω = x2T (t)dt = . . . = 2π −∞ −∞ −∞ LT (ω) naz´yv´ame (dvoustrann´ a) spektr´ aln´ı hustota energie 5
t
|X(jω)|2 LT (jω) = 2π • pokus´ıme se “nat´ahnout” T aˇz do ∞. V tomto pˇr´ıpadˇe by ale energie (i jej´ı hustota) rostly nade vˇsechny meze. Definujeme tedy (dvoustrannou) spektr´ aln´ı hustota v´ ykonu dˇelen´ım T (analogie P = E/T ): GT (jω) =
LT (jω) T
• nyn´ı uˇz m˚ uˇzeme “nataˇzen´ı” prov´est a spoˇc´ıtat spektr´aln´ı hustotu v´ykonu nejen pro u ´sek d´elky T , ale pro cel´y sign´al: |XT (jω)|2 G(jω) = lim GT (jω) = lim T →∞ T →∞ 2πT Vlastnosti G(jω) • Pro spektr´aln´ı funkci re´aln´eho sign´alu XT (jω) plat´ı: XT (jω) = XT⋆ (−jω) G(jω) je prakticky d´ana jej´ımi hodnotami na druhou, bude tedy ˇ cistˇ e re´ aln´ a a sud´ a. 6
• V´ykon n´ahodn´eho procesu v intervalu [ω1 , ω2 ] m˚ uˇzeme spoˇc´ıtat jako: Z ω2 Z −ω1 Z ω2 G(jω)dω. G(jω)dω = 2 G(jω)dω + P[ω1 , ω2 ] = ω1
−ω2
ω1
• Pro celkov´ y stˇredn´ı v´ykon n´ahodn´eho sign´alu plat´ı: Z +∞ G(ω)dω P = −∞
ˇ • Casto je ve sdˇelovac´ı technice a = 0. Pak je stˇredn´ı v´ykon roven rozptylu: P = D a efektivn´ı hodnota Xef = σ.
7
Wiener-Chinchinovy vztahy Spektr´aln´ı hustota v´ykonu je v´az´ana FT s autokorelaˇcn´ı funkc´ı R(τ ) (nˇekdy se tak dokonce definuje – je to jednoduˇsˇs´ı neˇz limitn´ı pˇrechod T → ∞): Z +∞ 1 G(jω) = R(τ )e−jωτ dτ 2π −∞ R(τ ) =
Z
+∞
G(jω)e+jωτ dω
−∞
8
Spektr´ aln´ı hustota v´ ykonu – diskr´ etn´ı ˇ cas PSD n´ahodn´eho procesu s diskr´etn´ım ˇcasem budeme definovat pˇr´ımo pomoc´ı autokorelaˇcn´ıch koeficient˚ u (vˇsimnˇete si terminologie: autokorelaˇcn´ı funkce R(τ ) pro spojit´y ˇcas, autokorelaˇcn´ı koeficienty R[k] pro diskr´etn´ı ˇcas): G(ejω ) =
∞ X
R[k]e−jωk
k=−∞
(jak´a je zde ω kruhov´a frekvence ?). G(ejω ) je Fourierovou transformac´ı s diskr´etn´ım ˇcasem (DTFT) autokorelaˇcn´ıch koeficient˚ u. Pokud tyto odhadneme (souborov´y odhad, ˇcasov´y odhad u ergodick´ych), m˚ uˇzeme G(ejω ) spoˇc´ıtat. Zpˇetn´y pˇrechod od G(ejω ) k autokorelaˇcn´ım koeficient˚ um (zpˇetn´a DTFT): Z π 1 G(ejω )e+jωk dω R[k] = 2π −π 9
Pro tekouc´ı vodu (R[k] odhadov´any z jedn´e realizace): normalized omega 0.5 0.4 0.3 0.2 0.1
−15
−10
−5
0
5
10
15
f 0.5 0.4 0.3 0.2 0.1
−4
−3
−2
−1
0
1
2
3
4 4
x 10
10
zoom od −F s/2 do F s/2:
normalized omega
0.7 0.6 0.5 0.4 0.3 0.2 0.1 −3
−2
−1
0
1
2
3
f 0.7 0.6 0.5 0.4 0.3 0.2 0.1 −8000
−6000
−4000
−2000
0
2000
4000
6000
⇒ tekouc´ı voda m´a silnou hodnˇe v´ykonu okolo 700 Hz, ˇze by rezonance trubky ? 11
8000
vlastnosti G(ejω ) jsou opˇet d´any standardn´ımi vlastnostmi obrazu DTFT: • autokorelaˇcn´ı koeficienty jsou re´aln´e, proto bude pro G(ejω ) platit G(ejω ) = G⋆ (e−jω ), • autokorelaˇcn´ı koeficienty jsou symetrick´e (sud´e), proto bude G(ejω ) vˇsude re´aln´a. • z toho vypl´yv´a, ˇze G(ejω ) bude re´aln´a a sud´a, podobnˇe jako G(jω). • je periodick´a (s periodou 2π, 1, Fs , 2πFs podle frekvence, kterou si vyberete) protoˇze sign´al je diskr´etn´ı. • V´ykon n´ahodn´eho procesu v intervalu [ω1 , ω2 ] m˚ uˇzeme spoˇc´ıtat jako: Z ω2 Z −ω1 Z ω2 1 1 1 G(ejω )dω + G(ejω )dω = G(ejω )dω. P[ω1 , ω2 ] = 2π ω1 2π −ω2 π ω1 • Pro celkov´ y stˇredn´ı v´ykon n´ahodn´eho sign´alu plat´ı: Z +π 1 P = G(ejω )dω 2π −π 12
to je ale hodnota zpˇetn´e DTFT pro k = 0: Z π Z +π 1 1 jω +jω0 R[0] = G(e )e dω = R[0] = G(ejω )dω 2π −π 2π −π takˇze jsme dostali vztah: R[0] = P
13
Odhad spektr´ aln´ı hustoty v´ ykonu G(ejω ) pomoc´ı DFT M´ame-li k disposici realizaci n´ahodn´eho procesu x[n] o N vzorc´ıch, m˚ uˇzeme PSD odhadnout pomoc´ı Diskr´etn´ı Fourierovy transformace (DFT – to je ta, co se jako jedin´a d´a sluˇsnˇe spoˇc´ıtat), pro pˇripomenut´ı: X[k] =
N −1 X
−j 2π N kn
x[n]e
n=0
G(ejω ) dostaneme pouze pro diskr´etn´ı frekvence: ωk =
2π N k:
1 jωk ˆ G(e ) = |X[k]|2 . N Tento odhad b´yv´a nˇekdy velice nespolehliv´y (zaˇsumˇen´y), proto se ˇcasto vyuˇz´ıv´a Welchova metoda – pr˚ umˇerov´an´ı odhadu PSD pˇres nˇekolik ˇcasov´ych u ´sek˚ u: • Sign´al rozdˇel´ıme na M u ´sek˚ u po N prvc´ıch a pro kaˇzd´y spoˇc´ıt´ame DFT: Xm [k] =
N −1 X
−j 2π N kn
xm [n]e
n=0
14
pro 0 ≤ m ≤ M − 1
• Odhad PSD spoˇc´ıt´ame: M −1 X 1 1 ˆ W (ejωk ) = |Xm [k]|2 . G M m=0 N
15
Uk´azka pro odhad PSD z jednoho 320-vzorkov´eho segmentu a z 1068 takov´ych segment˚ u: one realization 1.4 1.2 1 0.8 0.6 0.4 0.2 0
0
0.5
1
1.5
2
2.5
3
3.5
2.5
3
3.5
averaging − Welch 1 0.8 0.6 0.4 0.2 0
0
0.5
1
1.5
2
⇒ odhad z´ıskan´y pr˚ umˇerov´an´ım je mnohem hladˇs´ı. 16
Pr˚ uchod n´ ahodn´ ych sign´ al˚ u line´ arn´ımi syst´ emy • spojit´ yˇ cas: line´arn´ı syst´em m´a komplexn´ı kmitoˇctovou charakteristiku H(jω). Pro vstupn´ı sign´al se spektr´aln´ı hustotou v´ykonu Gx (jω) je v´ystupn´ı spektr´aln´ı hustota v´ykonu d´ana: Gy (jω) = |H(jω)|2 Gx (jω) • diskr´ etn´ı ˇ cas: line´arn´ı syst´em m´a komplexn´ı kmitoˇctovou charakteristiku H(ejω ). Pro vstupn´ı sign´al se spektr´aln´ı hustotou v´ykonu Gx (ejω ) je v´ystupn´ı spektr´aln´ı hustota v´ykonu d´ana: Gy (ejω ) = |H(ejω )|2 Gx (ejω ) V obou pˇr´ıpadech n´asob´ıme vstupn´ı PSD druhou mocninou modulu komplexn´ı kmitoˇctov´e charakteristiky. Tvar funkce hustoty rozdˇelen´ı pravdˇepodobnosti se pr˚ uchodem line´arn´ım syst´emem nemˇ en´ı – mˇen´ı se jen parametry. 17
Pˇr´ıklad: filtrov´an´ı jedn´e realizace teˇcen´ı vody filtrem H(z) = 1 − 0.9z −1 . Vstupn´ı sign´al a jeho PSD: 0.2
0
−0.2 50
100
150
200
250
300
0.5 0.4 0.3 0.2 0.1
0
0.5
1
1.5
18
2
2.5
3
Modul komplexn´ı kmitoˇctov´e charakteristiky a jeho druh´a mocnina: |H|
1.5
1
0.5
0
0.5
1
1.5
2
2.5
3
2
2.5
3
2
|H| 3.5 3 2.5 2 1.5 1 0.5 0
0.5
1
1.5
19
V´ystupn´ı sign´al a jeho PSD: 0.2
0
−0.2
50
100
150
200
250
300
0.1 0.08 0.06 0.04 0.02
0
0.5
1
1.5
20
2
2.5
3
Pˇr´ıklad n´ ahodn´ eho procesu – Gaussovsk´ y b´ıl´ y ˇsum V Matlabu funkce randn – jednotliv´e vzorky na sobˇe nez´avis´ı, funkce hustoty rozdˇelen´ı pravdˇepodobnosti je d´ana Gaussov´ym rozloˇzen´ım: [x−µ]2 1 − 2σ2 p(x) = N (x; µ, σ) = √ e σ 2π
µ=0, σ=1
µ=0, σ=1 1
0.9
0.4
0.8
0.7 0.3
p(x)
F(x)
0.6
0.5
0.2 0.4
0.3 0.1
0.2
0.1
0 −5
−4
−3
−2
−1
0 x
1
2
3
4
5
0 −5
−4
−3
Generov´an´ı v Matlabu: x = sigma*randn(1,N) + mu 21
−2
−1
0 x
1
2
3
4
5
3
2
1
0
−1
−2
−3
0
20
40
60
80
100
120
140
Autokorelaˇcn´ı koeficienty: R[0] = µ2 + D,
R[k] = µ2 pro k 6= 0
22
160
180
200
1
0.8
0.6
0.4
0.2
0
−150
−100
−50
0
50
100
150
Spektr´aln´ı hustota v´ykonu je u b´ıl´eho ˇsumu konstantn´ı (proto b´ıl´y): G(ejω ) = R[0]
23
5 4 3 2 1
0
0.5
1
1.5
2
2.5
3
2
2.5
3
averaging − Welch 1.06 1.04 1.02 1 0.98 0.96 0.94 0.92 0
0.5
1
1.5
Pro spojit´y ˇcas nejde ˇcistˇe b´ıl´y ˇsum vygenerovat: pokud by G(jω) byla nenulov´a pro vˇsechny ω, mˇel by nekoneˇcn´y v´ykon. . .
24
Kvantov´ an´ı Representace vzork˚ u diskr´etn´ıho sign´alu x[n] nen´ı moˇzn´a s libovolnou pˇresnost´ı ⇒ kvantov´ an´ı. Nejˇcastˇeji zaokrouhlujeme na fixn´ı poˇcet L kvantovac´ıch hladin, kter´e jsou oˇc´ıslovan´e od 0 do L − 1. Pokud m´ame na kvantov´an´ı k disposici b bit˚ u, L = 2b . Uniformn´ı kvantov´ an´ı m´a rovnomˇern´e rozloˇzen´ı kvantovac´ıch hladin q0 . . . qL−1 od minim´aln´ı hodnoty sign´alu xmin do maxim´aln´ı hodnoty xmax :
Kvantovac´ı krok ∆ je d´an:
xmax − xmin ∆= , L−1 25
pro velk´a L m˚ uˇzeme pouˇz´ıt pˇribliˇzn´y vztah: ∆=
xmax − xmin . L
Kvantov´ an´ı: pro hodnotu x[n] je index nejlepˇs´ı kvantovac´ı hladiny d´an: i[n] = arg
min
l=0...L−1
|x[n] − ql |,
a kvantovan´y sign´al je: xq [n] = qi[n] . Chyba kvantov´ an´ı: e[n] = x[n] − xq [n]. m˚ uˇze b´yt tak´e povaˇzov´ana za sign´ al.
26
Ilustrace na kvantov´an´ı kosinusovky x[n] = 4 cos(0.1n), L = 8: 4 2 0 −2 20
40
60
80
100
120
140
160
180
20
40
60
80
100
120
140
160
180
20
40
60
80
100
120
140
160
180
4 2 0 −2
0.5
0
−0.5
27
Abychom zjistili, jak je sign´al kvantov´an´ım naruˇsen, bude dobr´e spoˇc´ıtat v´ ykon chybov´eho sign´alu Pe a srovnat jej s v´ykonem uˇziteˇcn´eho sign´alu Ps : pomˇ er sign´ alu k ˇsumu – signal-to-noise ratio (SNR): SN R = 10 log10
Ps Pe
[dB].
Pro v´ypoˇcet v´ykonu chybov´eho sign´alu vyuˇzijeme teorie n´ ahodn´ ych proces˚ u: nezn´ame ∆ ze budou rovnomˇernˇe rozloˇzen´e. hodnoty e[n], ale v´ıme, ˇze budou v intervalu [− ∆ 2, +2] a ˇ Funkce hustoty rozdˇelen´ı pravdˇepodobnosti pro e[n] bude:
28
. . . v´yˇska
1 ∆
vych´az´ı z toho, ˇze plocha: ∆ 2
Z
−∆ 2
!
pe (g)dg = 1.
Tento proces m´a nulovou stˇredn´ı hodnotu (snadno bychom zjistili, ˇze v´ykon bude tedy roven rozptylu: P e = De =
Z
∞
−∞
2
g pe (g)dg =
Z
∆ 2
−∆ 2
·
3
1 g g pe (g)dg = ∆ 3 2
29
¸ ∆2
−∆ 2
1 = 3∆
R µ
∆ 2 −∆ 2
3
gpe (g)dg = 0),
3
∆ ∆ + 8 8
¶
∆2 = 12
V´ ypoˇ cet SNR pro kosinusovku amplituda A, kosinusovka m´a v´ykon Ps = xmin = −A, xmax = A, takˇze 2A ∆= L
A2 2
A2 ∆2 4A2 = . Pe = = 2 2 12 12L 3L
Pomˇer sign´alu k ˇsumu: SN R = 10 log10
Ps = Pe
A2 10 log10 A22 3L2
= 10 log10
3L2 . 2
Pokud m´ame k disposici b bit˚ u a poˇcet kvantovac´ıch hladin L = 2b : 3 3 SN R = 10 log10 (2b )2 = 10 log10 + 10 log10 22b = 1.76 + 20b log10 2 = 1.76 + 6 b dB. 2 2 Konstanta 1.76 z´avis´ı na charakteru sign´alu (cos, ˇsum), ale plat´ı, ˇze pˇrid´ an´ı/ubr´ an´ı jednoho bitu zlepˇsuje/zhorˇsuje SNR o 6 dB. 30
Pˇr´ıklad: kosinusovky s A = 4 kvantovan´e na L = 8 hladin´ach: SN Rteor = 1.76 + 3 × 6 = 19.76 dB
SN Rexp = 10 log10
N −1 1 X 2 s [n] N n=0
1 N
N −1 X
e2 [n]
n=0
Matlab: snr = 10*log10 (sum(x.^2) / sum(e.^2)) . . . celkem to vych´az´ı.
31
= 19.36 dB
