Ing. Ji ř í Kozumplík, CSc. Habilita č ní práce Vlnkové transformace a jejich využití pro filtraci signál ů EKG

VYSOKÉ UČENÍ TECHNICKÉ V BRNĚ

FAKULTA ELEKTROTECHNIKY A KOMUNIKAČNÍCH TECHNOLOGIÍ ÚSTAV BIOMEDICÍNSKÉHO INŽENÝRSTVÍ

Vlnkové transformace a jejich využití pro filtraci signálů EKG Habilitační práce

Ing. Jiří Kozumplík, CSc.

Listopad 2004

Obsah

Úvodní poznámky ...4

1 Vlnková transformace s diskrétním časem...5

1.1 Vlnkové transformace spojitého signálu ...5

1.2 Reálná dyadická vlnková transformace s diskrétním časem (DTWT)...6

1.2.1 Inverzní transformace (IDTWT) ...9

1.2.2 Půlpásmové filtry a jejich návrh...12

1.2.2.1 Nekauzální půlpásmové filtry s nulovou fázovou charakteristikou ...12

1.2.2.2 Kauzální půlpásmové filtry s lineární fázovou charakteristikou ...14

1.2.2.3 Návrh vycházející z Lagrangeova interpolačního vzorce...15

1.2.2.4 Návrh vycházející z rozložení nulových bodů v Z=(z+z

^-1

)/2...17

1.2.3 Odvození filtrů pro biortogonální DTWT ...18

1.2.4 Odvození filtrů pro ortogonální DTWT ...20

1.2.4.1 Filtry pro ortogonální DTWT podle Daubechiesové...23

1.2.4.2 Konstrukce filtrů pro ortogonální DTWT z rozložení nulových bodů...25

1.3 Komplexní dyadická vlnková transformace (DTCWT) ...27

1.4 Redundantní DTWT ...31

1.5 Paketová DTWT...32

2 Využití DTWT pro filtraci signálů EKG...33

2.1 Proč se může uplatnit DTWT při zpracování signálů EKG ...33

2.2 Filtrace signálů s využitím DTWT ...33

2.2.1 Výběr typu DTWT ...34

2.2.2 Prahování koeficientů DTWT ...37

2.2.2.1 Prahování koeficientů reálné DTWT...37

2.2.2.2 Prahování koeficientů komplexní DTWT ...38

2.2.3 Stanovení prahových hodnot pro vlnkovou filtraci...39

2.2.3.1 Rozptyly koeficientů DTWT bílého šumu ...40

2.2.3.2 Univerzální práh ...41

2.2.3.3 Empirický práh ...41

2.2.3.4 Práh vycházející ze zobecněného Gaussova rozložení koeficientů...41

2.2.4 Wienerovská filtrace v časově-měřítkové oblasti...42

2.2.4.1 Hybridní prahování...44

2.2.4.2 Metoda pilotního odhadu...45

2.3 Vlnková filtrace signálů EKG ...46

2.3.1 Testované signály a model rušení ...48

2.3.2 Shrnutí výsledků dřívějších experimentů ...49

2.3.2.1 Odhad rozptylu šumu...49

2.3.2.2 Experimenty s empirickými prahy ...49

2.3.2.3 Experimenty s hybridním prahováním ...52

2.3.3 Výsledky nových experimentů...53

2.3.3.1 Předzpracování testovacích signálů...53

2.3.3.2 Výběr testovacích signálů a metod vlnkové filtrace...55

2.3.3.3 Dosažené výsledky a jejich zhodnocení ...55

3 Závěr...68

4 Dodatky ...69

4.1 Podvzorkování a expanze diskrétního signálu ...69

4.2 Lineární fázové frekvenční charakteristiky u komplexních filtrů ...71

4.3 Úprava Matlabu pro realizaci komplexní DTWT ...72

Literatura ...73

Zkratky a symboly...76

Seznam obrázků a tabulek ...78

Poděkování

Tato práce byla podporována projekty

GAČR 102/02/0890

, GAČR 102/04/0472 a výzkumným záměrem MSM 262200011 (CZ 410011).

Úvodní poznámky

Tato práce je studií zaměřenou na filtraci elektrokardiografických signálů s využitím vlnkové transformace s diskrétním časem (DTWT).

Problematika vlnkové transformace a jejích aplikací se v posledních letech těší značné pozornosti. Svědčí o tom řada knižních publikací úzce zaměřených jen na tuto tematiku. Můžeme se setkat s monotematicky zaměřenými publikacemi čistě matematickými [8][17], s knihami psanými z užšího pohledu číslicového zpracování signálů [12][45][49] nebo s publikacemi, ve kterých se mísí oba pohledy [31][50]. Vedle monotematických publikací existuje celá řada dalších, které věnují problematice vlnkové transformace některou z kapitol – za všechny jmenujme dvě české knihy [20][51]. Teoreticky nebo aplikačně zaměřených časopiseckých publikací je obrovské množství. Pro zpracování úvodu do vlnkové transformace na několika stranách tohoto textu byl tedy k dispozici velký objem solidních pramenů.

