Užtriukšmintų signalų modeliai

Užtriukšmintų signalų modeliai yra labai populiarūs signalų apdorojimo literatūroje (žr., pvz., [1] - [6] ). Vienas iš paprasčiausių modelių yra užtriukšmintų sinusoidžių suma. Paprastai daroma prielaida, kad triukšmas yra baltas Gauso (normalinis) su nuliniu vidurkiu ir baigtine dispersija. Pagrindiniai parametrai, kuriuos reikia įvertinti yra dažniai. Sudėtingesnis modelis yra gaunamas, įvedus gesimą. Šiuo atveju pagrindiniai parametrai, kuriuos reikia įvertinti yra dažniai ir gesimo koeficientai. Šie parametrai įeina į modelį netiesiniu būdu.

Šiame darbe yra nagrinėjama modelių klasė, sudaryta iš gęstančių eksponenčių ir negęstančių / gęstančių sinusoidžių sumos. Daroma prielaida, kad ši suma yra užtriukšminta balto Gauso triukšmo seka.

Darbo tikslas yra ištirti Prony klasės metodų panaudojimo galimybę įvertinti nagrinėjamų modelių parametrus, bei išanalizuoti gautų įverčių tikslumą skaitinio modeliavimo pagalba. Darbo uždaviniai yra tokie:

1.      Išnagrinėti klasikinį Prony metodą.

2.      Išvesti formules, skirtas negęstančių ir gęstančių sinusoidžių parametrų įvertinimui klasikiniu Prony metodu.

3.      Išnagrinėti Prony metodą su priešfiltriu.

4.      Sukurti programinę įrangą MATLAB'o terpėje, skirtą nagrinėjamų modelių parametrų įvertinimui Prony klasės metodais.

5.      Naudojant sukurtą programinę įrangą atlikti skaitinius eksperimentus.

Prony metodas (Prony 1795; Hilderbrand 1956; Kay ir Marple 1981) yra metodas, kuris modeliuoja duomenis (vienodai nutolusias signalo atskaitas) tiesine eksponenčių kombinacija. Šis metodas aproksimuoja N duomenų taškų kreive, sudaryta iš p eksponenčių. Kadangi kiekvienas eksponentinis narys a_ke^Sk" turi du parametrus (amplitudę a_k ir rodiklį s_k ), Prony metodui reikia ne mažiau kaip 2p duomenų, t.y. N > 2p. Pradinis Prony metodas bando tiksliai pritaikyti kreivę, sudarytą iš p eksponentinių narių, prie N = 2 p elementų aibės. Kai N > 2 p, naudojamas mažiausių kvadratų metodas.

Tarkime, kad turime tam tikros funkcijos v(t) matavimo vienodai nutolusiuose taškuose rezultatus:

v_n = v(nAt), n = 0,1,...,N -1.                                                                         (2.1.1)

Mes norime aproksimuoti duomenis { v_n } sekančiu eksponentiniu modeliu:

Ža_ke^Skn ,                                                                                                  (2.1.2)

k=1

čia a_k, s_k g C, ir

(2.1.3)

