Cislicove_zpracovani_a_analyza_signalu_P

-n), což je reálné jen u filtrů FIR. Fázová
charakteristika takového kauzálního filtru s impulsní charakteristikou {h(0), h(1),…, h(N-1)}
je
() T
N
G ωω
2
1
arg
−
−= , (3.19)
Číslicové zpracování a analýza signálů 17
Filtr zavádí zpoždění (N-1)/2 vzorkovacích intervalů, což odpovídá polovině délky
impulsní charakteristiky.
Z rozložení nulových bodů a pólů lze posuzovat stabilitu filtru. Stabilní systém musí
mít všechny póly uvnitř jednotkové kružnice - rekursivní realizace v příkladu P.3. má sice 1
pól na kružnici, ale ten je eliminován nulovým bodem. Dodejme, že impulsní charakteristika
stabilního systému musí být absolutně sumovatelná (součet absolutních hodnot h(n) musí být
konečné číslo, filtry FIR jsou tedy vždy stabilní).
Z rozložení nulových bodů a pólů lze odhadnout přibližný tvar frekvenční
charakteristiky filtru. To může být užitečné nejen pro orientační odhad frekvenčních
vlastností filtru, ale i při návrhu filtru z předem zvolené konfigurace nulových bodů a pólů
(viz kap.3.4).
Příklad P. 4
Odhadněme tvar frekvenční charakteristiky filtru 2. řádu, který je popsán nulovými
body (n1, n2) a póly (p1, p2), jejichž polohy jsou naznačeny na obr. 3.5.

obr. 3.5. Rozložení nulových bodů a pólů pásmové propusti 2. řádu a naznačení principu odhadu tvaru
|G(ω)|, viz (3.20).

Dosadíme-li v (3.16) za komplexní proměnnou z výraz e
jωT
, znamená to, že vyčíslujeme
komplexní přenos H(z) pro hodnoty komplexní proměnné z ležící na jednotkové
kružnici. Každý úhlový kmitočet ω tedy koresponduje s nějakým bodem na jednotkové
kružnici (pro ω=0 je z=1, pro ω=ω
vz
/4 je z=j, ω=ω
vz
/2 je z=-1, atd.). Pak můžeme
vyjádřit modulovou frekvenční charakteristiku pomocí modulů rozdílů komplexních
výrazů,
()
21
21
21
21
cc
dd
K
pepe
nene
KG
TjTj
TjTj
=
−−
−−
=
ωω
ωω
ω ,
(3.20)
které můžeme vyjádřit jako rozdíly vektorů (na obr. 3.5 je tato vektorová operace
naznačena pro pól p
1
), které jsou představovány vzdálenostmi daného místa na
jednotkové kružnici od nulových bodů (d
1
, d
2
) a pólů (c
1
, c
2
). K (3.20) poznamenejme, že
18 Fakulta elektrotechniky a komunikačních technologií VUT v Brně
konstanta K může být libovolná - může být např. nastavena tak, aby byl maximální
modul přenosu roven jedné.
Leží-li na jednotkové kružnici nulový bod, má filtr na korespondujícím kmitočtu nulový
přenos. Leží-li na jednotkové kružnici pól, má filtr na korespondujícím kmitočtu
nekonečný přenos (filtr je na mezi stability), resp. čím blíže k jednotkové kružnici se pól
nachází, tím větší je zesílení filtru na odpovídající frekvenci.
Ukázka rozložení nulových bodů a pólů dvou pásmových propustí (které se liší jen
polohami nulových bodů) a korespondujících modulových frekvenčních charakteristik je
na obr. 3.6.

obr. 3.6. Rozložení nulových bodů a pólů dvou pásmových propustí a korespondující modulové frekvenční
charakteristiky.
null