První kapitola předložené práce je pokusem o stručný a srozumitelný úvod do DTWT. Mým záměrem bylo na velmi omezeném rozsahu (30 stran) textu vysvětlit problematiku transformace tak, aby byl výklad srozumitelný pro studenty, kteří absolvovali předmět Číslicové zpracování a analýza signálů na naší fakultě a orientují se v doporučené literatuře [20][27]. Ve snaze o srozumitelnost textu jsem se tedy neomezoval na pouhé uvádění „známých“ vzorců s příslušnými odkazy na literaturu, ale považoval jsem za žádoucí uvést všechna nutná odvození (některá jsem musel i přes obsáhlé zdroje informací někdy domýšlet) a souvislosti, s výjimkou definice vlnkové transformace (1.1) a vzorců (1.41) a (1.45). Kromě těchto tří vzorců jsem v textu upustil od odkazů na literaturu, proto alespoň na tomto místě konstatuji, že jsem hledal inspiraci ve všech publikacích citovaných v předchozím odstavci a navíc ještě v článcích [3][13][14][29], z nichž se poslední tři zabývají problematikou komplexní DTWT, kterou jsem ve výše uvedených knihách nenašel. Přiznávám, že snaha o relativně jednoduchý a výstižný text způsobila, že tato kapitola byla časově nejnáročnější. Záměrně jsem přitom vypustil úvahy o zasazení vlnkové transformace do širších souvislostí, zejména s Fourierovou transformací, protože jsou tyto úvahy velmi dobře zpracovány v knize [20], která je našim studentům dostupná. Tuto kapitolu zamýšlím použít jako část učebního textu mého předmětu Nové algoritmy zpracování signálů.

Druhá část práce je zaměřena na principy metod vlnkové filtrace signálů, tj. na způsoby upravování koeficientů DTWT. Teorie zabývající se touto problematikou je poměrně rozsáhlá a náročná (viz např. [9][10][21]), v zájmu stručnosti se tato část práce omezuje na známé metody a vzorce bez složitého odvozování a proto s potřebnými odkazy na literaturu. Jen snad nejzajímavější a v této práci preferovaná metoda vlnkové wienerovské filtrace je uvedena podrobněji včetně souvisejících odvození (která jsou v dostupné literatuře velmi strohá). Součástí této kapitoly jsou úvahy o možnosti použití této nelineární filtrace pro potlačení náhodných myoptenciálů, které bývají obsaženy v signálech EKG, v kap.2.3 jsou shrnuty výsledky našich dřívějších studií a v závěrečné části výsledky nové.

V minulosti jsem se zabýval problematikou komprese 1D i 2D dat s využitím DTWT, jedním z výstupů byla studie zaměřená na kompresi signálů EKG [38]. V současnosti existuje poměrně bohatý výběr časopiseckých publikací na téma ztrátové komprese EKG s využitím DTWT, přesto si dovoluji tvrdit, že s rozvojem informačních technologií pozbývá toto téma na významu, nemluvě o vcelku pochopitelné neochotě lékařské komunity smířit se s „poškozenými“ vstupními daty. Téma vlnkové filtrace signálů EKG jsem vybral proto, že je mu podle mého soudu věnována nezaslouženě malá pozornost, dostupné publikace jsou pouze stručné příspěvky z různých konferencí a většina publikovaných metod nemůže dávat dobré výsledky už proto, že je založena na „klasické“ DTWT s decimací (která je na druhé straně jediná použitelná pro ztrátovou kompresi dat). Navíc, s ohledem na zmíněný rozvoj komunikačních technologií, se jedná do značné míry o téma nadčasové.

V Brně, v listopadu 2004 Jiří Kozumplík

1 Vlnková transformace s diskrétním časem

1.1 Vlnkové transformace spojitého signálu

Vlnková transformace se spojitým časem (WT - Wavelet Transform) signálu x(t) je definována [8]

jako

∫

∞

∞

−

⎟ ⎠

⎜ ⎞

⎝

= ⎛ − dt

a b t t

a x b a

y 1 ( ) *

) ,

( ψ

. (1.1)

Jedná se o časově-frekvenční rozklad, který můžeme interpretovat jako korelaci signálu x(t) s funkcemi (vlnkami - z angl. wavelets) odvozenými z obecně komplexní mateřské vlnky ψ(t). Pro funkce ψ(t) se vžil název vlnky s ohledem na jejich tvary - ψ(t) musí mít nulovou střední hodnotu a tvarem často připomíná vlnku. Symbol * značí komplexně sdruženou funkci. Výsledná funkce y(a,b), stejně jako jednotlivé vlnky ψa,b(t), je popsána dvěma (spojitě proměnnými) parametry: časovým posunutím b a dilatací a, která určuje frekvenční spektrum příslušné vlnky. Konstanta a^-1/2 normalizuje energii jednotlivých vlnek.

Zvláštním případem transformace signálu se spojitým časem je diskrétní vlnková transformace (DWT) s parametry a=a₀^m a b=a₀^mkT , kde a₀>1, T>0 a m,k jsou celočíselné.

Nejčastější je dyadická DWT pro a=2^m , b=2^mkT, m>0. Koeficienty dyadické DWT jsou dt

kT t t

x k

m

y ^m

m

∫

^∞

∞

−

− −

= ( ) (2 )

2 ) 1 ,