^ak = ^Ak ^exp(M ^ ^sk ={^k + ^jwk )^At.

Įprastas parametrų { A_k ,ę_k ,X_k, w_k } ir p įverčių radimo būdas yra minimizuoti paklaidą:

N-1                   ²

ŽK -x„\ .                                                                                                         (2.1.4)

n=0

Tačiau, tai yra sudėtingas netiesinis mažiausių kvadratų uždavinys. Literatūroje, galima rasti iteratyvių metodų, sprendžiančių (2.1.4). Šiuose metoduose, pradinis nežinomų parametrų įvertis yra nuosekliai gerinamas. Prony metodas yra dažnai naudojamas modelio (2.1.2) parametrų įvertinimui. Nors šis metodas duoda suboptimalų sprendinį, kuris neminimizuoja (2.1.4), tačiau jo rezultatai yra gana patenkinami daugeliu atvejų. Pažymime

Zk = e^Sk     (k = 1,...,p).                                                                                    (2.1.5)

Dydžiai z_k yra vadinami modelio (2.1.2) poliais.                            Tada lygtis (2.1.2) gali būti užrašyta tokiu
pavidalu:

Xn =Ž a_kzl.                                                                                              (2.1.6)

k=1

v_n = į a_kzⁿ_k      (n = 0,1,...N -1).                                                          (2.1.7)

k=1

Perrašykime šį sąryšį kaip lygčių sistemą:

a ₁ + a ₂ + ... + a   = v ₀.

^ai ^Z1 + ^a2 ^Z2 +^... + ^ap^Zp = ^V1,

a, z,² + a₂ z ₂² +... + a_nz p = v₂
^{1122            p p     2}                                                                               (2.1.8)

a, z^N-1 + a2   ^-1 +... + a_pz^N-1 = v_N __v

Jei z₁,z₂,...,z_p būtų žinomi, tuomet lygtys (2.1.8) sudarytų N tiesinių lygčių sistemą su p nežinomųjų a₁, a₂,..., a_p. Jei N = p, tai ši sistema gali būti išspręsta tiksliai, o jei N > p, tai ji gali būti išspręsta apytiksliai mažiausių kvadratų prasme.

Tačiau, jei { z_n} yra nežinomi, mums reikia turėti mažiausiai 2p lygčių. Be to, { z_n} pasirodo

sistemoje (2.1.8) netiesiniu būdu, kas sukuria papildomų sunkumų. Šių sunkumų minimizavimas sudaro Prony metodo esmę.

Tegul z₁,z₂,...,z_p yra tokios algebrinės lygties šaknys

B(z) = z^p -fi, z^p-1 -fi₂ z^p-2 -...-P_p_i z-P_p =(z - zi )■...-(z - zp )= 0.                                             (2.1.9)

Mes turime nustatyti koeficientus fi_x,...,fi . Šiuo tikslu, padauginame pirmą sistemos (2.1.8) lygtį

iš fi_p, antrą iš fi_p-1, ... , p-ąją iš fi₁, o (p+/)-ąją iš (-1) ir sudedame rezultatus. Kadangi

z₁,z₂,...,z_p yra lygties (2.1.9) šaknys, mes gauname tokį sąryšį

Vp -fiVp-1 -... -fipVo = 0.                                                                                     (2.1.10)

Tą pačią procedūrą kartojame, pradedant nuo antros, trečios, ..., (N - p) -osios lygties. Tokiu būdu gausime N - p -1 papildomas lygtis. Taigi, naudojant (2.1.8) ir (2.1.9), mes gauname tokią sistemą, sudarytą iš N- p tiesinių lygčių

^Vp-1fi + ^Vp-2fi2 + ^... + ^V0fip = ^Vp ^{, v}pfi1 + ^Vp-fi2 + ^... + ^Vfip = ^Vp^

^VN-2^1 + ^VN-3^2 + ^... + ^VN - p-1^ßp = ^VN-1,

arba sutrumpinta forma

S v_n-k ßk = v_n,   n = p, p +1,..., N -1.                                                          (2.1.12)

k=1

čia f₀ = -1.

Kadangi v₁,v₂,...,v_N-1 yra žinomi, tiesinių lygčių sistema (2.1.11) gali būti išspręsta tiksliai, jei N = 2p, arba išspręsta apytiksliai mažiausių kvadratų metodu, jei N > 2p .

Suradus koeficientus ff,...,ff , nežinomieji z_1v..,z_v gali būti apskaičiuoti kaip lygčių sistemos

(2.1.8) šaknys. Jos yra realios arba kompleksinės. Tada (2.1.8) tampa sistema, sudaryta iš N tiesinių lygčių su p  nežinomųjų a₁,a₂,...,a_p. Šie nežinomieji gali būti nustatyti naudojant

pirmąsias sistemos (2.1.8) lygtis, arba priimtiniau, naudojant mažiausių kvadratų metodą visai sistemai (2.1.8).

Aprašytasis metodas yra vadinamas Prony metodu. Viena iš jo svarbiausių savybių yra ta, kad lygčių (2.1.8) netiesiškumas yra perkeliamas algebrinės lygties (2.1.9) šaknų suradimui.

Klasikinio Prony metodo santrauka:

Žingsnis 1. Sprendžiame lygčių sistemą (2.1.11) ir nustatom nežinomus polinomo (2.1.9) koeficientus.

Žingsnis 2. Apskaičiuojame polinomo (2.1.9) šaknis z₁,z₂,...,z_p .

Žingsnis 3. Sprendžiame lygčių sistemą (2.1.8) ir surandame kompleksines amplitudes

2.2. Prony metodas su priešfiltriu

Modeliavimo rezultatai rodo, kad Prony metodas yra labai jautrus duomenų matavimo paklaidoms, t.y. kuo paklaidos didesnės, tuo įverčiai blogesni - jie labai paslinkti nuo tikrųjų parametrų reikšmių. Norėdami pagerinti Prony metodo įverčių tikslumą, amerikiečių mokslininkai Kumaresan ir Feng [1] pasiūlė naują metodą - vadinamąjį Prony metodą su priešfiltriu. Metodo esmė yra panaudoti duomenų filtravimą prieš naudojant Prony metodą. Filtruojama tam, kad apriboti triukšmo įtaką. Jie pasiūlė metodą, kaip apskaičiuoti priešfiltrio

W (z )= 1 + w₁ z ^- +... + w_qz-^q,                                                                 (2.2.1)

(čia 0 < q < N - 2 p ) impulsinę reakciją 1,w_l,...,w_ą. Pasiūlytame metode priešfiltrio impulsinė

reakcija yra gaunama iteraciniu būdu, naudojant duomenis.

Nagrinėkime sekos x_n, apibrėžiamos sąryšiu (2.1.2),                          z -transformaciją. Ši transformacija yra
apibrėžiama tokiu būdu:

N-1

X(z) = £ x_nz ⁿ .                                                                                                (2.2.1.1)

n=0

Kadangi galioja sąryšis (2.1.2), tai

N-1   p

X(z) = YZ^a_k^eSknz ⁿ.                                                                                (2.2.1.2)

=0 k=1

Kadangi sąryšyje (2.2.1.2) sumos yra baigtinės, tai galima sukeisti sumavimo tvarką, t.y.

p     N-1

^p        N-1                        _F        N -i  ,              .

(2.2.1.3)

k=1       n=0

k=1       n=0

Pažymėkime

Y = e¹ z

(2.2.1.4)

Tada

N-1

n=0                         n=0

(2.2.1.5)

o tai yra geometrinė progresija. Todėl

N-1 n=0

Vadinasi,

1 -Y^N 1 -Y

(2.2.1.6)

■e^SkNz ^{- N}

h   Ili^1-e"^z-¹)

k=1

^w   *-> ^k    1 -e^Skz-¹

(2.2.1.8)

k=1

Vadinasi, X (z) galima užrašyti kaip dviejų polinomų santykį tokiu būdu:

x (z y.

(2.2.1.9)

čia

B(z ) = n(l - e^skz ^-1 )= 1+ Z bk

.-k

(2.2.1.10)

k=1

k=1

yra laipsnio p polinomas su šaknimis e*¹,..., e ^p, o

C (z ) = Ž

s_kN   -N

e ^k z

(2.2.1.11)

k=1

j=1

yra polinomas laipsnio N + p -1, todėl jį galima užrašyti tokiu pavidalu:

N+p-1

N + p-1

-k

C(z)= Ž ^c_k

(2.2.1.12)

k=0

Panagrinėsime C (z) koeficientus atskirais atvejais: Tegul p = 1, tada

eSlNz - N ) = co +    - + c2z -2 + k+Cn-i (z-))-1 + Cn (z^)