Na základě výše uvedeného postupu si může čtenář ověřit, že také tvar modulové
frekvenční charakteristiky na obr. 3.4 je docela pochopitelný.
Na obr. 3.7 jsou graficky (a snad i přehledně) shrnuty jednotlivé způsoby popisu
lineárních diskrétních systémů a naznačeny vzájemné souvislosti.
Číslicové zpracování a analýza signálů 19
diferenční rovnice
- rekursivní
- nerekursivní
přenosová funkce
H(z)
Nulové body a póly
Z- tr.
Z
-1
- tr.
kořeny
čitatele a
jmenovatele
intuitivní návrhy
filtrů
odezva na
δ(
n
)
impulsní
charakteristika
(FIR, IIR)
konvoluce
(FIR,
nerekursivní)
frekvenční
charakteristika
H(ω)
DTFT
návrh FIR
z=e
j
ω
T
Z - tr.,
Z
-1
- tr.
výpočet odezvy
posouzení
stability
návrhy na základě
analogových
předloh
odhad
diskretizovaná
frekvenční
charakteristika
H(kΩ)
ω
= k
Ω
DFT
návrh FIR (DFT
-1
)
výpočet odezvy
ve frekvenční
oblasti
posouzení
stability
vhodná
transformace

obr. 3.7. Souvislosti mezi jednotlivými způsoby popisu lineárních diskrétních systémů.