( ψ^* . (1.2)

Index m reprezentuje kmitočtové měřítko, index k časové měřítko. Konstanta T (která závisí na šířce pásma B mateřské vlnky, když T=1/(2B)) určuje hustotu vzorkování koeficientů na časové ose pro jednotlivé kmitočtové úrovně dané indexem m.

Fourierův obraz mateřské vlnky označme jako Ψ(ω)=F{ψ(t)}. Z Fourierova obrazu m-té vlnky (s normalizovanou energií)

( )

^j

⁽

^{x kT}

⁾

^{m m j kT}

m

m m m

m m

t j m m m m

m m

m m

m dx e

e x

dx dt

kT x t

kT x t dt kT e

t kT

F t

2 2

2 2 (2 )

2

2 2 22

2

2 2 2

1 2

2 2

1

ω ω

ω

ω ψ

ψ ψ

−

∞

∞

−

+

−

∞

∞

−

−

Ψ

=

=

=

= +

=

− =

⎟⎟ =

⎠

⎜⎜ ⎞

⎝

= ⎛ −

⎪⎭

⎪⎬

⎫

⎪⎩

⎪⎨

⎧

⎟⎟⎠

⎜⎜ ⎞

⎝

⎛ −

∫

∫

(1.3)

vyplývá, že časová expanze vlnky na 2^m-násobnou délku se projeví kompresí jejího spektra na 1/2^m- násobek výchozí šířky a jeho posunem k nižším frekvencím, se středním kmitočtem na 1/2^m-násobku výchozího. Dyadická DWT je tedy charakterizována oktávovou podobou spekter soustavy vlnek. S rostoucím m se krok posunutí a zvětšuje 2^m-krát. Výsledkem je (obecně nekonečná) množina koeficientů y(m,k) nerovnoměrně rozložených v časově-frekvenční rovině.

Zaveďme do (1.2) substituci

dt d kT t

kT

t

^{m m}

m

− =

⁻

⇒ = − =

−

2 τ τ 2 , τ

2

.

Rovnici (1.2) pak můžeme přepsat do podoby

∫

∫

^∞

∞

−

−

∗

∞

∞

−

−

∗ = −

+

= xτ kTψ τ dτ xτ ψ τ kT dτ

k m

y ^{m m}

m m

m

m ( ) (2 2 )

2 ) 1

2 ( ) 2 2 (

) 1 ,

( . (1.4)

Korelace signálu x(t) s jednotlivými vlnkami

) 2 ( ) 2 2 2 (

1 ^mt ^mkT t ^mkT

mψ^{∗ −} − =ψ^∗m − _(1.5)

může být realizována konvolucí s funkcemi časově reverzními,

( ) ( )

m,k xt (2 kT t) x

( )

t h (2 kT t)

y = ∗ψ_m^{∗ m} − = ∗ _m ^m − . (1.6)

Dyadická DWT pak může být vyjádřena jako

∫ ( )

∫

^∞

∞

−

∞

∞

−

−

=

−

= x τ h kT τ d τ h τ x kT τ d τ

k m

y ( , ) ( )

_m

( 2

^m

)

_m

( 2

^m

)

(1.7)

a lze ji realizovat rozkladem signálu bankou lineárních spojitých oktávových filtrů s impulsními charakteristikami h_m(t).

Má-li být možná přesná rekonstrukce signálu x(t)=DWT^-1 [y(m,k)], musí být jednotlivé koeficienty y(m,k) vzájemně nezávislé, což je splněno, tvoří-li vlnky ortogonální systém.

1.2 Reálná dyadická vlnková transformace s diskrétním časem (DTWT)

Dyadická vlnková transformace s diskrétním časem (DTWT) ym(n) diskrétního signálu x(n) je definována analogicky k (1.7) diskrétní konvolucí,

∑

∑

^∞

−∞

=

∞

−∞

=

−

=

−

=

i

m m i

m m

m n x i h n i h i x n i

y ( ) ( ) (2 ) ( ) (2 ) , (1.8)

tj. rozkladem signálu bankou diskrétních oktávových filtrů s impulsními charakteristikami hm(n).

Uvažujme nejdříve případ reálné transformace, tj. s filtry s reálnými impulsními charakteristikami.

Vzorkovací frekvence signálu ym(n) na výstupu m-tého filtru je 2^m-krát nižší než vzorkovací frekvence fvz vstupního signálu x(n).

Při použití kauzálních FIR filtrů s impulsními charakteristikami h_m(n), n=0,1,...,N_m-1 a při předem zvoleném stupni rozkladu M, kdy m = 1, 2, ... , M, můžeme psát

( ) ( ) ( 2 ) ,

, ,..., 2 , 1 ,

) 2 ( ) ( )

(

1 0

1 1

1

0

∑

1

∑

−

= +

+

−

=

+

−

=

=

−

=

M m

N

i

M M

M N

i

m m

m

i n x i h n

y

M m

i n x i h n

y

(1.9)

kde y_M+1(n) jsou koeficienty korespondující s nejnižším frekvenčním pásmem po M-stupňovém rozkladu.

Transformace s třístupňovým rozkladem (tj. pro M=3) a modulové frekvenční charakteristiky odpovídajících ideálních oktávových filtrů jsou uvedeny na obr.1.1.

H₁(z)

H₂(z)

H₃(z)

H₄(z) x(n)

2

4

8

8

y₁(n)

H₁ H₂

H₃ H₄

f/f_vz 1/2 1/4

1/8 1/16 0

2^1/2 2.2^1/2

2

|H_i|

y₂(n)

y₃(n)

y₄(n)

obr.1.1 Vlevo: realizace třístupňové dyadické DTWT bankou oktávových filtrů s podvzorkováním výstupů.

Výstupní posloupnosti jsou koeficienty dyadické DTWT. Bloky se symbolem ↓2^m zajišťují

podvzorkování, tj. výběr každého 2^m-tého vzorku signálu. Vpravo: modulové frekvenční charakteristiky ideálních oktávových filtrů.

Koeficienty dyadické DTWT jsou tedy tvořeny výstupními vzorky banky filtrů. Vzhledem k tomu, že jsou výstupy filtrů podvzorkovány, jak vyplývá z pravé strany (1.8), je počet koeficientů transformace shodný s počtem vzorků vstupního signálu x(n).

Předpokládejme dvojici zrcadlových filtrů - ideální dolní propust Hd a ideální horní propust Hh

s modulovými charakteristikami navzájem symetrickými okolo ω_vz/4=π/2,

( ) )