C (z ) = a₁ (l - e^SlNz ^{- N}) Iš šio sąryšio matome, kad c₀ = a₁,   c₁ = c₂ =... = c_N-1 = 0,   c_N = -a₁e^SiN. Tegul p = 2 , tada

C (z ) = a₁ (l - e^SlNz ^{- N} )(l - e^S2z ^- )+ a₂ (l - e^S2N z ^{- N} - e^Sl z ^- ) = a₁ - a₁e^S2z ^- --ae^SlNz^-N + ae^SlNe^S2 (z^-l f+^l + a₂ -a₂e^Slz^-l -a₂e^S2Nz^-N + a₂e^S2Ne^Sl (z^-l f+^l = (a_le^S2 + a₂e^Sl)z"^l -(a_le^SlN + a₂e^S2N)z^^N +(a_le^SlNe^S2 + a₂e'^2Ne^Sl)(z^- f

= a_l + a₂

(2.2.1.13)

(2.2.1.14)

Iš sąryšio matome, kad

c 1       e       e , Cn =-ae^S1N - ae^S2N,

(2.2.1.15)

sN   s₂   .              s₂N

e ¹ e ² + a₂e ² e

^CN+1      ^a1

Analogiškai galima gauti, kad bendru atveju

^Cp = ^Cp+1 = ^K = ^CN-1 = ^0.

Remiantis (2.2.1.9), turime

X(z)B(z ) = C (z).

Sąryšyje (2.2.1.17) paėmę abiejų pusių atvirkštinę z -transformaciją, gauname

čia * žymi tiesinės sąsūkos operaciją, t.y.

^p

^Xn * ^bn = Ž ^Xn-k^bk = ^Xn^b0 + ^Xn-l^bl + ^Xn-2^b2 + - + ^Xn-p^bp

k=0

Kadangi galioja sąryšis (2.2.1.16), tai

Ž ^Xn-k^bk = ^     ⁿ = A ^p + ^{1,kN - 1   (b}0 = ^1).

(2.2.1.16) (2.2.1.17) (2.2.1.18)

(2.2.1.19) (2.2.1.20)

k=0

2.2.2. Filtro tipo nustatymas (ar FIR, ar IIR)

Nagrinėjame baigtinės impulsinės charakteristikos filtrą su impulsine charakteristika:

1,w₁,...,w_ą.                                                                                  (2.2.2.1)

Jos z -transformaciją žymėsime W(z)

W(z) = 1 + w₁z ^- +... + w_qz^{- q}.                                                     (2.2.2.2)

Mes naudosime W(z) duomenų filtravimui. Paėmę abiejų sąryšio (2.2.1.18) pusių sąsūką su w_n, gauname:

Pažymėkime

Gauname

Kadangi

tai

w * x * b = w * c .

n       n       n               n       n

x = w * x

n       n            n

^dp+q = ^dp_+U+1 = ^... = ^dN-1 = 0'

S ^xn-k^bk = ⁰'     ⁿ = P + ^q,P + ^q + ^{1,...N - 1}

k=0

(2.2.2.3) (2.2.2.4)

(2.2.2.5)

(2.2.2.6) (2.2.2.7)

(2.2.2.7) yra N-(p + q) tiesinių lygčių sistema su p nežinomųjų b₁,b₂,...,b_p. Todėl mes turime turėti sąryšį:

N-(p + q)> p.                                                                     (2.2.2.8)

Atskliaudę, gauname:

N - p - q > p.                                                                      (2.2.2.9)

Sutvarkę šią nelygybę, gauname:

N - 2p > q, t.y. q < N - 2p.                                         (2.2.2.10)

Kadangi q > 0 , tai turime 0 < q < N - 2p . Reiškia filtras W(z)                               yra baigtinės impulsinės
charakteristikos (reakcijos) filtras (angl. FIR - filtras).

2.2.3. Filtro impulsinės charakteristikos apskaičiavimas

Lygybė sąryšyje (2.2.1.20) galioja tik kompleksiniams eksponentiniams signalams be triukšmo. Šį sąryšį (2.2.1.20) perrašome tokiu pavidalu:

x_n + xn-1b + x_n-₂b₂ +... + xn-p+bp-1 + xn-pbp = 0, kai n = p,p + 1,...,N-1.                                     (2.2.3.1)

Apibrėžkime postūmio operatorių z ¹ tokiu būdu:

z '¹x_n = x_n-1.                                                                                  (2.2.3.2)

Tada

x_n + z-^lx_nb_: + z-²x_nb₂ +... + z-^{p-\bp-1 + z-^px_nb_p = 0, kai n = p,p + 1,...,N-1.                              (2.2.3.3)
Iškėlę x_n už skliaustų, turime:

(1 + z-¹b₁ + z ²b₂ +... + z^Ap-1)b_p-1 + z^-pb_p)■ x_n = 0, kai n = p,p + 1,...,N-1.                     (2.2.3.4)

Remiantis polinomo B(z) apibrėžimu, gauname:

B(z)-x_n = 0,   n = p,p + 1,k ,N -1.                                                      (2.2.3.5)

Matome, kad polinomas B(z) anuliuoja signalą. Todėl polinomo B(z) šaknys yra artimos sričiai, kurioje koncentruojasi signalo energija, o 1/B(z) turi pikus tuose taškuose. Todėl W(z) yra parenkama taip, kad aproksimuotų 1/B(z) :

B(z )■ W (z ) = 1 + E(z ),                                                                   (2.2.3.6)

čia E (z) yra paklaidos sekos e_n z -transformacija. Laiko srityje sąryšį (2.2.3.6) galima užrašyti tokiu pavidalu:

b_n * w_n =S_n *e_n.                                                                                         (2.2.3.7)

n        n            n       n                                                                                            \               /

Priešfiltrio impulsinė reakcija w₁, w₂,..., w_q yra gaunama, minimizuojant paklaidų kvadratų sumą:

2                                       2

Z N =ZK - k * w_n\ ,                                                                   (2.2.3.8)

n=0                n=0

čia w₀ = 1, t.y. ieškomi tokie w₁, w₂,..., w_q, kurie minimizuoja šią paklaidų kvadratų sumą. Dydis 8_n yra vadinamas Kronekerio delta. Jis apibrėžiamas taip:

5_n = 1, kai n = 0 ir 5_n = 0, kai n * 0 .                                                           (2.2.3.9)

Apskaičiavus filtro impulsinę reakciją, duomenys filtruojami. Gaunami nufiltruoti duomenys {x_n} Su šiais nufiltruotais duomenimis sudaromos Prony lygtys:

Z x'A = 0,   n = p + q, p + q +1, k , N -1,                             (2.2.3.10)

k=0 kurios yra sprendžiamos tam, kad gauti naujus koeficientus Ą,...,b_p (b₀ = 1). Procedūra kartojama

tol, kol koeficientai {b_i} pradeda žymiai nesikeisti pereinant prie naujos iteracijos. Tuomet

apskaičiuojama polinomo B(z) šaknys. Apskaičiavus šaknis e*¹,...,e^Sp, koeficientai a_l,a₂,...,a_p

yra gaunami mažiausių kvadratų metodu. Prony metodo su priešfiltriu santrauka:

Žingsnis 1. Imkime w₀ = l ir w_n = 0, kai n = l,2,...,q tam tikram fiksuotam q (q < N - 2p).

Žingsnis 2. Naudodami w₁,w₂,...,w_q iš Žingsnis l, sudarome lygčių sistemą (2.2.3.10).

Priminsime, kad x'_n yra apibrėžiama kaip: x'_n = w_n * x_n. Žingsnis 3. Sprendžiame lygčių sistemą (2.2.3.10) mažiausių kvadratų metodu ir gauname

^bl^,b2^,K,bp ^.

Žingsnis 4. Naudodami gautus b_l,b₂,...,b_p , ieškome naujų w₁,w₂,...,w_q minimizuodami

paklaidų kvadratų sumą (2.2.3.8). Žingsnis 5. Grįžtame į Žingsnį 2, t.y. vėl sudarome lygčių sistemą (2.2.3.10), naudodami gautus w₁,w₂,...,w_q. Kartojame Žingsnį 3 ir gauname naujus b_l,b₂,...,b_p . Juos

naudojame Žingsnyje 4 naujų w₁,w₂,...,w_q gavimui. Žingsnį 5 kartojame tol,

kol b_l,b₂,. ,b_p pradeda reikšmingai nesikeisti.

Žingsnis 6. Naudodami gautus koeficientus b_l,b₂,...,b_p   sudarome polinomą B(z) ir

randame jo šaknis, kurios yra signalo poliai. Naudodami gautus polius, sprendžiame klasikinio Prony metodo (2.1.8) lygčių sistemą ir gauname

kompleksines signalo amplitudes a₁ , a₂ ,... , a _p .