Na závěr této kapitoly ještě připomeneme základní členění lineárních číslicových filtrů.
Členění filtrů podle realizace:
• Rekursivní filtry – přenosové funkce mají nejméně jeden pól mimo počátek.
• Nerekursivní filtry – přenosové funkce mají všechny póly v počátku. Filtry jsou
vždy stabilní.
Členění filtrů podle impulsní charakteristiky, resp. podle filozofie návrhu:
• Filtry IIR - mají nekonečnou impulsní charakteristiku. Původně vznikly jako
číslicové protějšky analogových filtrů (Butterworthových, Čebyševových,
eliptických). Realizovány jsou rekursivně, mají vždy nelineární fázovou
charakteristiku. (Princip intuitivního návrhu uvádíme v kap.3.4, princip návrhu
metodou bilineární transformace je uveden v kap.3.7).
• Filtry FIR - mají konečnou impulsní charakteristiku. Nemají analogové
protějšky. Realizovány bývají obvykle nerekursivně. Mohou být navrhovány
s lineární fázovou charakteristikou. (Princip intuitivního návrhu uvádíme
20 Fakulta elektrotechniky a komunikačních technologií VUT v Brně
v kap.3.4, principy dalších dvou způsobů návrhu jsou uvedeny v kap.3.5 a v
kap.3.6).
3.3.1 Příklady ke kapitole 3.3
Příklad P. 5
V prostředí Blockset zobrazte přechodovou charakteristiku systému:
85.0
1
)(
2
++
=
zz
zH.
Řešení

obr. 3.8. Jednoduché schéma pro příklad Příklad P. 5

Pro řešení potřebujeme zdroj jednotkového skoku, který najdeme v DSP Blockset/DSP
Sources/Constant. Zde necháme původní nastavení (Constant Value: 1). Dále musíme zadat
diskrétní systém (filtr). Použijeme blok Simulink/Discrete/Discrete Transfer Fcn. Zde zadáme
koeficienty polynomu v čitateli a ve jmenovateli. Pro zobrazení použijeme DSP Blockset/DSP
Sinks/Time Scope. Dvojklikem otevřeme okno pro zobrazení.
V menu Simulation/Simulation parameters/Stop Time nastavíme hodnotu 100.

Příklad P. 6
Zobrazte impulsovou, amplitudovou a fázovou frekvenční charakteristiku systému:
85.0
1
)(
2
++
=
zz
zH
Řešení
Pro řešení využijeme předchozího modelu a vztahu mezi impulsní charakteristikou a
diskretizovanou frekvenční charakteristikou LDS (který je dán diskrétní FT).
Číslicové zpracování a analýza signálů 21

obr. 3.9. Schéma pro zobrazení frekvenční charakteristiky systému v prostředí Blockset

Jako zdroj signálu vezmeme blok jednotkového impulsu DSP Blockset/DSP Sources/Discrete
Impulse. Dostaneme tak na výstupu impulsní charakteristiku systému (jejích prvních několik
vzorků – podle zvolené doby simulace). Pak provedeme DFT impulsní charakteristiky.
K tomu potřebujeme data „poskládat“ do zásobníku („buffer“), ze kterého pak tento vektor
dat použijeme pro výpočet frekvenční charakteristiky (DFT se musí vždy provádět z vektoru –
posloupnosti vzorků signálu - o konečné délce). Proto pro zobrazení spektra použijeme
nejprve blok DSP Blockset/Signal Management/Buffer a parametr Output buffer size
nastavíme na hodnotu 512. Tento blok dat potom přivedeme na vstup bloku, který realizuje
výpočet DFT (DSP Blockset/Transforms/FFT). Výsledek je samozřejmě komplexní a proto
bychom chtěli zobrazit zvlášť amplitudu a zvlášť fázi. K tomu použijeme blok Simulink/Math
Operation/Complex to Magnitude-Angle, který má dva výstupy reprezentující amplitudu a
fázi. Pro jejich zobrazení použijeme dva bloky DSP Blockset/DSP Sinks/Vector Scope. V
jejich parametrech nastavíme Input Domain jako Frequency a Frequency Range jako [0..Fs].
Říkáme tím, že se bude zobrazovat frekvenční charakteristika v rozsahu od nulového kmitočtu
do vzorkovací frekvence.
V menu Simulation/Simulation parameters/Stop Time nastavíme hodnotu 512.

Příklad P. 7
V prostředí Blockset realizujte systém (filtr), který je dán následující diferenční rovnicí
)4(
3
1
)2(
3
1
)(
3
1
)( −+−+= nxnxnxny
a vykreslete jeho impulsovou a amplitudovou frekvenční charakteristiku.

Řešení
Jedná se o nerekurzivní FIR systém, který pracuje pouze s pěti vstupními vzorky, přičemž
počítá průměr z aktuálního vzorku a vzorku zpožděného o dva a čtyři takty(vzorkovací
intervaly). Realizace v prostředí Blockset může být například následující (obr. 3.10).
22 Fakulta elektrotechniky a komunikačních technologií VUT v Brně

obr. 3.10. Realizace jednoduchého nerekurzivního systému pomocí elementárních bloků

Pro vykreslení impulsové charakteristiky a amplitudové frekvenční charakteristiky použijeme
bloky z minulého příkladu. Pro samotné sestavení modelu systému potřebujeme pouze blok
pro zpoždění (delay) signálu DSP Blockset/Signal Operations/Integer Delay a blok pro
násobení (gain) Simulink/Math Operations/Gain.
Pro vykreslení amplitudové frekvenční charakteristiky použijeme stejného přístupu jako v
předchozím případě (ukládání do zásobníku - bufferu). Jeho velikost nastavíme na 128
vzorků. Celkovou dobu simulace pak také na 128.

Příklad P. 8
V prostředí Blockset realizujte systém pomocí zpožďovacích členů (viz. předchozí
příklad):
85.0
1
)(
2
++
=
zz
zH
a vykreslete jeho impulsovou a amplitudovou frekvenční charakteristiku.

3.4 Intuitivní návrhy filtrů z nulových bodů a pólů v rovině "z".
Předem poznamenejme, že komplexní nulové body i póly se musí vyskytovat v
komplexně sdružených párech, jinak by přenosová funkce neobsahovala pouze reálné
konstanty a následkem toho by nebyla reálná ani impulsní charakteristika - na výstupu filtru s
komplexními konstantami by byl komplexní signál.


Návrhy filtrů FIR
Návrhy vychází z tzv. hřebenových filtrů, které jsou charakterizovány N rovnoměrně
rozloženými nulovými body na jednotkové kružnici (a N-násobným pólem v počátku pro
zajištění kauzality). Některé z nulových bodů jsou eliminovány rovnoměrně rozloženými
Číslicové zpracování a analýza signálů 23
póly; kmitočty, které korespondují s těmito póly, jsou středními kmitočty propustných pásem
výsledného filtru.

Hřebenové filtry – základní 2 typy mají přenosové funkce
()
() ...,3,2,1,)1(
2
1
...,3,2,1,)1(
2
1
=−=
=+=
−
−
NzzG
NzzH
N
N

(3.21)
Nejjednoduššími filtry tohoto typu jsou:
- dolní propust, H(z)=(1+z
-1
)/2 , max. přenos H(0)=1, diferenční rovnice y(n)=(x(n)+ x(n-
1))/2,
- horní propust, G(z)=(1-z
-1
)/2 , max. přenos H(f
vz
/2)=1, diferenční rovnice y(n)=(x(n)- x(n-
1))/2.

Lynnovy filtry - eliminace některých nulových bodů hřebenových filtrů póly vede k
filtrům, které mají charakter více či méně úzkopásmových (popř. vícepásmových) propustí 3
základních typů s obecnými přenosovými funkcemi
()
()
()
()
()
()
()
()Klichézzz
KzK
z
zF
Ksudézzz
KzK
z
zG
pKcelézzz
KzK
z
zH
ppKpp
p
pK
ppKpp
p
pK
ppKpp
p
pK
,...1
1
)1(
1
,...1
1
)1(
1
,,...1
1
)1(
1
2
2
2
−−−−
−
−
−−−−
−
−
−−−−
−
−
+−+−=
+
+
=
−−+−=
+
−
=
++++=
−
−
=
(3.22)
Konstanta K ve jmenovateli zajišťuje, aby byl maximální modul přenosu v propustných
pásmech jednotkový. Polynomy přenosových funkcí (3.22) jsou dělitelné (jinak by nešlo o
filtry FIR), výrazy vlevo představují rekursivní realizace, výrazy vpravo realizace
nerekursivní. Počet ekvidistantních propustných pásem udává konstanta p, šířka pásem bude
závislá na hodnotě konstanty K - čím bude vyšší, tím užší budou propustná pásma. Pro p=1
bude mít filtr H charakter dolní propusti, filtry G a F budou horní propusti.
Vyjdeme-li při návrhu z rozložení nulových bodů a pólů rekursivní dolní propusti
na dříve uvedeném obr. 3.4, snadno dospějeme k přenosové funkci (3.13), jen použijeme
opačný postup, než který jsme uvedli v příkladu Příklad P. 3. Z filtrů uvedených v (3.22) se
jedná o filtr H s konstantami K=5 a p=1. Vzorek výstupního signálu v tomto případě získáme
zprůměrněním posledních 5 vzorků vstupního signálu.
Pro všechny filtry v (3.22) je charakteristické poměrně velké zvlnění modulových
frekvenčních charakteristik v nepropustných pásmech (viz obr. 3.4). Výrazného zlepšení
dosáhneme použitím dvojice shodných filtrů v sérii.


Návrhy filtrů IIR
Systémy 1. řádu, např. s diferenční rovnicí typu y(n) = x(n) + Ky(n-1), stabilní pro -1 <
K < 1, nejsou zvlášť zajímavé a čtenář si snadno ověří jejich chování sám.
24 Fakulta elektrotechniky a komunikačních technologií VUT v Brně
Systémy 2. řádu nabízí dva zajímavé typy filtrů – pásmové zádrže a pásmové propusti.
Pásmové zádrže s obecnou přenosovou funkcí
() 10,
))((
))((
)1)(1(
)1)(1(
00
00
00
00
11
11
N, předpokládejme y(n)=0 pro n0 dostaneme shodnou posloupnost s
0
(n).
Číslicové zpracování a analýza signálů 61
V následující tabulce jsou uvedeny vlastnosti spekter a jejich výpočet pro jednotlivé
spojité a diskrétní signály (periodické i neperiodické). Vlastnosti jednotlivých transformací
jsou uvedeny v Dodatku 9.1.
Tab. 2. Vlastnosti spekter signálů
Spojité signály Diskrétní signály
Časová oblast
periodické neperiodické periodické neperiodické
Spektrum diskrétní,
neperiodické
spojité,
neperiodické
diskrétní,
periodické
spojité,
periodické
Výpočet spektra FŘ integrální FT DFT DTFT
Poznámka
spektrem jsou
koeficienty FŘ
používá se pojmu
spektrální hustota
spektrum je tvořeno
koef. DFŘ, které jsou
až na multiplikativní
konstantu shodné
s koef. DFT
prakticky počítáme jen
vzorky spektra pomocí
DFT

Uvažujme diskrétní, nekonečně dlouhý harmonický signál. Spektrum však počítáme
z konečného počtu (N) vzorků. Původně nekonečný signál je vynásoben oknem o N vzorcích.
V ideálním případě úsek periodického diskrétního signálu, který podrobíme DFT, obsahuje
celistvý počet period. Uvažme situaci na obr. 7.1, kde je harmonický signál s periodou 1/5
s vzorkovaný s periodou vzorkování T=1/32 s. Posloupnost vzorků je zřejmě periodická se 32
vzorky v periodě. Máme- li po DFT získat očekávané spektrum pouze se dvěma nenulovými
čarami, na 5 Hz a 27 Hz (tj. f
vz
-5 Hz), musíme transformovat celistvý násobek 32 vzorků
signálu. Jinak se v diskrétním spektru (po DFT) na uvedených kmitočtech nemůže žádná
spektrální čára objevit.


obr. 7.1. Diskrétní harmonický signál (1 s) o kmitočtu 5 Hz (perioda signálu 1/5 s, vzorkovací interval
T=1/32 s) a odpovídající 3 periody modulového a fázového spektra.

Častěji můžeme očekávat situaci, kdy úsek periodického signálu, který vstupuje do DFT,
obsahuje necelý počet period. Ukázka části takového spektra je na obr. 7.2a, obsahuje řadu
nenulových spektrálních čar v okolí kmitočtu 5Hz, na kterém bychom si přáli jedinou čáru.
Spektrální čára přesně na kmitočtu 5 Hz totiž ve spektru zastoupena není, spektrum se
„rozmaže“. Rozmazání spektra je způsobeno useknutím posloupnosti, tj. násobením
obdélníkovým oknem v časové oblasti. Ve frekvenční oblasti totiž dojde ke konvoluci spektra
62 Fakulta elektrotechniky a komunikačních technologií VUT v Brně
signálu (jedna spektrální čára) a spektra obdélníkového okna (funkce typu sin(x)/x).
Ekvidistantní vzorky modulového spektra leží na funkci |sin(x)/x|, která má střed na 5 Hz. Ke
konvoluci spektra harmonického signálu se spektrem obdélníkového okna dojde samozřejmě i
v předchozí situaci (viz bod 1.), vzorky modulového spektra také leží na funkci |sin(x)/x| se
středem na 5 Hz, ale leží právě na tomto středu a na průchodech nulou.
Možnost částečné nápravy spočívá ve volbě vhodnějšího okna, tedy takového, které bude mít
hlavní lalok frekvenční charakteristiky co nejužší a vedlejší laloky co nejnižší. Existují
standardní okna (váhovací funkce), které byly zmíněny u návrhu FIR filtrů v kapitole 3.6 a
které mají příznivější vlastnosti. To je také patrné z obr. 7.2, kde jsou kromě spektra
získaného váhováním obdélníkovým oknem (a), zobrazeny i spektra signálů váhovaných
trojúhelníkovým (b) a Hammingovým (c) oknem.



obr. 7.2. Část spektra harmonického signálu o kmitočtu 5 Hz, f
vz
=500 Hz získaného transformací
z necelého počtu period. (a) za použití obdélníkového okna, (b) trojúhelníkového okna, (c) Hammingova
okna.

7.2 Spektrální analýza stochastických signálů
Při spektrální analýze stacionárních stochastických signálů nám jde o výpočet
výkonového spektra signálu (nejčastěji pomocí DFT) a jeho následnou analýzu. Jedná se o
transformaci signálu z časové oblasti do oblasti kmitočtové, kde je obvykle tato analýza
snadněji proveditelná. Ztrácíme tak informaci o čase a proto tato analýza má smysl pouze pro
procesy stacionární popřípadě i ergodické. Pro výpočet výkonových spekter pomocí DFT
slouží dvě základní metody – metoda korelogramu a metoda periodogramu.
7.2.1 Odhad výkonového spektra metodou periodogramu
Nejsnadněji lze výkonové spektrum určit jako souborovou střední hodnotu kvadrátů
spekter jednotlivých realizací:
Číslicové zpracování a analýza signálů 63
() () () () ()
∑
=
≈
⎭
⎬
⎫
⎩
⎨
⎧
Ε=
⎭
⎬
⎫
⎩
⎨
⎧
Ε=
M
w
ww
wwwwxx
kX
NM
kX
N
kXkX
N
kS
1
22
*
1111
.
(7.3)
Máme tedy změřeno M realizací x
w
, každá o délce N, z nich určíme kvadrát spektra |X
w
(k)|
2
a
ta potom zprůměrujeme. Zde je na místě zdůraznit, že výkonové spektrum je funkce
deterministická, ale fázové spektrum jednotlivé realizace je náhodné a proto nemá smysl jej
průměrovat. Průměrováním tedy snižujeme rozptyl odhadu výkonového spektra.
Jestliže můžeme předpokládat, že analyzovaný signál je ergodický, lze výkonové spektrum
odhadnout pouze z jedné realizace. To provedeme tak, že změřený signál o délce N´
rozdělíme na M úseků o délce N=N´/M. Z každého takového segmentu je pak vypočteno
spektrum a výsledky pak průměrovány dle (7.3). Tímto postupem ovšem dojde ke zmenšení
rozlišovací schopnosti, protože pro výpočet individuálních spekter používáme menší počet
vzorků, než má celý měřený signál.
Na obr. 7.3a je jeden segment signálu, který obsahuje dvě harmonické složky o frekvencích
20 Hz a 40 Hz a náhodný signál se složkami od 0 do cca 150 Hz. Vzorkovací frekvence je
500 Hz. Spektrum tohoto jednoho segmentu obsahuje kromě výrazných složek na příslušných
kmitočtech také další složky díky přítomnosti šumu. Průměrováním několika takovýchto
úseků se složky náhodné složky potlačí.
Pozn. Před výpočtem individuálních spekter se standardně používá váhování některým oknem pro
zmenšení prosakování.


obr. 7.3. Výsledek získaný metodou periodogramu
7.2.2 Odhad výkonového spektra metodou korelogramu
Tato metoda využívá Wiener-Chinčinův teorém (viz také (6.6)) pro výpočet
výkonového spektra. Ten dává do souvislosti autokorelační funkci náhodného procesu r
xx
(τ)
s jeho výkonovým spektrem S
xx
(k)
64 Fakulta elektrotechniky a komunikačních technologií VUT v Brně
() (){} ()
∑
−
+−=
Ω−
==
1
1
N
N
Tjk
xxxxxx
erRDFTkS
τ
τ
ττ .
(7.4)
Autokorelační funkce r
xx
(τ) je funkce deterministická a symetrická kolem středu. Fourierova
transformace této funkce, tedy výkonové spektrum je také funkce deterministická a navíc díky
sudé symetrii i funkce reálná.
Pro stacionární náhodné procesy počítáme autokorelaci přes jednotlivé realizace:
() () ( )
∑
=
−≈
⎭
⎬
⎫
⎩
⎨
⎧
−Ε=
M
ii
w
ww
wwwwxx
nxnx
M
nxnx
N
r
1
)()(
11
τττ .
(7.5)
Pro ergodické procesy vypočteme autokorelační funkci z jedné změřené realizace dle vztahu
(6.3), který v sobě zahrnuje průměrování, takže výsledné spektrum má již menší rozptyl.
Nutno poznamenat, že odhad autokrelační funkce je pro velké hodnoty τ značně nepřesný,
protože do průměru vstupuje již málo hodnot. Proto se před vlastním výpočtem výkonového
spektra často provede odříznutí několika posledních hodnot r
xx
(τ).
Existuje ovšem řada modifikací této metody. Například je možné jako u metody
periodogramu rozdělit dlouhý úsek ergodického signálu na kratší segmenty a z každého určit
dílčí autokorelační funkci. Výslednou autokorelační funkci pak získáme jejich průměrováním.
Na výpočet výkonového spektra nám pak postačí aplikace jedné DFT. V případě metody
periodogramu jsme jich totiž potřebovali M.
Pozn. Vztah (7.4) je platný v případě, že autokorelační funkce byla počítána dle vychýleného odhadu
(6.3) resp. (7.5). Je