( ^{π π}

ω ω π ω π

ω

, pro

pro , pro

e H

_d ^j

0 2

2

1 2 0 2

1

∈

=

∈

=

,

( ) )

( ^{π π}

ω ω π ω π

ω

, pro

pro

, pro

e H

_h ^j

1 2

2

1 2 0 2

0

∈

=

∈

=

. (1.10)

Dále připomeňme jednoduchou frekvenční transformaci: provedeme-li v přenosové funkci H(z) substituci z→z^k, obdržíme systém H(z^k) s k-krát „stlačenou“ frekvenční charakteristikou. Na obr.1.2 jsou naznačeny frekvenční charakteristiky filtrů odvozených z výchozích ideálních dolních a horních propustí (1.10) a na obr.1.3 jsou oktávové filtry pro třístupňovou DTWT, jejichž konstrukce vychází z dvojice těchto zrcadlových filtrů. Blokové schéma realizace z obr.1.1 pak můžeme překreslit do podoby uvedené na obr.1.4.

f/f_vz 1/2 1/4

1/8 1/16 0

H_h(z)

|H_h1|

H_h(z²)

|H_h2|

f/f_vz 1/2 1/4

1/8 1/16 0

H_h(z⁴)

|H_h4|

f/f_vz 1/2 1/4

1/8 1/16 0

f/f_vz 1/2 1/4

1/8 1/16 0

H_d(z)

|H_d1|

H_d(z²)

|H_d2|

f/f_vz 1/2 1/4

1/8 1/16 0

H_d(z⁴)

|H_d4|

f/f_vz 1/2 1/4

1/8 1/16 0

obr.1.2 Filtry odvozené z výchozích dolních a horních propustí frekvenční transformací vycházející ze substituce z→z^k v přenosové funkci.

H_h(z)

|H|

f/f_vz

1/4 1/2 1/16 1/8

0

H_d(z)H_d(z²)H_h(z⁴)

H_d(z)H_h(z²) H_d(z)H_d(z²)H_d(z⁴)

obr.1.3 Oktávové filtry pro třístupňovou DTWT odvozené z transformovaných filtrů na obr.1.2.

H_h(z)

H_d(z) H_h(z²) H_d(z²) x(n)=c₀(n)

2 4 8 8

y₁(n)

H_h(z⁴) H_d(z⁴)

y₂(n) y₃(n) y₄(n) c₁(n)

c₂(n)

c₃(n)

obr.1.4 Realizace třístupňové dyadické DTWT s filtry odvozenými ze zrcadlových dolních propustí Hd a horních propustí Hh.

Poznamenejme, že nastavením přenosů v propustných pásmech zrcadlových filtrů (1.10) na hodnotu 2^1/2 získáme žádoucí přenosy v charakteristikách na obr.1.1.

Substituce z→z^k v přenosové funkci se projeví „zředěním“ impulsní charakteristiky filtru, které spočívá ve vložení k-1 nulových vzorků mezi jednotlivé vzorky impulsní charakteristiky výchozího filtru. Všimněme si např. filtru Hh(z²) z obr.1.4, který vznikl z Hh(z). Díky zředěné impulsní charakteristice počítá filtr Hh(z²) v každém taktu jen s každým druhým vstupním vzorkem c1(n) a výstup je následně podvzorkován s faktorem 4. Z toho vyplývá, že lze podvzorkování realizovat postupně tak, že výstup každého filtru v sérii podvzorkujeme s faktorem 2. Dosáhneme toho, že všechny použité filtry budou z dvojice zrcadlových filtrů Hd aHh, jejichž frekvenční charakteristiky podvzorkování s faktorem 2 umožňují. Schéma na obr.1.4 tedy můžeme překreslit do podoby na obr.1.5.

H_h(z)

H_d(z) H_h(z)

H_d(z) x(n)=d₀(n)

2

2

y₁(n)

H_d(z) d₂(n)

d₁(n) y₂(n)

y₃(n) d₃(n)=y₄(n) 2

2 H_h(z) 2

2

obr.1.5 Realizace třístupňové rychlé dyadické DTWT se zrcadlovými dolními propustmi Hd a horními propustmi Hh.

Realizace dyadické DTWT pomocí stromové struktury bank filtrů na obr.1.5 bývá nazývána rychlou DTWT.

Při použití kauzálních FIR dolních a horních propustí s impulsními charakteristikami {hd(n), n=0,1,...,Nd-1} a {hh(n), n=0,1,...,Nh-1} a při zvoleném stupni rozkladu M můžeme podle obr.1.5 psát

( ) n h ( i ) d ( n i ) , y

M ,..., , m , ) i n ( d ) i ( h ) n ( y

, M ,..., , m , ) i n ( d ) i ( h ) n ( d

, ) n ( x ) n ( d

d h d

N

i d M

M

N

i

m h m

N

i

m d m

∑

∑

∑

−

= +

−

= −

−

= −

−

=

=

−

=

=

−

=

=

1 0 1

1 0

1 1

0

1 0

2

2 1 2

2 1 2

(1.11)

kde dm(n) je podvzorkovaný výstup m-té dolní propusti, ym(n) je vektor koeficientů DTWT v m-tém stupni rozkladu, když yM+1(n) je vektor koeficientů korespondujících s nejnižším frekvenčním pásmem po M-stupňovém rozkladu.

Výsledkem dyadické transformace při M-stupňovém rozkladu jsou koeficienty ym(n), m=1,2,

…,M+1, které jsou nerovnoměrně rozloženy v časově-frekvenční rovině, jak je naznačeno (pro M=3) na obr.1.6. Z obr.1.6 je zřejmé, že s rostoucím m klesá (díky podvzorkování) časová rozlišovací schopnost koeficientů. Naopak frekvenční rozlišovací schopnost s rostoucím m (tedy směrem k nižším kmitočtům) roste, což vyplývá z oktávové podoby banky rozkladových filtrů (viz obr.1.2).

y₁(0) y₁(2) y₁(4) y₁(6) y₁(8) y₁(10) y₁(12) y₁(14)

y₂(0) y₂(1) y₂(2) y₂(3) y₂(4) y₂(5) y₂(6) y₂(7)

y₃(0) y₃(1) y₃(2) y₃(3)

y₄(0) y₄(1) y₄(2) y₄(3)

n

m

obr.1.6 Nerovnoměrné rozložení koeficientů ym(n) dyadické DTWT v časově-frekvenční rovině při třístupňovém rozkladu.

1.2.1 Inverzní transformace (IDTWT)

Princip inverzní transformace (pro třístupňovou dyadickou DTWT na obr.1.5) je zachycen na obr.1.7.

Podvzorkované posloupnosti (koeficienty transformace) je nutné interpolovat, každý interpolátor je tvořen expanderem (který vkládá nulové vzorky mezi sousední vzorky posloupnosti) a interpolačním (resp. rekonstrukčním) filtrem, kterým je buď dolní nebo horní propust. Příslušný rekonstrukční filtr musí být vhodným protějškem korespondujícího filtru rozkladového. Uvažujeme-li pouze kauzální filtry, je nutné použít zpožďovací členy, jak je ukázáno na obr.1.7.

F_d(z) y₁(n)

y₄(n)

F_h(z) 2

2

2 z^-τ1

F_d(z) F_h(z) y₃(n)

2 2

2 z^-τ2

y₂(n)

F_d(z) F_h(z)

x'(n)=x(n−τ)

obr.1.7 Princip inverzní transformace (IDTWT) pro třístupňovou DTWT. Blok ↑2 realizuje expanzi posloupnosti (vložením jednoho nulového vzorku mezi sousední vzorky signálu), Fd (resp. Fh) je rekonstrukční dolní (resp. horní) propust. Vstupem jsou koeficienty DTWT z obr.1.5.

Z transformace na obr.1.5 a inverzní transformace na obr.1.7 je zřejmé, že základem je dvoukanálová banka rozkladových (H_h, H_d) a rekonstrukčních (F_h, F_d) filtrů na obr.1.8. Zaměřme se nyní na podmínky, jaké musí splňovat čtveřice použitých filtrů na obr.1.8, abychom na výstupu získali signál totožný se zpožděným (díky kauzalitě filtrů) vstupním signálem, x´(n)=x(n-τ). Tím stanovíme podmínky inverzibility transformace a zároveň potvrdíme správnost schématu pro IDTWT na obr.1.7.

H_d(z)

H_h(z) 2

2

2

2

F_h(z)

F_d(z)

x(n) c_h(n)

c_d(n) y_d(n)

y_h(n) q_h(n)

q_d(n)

x'(n)=x(n−τ)

obr.1.8 Dvoukanálová banka rozkladových (Hh, Hd) a rekonstrukčních (Fh, Fd) filtrů.

Zajímají nás podmínky věrné rekonstrukce vstupního signálu x(n) pro skutečné filtry, takže upustíme od původní představy ideálních zrcadlových fitrů H_h, H_d podle (1.10). Použití filtrů s neideálními frekvenčními charakteristikami povede po podvzorkování samozřejmě ke vzniku aliasingu.

Obrazy výstupních signálů rozkladových filtrů jsou

( ) z H ( z ) X ( z ) , i h , d

C

_i

=

_i

=

, (1.12)

po podvzorkování obdržíme signály, jejichž obrazy jsou podle (4.8) (viz Dodatek 4.1)

( ) z C z C z , i h , d

Y

_{i i i}

⎥⎦ ⎤ =

⎢⎣ ⎡ ⎟

⎠ ⎞

⎜ ⎝

⎛−

⎟ +

⎠ ⎞

⎜ ⎝

= ⎛

^{12 12}

2

1

. (1.13)

Obrazy signálů po expanzi jsou

( ) z Y ( ) z [ C ( ) z C ( ) z ] [ H ( z ) X ( z ) H ( z ) X ( z ) ] , i h , d

Q

_i

=

_i

=

_i

+

_i

− =

_i

+

_i

− − =

2 1 2

2

1

, (1.14)

obraz výstupního (rekonstruovaného) signálu je

( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( )

[

F z H z F z H z

]

X(z)

[

F

( ) ( )

z H z F

( ) ( )

z H z

]

X( z), z

Q z F z Q z F z ' X

h h d

d h

h d

d

h h d

d

−

− +

− +

+

=

= +

=

2 1 2

1 ^(1.15)

kde druhý člen představuje vliv aliasingu; abychom vliv aliasingu potlačili, musí být tento člen nulový. První člen v (1.15) by měl korespondovat se zpožděným vstupním signálem. Podmínky věrné rekonstrukce vstupního signálu, kdy x´(n)=x(n-τ), jsou

( ) ( ) z H z + F ( ) ( ) z H z = z

^−τ

F

_{d d h h}

2

(1.16)

a

( ) ( ) z H − z + F ( ) ( ) z H − z = 0

F

_{d d h h} , (1.17)

když τ je fázové zpoždění filtrů Hd(z)Fd(z) a Hh(z)Fh(z). Podmínka (1.17) vede k výběru antialiasingových filtrů

) ( )

(z H z

F_d = _h − a F_h(z)=−H_d(−z) (1.18) nebo

) ( )

( z H z

F

_d

= −

_h

−

a

F

_h

( z ) = H

_d

( − z )

(1.19)

Podmínku (1.16) lze po zavedení antialiasingových filtrů (1.18) nebo (1.19) upravit do tvaru

τ

= −

−

−

=

−

−

−F z H z P z P z z

z H z

F_d( ) _d( ) _d( ) _d( ) _d( ) _d( ) 2 . (1.20)

Filtry Pd(z) a Pd(-z) jsou zrcadlovými filtry, když Pd(z) je dolní propust a Pd(-z)=Ph(z) je horní propust.

Podívejme se nyní na vlastnosti a vhodný výběr těchto důležitých filtrů.

Předpokládejme kauzální dolní propust Pd(z) s impulsní charakteristikou

( ) n { p ( ) ( ) ( ) ( ) ( ) ( ) ( ) , p , p , p , p , p , p } .

p

_d

=

_d

0

_d

1

_d

2

_d

3

_d

4

_d

5

_d

6

(1.21)

Zrcadlová horní propust Pd(-z)=Ph(z) musí mít samozřejmě impulsní charakteristiku s opačnými znaménky u vzorků s lichými indexy,

( ) n { p ( ) , p ( ) ( ) , p , p ( ) ( ) , p , p ( ) ( ) , p } .

p

_h

=

_d

0 −

_d

1

_d

2 −

_d

3

_d

4 −

_d

5

_d

6

(1.22)

Má-li být splněna podmínka (1.20), musí být vzorky s lichými indexy s vyjímkou prostředního vzorku pd(3) nulové, abychom po odečtení přenosových funkcí Pd(z) a Pd(-z) obdrželi výraz 2z^-3. Prostřední vzorek by měl být pd(3)=1. Má-li být fázové zpoždění obou zrcadlových filtrů konstantní (v našem případě 3), musí být jejich impulsní charakteristiky symetrické. Filtry, které uvedeným podmínkám vyhovují, nazvěme půlpásmovými filtry (z angl. halfband filters). Podrobněji se na tyto filtry zaměříme v následující kapitole.

Návrh dvoukanálové banky filtrů (na obr.1.8) s možností přesné rekonstrukce signálu pak může být následující.

• Vybereme vhodnou dolní propust Pd(z), která vyhovuje podmínce (1.20). Systém Pd(z) je dolní propust a systém Pd(-z) horní propust, oba filtry mají zrcadlově symetrické frekvenční charakteristiky kolem relativní úhlové frekvence ω=π/2. Filtr Pd(z) je půlpásmový filtr, který má na ω=π/2 přenos poloviční oproti (maximálnímu) přenosu na ω=0.

• Systém P_d(z) převedeme na sériové spojení dvou filtrů H_d(z)F_d(z), přenosové funkce H_h(z) a F_h(z) odvodíme z (1.18) nebo (1.19).

1.2.2 Půlpásmové filtry a jejich návrh

Ideální půlpásmovou dolní propust jsme již popsali v (1.10). Nekauzální varianta s nulovou fázovou charakteristikou má reálnou frekvenční charakteristiku

( ) )

( ^{π π}

ω ω π ω π

ω

, pro

pro , pro

e H

_d ^j

0 2

2

1 2 0 2

1

0

∈

=

∈

=

. (1.23)

Díky periodicitě frekvenční charakteristiky představují vzorky impulsní charakteristiky koeficienty Fourierovy řady,

∞

−

−∞

=

=

= ∫

−

,..., , , ,..., n

n , sin n d

e )

n ( h

vz

vz

j vz

d

1 0 1

2 2 2 1

1

⁴

4 0

ω ω

ω

π

π

ω ω

^{. (1.24)}

Část nekonečné impulsní charakteristiky je zachycena na obr.1.9.

obr.1.9 Část nekonečné impulsní charakteristiky ideální půlpásmové dolní propusti.

1.2.2.1 Nekauzální půlpásmové filtry s nulovou fázovou charakteristikou

Abychom získali realizovatelnou nerekursivní půlpásmovou dolní propust, musíme impulsní charakteristiku hd0(n)vynásobit vhodnou konečnou symetrickou váhovou posloupností ρ0(n) k získání výsledné impulsní charakteristiky

2 ,..., 1 2 , 1

) ( ) ( )

( _{0 0}

0

−

− −

=

= N N

n n n h n

p_{d d} ρ ,

kde N je liché číslo představující počet jejích vzorků (má-li mít výsledný filtr nulovou fázovou frekvenční charakteristiku, nemůže být N sudé číslo).

Pro nekauzální dolní propust Pd0 s nulovou fázovou frekvenční charakteristikou (tj. se impulsní charakteristikou symetrickou kolem n=0), jejíž frekvenční charakteristika je

( )

^ω

( )

^ω

( )

^N

( ) [

^{jnω jnω}

]

n d d

N

n N

jn d

j

d

e p n e p p n e e

P = = +

⁻

+

−

=

−

− −

=

−

∑

∑

²

1

1 0 0

2 1

2 1 0

0

0

, (1.25)

požadujeme, aby

pro

0 ( ) ( ) 0 2 ( ) ( ) 0 ( ) 0 2

²

( ) 1

1

1 0 0

2 1

1 0 0

0

0

= + = + =

= ∑ ∑

−

=

−

=

N

n d d

N

n d d

j

d

e p p n cos p p n

P

ω :

, (1.26)

pro

( ) ( ) 0 2 ( ) ( ) ( ) 0 2

²

( ) 0

1

1 0 0

2 1

1 0 0

0

= + = − =

= ∑ ∑

−

=

−

=

N

n d d

N

n d d

j

d

e p p n cos n p p n

P

: π

π

ω

^π . (1.27)

K jednoduchým výsledným tvarům obou rovnic připomeňme, že pro sudé indexy n platí, že pd0(n)=0 (viz obr.1.9). Sečtením, resp. odečtením, obou předchozích rovnic obdržíme

2 ) 1 0

0( =

pd , resp.

∑

−

= 2 =

1

1

0 2

) 1 ( 2

N

n

d n

p . (1.28)

Dále požadujeme, aby

pro

( )

2 1 2

2

0

=

= π : P

_d

e

^jπ/

ω

. (1.29)

Podmínka (1.29) je splněna, protože platí p_d0(n)cos(nπ/2)=0. Skutečně, protože pro sudé n je p_d0(n)=0 a pro liché indexy n je cos(nπ/2)=0.

Frekvenční charakteristiky půlpásmové dolní propusti P_d0 (resp. horní propusti P_h0) jsou

( ) ( ) ∑ ( ) ( )

−

=

+

=

²

1

1 0 0

0

0 2 cos

N

n d d

j

d

e p p n n

P

^ω

π

, (1.30)

resp.

( ) ( ) ∑ ( ) ( )

−

=

−

=

²

1

1 0

0

0

0 2 cos

N

n d

d j

h

e p p n n

P

^ω

π

. (1.31)

Na obr.1.10 jsou zachyceny impulsní charakteristiky (nahoře) a frekvenční charakteristiky půlpásmové dolní propusti (vlevo) a horní propusti (vpravo). Součet frekvenčních charakteristik obou půlpásmových filtrů je

( )

₀

( ) 2

₀

( ) 0 1

0 ^j

+

_h ^j

=

_d

=

d

e P e p

P

^{ω ω} (1.32)

nebo jinak

( )

0

(

^{( )}

) ¹

0 ^jω

+

_d ^jω−π

=

d

e P e

P

. (1.33)

Pro přenosové funkce zřejmě platí

( )

₀

( ) 1

0

z + P z =

P

_{d h} , (1.34)

resp.

( )

₀

( ) 1

0

z + P − z =

P

_{d d} . (1.35)

obr.1.10 Impulsní charakteristiky (nahoře) a frekvenční charakteristiky (dole) půlpásmové dolní propusti (vlevo) a horní propusti (vpravo). Počet vzorků impulsní charakteristiky N=7.

1.2.2.2 Kauzální půlpásmové filtry s lineární fázovou charakteristikou

Kauzální půlpásmová dolní propust má frekvenční charakteristiku

( ) ( )

⎥ ⎥

⎥

⎦

⎤

⎢ ⎢

⎢

⎣

⎡

⎥ ⎥

⎦

⎤

⎢ ⎢

⎣

⎡ +

⎟ +

⎠

⎜ ⎞

⎝

= ⎛ −

^⎟⎠

⎜ ⎞

⎝

⎛ − −

⎟⎠

⎜ ⎞

⎝

⎛ − −

− −

−

=

− −

∑

^{2 1 21 21}

1

0 2

1

2

1

^{j n N j n N}

N

n d d

j N j

d

N p n e e

p e

e

P

^{ω ω ω ω} . (1.36)

Transformace z→ -z (pro realizaci zrcadlové horní propusti) vede ke změně znaménka prvního činitele,

( ^{− e}

^jω

)

^−N2−1

^{= − e}

^jωN2−1,

protože N=3, 7, 11, …; výraz (N-1)/2 má tedy vždy lichou hodnotu. Výsledná kauzální půlpásmová horní propust má frekvenční charakteristiku

( ) ( )

⎥ ⎥

⎥

⎦

⎤

⎢ ⎢

⎢

⎣

⎡

⎥ ⎥

⎦

⎤

⎢ ⎢

⎣

⎡ +

⎟ −

⎠

⎜ ⎞

⎝

− ⎛ −

=

^⎟⎠

⎜ ⎞

⎝

⎛ − −

⎟⎠

⎜ ⎞

⎝

⎛ − −

− −

−

=

− −

∑

^{2 1 21 21}

1

0 2

1

2

1

^{j n N j n N}

N

n d d

j N j

h

N p n e e

p e

e

P

^{ω ω ω ω} . (1.37)

Odečteme-li charakteristiku horní propusti (1.37) od charakteristiky dolní propusti (1.36) dostaneme

( ) ( )

^{21 21}

2 2 1

− −

− −

⎟ =

⎠

⎜ ⎞

⎝

= ⎛ −

−

_h ^{j j N} _d ^{j N}

j

d

N e

p e

e P e

P

^{ω ω ω ω} , (1.38)

resp.

( ) ( ^{− −} ) ( ) ^{= −}

d

(

^{j( + )}

) ⁼

^{−j N2−1}

j d j d j

d

e P e P e P e e

P

^{ω ω ω ω π ω} . (1.39)

Pro přenosové funkce platí

( ) −

_h

( ) =

_d

( ) −

_d

( ) − =

^−N2−1

d

z P z P z P z z

P

. (1.40)

Připomeneme-li podmínky pro věrnou rekonstrukci (1.20)

τ

= −

−

−

=

−

−

−F z H z P z P z z

z H z

F_d( ) _d( ) _d( ) _d( ) _d( ) _d( ) 2 ,

vidíme, že jim půlpásmové filtry vyhovují: jejich maximální přenos musí být 2 a fázové zpoždění τ=(N-1)/2. Poznamenejme, že rozkladové a rekonstrukční filtry musí splňovat podmínky

) ( ) ( ) ( )

( ) ( )

(z F z H z a P z F z H z

P_d = _{d d h} = _{h h} .

Podívejme se nyní na návrhy půlpásmových filtrů, odtud se pak dostaneme k návrhům různých typů rozkladových a rekonstrukčních filtrů.

1.2.2.3 Návrh vycházející z Lagrangeova interpolačního vzorce Z interpolačního vzorce podle Lagrangea, který lze nalézt např. v [17],

2 ..., 1 2 ,

; 1 2 2

1

,

−

− −

− = Π −

⎟=

⎠

⎜ ⎞

⎝

⎛ −

+ =− ≠

τ τ τ

α τ

^τ

τ k

i k

i

k i i k

(1.41)

(pro liché τ), vyplývají pro výpočet interpolované hodnoty následující váhy pro (τ+1)/2 vzorků z každé strany od hledané hodnoty:

a) pro τ=1: {1/2, 1/2} … interpolovaná hodnota jako průměr dvou sousedních hodnot;

b) pro τ=3: {-1/16, 9/16, 9/16, -1/16} … výpočet ze čtyř nejbližších hodnot;

