Pozn´amky k pˇredmˇetu Biologick´e a akustick´e sign´aly Zbynˇek Koldovsk´y March 8, 2011
1
Z´akladn´ı fakta k pˇredmˇetu: • K udˇelen´ı z´apoˇctu za cviˇcen´ı je nutn´a u´ cˇ ast (max. 2 neomluven´e absence) a odevzd´an´ı u´ loh, jsou-li zadan´e. Zkouˇska z pˇredmˇetu je u´ stn´ı. • Obsah pˇredmˇetu je orientov´an na v´yhradnˇe technickou str´anku problematiky zpracov´av´an´ı biosign´al˚u a akustick´ych sign´al˚u. Medic´ınsk´a problematika zde nen´ı pˇredmˇetem.
2
Pˇr´ıklady biologick´ych a akustick´ych sign´alu˚
1
Rozdˇelen´ı sign´alu˚ podle druhu
1.1
Podle rozmˇeru • jednorozmˇern´e - sledov´an´ı jedn´e veliˇciny, vˇetˇsinou v z´avislosti na cˇ ase • dvourozmˇern´e - obraz • v´ıcekan´alov´e - EEG, EKG,. . . • v´ıcerozmˇern´e - obrazov´y z´aznam v cˇ ase Podle p˚uvodu • elektrick´e • mechanick´e • magnetick´e • chemick´e a dalˇs´ı Podle pˇr´ıcˇ iny vzniku • spont´ann´ı - EEG, EKG • vznikl´e odezvou na dan´e buzen´ı Podle pravidelnosti • opakuj´ıc´ı se - EKG • nepravideln´e - EEG
1.2
Pˇr´ıklady biologick´ych sign´alu˚
ˇ Akˇcn´ı potenci´al bunky • Akˇcn´ı potenci´al je elektrick´y sign´al, kter´y doprov´az´ı mechanickou kontrakci buˇnky stimulovanou elektrick´ym proudem. • Z´akladn´ı komponenta bioelektrick´ych sign´al˚u. • Klidov´y potenci´al → depolarizace → repolarizace. • Akˇcn´ı potenci´al se sˇ´ıˇr´ı po axonu smˇerem k synaps´ım (nervov´e buˇnky) nebo po svalov´em vl´aknu. • Napˇr. nervov´e buˇnky generuj´ı max. asi 1000 impuls˚u za sekundu, srdeˇcn´ı mohou generovat dalˇs´ı impuls aˇz po 150-300ms • Akˇcn´ı potenci´al dan´e buˇnky je vˇzdy stejn´y, stejn´a napˇet´ı (“vˇse, nebo nic”). • (obr. 1) Elektroneurogram (ENG) • Elektroneurogram je grafick´e vyj´adˇren´ı stimulu a n´asledn´eho souhrnn´eho akˇcn´ıho potenci´alu v nervu. • Amplituda do 10 mV, frekvence do 1 kHz • M˚uzˇ e b´yt pouˇzit k mˇeˇren´ı rychlosti sˇ´ıˇren´ı stimulu nebo akˇcn´ıho potenci´alu v nervu pomoc´ı napˇr. jehlov´ych elektrod. 3
• Typick´e rychlosti sˇ´ıˇren´ı – 45-70 m/s v nervov´ych vl´aknech – 0.2-0.4 m/s v srdeˇcn´ım svalu ˇ – 0.03-0.05 m/s ve zpoˇzdovac´ ıch vl´aknech mezi srdeˇcn´ı pˇreds´ın´ı a komorou • Nˇekter´a nervov´a onemocnˇen´ı mohou b´yt prov´azena zmˇenou rychlosti sˇ´ıˇren´ı. Elektromyogram (EMG) • Z´aznam elektrick´e aktivity, kter´a vych´az´ı ze svalov´ych vl´aken. • Sloˇzen´ı svalu: Sval → pohybov´e jednotky → kaˇzd´a jednotka se skl´ad´a z: nervov´a buˇnka motoneuron, jej´ı nervov´y v´ybˇezˇ ek a svalov´a vl´akna inervovan´a t´ımto v´ybˇezˇ kem. • Pohybov´e jednotky velk´ych sval˚u pro hrubˇs´ı pohyby obsahuj´ı stovky vl´aken, zat´ımco svaly pro pˇresn´e pohyby m´enˇe. • Pohybov´a jednotka stimulovan´a nervov´ym sign´alem se stahuje a generuje elektrick´y sign´al, kter´y je souhrnem akˇcn´ıch potenci´al˚u bunˇek (SMUAP - single-motor-unit action potential). • SMUAP b´yv´a dvou aˇz tˇr´ıf´azov´y, 3-15ms, 100-300µV. • EMG je souˇcet cˇ i ˇretˇezec SMUAP˚u. • Tvar EMG (SMUAP) z´avis´ı na typu elektrod, jejich pozici a projekci elektrick´eho pole na elektrody. M˚uzˇ e b´yt ovlivnˇen r˚uzn´ymi nemocemi. • (obr. 2) Elektrokardiogram (EKG) • Elektrick´y projev srdeˇcn´ı aktivity z´ıskan´y povrchov´ymi elektrodami. • Amplituda do 1 mV, frekvence do 150 aˇz 500 Hz • Srdce → dvˇe pˇreds´ınˇe a dvˇe komory. • Odpoˇcinkov´a a pln´ıc´ı f´aze → diastola, kontrakce → systola. • Srdeˇcn´ı tep (norm´alnˇe 70 bpm, max. aˇz 200 bpm) je kontrolov´an SA (sinoatri´aln´ım) uzlem. • Sekvence ud´alost´ı: 1. Sign´al z SA uzlu. 2. Pomal´e sˇ´ıˇren´ı elektrick´e aktivity srdeˇcn´ım svalem zp˚usobuj´ıc´ı pomal´y stah pˇreds´ınˇe (P vlna). 3. PQ segment, zpoˇzdˇen´ı ˇr´ızen´e AV uzlem. Pauza pom´ah´a k dokonˇcen´ı toku krve z pˇreds´ın´ı do komor. 4. His˚uv svazek a Purkyˇnova vl´akna sˇ´ıˇr´ı vysokou rychlost´ı stimul do komor. 5. Rychl´a depolarizace (kontrakce) komor, QRS komplex, dvou nebo tˇr´ıf´azov´y, amplituda 1 mV, trv´an´ı kolem 80ms. 6. Relativnˇe dlouh´y akˇcn´ı potenci´al svalov´ych bunˇek srdeˇcn´ıch komor (300-350 ms) zp˚usobuje isoelektrick´y ST segment d´elky 100-120 ms po QRS. 7. Repolarizace (uvolnˇen´ı) komor zp˚usobuje pomalou T vlnu, 120-160 ms. • (obr. 3) • Zmˇeny v˚ucˇ i norm´aln´ımu srdeˇcn´ımu rytmu se naz´yvaj´ı arytmie, existuje jich velk´a ˇrada. • Mˇeˇren´ı EKG:
4
– Standardn´ı 12ti svodov´e EKG m´a 4 konˇcetinov´e a 6 hrudn´ıch svod˚u. – Referenˇcn´ı elektroda: prav´a noha. – Bipol´arn´ı svody na lev´e a prav´e ruce a lev´e noze se postupnˇe znaˇc´ı I, II a III. – Wilsonova centr´aln´ı svorka: kombinace svod˚u na lev´e a prav´e ruce a lev´e noze a je pouˇzita jako referenˇcn´ı pro hrudn´ı svody. – Zes´ılen´e unipol´arn´ı svody aVR a aVL (L a R ruka), aVF (L noha) jsou z´ısk´any odpov´ıdaj´ıc´ı elektrodou vztaˇzen´e k Wilsonovˇe centr´aln´ı svorce – Tzv. Einthoven˚uv troj´uheln´ık. – Hrudn´ı svody: V1 , . . . , V6 . • (obr. 4 a 5) • Fet´aln´ı EKG: EKG plodu v tˇele matky, amplituda do 100µV, frekvence do 150 Hz Pulzn´ı saturace (SpO2) • Pulsn´ı oxymetrie je neinvazivn´ı, opticko-elektronick´a metoda, mˇeˇr´ıc´ı v pulsuj´ıc´ı cˇ a´ sti arteri´aln´ı krve procentu´aln´ı saturaci arteri´aln´ı krve kysl´ıkem. • SpO2 ud´av´a procentu´aln´ı pomˇer mezi mnoˇzstv´ım hemoglobinu pr´avˇe nav´azan´eho s O2 a hemoglobinu schopn´eho v´azat O2 . ˇ y zp˚usob mˇeˇren´ı: prosvˇecov´an´ı prstu pacienta. Intenzita pr˚uchoz´ıho svˇetla se mˇen´ı s hustotou arteri´aln´ı • Cast´ krve. Elektroencefalogram (EEG) • Sn´ım´a elektrickou aktivitu mozku standardnˇe 10-20 elektrodami na hlavˇe. • Velmi slab´y sign´al, ˇra´ dovˇe des´ıtky µV. • Velmi sloˇzit´a data → akˇcn´ı potenci´aly milion˚u neuron˚u. Souˇcet rozliˇcn´e aktivity mnoha mal´ych z´on kortik´aln´ı vrstvy pod elektrodou. • Data zpracov´any 75 Hz low-pass filter. • Nejzn´amˇejˇs´ı periodick´e aktivity v EEG (rytmy): – Delta: 0.5 ≤ f < 4Hz a Theta: 4 ≤ f < 8Hz → identifik´atory r˚uzn´ych sp´ankov´ych stadi´ı. V bdˇel´em stavu u dospˇel´ych jsou patologick´ym pˇr´ıznakem. – Alfa: 8 ≤ f < 13Hz → u dospˇel´ych v occipit´aln´ı oblasti hlavy pˇri zavˇren´ych oˇc´ıch. – Beta: f > 13Hz → bˇezˇ n´a aktivita v napjat´em cˇ i nˇejak uzkostliv´em stavu. Chyb´ı-li → abnormalita. ERP - Event-Related Potencial, Evokovan´e potenci´aly • Technika pˇri kter´e se pouˇz´ıv´a pˇresnˇe definovan´a stimulace pro excitaci pˇr´ısluˇsn´eho mozkov´eho (nervov´eho) analyz´atoru a n´aslednˇe se mˇeˇr´ı EEG nebo ENG reakce. • Stimulace: zrakov´e, sluchov´e, elektrick´e, somatosenzorick´e. • ERP mohou m´ıt kr´atkou nebo dlouhou latenci (zpoˇzdˇen´ı reakce). Prvn´ı mohou z´aviset na vlastnostech stimulace, druh´e zase v´ıce na podm´ınk´ach, za kter´ych jsou poˇrizov´any. PCG - akustick´y srdeˇcn´ı pulz, t´ezˇ fonokardiogram • Nejzn´amˇejˇs´ı biosign´al - poslouch´an fonendoskopem (stetoskopem). • Akustick´y cˇ ili mechanick´y sign´al.
5
• Zvukov´e projevy tlaku v cel´em kardiovaskul´arn´ım syst´emu. Nikoliv pouze pohyb srdeˇcn´ı chlopnˇe jak se p˚uvodnˇe myslelo. Chov´a se tedy jako nafouknut´y bal´onek. • Mˇeˇren´ı mikrofony, tlakov´e senzory, cˇ idla zrychlen´ı apod. • Diagnostika: sˇelesty. • Hlavn´ı cˇ a´ sti: S1 a S2. S1: kontrakce komor a) pohyb krve smˇerem ke komor´am utˇesnˇen´e AV chlopn´ı b) n´ahl´e napˇet´ı uzavˇren´e AV chlopnˇe zadrˇzuj´ıc´ı krev c) otevˇren´ı chlopn´ı a vypumpov´an´ı krve z komor, S2 po systolick´e pauze: uzavˇren´ı p˚ulmˇes´ıcov´ych chlopn´ı a) srdeˇcn´ı b) plicn´ı. Nˇekdy m˚uzˇ e b´yt slyˇset i S3 odpov´ıdaj´ıc´ı n´ahl´emu ukonˇcen´ı pln´ıc´ı f´aze komory (hlubok´e frekvence). • Frekvenˇcn´ı rozsah S1: 15-800 Hz, S2: 25-800 Hz, S3 a S4: 10-40 Hz • (obr. 6) CP - karotick´y pulz • Sign´al vznik´a mˇeˇren´ım tlaku v krkavici na krku v m´ıstˇe, kde proch´az´ı bl´ızko povrchu tˇela. • Sign´al je vhodn´ym doplˇnkem PCG. • Vznik´a vypuzen´ım krve z lev´e komory do aorty. Vrcholem vznik´a vlna P, n´asledn´a rovina kˇrivky T zp˚usoben´a odrazem pulzu z horn´ı cˇ a´ sti tˇela, uz´avˇer srdeˇcn´ı chlopnˇe - stup´ınek D a odraz pulzu z doln´ı cˇ a´ sti tˇela DW. • (obr. 6) Dalˇs´ı • EOG → Elektrookulogram, amplituda do 10 mV, frekvence do 100 Hz • EGG → elektrogastrogram. Depolarizace a repolarizace hladk´e svaloviny ze stˇredn´ı cˇ a´ sti zˇ aludku. Ne kaˇzd´a aktivita vˇsak odpov´ıd´a kontrakc´ım zˇ aludku. Diagnostick´y potenci´al EGG nen´ı zcela potvrzen. • VMG → vibromyogram. Mechanick´y projev kontrakce sval˚u (vibrace) doprov´azej´ıc´ı EMG. • VAG → vibroarthrogram. Vibrace kolenn´ıho kloubu bˇehem pohybu. • OAE → otoakustick´a emise. Akustick´a energie vyd´avan´a Cortiho org´anem (hlem´yzˇ dˇ stˇredn´ıho ucha): spont´ann´ı nebo odezva na akustickou stimulaci.
1.3
C´ıle anal´yzy biologick´ych sign´alu˚ • sbˇer informac´ı • diagnostika • l´ecˇ ba, ˇr´ızen´ı • vyhodnocov´an´ı (napˇr. u´ cˇ innosti l´ecˇ by)
˚ Zpusoby z´ısk´av´an´ı dat • invazivn´ı - neinvazivn´ı (um´ıstˇen´ı cˇ idel uvnitˇr nebo vnˇe tˇela) • aktivn´ı - pasivn´ı (stimulace nebo bez) ˚ zit´e vlastnosti mˇerˇ´ıc´ı techniky Duleˇ • izolace subjektu - odstranˇen´ı nebezpeˇc´ı u´ razu elektˇrinou • rozsah operac´ı, parametr˚u, sign´al˚u 6
• citlivost • linearita • hystereze - z´avislost stavu na stavech pˇredchoz´ıch m˚uzˇ e zp˚usobovat odchylku v mˇeˇren´ı • frekvenˇcn´ı odezva - r˚uzn´a citlivost na frekvence • stabilita - opakovatelnost pokus˚u • pomˇer sign´al/ˇsum → signal-to-noise ratio (SNR) • pˇresnost
1.4
Hudebn´ı sign´aly
Zdroje sign´al˚u: • senzory – mikrofony - dynamick´e, kondenz´atorov´e, lampov´e – sn´ımaˇce - magnetick´e, piezoelektrick´e • gener´atory - syntez´ator
7
Frekvenˇcn´ı anal´yza sign´alu˚ a filtry - opakov´an´ı
2 2.1
Diskr´etn´ı Fourierovy transformace • DTFT sign´alu x[n], n ∈
Z X F (θ) =
+∞ X
x[n]e−inθ ,
θ ∈ (−π, π]
n=−∞
• Je-li x[n] koneˇcn´y d´elky N , n = 0, . . . , N − 1, pak X F (θ) =
N −1 X
x[n]e−inθ ,
θ ∈ (−π, π].
n=0
• DFT je definovan´a pro koneˇcn´y sign´al N hodnotami X F (θ) rovnomˇernˇe rozm´ıstˇen´ymi na intervalu θ ∈ [0, 2π), tedy v bodech θk = k 2π N , k = 0, . . . , N − 1, tedy X[k] =
N −1 X
2π
x[n]e−ink N ,
k = 0, . . . , N − 1.
n=0
• Rychl´y v´ypoˇcet DFT: FFT • Prodlouˇzen´ı sign´alu x[n] o M nul: DTFT X F (θ) z˚ust´av´a stejn´a, ale DFT se mˇen´ı t´ım, zˇ e vyhodnocuje DTFT v N + M bodech, tedy dojde k jemnˇejˇs´ımu dˇelen´ı intervalu [0, 2π). • DFT vyuˇz´ıv´ame k neparametrick´emu odhadu spektra, nejˇcastˇeji metodou ok´enek (Welchova metoda).
2.2
Parametrick´e odhady spektra • Pˇredpokl´ad´ame, zˇ e sign´al m´a vlastnosti uˇcit´eho modelu a odhadujeme parametry modelu, kter´e potom urˇcuj´ı spektrum.
• Autoregresn´ı model ˇra´ du p: sign´al x[n] splˇnuje x[n] = −
p X
ak x[n − k] + u[n],
k=1
kde u[n] je b´ıl´y sˇum s rozptylem (energi´ı) σ 2 . Sign´al x[n] tedy “ch´apeme” jako b´ıl´y sˇum filtrovan´y allpole filtrem. Spektr´aln´ı charakteristika filtru (urˇcen´eho koeficienty a1 , . . . , ap ) je tedy stejn´a jako spektrum sign´alu. Spektrum sign´alu je tedy rovno σ2 p SAR (ω) = Pp , | k=0 ak e−iωk |2 kde a0 = 1. Koeficienty a1 , . . . , ap lze odhadovat mnoha zp˚usoby (nejˇcastˇeji autokorelaˇcn´ı metodou s pouˇzit´ım Levinson-Durbinova algoritmu). • Dalˇs´ı modely: MA (moving averages), ARMA • V´yhodou parametrick´ych odhad˚u je vˇetˇs´ı pˇresnost a schopnost odhadovat spektrum i z kr´atk´eho z´aznamu sign´alu (odhadujeme mal´y poˇcet parametr˚u). Nev´yhodou je z´avislost na pˇredpokladu, zˇ e sign´al dan´y model splˇnuje.
8
2.3
Obecn´e vlastnosti LTI syst´emu˚ a filtru˚
LTI syst´em (diskr´etn´ı) se vstupn´ım sign´alem x[n] a v´ystupn´ım sign´alem y[n]: P+∞ • y[n] = k=−∞ h[k]x[n − k], kde h[n] je impulzn´ı odezva syst´emu • Pˇrenosov´a funkce: Z-transformace
+∞ X
H(z) =
h[k]z −k ,
k=−∞
kde z ∈
C, konkr´etnˇeji z ∈ ROC, kde ROC (Region-of-convergence) je oblast konvergence ˇrady
• DTFT transformace: H F (θ) = H(eiθ ) =
P+∞
k=−∞
h[k]e−ikθ , θ ∈ (−π, π]
• Vzorkovac´ı teor´em: θ = ωTs = 2π ffs Obecn´e vlastnosti: • stabiln´ı (ROC obsahuje jednotkovou kruˇznici - existuje DTFT) × nestabiln´ı • kauz´aln´ı (ROC obsahuje ∞) × nekauz´aln´ı • h[n] koneˇcn´a (FIR) × nekoneˇcn´a (IIR) Filtr: Stabiln´ı LTI syst´em naz´yv´ame filtr.
2.4
Specifikace filtru H F (θ) = |H F (θ)|eiψ(θ) ,
kde |H F (θ)| je magnituda filtru (v decibelech 20 log10 (|H F (θ)|)) a ψ(θ) je f´aze ( arctan2(=[H F (θ)], <[H F (θ)]) H F (θ) 6= 0 ψ(θ) = nedefinovan´a H F (θ) = 0, kde arctan2(y, x) = α, α ∈ (−π, π], kde plat´ı cos(α) = √
2.5
x x2 +y 2
a sin(α) = √
y . x2 +y 2
Kauz´aln´ı filtry s re´aln´ymi koeficienty a racion´aln´ı pˇrenosovou funkc´ı • Zkratka: RCSR filtr (real, causal, stable, rational) • Vztah mezi vstupn´ım sign´alem x[n] a v´ystupn´ım sign´alem y[n] pops´an diferenˇcn´ı rovnic´ı y[n] + a1 y[n − 1] + · · · + ap y[n − p] = b0 x[n] + b1 x[n − 1] + · · · + bq x[n − q] V Matlabu pˇr´ıkaz filter • Pˇrenosov´a funkce je racion´aln´ı v promˇenn´e z −1 H(z) =
b0 + b1 z −1 + · · · + bq z −q b(z −1 ) = , −1 −p 1 + a1 z + · · · + ap z a(z −1 )
kde ap 6= 0 a bq 6= 0. Jelikoˇz je filtr kauz´aln´ı (vypl´yv´a z tvaru diferenˇcn´ı rovnice), ROC je urˇcena nejvˇetˇs´ım p´olem v absolutn´ı hodnotˇe (zaˇca´ tek ROC) a nekoneˇcnem. Proto chceme-li, aby byl z´aroveˇn i stabiln´ı, mus´ı ROC obsahovat jednotkov´y kruh a tedy nejvˇetˇs´ı p´ol mus´ı m´ıt menˇs´ı velikost neˇz jedna. Jin´ymi slovy vˇsechny p´oly mus´ı leˇzet uvnitˇr jednotkov´eho kruhu. • P´oly pˇrenosov´e funkce H(z) jsou koˇreny polynomu a(z −1 ), kter´e z´aroveˇn nejsou koˇrenem b(z −1 ), anebo jsou-li koˇrenem b(z −1 ) potom niˇzsˇ´ıho ˇra´ du neˇz v a(z −1 ). • Analogicky definujeme nuly H(z) jako koˇreny b(z −1 ). 9
• Plat´ı: DTFT transformaci RCSR filtru lze vyj´adˇrit ve tvaru H F (θ) = A(θ)eiφ(θ) , kde A(θ) je re´aln´a (kladn´a i z´aporn´a) funkce, kterou naz´yv´ame amplitudou a φ(θ) je f´aze spojit´a v promˇenn´e θ
Specifikace filtru˚ podle magnitudy
2.6
• Ide´aln´ı charakteristiky nelze dos´ahnout pomoc´ı RCSR filtru • Pˇr. Low-pass (LP) filtr: rozsah v propustn´em p´asmu 1 ± δ, hranice propustn´eho p´asma θP , hranice nepropustn´eho p´asma θS , tolerance v nepropustn´em p´asmu δS • Analogicky definujeme high-pass, band-pass, band-stop, multi-band
2.7
Specifikace filtru˚ podle f´aze
2.7.1
Line´arn´ı f´aze
• Pˇr. h[n] = δ[n − L], H(z) = z −L , H F (θ) = e−iθL =⇒ φ(θ) = −Lθ • Myˇslenka: V p´asmu, kde je φ(θ) line´arn´ı je vstupn´ı sign´al pouze zpoˇzdˇen o L vzork˚u, tedy nemˇen´ı se jeho ”tvar” • Pro obecn´y filtr definujeme f´azov´e zpoˇzdˇen´ı τp (θ) = −
φ(θ) θ
• Filtr s line´arn´ı f´az´ı m´a konstantn´ı f´azov´e zpoˇzdˇen´ı 2.7.2
Zobecnˇen´a line´arn´ı f´aze
• Zobecnˇen´a line´arn´ı f´aze: φ(θ) = φ0 − θτg , kde φ0 je konstanta a τg naz´yv´ame skupinov´e zpoˇzdˇen´ı • Pro obecn´e filtry definujeme skupinov´e zpoˇzdˇen´ı jako funkci θ (je-li φ(θ) diferencovateln´a) τg (θ) = −
dφ(θ) dθ
• Filtr se zobecnˇenou line´arn´ı f´az´ı m´a konstantn´ı skupinov´e zpoˇzdˇen´ı • Pro RCSR filtr s line´arn´ı (zobecnˇenou) f´az´ı plat´ı – M´a koneˇcnou impulsn´ı odezvu (je FIR) – F´azov´e nebo skupinov´e zpoˇzdˇen´ı je rovno polovinˇe jeho ˇra´ du N – Filtr je budˇ ∗ symetrick´y, tedy h[n] = h[N − n] a plat´ı φ0 = 0 ∗ antisymetrick´y, tedy h[n] = −h[N − n] a plat´ı φ0 =
10
π 2
2.7.3
Filtry s minim´aln´ı f´az´ı
• Je-li β nulou RCSR filtru, pak lze filtr upravit tak, zˇ e nulu β zamˇen´ıme za β tudu filtru
−1
aniˇz bychom zmˇenili magni-
• Plat´ı: Filtr, kter´y m´a vˇsechny nuly uvnitˇr jednotkov´eho kruhu, m´a menˇs´ı skupinov´e zpoˇzdˇen´ı neˇz filtry, kter´e maj´ı stejnou magnitudu ale nˇekter´e nuly vnˇe jednotkov´eho kruhu • Takov´y filtr m´a tedy minim´aln´ı skupinov´e zpoˇzdˇen´ı, avˇsak naz´yv´a se filtr s minim´aln´ı f´az´ı • Vytvoˇren´ı filtru s minim´aln´ı f´az´ı 1. Navrhneme filtr 2. Nalezneme jeho nuly 3. Zamˇen´ıme jeho nuly tak, aby byly vˇsechny uvnitˇr jednotkov´eho kruhu a adekv´atnˇe nastav´ıme zes´ılen´ı • M´a-li FIR filtr (zobecnˇenou) line´arn´ı f´azi, pak plat´ı H(z −1 ) = ±z N H(z) =⇒ je-li β nulou tohoto filtru, pak nutnˇe tak´e β −1 je jeho nulou =⇒ m´a-li m´ıt tento filtr z´aroveˇn minim´aln´ı f´azi, pak nutnˇe leˇz´ı vˇsechny jeho nuly na jednotkov´em kruhu 2.7.4
All-pass filtry
• All-pass filtr splˇnuje |H F (θ)| = 1, tedy sign´al po pr˚uchodu t´ımto filtrem m´a zmˇenˇenou pouze f´azi • Pˇr´ıklad all-pass filtru prvn´ıho ˇra´ du: H(z) =
a − z −1 , 1 − az −1
|a| < 1
Podm´ınka |a| < 1 je nutn´a, aby mˇel filtr p´ol a uvnitˇr jednotkov´eho kruhu (kauz´aln´ı a stabiln´ı). M´a-li b´yt z´aroveˇn all-pass, pak k dan´emu p´olu mus´ı m´ıt nulu a−1 , kter´a je t´ım p´adem vnˇe jednotkov´eho kruhu. All-pass filtr m´a tedy maxim´aln´ı f´azi (vˇsechny nuly vnˇe jednotkov´eho kruhu) • RSCR all-pass filtr se nutnˇe skl´ad´a z nˇekolika all-pass filtr˚u prvn´ıho ˇra´ du H(z) =
p Y αk − z −1 1 − αk z −1
k=1
11
Metody navrhov´an´ı FIR filtru˚ s line´arn´ı f´az´ı
3
Typy FIR filtr˚u s line´arn´ı (zobecnˇenou) f´az´ı (jin´e moˇznosti FIR filtr˚u nejsou): • Typ I. τg = M , φ0 = 0, ˇra´ d N sud´y, h[n] symetrick´a, hod´ı se na LP, HP, BP, BS a multiband • Typ II. τg = M + 0.5, φ0 = 0, ˇra´ d N lich´y, h[n] symetrick´a, hod´ı se na LP, BP • Typ III. τg = M , φ0 = π2 , ˇra´ d N sud´y, h[n] antisymetrick´a, hod´ı se na aproximov´an´ı diferenci´ator˚u a Hilbertov´ych transform´ator˚u • Typ IV. τg = M + 0.5, φ0 = Hilbertov´ych transform´ator˚u
π 2,
ˇra´ d N lich´y, h[n] antisymetrick´a, hod´ı se na aproximov´an´ı diferenci´ator˚u a
Filtry s minim´aln´ı f´az´ı: • nuly pouze na jednotkov´em kruhu: tzv. NOTCH a COMB filtry
3.1
Filtr s ide´aln´ı charakteristikou • Amplitudov´a charakteristika ( A(θ) =
±1 0
|θ| je v propustn´em p´asmu jinak
• F´azi si m˚uzˇ eme zvolit → chceme (zobecnˇenou) line´arn´ı, tedy obecnˇe π
1
eiφ(θ) = ei(µ 2 − 2 N θ) •
– µ = 0 ⇒ typ I. nebo II. – µ = 1 ⇒ typ III. nebo IV.
ˇ nulovou (ei0 = 1) nebo π (ei π2 = i). V • Pozn´amka: Line´arn´ı f´azi ide´aln´ıho filtru v podstatˇe vol´ıme budto 2 prvn´ım pˇr´ıpadˇe jsou A(θ) a h[n] symetrick´e a v druh´em antisymetrick´e (ovˇsem h[n] kolem n = 0). Pˇr´ıd´an´ım N e−i 2 θ , tj. vyn´asoben´ım H F (θ), provedeme de facto posun stˇredu symetrie h[n] o N2 vzork˚u doprava (aby zkr´acen´y filtr byl kauz´aln´ı)
3.2
N´avrhy filtru˚ metodou zkracov´an´ı ide´aln´ı impulzn´ı odezvy
1. Zvol´ıme A(θ) 2. Zvol´ıme typ filtru (I., II.,. . . ) 3. Zvol´ıme ˇra´ d filtru. Pak
π
N
H F (θ) = A(θ)ei(µ 2 − 2 θ) 4. Spoˇcteme ide´aln´ı impulsn´ı odezvu inverzn´ı DTFT Z π N π 1 hi [n] = A(θ)ei(µ 2 − 2 θ) eiθ dθ 2π −π 5. Zkr´at´ıme odezvu
( hi [n] 0 ≤ n ≤ N h[n] = 0 jinak
• Pozn´amka: Vˇetˇsinou je A(θ) po cˇ a´ stech konstantn´ı, takˇze integr´al v inverzn´ı DTFT je jednoduch´y (integrace funkce ex ) 12
• Pˇr´ıklad: LP, HP nebo BP filtr, ( ±1 A(θ) = 0
θ1 ≤ |θ| ≤ θ2 jinak
Zvol´ıme typ I., takˇze N je sud´e, ( F
H (θ) =
N
±e−i 2 θ 0
θ1 ≤ |θ| ≤ θ2 jinak
Ide´aln´ı impulzn´ı odezva tedy je 1 hi [n] = 2π
Z
−θ1
e
θi(n− N 2 )
−θ2
Z
!
θ2
dθ +
e
θi(n− N 2 )
dθ
θ1
N N N N 1 e−θ1 i(n− 2 ) − e−θ2 i(n− 2 ) + eθ2 i(n− 2 ) − eθ1 i(n− 2 ) N 2πi(n − 2 ) 1 = 2i sin(θ2 (n − N/2)) − 2i sin(θ1 (n − N/2)) N 2πi(n − 2 ) 1 sin(θ2 (n − N/2)) − sin(θ1 (n − N/2)) = π(n − N/2)
=
Pˇr: LP filtr =⇒ θ1 = 0 ( h[n] =
sin(θ2 (n−N/2)) π(n−N/2) θ2 π
0 ≤ n ≤ N, n 6= n=
N 2
N 2
• Roste-li ˇra´ d filtru N → +∞, pak v´ıme, zˇ e se frekvenˇcn´ı charakteristika zkr´acen´eho filtru bl´ızˇ´ı ide´alu. Jelikoˇz vˇsak ide´aln´ı frekvenˇcn´ı charakteristika nen´ı spojit´a, nast´av´a v bodech jej´ı nespojitosti tzv. Gibbs˚uv jev. Jak je vidˇet v n´asleduj´ıc´ım grafu, kde θ1 = 0 a θ2 = π/4, rozkmit amplitudy filtru z˚ust´av´a v bodˇe nespojitosti st´ale “velk´y” i kdyˇz ˇra´ d filtru zvˇetˇsujeme. 1.4 N=10 N=100 N=1000
1.2
A(θ)
1 0.8 Gibbsuv jev 0.6 0.4 0.2 0 0
3.3
0.5
1
1.5
θ
2
2.5
3
N´avrhy filtru˚ metodou ok´enek
1. Navrhneme toleranˇcn´ı rozsahy θp a θs pro kaˇzd´y pˇrechod z propustn´eho do nepropustn´eho p´asma a ide´aln´ı odezvu definujeme tak, zˇ e hranou pˇrechodu je stˇred (θp + θs )/2 2. Navrhneme impulsn´ı odezvu h[n] metodou zkracov´an´ı 13
3. Uprav´ıme impulsn´ı odezvu ( h[n] =
hi [n]w[n] 0 ≤ n ≤ N 0 jinak
Vlastnosti ok´enkovac´ı funkce w[n]: • symetrick´a a kladn´a =⇒ zachov´av´a (zobecnˇenou) line´arn´ı f´azi filtru • n´asoben´ı w[n] v cˇ asov´e oblasti (ok´enkov´an´ı) znamen´a konvoluci ve frekvenˇcn´ı oblasti: H F (θ) ←−
1 1 {W F ? H F }(θ) = 2π 2π
Z
+∞
W F (λ)H F (θ − λ)dλ
−∞
=⇒ “low-pass filtrace frekvenˇcn´ı charakteristiky H F (θ)” −→ vyhlazen´ı charakteristiky zp˚usoben´e Gibbsov´ym jevem (viz pˇredeˇsl´y a n´asleduj´ıc´ı obr´azek po pouˇzit´ı Hanningova ok´enka) 1.4 N=10 N=100 N=1000
1.2
A(θ)
1 0.8 0.6 0.4 0.2 0 0
0.5
1
1.5
θ
2
2.5
3
• Odhad ˇra´ du filtru pro splnˇen´ı jeho detailn´ıch vlastnost´ı dle poˇzadavk˚u (rozsahy tolerance v propustn´ych a nepropustn´ych p´asmech) vyˇzaduje detailn´ı anal´yzu vlastnost´ı dan´e ok´enkovac´ı funkce (viz. napˇr. v Matlabu spoleˇcn´e pouˇzit´ı pˇr´ıkaz˚u fir1, ok´enka kaiser a odhadu ˇra´ du kaiserord)
14
N´avrhy IIR filtru˚
4
Postup n´avrhu IIR filtr˚u: • N´avrh analogov´eho LP filtru • Transformace na LP, HP, BP, BS analogov´y filtr • Transformace na digit´aln´ı filtr Analogov´e filtry: • Dobˇre prozkouman´a oblast • Zaj´ım´a n´as pˇredevˇs´ım magnituda filtru, f´azov´a charakteristika je u IIR filtr˚u aˇz druhoˇrad´a (f´aze RCSR IIR filtr˚u je neline´arn´ı) ˇ • Vˇetˇsina n´avrh˚u (aˇz na Cebyˇ sev˚uv II. typu) m´a tvar |H F (ω)|2 =
1 , 1 + Λ ωω0
kde Λ(·) > 0 je “mal´a” v propustn´em p´asmu a “velk´a” v nepropustn´em p´asmu • Pozn. Urˇc´ıme-li |H F (ω)|2 , je t´ım jiˇz urˇcen souˇcin Laplaceov´ych transformac´ı filtru H L (s)H L (−s), ale nen´ı jeˇstˇe jednoznaˇcnˇe urˇcena H L (s) (Laplaceova transformace analogov´eho filtru). P´oly souˇcinu H L (s)H L (−s) jsou p´oly H L (s) a jejich obrazy pˇres imagin´arn´ı osu. M´a-li b´yt H L (s) stabiln´ı, vˇsechny jeho p´oly mus´ı b´yt na lev´e polorovinˇe (z´aporn´a re´aln´a cˇ a´ st). T´ımto je potom H L (s) urˇcena jednoznaˇcnˇe.
4.1
˚ filtr Butterworthuv 1
|H F (ω)|2 = 1+
ω ω0
2N ,
kde N je ˇra´ d. •
ω ω0
je monotonn´ı funkce, takˇze i |H F (ω)|2 je monotonn´ı (nem´a zˇ a´ dn´e kmity) - sˇikovn´a vlastnost Butterworthova filtru
• Jelikoˇz s = iω, ω = s/i, dosazen´ım dost´av´ame 1
H L (s)H L (−s) = 1+
s iω0
2N =
1 1 − (−1)N
s ω0
2N
• P´oly filtru jsou sk = ω0 ei
(N +1+2k)π 2N
0 ≤ k ≤ 2N − 1,
,
tedy symetricky kolem imagin´an´ı osy. Prvn´ıch N p´ol˚u je na lev´e polorovinˇe. M´a-li b´yt filtr stabiln´ı, pak H L (s) =
N −1 Y k=1
−sk s − sk
• Filtr nem´a zˇ a´ dn´e nuly ˚ N´avrh Butterworthova analogov´eho LP filtru podle specifick´ych poˇzadavku: • M´ame uˇzivatelem zvolen´y konec propustn´eho p´asma ωp , zaˇca´ tek nepropustn´eho p´asma ωs , ωs > ωp , toleranci v propustn´em p´asmu δp a v nepropustn´em p´asmu δs . Naˇs´ım u´ kolem je zvolit nejniˇzsˇ´ı moˇzn´y ˇra´ d filtru N a frekvenci ω0 , ωp ≤ ω0 ≤ ωs , tak, aby tyto poˇzadavky byly splnˇeny. 15
• Definujeme: s d=
(1 − δp )−2 − 1 δs−2 − 1
ωp ωs
k=
• V´yznam definovan´e hodnoty d lze l´epe vysvˇetlit n´asledovnˇe. Definujeme-li tolerance v decibelech Ap = −10 log10 (1 − δp )2 As = −10 log10 δs2 , pak plat´ı s d=
100.1Ap − 1 . 100.1As − 1
• Z v´yznamu zadan´ych parametr˚u mus´ı platit nerovnosti 1 1+
ωp ω0
2 2N ≥ (1 − δp )
1 1+
ωs ω0
2 2N ≤ δs
´ • Upravami dost´av´ame
2N ωp ≤ (1 − δp )−2 − 1 ω0 2N ωs ≥ δs−2 − 1 ω0
• Slouˇcen´ım tˇechto nerovnost´ı dost´av´ame 2N 2N 1 1 ωs δs−2 − 1 ⇐⇒ = 2 ≥ ωp (1 − δp )−2 − 1 k d • Tedy mus´ı platit ln(1/d) . ln(1/k) N mus´ı b´yt cel´e a proto zvol´ıme nejmenˇs´ı takov´e splˇnuj´ıc´ı tuto podm´ınku. Zb´yv´a zvolit frekvenci ω0 . Pro tuto volbu N lze uk´azat, zˇ e mus´ı platit −1/(2N ) −1/(2N ) ωp (1 − δp )−2 − 1 ≤ ω0 ≤ ωs δs−2 − 1 N≥
Postup n´avrhu Butterworthova analogov´eho LP filtru 1. Spoˇcteme d a k podle zadan´ych parametr˚u l m ln(1/d) 2. N = ln(1/k) (tyto z´avorky znaˇc´ı zaokrouhlov´an´ı nahoru) 3. Zvol´ıme ω0 tak, aby splˇnovala v´ysˇe uvedenou nerovnost. 4. Spoˇcteme p´oly sk = ω0 ei
(N +1+2k)π 2N
,
0≤k ≤N −1
5. Filtr je urˇcen Laplaceovou transformac´ı H L (s) =
N −1 Y k=1
16
−sk . s − sk
4.2
Dalˇs´ı analogov´e filtry ˇ • Cebyˇ sev˚uv I. typu: |H F (ω)|2 =
1
1+
, 2 TN2 ( ωω0 )
ˇ kde TN (x) je Cebyˇ sev˚uv polynom N -t´eho ˇra´ du. ˇ • Cebyˇ sev˚uv II. typu: F
2 TN2 ( ωω0 )
2
|H (ω)| =
1 + 2 TN2 ( ωω0 )
• Eliptick´y: |H F (ω)|2 =
1
1+
2 ( ω ), 2 R N ω0
ˇ kde RN (x) je Cebyˇ sevova racion´aln´ı funkce N -t´eho stupnˇe. Vlastnosti: Butterworth monotonn´ı monotonn´ı zˇ a´ dn´e
propustn´e p. nepropustn´e p. nuly
N≥
ˇra´ d
ln(1/d) ln(1/k)
ˇ Cebyˇ sev I. RZ monotonn´ı zˇ a´ dn´e arccosh (1/d) N ≥ arccosh(1/k)
ˇ Cebyˇ sev II. monotonn´ı RZ arccosh(1/d) N ≥ arccosh (1/k)
Eliptick´y RZ RZ N≥
K(k2 )K(1−d2 ) K(1−k2 )K(d2 )
• RZ = rovnomˇernˇe zvlnˇen´y R π/2 • K(m) = 0 (1 − m sin2 (x))−1/2 dx, v Matlabu funkce ellipke
4.3
Transformace LP filtru na LP, HP, BP a BS • Mˇejme racion´aln´ı funkci f (·) slouˇz´ıc´ı k pˇrechodu od promˇenn´e s k s˜ s = f (˜ s) • Ve frekvenc´ıch s = iω, s˜ = i˜ ω , tzn. ω = −if (i˜ ω) • Laplaceova a Fourierova transformace transformovan´eho filtru ˜ L (˜ H s) = H L (s)|s=f (˜s) ˜ F (˜ H ω ) = H F (ω)|ω=−if (i˜ω) • Podstata myˇslenky: je-li nov´a promˇenn´a ω ˜ v propustn´em p´asmu ztransformovan´eho filtru, pak transformace f (·) mus´ı fungovat tak, aby ω byla v propustn´em p´asmu p˚uvodn´ıho filtru a naopak. • f (·) mus´ı zachov´avat stabilitu a b´yt co nejniˇzsˇ´ıho ˇra´ du, aby nezvyˇsovala ˇra´ d ztransformovan´eho filtru.
Pouˇz´ıvan´e transformace: • LP→LP: s =
s˜ ωc ,
ω=
ω ˜ ωc
• LP→HP: s =
ωc s˜ ,
ω = − ωω˜c
• LP→BP: s =
s˜2 +ωl ωh s˜(ωh −ωl ) ,
(roztahuje cˇ i ztahuje frekv. osu podle ωc )
ω=
ω ˜ 2 −ωl ωh ω ˜ (ωh −ωl ) ,
kde ωh > ωl
17
• LP→BS: s =
s˜(ωh −ωl ) s˜2 +ωl ωh ,
ω=
ω ˜ (ωh −ωl ) ωl ωh −˜ ω2
Pˇr´ıklad postupu n´avrhu HP filtru z LP filtru: • Volba ωc > 0, napˇr. ωc = 1 • Specifikace HP filtru ω ˜p, ω ˜ s , δ˜p , δ˜s pˇrevedeme na δp = δ˜p , δs = δ˜s , ωp =
ωc ω ˜p ,
ωs =
ωc ω ˜s
• Navrhneme LP filtr podle ωp , ωs , δp , δs ˜ L (˜ • Z H L (s) vytvoˇr´ıme H s)
4.4
Digitalizace analogov´eho filtru
Uvedeme zde pro u´ plnost a pochopen´ı podstaty tˇri metody digitalizace, z nichˇz prvn´ı dvˇe se prakticky nepouˇz´ıvaj´ı. 4.4.1
Metoda vzorkov´an´ı impulzn´ı odezvy
1. Vypoˇcteme inverzn´ı Laplaceovu transformaci navrˇzen´e pˇrenosov´e funkce H L (s), cˇ´ımˇz z´ısk´ame impulzn´ı odezvu analogov´eho filtru h(t) 2. Vzorkujeme h(t), h[n] = Ts h(nTs ) 3. Spoˇcteme Z-transformaci h[n], cˇ´ımˇz odvod´ıme koeficienty ak a bk digit´aln´ıho filtru pro jeho realizaci Vlastnosti: • Zachov´av´a racionalitu pˇrenosov´e funkce, tedy realizovatelnost filtru. • M´a-li jmenovatel H L (s) vyˇssˇ´ı ˇra´ d neˇz cˇ itatel, pak metoda zachov´av´a ˇra´ d i stabilitu. • Nezachov´av´a obecnˇe poˇzadovan´e vlastnosti filtru. • Aproximace je z´avisl´a na Ts • Analogick´e vlastnosti souvisej´ıc´ı se vzorkovac´ım teor´emem → aliasing 4.4.2
Metoda zpˇetn´e diference
Myˇslenka: H L (s) = s m´a v Laplaceovˇe transformaci v´yznam derivace. Pro pˇrevod do Z-transformace ji nahrad´ıme diferenc´ı, tedy provedeme substituci s←
1 − z −1 Ts
Vlastnosti: • Zachov´av´a ˇra´ d filtru i stabilitu. • Aby metoda pˇren´asˇela poˇzadovan´e vlastnosti analogov´eho filtru, kter´ych jsme dos´ahli jeho n´avrhem, na filtr digitalizovan´y, potˇrebujeme, aby se transformac´ı mapovala imagin´arn´ı osa s-roviny (H L (s) na imagin´arn´ı ose je rovna Fourierovˇe transformaci) na jednotkov´y kruh z-roviny. Sledujme tedy co dˇel´a metoda zpˇetn´e diference: Komplexn´ı cˇ´ıslo s lze ps´at jako s = σ + iω. Pro σ = 0 d´av´a transformace 1 1 + iωTs 1 1 1 ⇒z− = ⇒ z − = z= 1 − iωTs 2 2(1 − iωTs ) 2 2 Imagin´arn´ı osa se tedy mapuje na kruh o stˇredu poˇzadovan´e vlastnosti filtru.
1 2
a polomˇeru
1 2.
Metoda tedy obecnˇe nezachov´av´a
• Pˇr.: Po transformaci HP analogov´eho filtru vznikne filtr, kter´y m´a vˇsechny p´oly v lev´e polorovinˇe Z-roviny. Digitalizovan´y filtr tedy nebude HP. 18
4.4.3
Biline´arn´ı transformace
Definovan´a substituc´ı s←
2 z−1 Ts z + 1
Vlastnosti: • Metoda zachov´av´a stabilitu a poˇcet p´ol˚u. • Transformace mapuje levou polorovinu s-roviny na jednotkov´y kruh z-roviny, pˇriˇcemˇz imagin´arn´ı osu mapuje na |z| = 1 ⇒ zachov´av´a vlastnosti a hod´ı se na vˇsechny typy filtr˚u. Lze odvodit, zˇ e pro s = σ + iω, σ=0 1 + i 21 ωTs = eiθ , z= 1 − i 12 ωTs kde θ = 2 arctan(0.5ωTs ), coˇz je transformace, kter´a mapuje ω ∈ (−∞, +∞) ↔ θ ∈ (−π, π) (viz vlastnosti funkce arctan) • Ze znalosti tohoto vztahu mezi θ a ω se nab´ız´ı moˇznost navrhovat nejprve hraniˇcn´ı frekvence podle navrhovan´eho digit´aln´ıho filtru. Tedy napˇr. podle specifikace θp vypoˇc´ıtat 2 θp ωp = tan Ts 2 Shrnut´ı postupu n´avrhu digit´aln´ıho IIR filtru pomoc´ı biline´arn´ı transformace: 1. Hraniˇcn´ı frekvence p´asem poˇzadovan´eho digit´aln´ıho filtru pˇrepoˇcteme podle vztahu ω =
2 Ts
tan
θ 2
.
2. Navrhneme analogov´y filtr podle pˇrepoˇcten´ych specifikac´ı, resp. jeho pˇrenosovou funkci H L (s). 3. Transformujeme H L (s) pomoc´ı biline´arn´ı transformace. 4.4.4
F´azov´a charakteristika IIR filtru˚
• RCSR IIR filtr m´a neline´arn´ı f´azi. • Skupinov´e zpoˇzdˇen´ı m˚uzˇ e m´ıt nejv´ıce dramatick´y pr˚ubˇeh pr´avˇe v propustn´em p´asmu → kaˇzd´a propuˇstˇen´a frekvence m´a jin´e zpoˇzdˇen´ı. To zp˚usobuje zkreslen´ı tvaru zpracovan´eho sign´alu, coˇz m˚uzˇ e b´yt neˇza´ douc´ı napˇr. pˇri zpracov´an´ı biologick´ych sign´al˚u (EKG - zmˇena tvaru QRS komplexu) • Zkreslen´ı zpracovan´eho sign´alu je nejmenˇs´ı m´a-li filtr co nejv´ıce konstantn´ı zpoˇzdˇen´ı v propustn´em p´asmu (v Matlabu pˇr´ıkaz grpdelay).
19
5
Z´akladn´ı krit´eria pro porovn´av´an´ı sign´alu˚
Zde budeme uvaˇzovat pouze re´aln´e sign´aly. Zobecnˇen´ı na sign´aly komplexn´ı je jednoduch´e. Uvaˇzujme dva sign´aly x a y, jejich namˇeˇren´e vzorky x[n] a y[n], n = 1, . . . , N . Smyslem kapitoly je vysvˇetlit si souvislosti a v´yznamy hodnot jako jsou N 1 X bxx = σ C bx2 = x[n]2 N n=1 N X bxy = 1 x[n]y[n]. C N n=1
Prvn´ı veliˇcina m´a v´yznam energie sign´alu x, druh´a je tzv. odhadem korelace sign´al˚u x a y. Zp˚usob znaˇcen´ı tˇechto hodnot je standardn´ı ve statistice a jeho logika postupnˇe vyplyne. Napˇr. znak se stˇr´ısˇkou ve statistice znaˇc´ı odhad.
5.1
Sign´aly jako vektory
Kaˇzd´y sign´al si lze pˇredstavit jako vektor, jehoˇz prvky jsou vzorky sign´alu. Vektory znaˇc´ıme mal´ym tuˇcn´ym p´ısmenem, vektor je standardnˇe sloupcov´y. Sign´aly x a y m˚uzˇ eme reprezentovat jako x[1] y[1] x[2] y[2] x= . y = . . .. .. x[N ] y[N ] Standardn´ı skal´arn´ı souˇcin vektor˚u, znaˇc´ıme hx, yi, je dle definice roven hx, yi =
N X
x[n]y[n],
n=1
z´aroveˇn pomoc´ı maticov´eho n´asoben´ı plat´ı hx, yi = xT y = yT x. Norma vektoru (velikost) je definovan´a kxk2 = hx, xi =
N X
x[n]2 = xT x.
n=1
bxx a C bxy se liˇs´ı pouze n´asoben´ım konstantou 1/N Vid´ıme, zˇ e hodnoty C bxx = 1 kxk2 C N bxy = 1 hx, yi, C N bxx a C bxy m˚uzˇ eme vysvˇetlovat podle v´yznamu a vlastnost´ı skal´arn´ıho souˇcinu vektor˚u. Plat´ı tedy v´yznam C • x a y jsou line´arnˇe nez´avisl´e ⇔ x = αy tedy x = αy tedy jeden sign´al je roven α-n´asobku druh´eho • Schwarzova nerovnost ˇr´ık´a, zˇ e |hx, yi| ≤ kxk · kyk, kde rovnost nast´av´a tehdy a jenom tehdy jsou-li x a y line´arnˇe z´avisl´e • x a y jsou kolm´e ⇔ hx, yi = 0. 20
Naˇs´ım z´avˇerem tedy m˚uzˇ e b´yt, zˇ e Cxy je v absolutn´ı hodnotˇe o to vˇetˇs´ı, cˇ´ım v´ıce jsou x a y line´arnˇe z´avisl´e. Pˇr. Nechˇt x a y jsou n´ahodnˇe a nez´avisle vygenerovan´e vektroy podle nˇejak´eho spojit´eho pravdˇepodobnostn´ıho rozloˇzen´ı (napˇr. Gaussova) s nulovou stˇredn´ı hodnotou. Pak • x a y jsou line´arnˇe z´avisl´e s nulovou pravdˇepodobnost´ı (pravdˇepodobnost, zˇ e oba vektory leˇz´ı na jedn´e pˇr´ımce) ˇ ım delˇs´ı vektory generujeme (ˇc´ım vˇetˇs´ı N ), t´ım je tato pravdˇepodobnost “menˇs´ı”, protoˇze dimenze prostoru • C´ se zvˇetˇsuje.
5.2
N´ahodn´e sign´aly jako n´ahodn´e veliˇciny
ˇ hustotami rozloˇzen´ı fX a fY a sdruˇzenou hustotou fXY . (Hodnota Nechˇt X a Y jsou n´ahodn´e veliˇciny s pravd. fX (x) je u´ mˇern´a koncentraci pravdˇepodobnosti, zˇ e X je rovn´e x. Podobnˇe fXY (x, y) odpov´ıd´a pravdˇepodobnosti, zˇ e X je rovn´e x a z´aroveˇn Y je rovn´e y.) Z´akladn´ı charakteristiky n´ahodn´ych veliˇcin jsou Z stˇredn´ı hodnota X: E[X] = xfX (x)dx R Z 2 rozptyl X: E[(X − E[X])2 ] = (x − E[X])2 fX (x)dx = σX = CXX R Z Z korelace X a Y : E[XY ] = xyfXY (x, y)dxdy R R Z Z kovariance X a Y : E[(X − E[X])(Y − E[Y ])] = (x − E[X])(y − E[Y ])fXY (x, y)dxdy = CXY
R R
Vid´ıme, zˇ e kovariance odpov´ıd´a korelaci X a Y s odeˇcten´ymi stˇredn´ımi hodnotami. Budeme od nyn´ı proto pro jednoduchost pˇredpokl´adat, zˇ e E[X] = E[Y ] = 0. Plat´ı: Jsou-li X a Y nez´avisl´e, pak je jejich korelace (kovariance) nulov´a, tj. E[XY ] = 0. Tento fakt lze snadno odvodit na z´akladˇe definice nez´avislosti X a Y tj., zˇ e fXY (x, y) = fX (x)fY (y). Pˇripomeˇnme, zˇ e nez´avislost jakoˇzto pojem z teorie pravdˇepodobnosti odpov´ıd´a naˇs´ı pˇredstavˇe o tom, zˇ e nˇejak´e veliˇciny cˇ i jevy (sign´aly) mezi sebou nesouvis´ı. Line´arn´ı nez´avislost je v tomto ohledu pˇr´ıliˇs konkr´etn´ı a nedostaˇcuj´ıc´ı. Fakt, zˇ e jsou sign´aly (jako vektory) line´arnˇe nez´avisl´e (jeden nen´ı n´asobkem druh´eho) jeˇstˇe nem˚uzˇ e vypov´ıdat o tom, zˇ e spolu “nijak nesouvis´ı”. V praxi ovˇsem skuteˇcn´e hodnoty E[X], Cxy , atd. nezn´ame, protoˇze nezn´ame hustoty fX , fY a fXY . K dispozici m´ame pouze namˇeˇren´e vzorky sign´al˚u X a Y , kter´e ch´apeme jako realizace X a Y (tedy to je n´asˇ model sign´al˚u). Hodnoty odhadujeme na z´akladˇe z´akona velk´ych cˇ´ısel. XN =
N 1 X x[n] N n=1
(odhad E[X])
N 1 X 2 b Cxx = σ bx = x[n]2 N n=1 N X bxy = 1 C x[n]y[n], N n=1
coˇz jsou stejn´e hodnoty, kter´e uvaˇzujeme od zaˇca´ tku kapitoly. bxy , kter´a je odhadem Cxy , je t´ım vˇetˇs´ı, cˇ´ım v´ıce jsou X a Y z´avisl´e. D˚uleˇzit´ym z´avˇerem je, zˇ e hodnota C bxy je pˇr´ımo u´ mˇern´a hodnotˇe skal´arn´ıho Zde vid´ıme jednoduchou souvislost s pˇredchoz´ı podkapitolou. Hodnota C souˇcinu x a y. Z line´arn´ı algebry zn´ame, zˇ e skal´arn´ı souˇcin vektor˚u funguje jako krit´erium line´arn´ı (ne)z´avislosti vektor˚u. Je-li skoro roven nule, pak jsou vektory skoro kolm´e a tedy “dokonale line´arnˇe nez´avisl´e” a naopak. Odtud m˚uzˇ eme vidˇet jak´a je souvislost mezi pojmy “line´arn´ı nez´avislost vektor˚u” a “nez´avislost n´ahodn´ych veliˇcin”.
21
bxy bl´ızk´a nule, a tedy indikuje, zˇ e Cxy = 0, pak z toho neplyne, zˇ e X Pˇripomˇenˇ me, zˇ e pokud vych´az´ı C ˇ ık´ame pouze, zˇ e jsou nekorelovan´e. Skal´arn´ı souˇcin x a y je tedy “pouze” krit´erium a Y jsou nez´avisl´e. R´ (ne)korelovanosti. V praxi, chceme-li zjistit jsou-li X a Y nez´avisl´e, pak n´am cˇ asto staˇc´ı posouzen´ı hodnoty bxy . V anal´yze nez´avisl´ych komponent, kter´a bude l´atkou v jedn´e z dalˇs´ıch kapitol, je vˇsak podm´ınka C bxy ≈ 0 C pouze nutnou podm´ınkou nez´avislosti. Pravidla tedy zn´ı X a Y jsou korelovan´e ⇒ X a Y jsou z´avisl´e X a Y jsou nez´avisl´e ⇒ X a Y jsou nekorelovan´e avˇsak X a Y jsou nekorelovan´e ; X a Y jsou nez´avisl´e
5.3
Stacion´arn´ı n´ahodn´e procesy
Zde budeme uvaˇzovat tento obecnˇejˇs´ı model n´ahodn´eho sign´alu, kter´y jiˇz nen´ı charakterizov´an jedinou n´ahodnou veliˇcinou, ale posloupnost´ı n´ahodn´ych veliˇcin. Charakterizuje ho dalˇs´ı d˚uleˇzit´a veliˇcina autokovarianˇcn´ı funkce, kter´a je pˇr´ımo souvis´ı s jeho spektrem. Autokovarianˇcn´ı funkce procesu (sign´alu) x, nebo zkr´acenˇe autokovariance, je definovan´a jako Cxx [τ ] = E[(x[n] − E[x[n]])(x[n + τ ] − E[x[n + τ ]])]. Opˇet budeme pro jednoduchost pˇredpokl´adat, zˇ e stˇredn´ı hodnota x je nula, tedy E[x[n]] = E[x[n + τ ]] = 0. Potom Cxx [τ ] = E[x[n]x[n + τ ]]. Odhady tˇechto hodnot prov´ad´ıme prakticky obdobnˇe jako v pˇredchoz´ım. Teoretick´ym pˇredpokladem je, zˇ e proces je tzv. ergodick´y. Odhad autokovariance je N X bxx [τ ] = 1 C x[n]x[n + τ ], N n=1
kde hodnoty x[n], kdy n pˇresahuje hodnoty 1 aˇz N , nahrazujeme nulama. Pro dva procesy x a y definujeme tzv. vz´ajemnou kovarianci (cross-kovarianci) a jej´ı odhad jako Cxy [τ ] = E[x[n]y[n + τ ]] N X bxy [τ ] = 1 C x[n]y[n + τ ]. N n=1
Obdobn´e souvislosti mezi pojmy skal´arn´ı souˇcin a kovariance jako v pˇredchoz´ıch podkapitol´ach lze odvodit pomoc´ı definice vektoru 0 .. k . 0 . x[1] x[k] = x[2] . .. x[N − k] Pak plat´ı bxx [τ ] = 1 x[0]T x[τ ] C N 1 b Cxy [τ ] = x[0]T y[τ ]. N 22
5.4
˚ erov´an´ı Metoda synchronizovan´eho prumˇ
Uvaˇzujme pˇr´ıpad, kdy m´ame k dispozici nˇekolik mˇeˇren´ı stejn´eho sign´alu x. Kaˇzd´e mˇerˇen´ı obsahuje jin´e zaruˇsen´ı a mˇeˇren´ı nejsou cˇ asovˇe synchronizovan´a. Sign´al z´ıskan´y kt´ym mˇeˇren´ım lze popsat jako yk [n] = x[n − dk ] + ξk [n], ´ kde dk je nezn´am´e zpoˇzdˇen´ı x vlivem nesynchronizovan´eho mˇeˇren´ı a ξk je zaruˇsen´ı. Ukolem je z´ıskat co nejlepˇs´ı odhad x ze vˇsech mˇeˇren´ı yk . Jednoduchou u´ vahou dojdeme k pˇredpoklad˚um, zˇ e x a yk jsou z´avisl´e, protoˇze nesou spoleˇcnou informaci. Naopak x a ξk jsou nez´avisl´e, protoˇze procesy zp˚usobuj´ıc´ı ruˇsen´ı vˇetˇsinou nijak nesouvis´ı s procesy, kter´e generuj´ı sledovan´y sign´al (napˇr. EKG a sˇum aparatury). Jako krit´erium z´avislosti x a yk pouˇzijeme odhad jejich vz´ajemn´e bxy [τ ], jej´ızˇ hodnota by mˇela b´yt nejvˇetˇs´ı kdyˇz τ = dk . Dosazen´ım dostaneme kovariance C N N X 1 X bxy [τ ] = 1 C x[n]y[n + τ ] = x[n](x[n − dk + τ ] + ξk [n + τ ]) = N n=1 N n=1
=
N N 1 X 1 X x[n]x[n − dk + τ ] + x[n]ξk [n + τ ], N n=1 N n=1 | {z } | {z } bxx [τ −dk ] C
bxξ [τ ]≈0 C k
tedy vid´ıme, zˇ e hodnota je rovna souˇctu odhad˚u autokovarianˇcn´ı funkce x, jej´ızˇ hodnota je maxim´aln´ı kdyˇz τ = dk , a vz´ajemn´e kovariance x a ξk , kter´a je rovna nule, protoˇze x a ξk jsou nez´avisl´e. Metodu synchronizovan´eho pr˚umˇerov´an´ı prov´ad´ıme v n´asleduj´ıc´ıch kroc´ıch: 1. Zvol´ıme referenˇcn´ı mˇeˇren´ı, kter´ym de facto nahrazujeme n´am nezn´am´y sign´al x, napˇr. zvol´ıme y1 . 2. Pomoc´ı vz´ajemn´e kovariance postupnˇe odhadneme vz´ajemn´a zpoˇzdˇen´ı sign´al˚u y1 a yk , k = 2, . . . . 3. Pomoc´ı odhadnut´ych zpoˇzdˇen´ı sign´aly posuneme (synchronizujeme) a zpr˚umˇerujeme, cˇ´ımˇz z´ısk´ame odhad x. Zd˚uvodnˇen´ı posledn´ıho kroku je n´asleduj´ıc´ı. Pˇredpokl´adejme, zˇ e y1 a y2 jsou jiˇz synchronizovan´e, takˇze y1 [n] = x[n] + ξ1 [n] y2 [n] = x[n] + ξ2 [n]. Jejich pr˚umˇerem vznikne sign´al y[n] =
1 1 1 (y1 [n] + y2 [n]) = x[n] + ξ1 [n] + ξ2 [n]. 2 2 2
Energie sign´alu y je rovna (st´ale pˇredpokl´ad´ame nulovou stˇredn´ı hodnotu vˇsech sign´al˚u) 1 1 E[y[n]2 ] = E[x[n]2 ] + E[ξ1 [n]2 ] + E[ξ2 [n]2 ] + E[x[n]ξ1 [n]] + E[x[n]ξ2 [n]] + E[ξ1 [n]ξ2 [n]], {z } | 4 4 vz´ajemn´e kovariance nez´avisl´ych proces˚u=0
pˇriˇcemˇz prvn´ı cˇ len odpov´ıd´a energii sign´alu x a dalˇs´ı dva energi´ım ruˇsen´ı. Oznaˇcme σx2 = E[x[n]2 ] a pˇredpokl´adejme, zˇ e energie vˇsech ruˇsen´ı je stejn´a tj., zˇ e E[ξ1 [n]2 ] = E[ξ2 [n]2 ] = σξ2 . Pomˇer sign´al/ˇsum (Signalto-Noise Ratio - SNR) sign´alu y je roven σ2 SNRy = 2 x2 . σξ Stejn´y v´ypoˇcet SNRy1 , tj. sign´alu y1 , d´av´a poloviˇcn´ı hodnotu SNRy . Obecnˇe pokud bychom zpr˚umˇerovali M sign´al˚u y1 , . . . , yM , pak SNR v´ysledn´eho sign´alu bude M σx2 /σξ2 . Vyj´adˇreno v decibelech ! σ2 σx2 10 log10 M 2 = 10 log10 M + 10 log10 x2 , σξ σξ 23
takˇze pr˚umˇerov´an´ım M synchronizovan´ych sign´al˚u zlepˇs´ıme SNR o 10 log10 M decibel˚u. Pˇr´ıkladn´ym pouˇzit´ım metody synchronizovan´eho pr˚umˇerov´an´ı je pˇri zpracov´an´ı EKG. Z kr´atk´eho z´aznamu extrahujeme nˇekolik QRS komplex˚u, provedeme jejich synchronizaci a zpr˚umˇerujeme. Z´ısk´ame tak m´enˇe zaˇsumˇen´y nebo pˇresnˇejˇs´ı pr˚ubˇeh QRS. Podobnˇe se pouˇz´ıv´a metoda pˇri zpracov´av´an´ı ERP.
24
6
Filtry minimalizuj´ıc´ı kvadratick´a krit´eria
Pˇredchoz´ı kapitoly se zab´yvaly n´avrhem filtr˚u, pˇriˇcemˇz c´ılem bylo navrhnout filtr co nejl´epe s ohledem na ide´aln´ı frekvenˇcn´ı a f´azovou charakteristiku. V t´eto kapitole budeme filtr navrhovat podle jeho v´ystupu, je-li vstupem konkr´etn´ı sign´al. Nav´ıc se posuneme za hranice LTI syst´em˚u, tedy filtr˚u, kter´e se v cˇ ase nemˇen´ı a budeme se zab´yvat n´avrhem tzv. adaptivn´ıch filtr˚u. Vstupn´ı sign´al, kter´y chceme zpracovat, budeme oznaˇcovat x[n], n = 1, . . . , N . V´ystup, cˇ ili sign´al zpracovan´y, bude y[n]. D´ale budeme d[n] oznaˇcovat sign´al, kter´y chceme zpracov´an´ım z´ıskat, a tedy y[n] by mu mˇel b´yt co nejv´ıce podobn´y, ide´alnˇe roven. Rozd´ıl mezi y[n] a d[n] m´a v´yznam chyby, kter´e se v dan´y okamˇzik dopouˇst´ıme a budeme j´ı definovat sign´al e[n] = d[n] − y[n]. Filtr, kter´y budeme hledat, bude FIR d´elky L a jeho koeficienty (impulzn´ı odezvu) budeme oznaˇcovat w[k], k = 0 . . . , L − 1. Vztah mezi zpracov´avan´ym sign´alem x a v´ystupn´ım sign´alem y je tedy pops´an y[n] = w[0]x[n] + w[1]x[n − 1] + · · · + w[L − 1]x[n − L] =
L−1 X
w[k]x[n − k].
k=0
Zavedeme-li vektory xn =
x[n] x[n − 1] .. .
w=
x[n − L]
w[0] w[1] .. .
,
w[L]
m˚uzˇ eme vztah x a y popsat pomoc´ı maticov´eho souˇcinu y[n] = wT xn . Vˇsechny filtry, kter´e zde budeme uvaˇzovat, jsou navrˇzeny tak, aby minimalizovaly konkr´etn´ım zp˚usobem kvadr´at chybov´eho sign´alu e[n]. Zavedeme krit´erium, kter´e je funkc´ı w (tedy funkc´ı L promˇenn´ych w[k], k = 0 . . . , L − 1, tj. koeficient˚u hledan´eho filtru) a je rovno kvadr´atu chyby v cˇ asov´y okamˇzik n Jn (w) = e[n]2 . Gradient Jn (w), tj. vektor parci´aln´ıch derivac´ı podle jednotliv´ych sloˇzek w, lze zapsat vektorovˇe 4Jn (w) = −2xn d(n) + 2xn xTn w. Zavedeme oznaˇcen´ı Rn = xn xTn pn = xn d[n], pak gradient m˚uzˇ eme zapsat jako 4Jn (w) = −2pn + 2Rn w.
6.1
Least Mean Square (LMS) filtr
T´ımto filtrem (neadaptivn´ım) minimalizujeme pr˚umˇern´y kvadr´at chyby na dan´em cˇ asov´em u´ seku n = n1 , . . . , n2 . Minimalizujeme tedy pr˚umˇernou hodnotu krit´eria Jn (w) na tomto intervalu, tedy minimalizujeme Jn1 ,n2 (w) =
n2 n2 X X 1 1 Jn (w) = e[n]2 . n2 − n1 + 1 n=n n2 − n1 + 1 n=n 1
1
Pro jednoduchost budeme uvaˇzovat cel´y u´ sek dostupn´ych dat. Tedy n1 = 1 a n2 = N a krit´erium m˚uzˇ eme znaˇcit bez doln´ıho indexu N 1 X J(w) = e[n]2 . N n=1 25
Derivace je line´arn´ı operace, takˇze gradient J(w) je roven pr˚umˇeru gradient˚u Jn (w) a plat´ı 4J(w) = −2p + 2Rw, kde R=
p=
N 1 X xn xTn N n=1 N 1 X xn d[n]. N n=1
Minimum krit´eria J(w) vzhledem k w najdeme, kdyˇz poloˇz´ıme gradient rovn´y nule, tj. 4J(w) = 0. Dostaneme ˇreˇsen´ı w = R−1 p. Aby mˇela pˇredchoz´ı rovnice smysl, mus´ı inverze matice R existovat. Ta napˇr´ıklad nemus´ı existovat, pokud uvaˇzujeme pˇr´ıliˇs kr´atk´y interval, kde filtr poˇc´ıt´ame. Inverze napˇr. neexistuje pokud filtr poˇc´ıt´ame na intervalu d´elky jedna, protoˇze pak R = Rn a m´a hodnost 1. LMS filtr m˚uzˇ eme uˇcinit adaptivn´ım napˇr. tak, zˇ e ho budeme poˇc´ıtat na kaˇzd´em bloku dat zvl´asˇˇt. Z v´ysˇe uveden´ych d˚uvod˚u vˇsak mus´ı b´yt kaˇzd´y blok dostateˇcnˇe dlouh´y. Nav´ıc se budeme muset vypoˇra´ dat s nespojit´ymi pˇrechody mezi bloky, protoˇze filtr se m˚uzˇ e mˇenit mezi bloky pˇrekotnˇe.
6.2
˚ filtr Wieneruv
LMS filtr b´yv´a nˇekdy nespr´avnˇe (nebo pˇripusˇtme alespoˇn, zˇ e nepˇresnˇe) zamˇenˇ ov´an s Wienerov´ym filtrem. Zde si vysvˇetl´ıme jak pˇresnˇe je Wiener˚uv filtr definov´an a jak s LMS filtrem souvis´ı. Pˇredpokladem pro definici Wienerova filtru je, zˇ e sign´aly x a d jsou slabˇe stacion´arn´ı procesy, pro jednoduchost, s nulovou stˇredn´ı hodnotou. Z toho plyne, zˇ e i v´ystupn´ı sign´al y a chybov´y sign´al e jsou slabˇe stacion´arn´ı. V krit´eriu, kter´e minimalizujeme, nahrazujeme aktu´aln´ı hodnotu kvadr´atu chyby e[n] jej´ı teoretickou stˇredn´ı hodnotou, takˇze krit´erium definujeme jako J(w) = E[e[n]2 ]. Jelikoˇz je e slabˇe stacion´arn´ı, je toto krit´erium nez´avisl´e na cˇ ase n. Z toho jiˇz m˚uzˇ eme ch´apat, zˇ e Wiener˚uv filtr je filtr “teoretick´y”, zat´ımco LMS filtr je spoˇcten´y z aktu´alnˇe namˇeˇren´ych dat. Odvozen´ı Wienerova filtru je analogick´e odvozen´ı LMS filtru. Sumy poˇc´ıtaj´ıc´ı pr˚umˇern´e hodnoty nahrad´ıme oper´atorem stˇredn´ı hodnoty E[·] a vyjde w = R−1 p, kde R = E xn xTn p = E xn d[n] . Opˇet zde R a p nez´avis´ı na cˇ ase, protoˇze x a d jsou slabˇe stacion´arn´ı. Hlubˇs´ım n´ahledem m˚uzˇ eme vidˇet, zˇ e prvky R obsahuj´ı hodnoty autokovariance x, Cxx [τ ], a d obsahuje hodnoty vz´ajemn´e kovariance x a d, Cxd [τ ]. Kdybychom nyn´ı postupovali opaˇcnˇe a odvozovali Wiener˚uv filtr ze vzork˚u namˇeˇren´ych sign´al˚u, nahradili bychom stˇredn´ı hodnoty jejich odhady, a tedy oper´ator stˇredn´ı hodnoty aritmetick´ymi pr˚umˇery. T´ım bychom dostali opˇet LMS filtr. Spr´avn´y z´avˇer tedy je, zˇ e LMS filtr je odhadem Wienerova filtru za pˇredpokladu, zˇ e x a d jsou slabˇe stacion´arn´ı. Zvˇetˇsujeme-li d´elku dostupn´ych dat, pak LMS filtr konverguje k Wienerovu filtru. Pokud jsou x a d nestacion´arn´ı (ˇreˇc, EEG), pak nem˚uzˇ eme o Wienerovu filtru hovoˇrit a LMS filtr obecnˇe nekonverguje. ˇ Casto lze ale pˇredpokl´adat, zˇ e data jsou stacion´arn´ı alepsoˇn na kr´atk´ych intervalech (ˇreˇc na intervalech max. d´elky ˇ ım je vˇsak interval kratˇs´ı, 20-40ms). Tam tedy Wiener˚uv filtr teoreticky existuje a LMS filtr je jeho odhadem. C´ t´ım je pochopitelnˇe odhad teoreticky m´enˇe pˇresn´y. To se projevuje pr´avˇe t´ım, zˇ e se LMS filtr vlivem odchylky pˇrekotnˇe mˇen´ı.
26
6.3
Adaptivn´ı LMS filtr
V pˇredchoz´ım jsme jiˇz zm´ınili moˇznost poˇc´ıtat LMS filtr po bloc´ıch a t´ım ho adaptovat na aktu´aln´ı data. Probl´emem jsou vˇsak pˇrekotn´e zmˇeny, kter´e se zvˇetˇsuj´ı cˇ´ım v´ıce je interval kr´atk´y. Nav´ıc nelze interval zkr´atit na minim´aln´ı d´elku kv˚uli invertibilitˇe matice R. Adaptivn´ı LMS filtr je proto navrˇzen jinak. Jeho c´ılem je upravit w tak, aby byla minimalizovan´a aktu´aln´ı hodnota chyby v cˇ ase n dan´a krit´eriem Jn (w). Jelikoˇz se tedy hodnota w mˇen´ı v cˇ ase, zavedeme doln´ı index oznaˇcuj´ıc´ı cˇ as, tedy wn . Minimalizaci Jn (w) nelze vyˇreˇsit poloˇzen´ım gradientu rovn´y nule, protoˇze, jak jiˇz bylo ˇreˇceno, takov´a soustava nem´a ˇreˇsen´ı, protoˇze Rn nem´a inverzi. Adaptivn´ı LMS algoritmus proto prov´ad´ı u´ pravu metodou nejvˇetˇs´ıho sp´adu wn+1 = wn − µ4Jn (w), kde µ je d´elka kroku. Dosazen´ım dostaneme krok wn+1 = wn − µxn e[n], Metoda nejvˇetˇs´ıho sp´adu je iteraˇcn´ı metoda pro hled´an´ı minima funkce. Proveden´ım jedn´e jej´ı iterace tedy teoreticky zmenˇs´ıme chybu Jn (w). V n´asleduj´ıc´ım cˇ ase je jiˇz hodnota chyby Jn+1 (w) jin´a. Proveden´ım jedn´e iterace tuto chybu opˇet zmenˇs´ıme atd. Filtr se takto neust´ale adaptuje s ohledem na aktu´aln´ı chybu. D˚uleˇzit´a je ot´azka konvergence. Zde mus´ı b´yt opˇet pˇredpokl´adan´e, zˇ e sign´aly x a d jsou stacion´arn´ı. Pak existuje Wiener˚uv filtr a za urˇcit´ych podm´ınek plat´ı, zˇ e adaptivn´ı LMS filtr k nˇemu konverguje. Vˇecˇ nou ot´azkou z˚ust´av´a volba d´elky kroku µ. Je zˇrejm´e, zˇ e pokud je µ pˇr´ıliˇs mal´e, pak je konvergence pomal´a a filtr se m˚uzˇ e adaptovat nedostateˇcnˇe rychle. Je-li naopak µ pˇr´ıliˇs velk´e, m˚uzˇ e to zp˚usobit divergenci a filtr ˇ nefunguje spr´avnˇe. Ot´azka optim´aln´ı volby µ je sloˇzit´a a z´avis´ı na konkr´etn´ıch pˇredpokladech. Casto b´yv´a jej´ı fixn´ı hodnota nebo funkˇcn´ı z´avislost urˇcov´ana experiment´alnˇe na z´akladˇe zkuˇsenost´ı. Obecnˇe je podm´ınkou konvergence metody nejvˇetˇs´ıho sp´adu 0<µ<
2 λmax
,
kde λmax je nejvˇetˇs´ı vlastn´ı cˇ´ıslo matice Rn . Toto je ovˇsem podm´ınka pro pevnou Rn , kter´a nez´avis´ı na cˇ ase, tedy i pro pevn´e krit´erium Jn (w), kter´e minimalizujeme. Tyto veliˇciny se pˇri pouˇzit´ı adaptivn´ıho filtru v kaˇzd´em kroce mˇen´ı.
6.4
Normalizovan´y adaptivn´ı LMS filtr
Podle definice lze chybov´y sign´al e[n] ps´at ve tvaru e[n] = d[n] − wnT x[n]. Ten odpov´ıd´a chybˇe, kter´e se dopouˇst´ıme v aktu´aln´ı okamˇzik n. Uvaˇzujme jin´y chybov´y sign´al T [n] = d[n] − wn+1 x[n],
kter´y odpov´ıd´a chybˇe ve stejn´y cˇ asov´y okamˇzik n, avˇsak zp˚usoben´e jiˇz adaptovan´ym filtrem, kter´y aplikujeme aˇz v n´asleduj´ıc´ım okamˇziku n + 1. Normalizovan´y adaptivn´ı LMS filtr je definovan´y jako takov´y, kter´y prov´ad´ı minim´aln´ı zmˇenu filtru definovanou jako δwn+1 = wn+1 − wn , (tedy minimalizujeme kδwn+1 k2 ) za podm´ınky, zˇ e [n] = 0. Vyˇreˇsen´ım t´eto u´ lohy vych´az´ı adaptace normalizovan´eho filtru jako xn e[n]. wn+1 = wn − µ kxn k2 Jak vid´ıme, od p˚uvodn´ıho adaptivn´ıho LMS filtru se jeho normalizovan´a verze v podstatˇe liˇs´ı pouze v adaptivn´ı volbˇe d´elky kroku, tedy µ je nahrazeno µ/kxn k2 . 27
6.5
Adaptivn´ı RLS filtr
Zkratka RLS znamen´a Recursive Least Square. Tento filtr v kaˇzd´em cˇ asov´em okamˇziku minimalizuje krit´erium JnRLS (wn ) =
n X
λn−k en [k]2 ,
k=1
kde 0 < λ ≤ 1 a chybov´y sign´al je definovan´y jako en [k] = d[k] − wnT xk . Toto krit´erium se liˇs´ı od pˇredchoz´ıch v n´asleduj´ıc´ıch tˇrech bodech. 1. Horn´ı mez sumy ve vzorci nen´ı pevn´a, je rovna aktu´aln´ımu cˇ asu n. Krit´erium tedy sˇc´ıt´a vˇsechny chyby od zaˇca´ tku aˇz po aktu´aln´ı okamˇzik. 2. Chybov´y sign´al je zde nutn´e definovat s doln´ım indexem n. en [k] odpov´ıd´a hypotetick´e chybˇe zp˚usoben´e filtrem wn v cˇ ase k. V krit´eriu tedy sˇc´ıt´ame hypotetick´e chyby zp˚usoben´e aktu´aln´ım filtrem wn . Vztah s chybov´ym sign´alem uvaˇzovan´ym v pˇredeˇsl´ych kapitol´ach je e[n] = en [n]. 3. Kaˇzd´a chyba je v´azˇ ena hodnotou, kter´a je rovna mocninˇe cˇ´ısla λ. ˇ ım vˇetˇs´ı je mocnina λ, t´ım menˇs´ı je v´aha chyby. Aktu´aln´ı chyba e[n] m´a nejvˇetˇs´ı v´ahu, a tedy nejv´ıce ovlivˇnuje C´ hodnotu krit´eria. V´ahy pˇredchoz´ıch chyb se exponencielnˇe zmenˇsuj´ı, takˇze pˇredchoz´ı chyby ovlivˇnuj´ı krit´erium ˇ ım menˇs´ı je λ, t´ım menˇs´ı maj´ı na krit´erium vliv chyby minul´e. V extr´emn´ım pˇr´ıpadˇe λ = 0 je krit´erium m´enˇe. C´ rovno pouze kvadr´atu aktu´aln´ı chyby en [n]. Na druh´e stranˇe pokud λ = 1, kaˇzd´a chyba m´a v´ahu 1. λ tedy ovlivˇnuje “pamˇeˇt” krit´eria. Pokud je λ = 0, pak JnRLS (wn ) = Jn (wn ), tedy je stejn´e jako krit´erium minimalizovan´e adaptivn´ım LMS filtrem. Pokud je λ = 1, pak je adaptivn´ı RLS filtr v kaˇzd´em okamˇziku roven LMS filtru (neadaptivn´ımu) spoˇcten´em na intervalu od poˇca´ tku do aktu´aln´ıho okamˇziku. Minimalizaci JnRLS (wn ) lze prov´est stejn´ym zp˚usobem jako v odvozen´ı LMS filtru, tj. poloˇzen´ım gradientu RLS Jn (wn ) rovno nule. V´ysledkem je tzv. norm´alov´a rovnice Φn wn = zn , kde Φn = zn =
n X k=1 n X
λn−i xk xTk λn−i xk d[k].
k=1
ˇ sen´ı norm´alov´e rovnice pˇr´ım´ym zp˚usobem v kaˇzd´y Vid´ıme, zˇ e matice Φn a vektor zn jsou obdobou Rn a pn . Reˇ cˇ asov´y okamˇzik obn´asˇ´ı v´ypoˇcet matice Φn a jej´ı inverze a v´ypoˇcet vektoru zn . To je v´ypoˇcetnˇe velmi n´aroˇcn´e. D´ıky tvaru krit´eria JnRLS (wn ) vˇsak lze dos´ahnout v´yrazn´ych u´ spor, kter´e spoˇc´ıvaj´ı v tom, zˇ e pˇri v´ypoˇctech vyuˇz´ıv´ame hodnoty, kter´e jiˇz byly vypoˇcteny v pˇredchoz´ıch kroc´ıch. Konkr´etnˇe lze vyj´ıt z rekursivn´ıch vztah˚u, kter´e je snadn´e odvodit, Φn = λΦn−1 + xn xTn zn = λzn−1 + xn d[n]. Odtud tak´e n´azev filtru “Recursive”. Odvozen´ı je n´asleduj´ıc´ı. K v´ypoˇctu inverze Φn pouˇzijme zn´am´e inverzn´ı lemma, ze kter´eho plyne, zˇ e inverze matice A, takov´e kterou lze napsat ve tvaru A = B−1 + CD−1 CT , je rovna A−1 = B − BC(D + CT BC)−1 CT B.
28
Pr´avˇe v tomto tvaru je matice Φn zapsan´a v´ysˇe uveden´ym rekursivn´ım vztahem. Do inverzn´ıho lemma tedy m˚uzˇ eme dosadit n´asledovnˇe A = Φn −1
B
= λΦn−1
C = xn D = 1, cˇ´ımˇz dostaneme −1 −1 Φ−1 Φn−1 − n =λ
T −1 λ−2 Φ−1 n−1 xn xn Φn−1
, 1 + λ−1 xTn Φ−1 n−1 xn coˇz je vzorec pro rekursivn´ı v´ypoˇcet inverze Φn . Nyn´ı nadefinujeme nˇekolik pomocn´ych promˇenn´ych a provedeme substituce. Definujeme Pn = Φ−1 n kn =
λ−1 Pn−1 xn . 1 + λ−1 xTn Pn−1 xn
Snadno lze ovˇeˇrit, zˇ e kn = Pn xn a rekurzivn´ı vztah Pn = λ−1 Pn−1 − λ−1 kn xTn Pn−1 . Dosad´ıme nyn´ı do rˇeˇsen´ı norm´alov´e rovnice, pouˇzijeme rekurzivn´ı vztah pro zn a Pn a odvod´ıme tak vzorec pro rekurzivn´ı v´ypoˇcet filtru. wn = Φ−1 x d[n] n zn = Pn zn = λPn zn−1 + P | n{z n} kn
= λ(λ−1 Pn−1 − λ−1 kn xTn Pn−1 )zn−1 + kn d[n] = Pn−1 zn−1 −kn xTn Pn−1 zn−1 +kn d[n] | {z } | {z } wn−1
= wn−1 + = wn−1 +
wn−1
kn d[n] − xTn wn−1 T kn d[n] − wn−1 xn .
|
{z
apriorn´ı chyba ξ[n]
}
V´yraz v z´avorce naz´yv´ame apriorn´ı chybou, protoˇze je rovna hypotetick´e chybˇe filtru wn−1 v cˇ ase n, tj. ξ[n] = en−1 [n]. K proveden´ı adaptace filtru tedy nen´ı pokaˇzd´e zapotˇreb´ı poˇc´ıtat inverzi matice Φn . Probl´emem m˚uzˇ e ovˇsem b´yt jej´ı existence. V kapitole o LMS filtru bylo zm´ınˇeno, zˇ e matice R m˚uzˇ e b´yt singul´arn´ı, je-li spoˇctena z pˇr´ıliˇs kr´atk´eho u´ seku dat. Toto pravidlo plat´ı analogicky pro matici Φn . Ta je poˇc´ıt´ana z kr´atk´eho u´ seku je-li λ bl´ızk´a nule, proto mus´ıme λ volit dostateˇcnˇe velkou. Φn je tak´e poˇc´ıt´ana z velmi kr´atk´eho u´ seku pˇri spuˇstˇen´ı algoritmu. V´yhodou rekursivn´ıho v´ypoˇctu ale je, zˇ e m˚uzˇ e b´yt inicializovan´y libovolnou regul´arn´ı matic´ı, cˇ´ımˇz je probl´em v mnoha pˇr´ıpadech vyˇreˇsen. V´ypoˇcet adaptivn´ıho RLS filtru shrneme v n´asleduj´ıc´ıch kroc´ıch. Zavedeme v nich jeˇstˇe pomocn´y vektor hn , d´ıky nˇemuˇz uˇsetˇr´ıme dalˇs´ı operace. Filtr nejˇcastˇeji inicializujeme 1 0 w0 = . P0 = δI, .. 0 kde I je jednotkov´a matice a δ je voliteln´a konstanta. Kaˇzdou n´asleduj´ıc´ı adaptaci filtru poˇc´ıt´ame v kroc´ıch hn = Pn−1 xn hn kn = λ + xTn hn T ξ[n] = d[n] − wn−1 xn
wn = wn−1 + kn ξ[n] Pn = λ−1 Pn−1 − λ−1 kn xTn Pn−1 . 29
´ 7 Uvod do slep´e separace sign´alu˚ a anal´yzy nez´avisl´ych komponent ´ Ulohou slep´e separace (Blind Source Separation - BSS) je z´ıskat p˚uvodn´ı sign´aly ze sign´al˚u namˇeˇren´ych, kter´e jsou jejich smˇes´ı. Ani p˚uvodn´ı sign´aly a syst´em, kter´ym byly zpracov´any neˇz byly namˇeˇreny, nejsou zn´amy. Pˇredpokl´adaj´ı se pouze co moˇzn´a nejobecnˇejˇs´ı vlastnosti sign´al˚u a parametrick´y model syst´emu. Modely syst´emu rozliˇsujeme na • line´arn´ı, kde plat´ı princip superpozice a kter´e rozliˇsujeme na – okamˇzit´y a – konvolutorn´ı, a • neline´arn´ı. Za dostateˇcnˇe obecn´e vlastnosti p˚uvodn´ıch sign´al˚u povaˇzujeme napˇr´ıklad • vz´ajemnou nez´avislost, • ˇr´ıdkost, • nez´apornost. Tyto vlastnosti, d´ıky jejich obecnosti, m˚uzˇ eme uplatnit v nejr˚uznˇejˇs´ıch oblastech zpracov´an´ı sign´al˚u. Nez´avislost lze pˇredpokl´adat u sign´al˚u, kter´e spolu nesouvis´ı, napˇr. srdeˇcn´ı aktivita (EKG), svalov´a cˇ innost (EMG), cˇ innost mozku (EEG) a ruˇsen´ı (ˇsum aparatury, 60Hz ze s´ıtˇe). Lze-li sˇ´ıˇren´ı sign´al˚u v prostˇred´ı povaˇzovat za nekoneˇcnˇe rychl´e (napˇr. je-li vzorkovac´ı frekvence v porovnan´ı s rychlost´ı sˇ´ıˇren´ı sign´alu velmi mal´a) a na senzorech se superponuj´ı, pak lze syst´em povaˇzovat za line´arn´ı okamˇzit´y (EKG, EEG, EMG). Pokud je nutn´e povaˇzovat rychlost sˇ´ıˇren´ı za koneˇcnou, sign´aly nav´ıc interferuj´ı s vlastn´ımi odrazy v prostˇred´ı, pak je syst´em line´arn´ı konvolutorn´ı (akustick´e sign´aly, sign´aly v telekomunikaˇcn´ı oblasti). S neline´arn´ım syst´emem se setk´av´ame v BSS zˇr´ıdka. Pˇr´ıkladem je situace, kdy sign´aly maj´ı vˇetˇs´ı dynamick´y rozsah neˇz A/D pˇrevodn´ıky a doch´az´ı pˇri jejich mˇeˇren´ı k pˇreteˇcen´ı, limitaci nebo jin´e neline´arn´ı operaci. ˇ ıdkost´ı se rozum´ı vlastnost, zˇ e sign´al m´a mnoho koeficient˚u (vzork˚u) v nˇejak´e dom´enˇe nulov´e nebo skoro R´ nulov´e. Lidsk´a ˇreˇc je napˇr´ıklad ˇr´ıdk´a v cˇ asovˇe-frekvenˇcn´ı oblasti (spektrogram), fotografie pˇrirozen´eho prostˇred´ı pak ve waveletov´e oblasti. T´eto vlastnosti se vyuˇz´ıv´a pˇri komprimaci, protoˇze mal´e koeficienty sign´alu lze zanedbat. V BSS lze t´eto vlastnosti vyuˇz´ıt ke vz´ajemn´emu oddˇelen´ı sign´al˚u s pˇredpokladem, zˇ e se v dan´e oblasti sign´aly nepˇrekr´yvaj´ı (protoˇze jsou ˇr´ıdk´e). Nez´apornost lze pˇredpokl´adat u sign´al˚u, kter´e odpov´ıdaj´ı nˇejak´e nez´aporn´e fyzik´aln´ı veliˇcinˇe, napˇr. amplitudov´e spektrum sign´alu.
7.1
Line´arn´ı okamˇzit´y model
Nejprve zavedeme standardn´ı znaˇcen´ı. Poˇcet pozorovan´ych (mˇerˇen´e) sign´al˚u bude m a budeme je znaˇcit x1 [n], . . . , xm [n]. Matici X definujeme tak, zˇ e kaˇzd´y jej´ı ˇra´ dek obsahuje vzorky odpov´ıdaj´ıc´ıho pozorovan´eho sign´alu. Jej´ı rozmˇery jsou tedy m × N , kde N je celkov´y poˇcet vzork˚u. Kaˇzd´y jej´ı sloupec odpov´ıd´a vzork˚um od vˇsech sign´al˚u v jeden stejn´y cˇ asov´y okamˇzik. Obdobnˇe definujeme p˚uvodn´ı sign´aly (nezn´am´e) s1 [n], . . . , sd [n], jejich poˇcet d a matici S o rozmˇerech d × N . Line´arn´ı okamˇzit´y model pˇredpokl´ad´a nekoneˇcnˇe rychl´e sˇ´ıˇren´ı sign´al˚u a jejich superpozici. Kaˇzd´y vzorek pozorovan´ych sign´al˚u tak m˚uzˇ eme zapsat jako line´arn´ı kombinaci sign´al˚u p˚uvodn´ıch bez zpoˇzdˇen´ı, tedy napˇr. xk [n] = ak1 s1 [n] + ak2 s2 [n] + · · · + akd sd [n]. Vid´ıme, zˇ e nt´y vzorek pozorovan´eho sign´alu xk [n] je z´avisl´y pouze na vzorc´ıch p˚uvodn´ıch sign´al˚u z odpov´ıdaj´ıc´ıho cˇ asov´eho okamˇziku. D´ale ho ovlivˇnuj´ı koeficienty ak1 , . . . , akd , z nichˇz kaˇzd´y urˇcuje intenzitu, se kterou se p˚uvodn´ı sign´al pˇriˇc´ıt´a k ostatn´ım. Tyto koeficienty se v cˇ ase nemˇen´ı. Definujeme proto tzv. mixovac´ı matici A, jej´ızˇ ijt´y prvek je roven Aij = aij a jej´ı rozmˇery jsou m × d. Cel´y model je pops´an soustavou m line´arn´ıch rovnic, kterou m˚uzˇ eme vyj´adˇrit pomoc´ı matic jako X = AS. Charakteristick´e pro u´ lohu BSS je, zˇ e celou pravou stranu t´eto soustavy nezn´ame. Rozliˇsujeme n´asleduj´ıc´ı tˇri pˇr´ıpady podle poˇctu sign´al˚u na soustavy 30
• pˇreurˇcen´e (overdetermined), kdy m > d, • urˇcen´e cˇ i regul´arn´ı (determined), kdy m = d a • nedourˇcen´e (underdetermined), kdy m < d. V tomto rozdˇelen´ı z´aroveˇn pˇredpokl´ad´ame, zˇ e hodnost matice A je maxim´aln´ı moˇzn´a tj. rovna menˇs´ımu z cˇ´ısel m a d. V opaˇcn´em pˇr´ıpadˇe lze takovou soustavu redukovat tak, aby tato podm´ınka platila. V n´asleduj´ıc´ım se omez´ıme pouze na pˇr´ıpad regul´arn´ı soustavy. V takov´em pˇr´ıpadˇe existuje inverze matice A, tedy A−1 , a p˚uvodn´ı sign´aly m˚uzˇ eme z´ıskat S = A−1 X. Znamen´a to, zˇ e odhad p˚uvodn´ıch sign´al˚u a matice A jsou ekvivalentn´ı u´ lohy. V pˇr´ıpadˇe nedourˇcen´e soustavy tomu tak nen´ı.
7.2
Anal´yza nez´avisl´ych komponent
Budeme pˇredpokl´adat line´arn´ı okamˇzit´y regul´arn´ı syst´em. Z´asadn´ım pˇredpokladem anal´yzy nez´avisl´ych komponent (Independent Component Analysis - ICA) je pˇredpoklad, zˇ e p˚uvodn´ı sign´aly s1 [n], . . . , sd [n] jsou nez´avisl´e. ´ Vlivem syst´emu m´ısˇen´ı sign´al˚u pozorujeme sign´aly x1 [n], . . . , xd [n], kter´e jsou z´avisl´e. Ulohu ICA lze tedy formulovat tak, zˇ e hled´ame takovou matici W, aby sign´aly WX byly nez´avisl´e. 7.2.1
Neurˇcitosti ICA
Je zˇrejm´e, zˇ e v ide´aln´ım pˇr´ıpadˇe bude nalezen´a matice W rovna matici A−1 . Vzhledem k principu nez´avislosti, na jehoˇz z´akladˇe W hled´ame, lze jiˇz nyn´ı konstatovat, zˇ e ˇreˇsen´ı je nejednoznaˇcn´e v tom, zˇ e jiˇz nelze zpˇetnˇe urˇcit • poˇrad´ı p˚uvodn´ıch sign´al˚u, • rozptyl (energii) sign´al˚u a • znam´enko. Pokud totiˇz tyto vlastnosti p˚uvodn´ıch sign´al˚u zmˇen´ıme, na jejich nez´avislost to nem´a vliv. Dalˇs´ı neurˇcitosti plynou aˇz z konkr´etn´ıch model˚u sign´al˚u. 7.2.2
Modely sign´alu˚
Model p˚uvodn´ıch sign´al˚u urˇcuje to, jak´e m˚uzˇ eme pouˇz´ıt krit´erium jejich nez´avislosti. Z´akladn´ı uvaˇzovan´e modely sign´al˚u pˇredpokl´adaj´ı, zˇ e sign´al je • posloupnost nez´avisl´ych realizac´ı negaussovsk´e n´ahodn´e veliˇciny (princip negaussovskosti), • posloupnost nez´avisl´ych realizac´ı gaussovsk´ych n´ahodn´ych veliˇcin s jejichˇz rozptyl z´avis´ı na cˇ ase (princip nestacionarity), • gaussovsk´y n´ahodn´y slabˇe stacion´arn´ı proces (princip spektr´aln´ı diversity). D´ale se pouˇz´ıvaj´ı kombinovan´e modely tˇechto tˇr´ı pˇredchoz´ıch. Zde se budeme vˇenovat popisu pouze prvn´ıho z nich. 7.2.3
Princip negaussovskosti
Jelikoˇz pˇredpokl´ad´ame, zˇ e sign´al je i.i.d. posloupnost n´ahodn´e veliˇciny, jeho vlastnosti jsou urˇcen´e hustotou pravdˇepodobnosti t´eto veliˇciny. Hustoty p˚uvodn´ıch sign´al˚u s1 [n], . . . , sd [n] oznaˇc´ıme fs1 , . . . , fsd a jejich sdruˇzenou hustotu oznaˇc´ıme fs1 ,...,sd . Jelikoˇz pˇredpokl´ad´ame, zˇ e p˚uvodn´ı sign´aly jsou nez´avisl´e, plat´ı, zˇ e fs1 ,...,sd = fs1 · fs2 · · · fsd . Nyn´ı hled´ame matici W takovou, aby sign´aly Y = WX byly nez´avisl´e, a tedy jejich hustoty a sdruˇzen´a hustota splˇnovaly stejnou podm´ınku. Za t´ımto u´ cˇ elem potˇrebujeme krit´erium, kter´e vyhodnocuje do jak´e m´ıry je podm´ınka splnˇena, resp. “jak hodnˇe” jsou sign´aly Y nez´avisl´e. Toto krit´erium naz´yv´ame objektivn´ı funkce.
31
Vhodn´ym krit´eriem nez´avislosti je tzv. vz´ajemn´a informace Z Z fy1 ,...,yd (y1 , . . . , yd ) I(Y) = dy1 , . . . , dyd , ··· fy1 ,...,yd (y1 , . . . , yd ) ln f (y y1 1 ) · fy2 (y2 ) · · · fyd (yd ) R R kde fy1 , . . . , fyd a fy1 ,...,yd jsou hustoty sign´al˚u Y. Vz´ajemn´a informace je rovna nule pouze pokud jsou Y nez´avisl´e, coˇz je vidˇet z toho, zˇ e cˇ itatel i jmenovatel ve zlomku si budou rovny, hodnota zlomku tedy bude 1 a jeho logaritmu 0. Vlastnost´ı vz´ajemn´e informace je, zˇ e jsou-li sign´aly Y nekorelovan´e (nutn´a podm´ınka nez´avislosti) a normovan´e na jednotkov´y rozptyl, pak d X I(Y) = H(yi ) + const., i=1
kde H(yi ) je entropie sign´alu yi definovan´a Z H(yi ) = −
R
fyi (y) ln fyi (y)dy.
Minimalizaci vz´ajemn´e informace sign´al˚u Y tedy m˚uzˇ eme formulovat tak, zˇ e hled´ame nekorelovan´e sign´aly s minim´aln´ım souˇctem entropi´ı. Z v´yznamu a vlastnost´ı entropie z´aroveˇn m˚uzˇ eme odvodit dalˇs´ı interpretace tohoto principu ICA. • Sign´al, kter´y m´a Gaussovo rozloˇzen´ı a jednotkov´y rozptyl, m´a ze vˇsech sign´al˚u s jednotkov´ym rozptylem nejvˇetˇs´ı entropii. Jestliˇze se snaˇz´ıme entropii minimalizovat, pak se de facto snaˇz´ıme minimalizovat “Gaussovost” sign´alu. • Centr´aln´ı limitn´ı teor´em, jedno z nejz´asadnˇejˇs´ıch tvrzen´ı v teorii pravdˇepodobnosti, m˚uzˇ eme neform´alnˇe interpretovat tak, zˇ e line´arn´ı kombinac´ı n´ahodn´ych veliˇcin (sign´al˚u) dost´av´ame veliˇcinu “v´ıce gaussovskou”. Zam´ıchan´e sign´aly X jsou tedy “v´ıce gaussovsk´e” neˇz sign´aly p˚uvodn´ı. To odpov´ıd´a naˇsemu snaˇzen´ı z´ıskat pomoc´ı ICA sign´aly “co nejm´enˇe gaussovsk´e”. Vzorec entropie H(yi ) lze ps´at ve tvaru H(yi ) = −E[ln fyi ], tedy stˇredn´ı hodnoty logaritmu hustoty sign´alu. Tuto hustotu v praxi nezn´ame a nahrazujeme ji vhodnou neline´arn´ı funkc´ı. V mnoha algoritmech se pouˇz´ıv´a jako n´ahrada logaritmu hustoty funkce x3 nebo tanh(x). Je ovˇsem moˇzn´e pouˇz´ıt parametrick´e odhady hustoty, kter´e vedou k m´ırnˇe vˇetˇs´ı v´ypoˇcetn´ı n´aroˇcnosti, avˇsak vˇetˇs´ı pˇresnosti. Neparametrick´e odhady hustot jsou nejuniverz´alnˇejˇs´ı avˇsak za cenu velk´e v´ypoˇcetn´ı n´aroˇcnosti a cˇ asto i nestability.
7.3
Zpracov´an´ı v´ıcekan´alov´ych z´aznamu˚ pomoc´ı ICA a rekonstrukce
Tuto metodu zpracov´an´ı lze aplikovat na v´ıcekan´alov´e sign´aly, jako jsou napˇr. EEG nebo EKG. Nejprve sign´aly z elektrod uloˇzen´e v ˇra´ dc´ıch matice X analyzujeme pomoc´ı ICA algoritmu, cˇ´ımˇz z´ısk´ame matici W a “p˚uvodn´ı sign´aly” Y = WX. Jelikoˇz nejˇcastˇeji nedok´azˇ eme konkretizovat co jsou v tomto pˇr´ıpadˇe “p˚uvodn´ı sign´aly”, v´ıme pouze, zˇ e jsou nez´avisl´e, naz´yv´ame Y nez´avisl´e komponenty nebo zkr´acenˇe komponenty. T´ım, zˇ e jsou Y nez´avisl´e, mohou nˇekter´e z nich odpov´ıdat izolovan´ym sign´al˚um od proces˚u, kter´e nesouvis´ı s c´ılov´ym sign´alem, j´ımˇz je EKG, EEG, EMG, nebo jejich cˇ a´ st. Chceme-li vybran´e komponenty z p˚uvodn´ıho z´aznamu odstranit, m˚uzˇ eme toho dos´ahnout pomoc´ı rekonstrukce p˚uvodn´ıch sign´al˚u pouze z komponent, kter´e zachov´ame. Postup lze shrnout do n´asleduj´ıc´ıch bod˚u. 1. Ze z´aznamu vytvoˇr´ıme matici X. 2. Pomoc´ı ICA algoritmu nalezneme nez´avisl´e komponenty X, tedy matici W a komponenty Y = WX. 3. Vybereme komponenty, kter´e chceme z dat odstranit (poloˇz´ıme rovny nule) nebo jinak zpracovat (napˇr. e filtrace). Z´ısk´ame tak upravenou matici Y. e = W−1 Y. e 4. P˚uvodn´ı z´aznam rekonstruujeme X 32
e = Y, Posledn´ı krok naz´yv´ame rekonstrukce, protoˇze pokud bychom nechali komponenty Y beze zmˇen, tedy Y pak by byla i p˚uvodn´ı zaznamenan´a data beze zmˇen. Konkr´etnˇe by platilo e = W−1 Y e = W−1 Y = W−1 WX = X. X
References ˇ [1] J. Svatoˇs, Biologick´e sign´aly I., Nakladatelstv´ı CVUT, 1992. ˇ [2] J. Holˇc´ık, Modelov´an´ı a simulace biologick´ych syst´em˚u, Nakladatelstv´ı CVUT, 2006. [3] B. Porat, A Course in Digital Signal Processing, John Wiley & Sons, New York, 1997. [4] R. M. Rangayyan, Biomedical Signal Analysis: A Case Study, John Wiley and Sons, 2002.
33