c) pro τ=5: {3/256, -25/256, 150/256, 150/256, -25/256, 3/256} … atd.

S uvedenými vahami korespondují první tři (pro τ=1, 3, 5) nejjednodušší půlpásmové dolní propusti

τP_d, které by se aplikovaly na posloupnost expandovanou s faktorem 2 (tj. na diskrétní signál se zdvojnásobenou vzorkovací frekvencí proložením vzorky nulových hodnot):

( ) (

^{1 2}

)

1

1 2

2

1 +

₋

+

₋

= z z

z

P

_d , (1.42)

( ) (

^{2 3 4 6}

)

3

1 9 16 9

16

1

_{− − − −}

− + +

+

−

= z z z z

z

P

_d , (1.43)

( ) (

^{2 4 5 6 8 10}

)

5

3 25 150 256 150 25 3

256

1 −

₋

+

₋

+

₋

+

₋

−

₋

+

₋

= z z z z z z

z

P

_d . (1.44)

K významu indexu τ u přenosové funkce ^τP_d(z) půlpásmové dolní propusti můžeme také uvést, že udává velikost fázového zpoždění τ=(N-1)/2, kde N je počet vzorků impulsní charakteristiky.

Obecný vzorec pro vyjádření přenosové funkce půlpásmové dolní propusti, který lze nalézt např. v [8] [45], má tvar

( )

^{p k k}

k p p

d

z z

k k z p

z z

P ⎟⎟

⎠

⎜⎜ ⎞

⎝

⎟ ⎛ −

⎠

⎜ ⎞

⎝

⎟⎟ ⎛ −

⎠

⎜⎜ ⎞

⎝

⎛ + −

⎟⎟ ⎠

⎜⎜ ⎞

⎝

⎟ ⎛ +

⎠

⎜ ⎞

⎝

= ⎛ +

^{− −}

=

−

∑ ^{1 1} ₂ ¹ ₂

2 1 2

2 1

¹

1 0 1

, (1.45)

odkud je zřejmé, že počet nulových bodů v z=-1 je 2p.

Na obr.1.11 jsou zobrazeny modulové frekvenční charakteristiky uvedených půlpásmových dolních propustí (1.42), (1.43) a (1.44). Se strmostí přechodové části charakteristiky samozřejmě roste i délka impulsní charakteristiky filtru. Za zmínku stojí, že jsou uvedené charakteristiky monotónní.

obr.1.11 Modulové frekvenční charakteristiky půlpásmových dolních propustí ¹Pd, ²Pd a ³Pd.

Rozložení nulových bodů těchto půlpásmových dolních propustí v rovině „z“ (všechny póly filtrů FIR jsou vždy v počátku) je zachyceno na obr.1.12. Plochá modulová charakteristika v oblasti kmitočtů pod ω=π je zajištěna násobným nulovým bodem v z=-1, monotónní tvar charakteristiky nepřipouští žádné jiné nulové body na jednotkové kružnici, všechny zbylé nulové body se nacházejí v pravé polorovině „z“ a jejich poloha zajišťuje plochou modulovou charakteristiku v oblasti kmitočtů těsně nad ω=0.

obr.1.12 Rozložení nulových bodů půlpásmových dolních propustí ¹Pd, ³Pd a ⁵Pd.