MOKSLASplius.lt

Skaitiniai filtrai

Ankstesniame eksperimente interferogramai iššifruoti naudojome diskrečiąją Fourier transformaciją, o dabar susipažinsime su Laplace'o transformacijos diskretiniu analogu — \twcal Z transformacija. Abi jos yra labai svarbios apdorojant diskretinius signalus ir vaizdus. \twcal Z transformacija atsirado palyginti neseniai, kartu su kompiuteriais. Ji tampriai siejasi su rekurentinių lygčių teorija. Skyriuje susipažinsime su dviejų tipų skaitiniais filtrais: rekursiniais ir nerekursiniais. Parodysime, kaip jie filtruoja realiu laiku siunčiamą signalą. Panagrinėsime pereinamąsias tokių filtrų charakteristikas bei atskleisime tamprų ryšį tarp \twcal Z transformacijos ir filtro dažninės charakteristikos. Kadangi matavimų duomenys labai dažnai apdorojami kompiuteriais, skaitiniai filtrai plačiai taikomi sudiskretintų duomenų apdorojimui.

\bf \twcal Z transformacija

Tarkime, kad laiko momentais $t=n\,T$, kur $n=0,1,2,\dotsc$, o $T$ yra signalo diskretizavimo žingsnis ar periodas, signalas įgyja $x(n)$ reikšmes. Diskretinio signalo arba diskretinės funkcijos $x(n)$ \twcal Z vaizdu (atvaizdu) vadinsime [Krivickas84,Doetsch67] kompleksinio kintamojo $z$ funkciją $X(z)$, gaunamą tokiu būdu: X(z)=\sum_{n=0}^{n=\infty}x(n)\, z^{-n}.\tag{1} Visada egzistuoja ir atvirkštinė \twcal Z ransformacija. Iš vaizdo $X(z)$ ji atstato pirminę funkciją $x(n)$: x(n)=\frac{1}{2\pi\,\mathrm{i}}{\int_C} X(z){z^{n-1}}\,\mathrm{d} z,\quad n=0,1,2,\dots\qquad \textrm{\v cia}\quad \mathrm{i}={\sqrt{-1}}.\tag{2} Integravimo kontūras $C$ apima koordinačių pradžią kompleksinėje $z$ plokštumoje.

Apskaičiuosime kelių elementarių funkcijų \twcal Z atvaizdus. Diskretinis Heavyside'o funkcijos analogas yra žemiau pavaizduotas diskretinis laiptelis $x(n)$, tenkinantis sąlygas $x(n)=0$, kai $n<0$ ir $x(n)=1$, kai $n\geq 0$. SF01<\![CDATA[ \boldmath$ \mathrm{ListPlot}[\mathrm{Table}[\{n,$ =0,1,0]"]" > \boldmath$\},\{n,\fcolorbox[rgb]{1,0,0}{1,1,1}{$\vphantom{\vbox to 6pt{}}\smash{n_{start}}$},\fcolorbox[rgb]{1,0,0}{1,1,1}{$\vphantom{\vbox to 6pt{}}\smash{n_{end}}$}\}],\mathrm{AxesLabel}\to\{"n","x"\},\fcolorbox[rgb]{1,0,0}{1,1,1}{$\vphantom{\vbox to 4pt{}}\smash{parinktys}$}]$

Duomenys nuo pradinio \boldmath$n_{start}$=" > iki \boldmath$n_{end}$=" >
Piešimo parinktys: All, ImageSize -> 360,AxesOrigin->{0,0},PlotStyle->{RGBColor[0,0,1],PointSize[0.03]}}"]">

MSPBlock[{$$r1,$$r2,$$funcv,$$parinktys}, duomenys = Table[{n,$$funcv}, {n, $$r1, $$r2}]; MSPShow[ListPlot[duomenys,AxesLabel->{"n","x"},Evaluate[Sequence@@$$parinktys]]] ]
]]>
20070824181129
11 submitas Vykdyti

Jo \twcal Z vaizdas yra geometrinės progresijos narių suma, $X(z)=\smash{\sum\limits_{n=0}^\infty} z^{-n}$, kurią visai nesunku apskaičiuoti analiziškai: SF02<\![CDATA[ \boldmath\begin{eqnarray*} &&\mathrm{Sum}[z^{-n},\{n,0,\infty\}] \end{eqnarray*}

MSPFormat[Sum[z^(-n),{n,0,Infinity}],StandardForm]
]]>
20070824181657
20 submitas Vykdyti

Jei diskretinė funkcija $x(n)$ turi gęstančios eksponentės pavidalą, pavyzdžiui, $\mathrm{e}^{-n/5}$, SF03<\![CDATA[ \boldmath$ \mathrm{ListPlot}[\mathrm{Table}[\{n,$ =0,Exp[-0.2*n],0]"]" > \boldmath$\},\{n,\fcolorbox[rgb]{1,0,0}{1,1,1}{$\vphantom{\vbox to 6pt{}}\smash{n_{start}}$},\fcolorbox[rgb]{1,0,0}{1,1,1}{$\vphantom{\vbox to 6pt{}}\smash{n_{end}}$}\}],\mathrm{AxesLabel}\to\{"n","x"\},\fcolorbox[rgb]{1,0,0}{1,1,1}{$\vphantom{\vbox to 4pt{}}\smash{parinktys}$}]$

Duomenys nuo pradinio \boldmath$n_{start}$=" > iki \boldmath$n_{end}$=" >
Piešimo parinktys: All, ImageSize -> 360,AxesOrigin->{0,0},PlotStyle->{RGBColor[0,0,1],PointSize[0.015]}}"]">

MSPBlock[{$$r1,$$r2,$$funcv,$$parinktys}, duomenys1 = Table[{n,$$funcv}, {n, $$r1, $$r2}]; MSPShow[ListPlot[duomenys1,AxesLabel->{"n","x"},Evaluate[Sequence@@$$parinktys]]] ]
]]>
20070825095243
31 submitas Vykdyti

tai funkciją $X(z)=\mathop{\sum}\limits_{n=0}^\infty \mathrm{e}^{-\frac15 n } z^{-n}$ taip pat galima tiksliai apskaičiuoti: $X(z)=z\big/\big(z-\mathrm{e}^{-\frac15 }\big)$: SF04<\![CDATA[ \boldmath\begin{eqnarray*} &&\mathrm{Sum}[\ee^{-\frac15 n} z^{-n},\{n,0,\infty\}] \end{eqnarray*}

MSPFormat[Sum[z^(-n)*Exp[-(1/5)*n],{n,0,Infinity}],StandardForm]
]]>
20070824182855
40 submitas Vykdyti

Begalybės simbolį $\infty$ surinksite klavišais Esc inf Esc. Kaip matome, \boldmath$\mathrm{Sum}[~]$ komanda gali sumuoti ne tik baigtines, bet ir begalines simbolines eilutes. Simboliniai begalinių eilučių sumavimo algoritmai neturi nieko bendra su įprasto baigtinio sumavimo veiksmais. Paprastai šie algoritmai stengiasi surasti rekurentinę lygtį, kurią tenkina ieškoma begalinė suma. Jei rastą rekurentinę lygtį pavyksta išspręsti, pateikiamas atsakymas. Rekurentinėms lygtims spręsti dažnai panaudojama ir \twcal Z tranformacija. Mathematica sistemos galimybes spręsti rekurentines lygtis išplečia standartinis \boldmath$\mathrm{DiscreteMath}{}^\backprime}\mathrm{RSolve}{}^\backprime paketas bei šioje srityje dirbančių mokslininkų paruoštos programos. Pastarąsias rasite interneto svetainėje http://www.risc.uni-linz.ac.at.

Lengva įsitikinti, kad bendresniu atveju diskretinio eksponentinio pereinamojo proceso $\mathrm{e}^{-b n}$ \twcal Z vaizdas yra $z/(z-\mathrm{e}^{-b})$, kur $b$ — teigiamas skaičius.

Kadangi \twcal Z transformacija plačiai naudojama, \twcal Z vaizdams skaičiuoti Mathematica turi specialią komandą \boldmath$\mathrm{ZTransform}[x(n), n, z]$. Ankstesnėse versijose \twcal Z transformacijai atlikti reikia iškviesti \boldmath$\mathrm{DiscreteMath}{}^\backprime}\mathrm{ZTransform}{}^\backprime paketą. Taigi, norėdami rasti diskretinės funkcijos $x(n)=\mathrm{e}^{-b n}$ \twcal Z vaizdą, rašome: SF05<\![CDATA[ \boldmath\begin{eqnarray*} &&\mathrm{ZTransform}[\mathrm{Exp}[-b*n] ,n,z] \end{eqnarray*}

MSPFormat[ZTransform[Exp[-b n] ,n,z],StandardForm]
]]>
20070825094507
50 submitas Vykdyti

Prisiminkime, kad Laplace'o transformacijoje (žr. Pereinamieji virpesiai LC konture eksperimentą) laikas buvo atskaitomas nuo $t=0$ momento. Panašiai \twcal Z transformacijoje diskretinio kintamojo numeracija prasideda nuo $n=0$.

Paimkime dar vieną diskretinį signalą — diskretinio kosinuso funkciją $\cos (n \omega T)$: SF06<\![CDATA[ \boldmath$ \mathrm{ListPlot}[\mathrm{Table}[\{n,$ =0,Cos[n*w*T]/.{w->0.3,T->1},0]"]" > \boldmath$\},\{n,\fcolorbox[rgb]{1,0,0}{1,1,1}{$\vphantom{\vbox to 6pt{}}\smash{n_{start}}$},\fcolorbox[rgb]{1,0,0}{1,1,1}{$\vphantom{\vbox to 6pt{}}\smash{n_{end}}$}\}],\mathrm{AxesLabel}\to\{"n","x"\},\fcolorbox[rgb]{1,0,0}{1,1,1}{$\vphantom{\vbox to 4pt{}}\smash{parinktys}$}]$

Duomenys nuo pradinio \boldmath$n_{start}$=" > iki \boldmath$n_{end}$=" >
Piešimo parinktys: All, ImageSize -> 360,AxesOrigin->{0,0},PlotStyle->{RGBColor[0,0,1],PointSize[0.02]}}"]">

MSPBlock[{$$r1,$$r2,$$funcv,$$parinktys}, duomenys1 = Table[{n,$$funcv}, {n, $$r1, $$r2}]; MSPShow[ListPlot[duomenys1,AxesLabel->{"n","x"},Evaluate[Sequence@@$$parinktys]]] ]
]]>
20070825100033
61 submitas Vykdyti

Esant neigiamoms $n<0$ vertėms paprastumo dėlei laikysime (iš tiesų šioje srityje $x(n)$ nėra apibrėžta), kad $x(n)=0$. Pavaizduoto diskretinio signalo \twcal Z vaizdas yra SF07<\![CDATA[ \boldmath\begin{eqnarray*} &&zCos=\mathrm{ZTransform}[\mathrm{Cos}[-n*\omega*T] ,n,z]//\mathrm{Simplify} \end{eqnarray*}

MSPSessionVariable[zCos,Null]; MSPFormat[(zCos=Simplify[ZTransform[Cos[-n*\[Omega]*T] ,n,z]]),StandardForm]
]]>
20070825125138
70 submitas Vykdyti

Atvirkštinę \twcal Z transformaciją atlieka komanda \boldmath$\mathrm{InverseZTransform}[X(z),z,n], kuria iš gauto \twcal Z vaizdo atstatysime pirminę diskretinę funkciją: SF08<\![CDATA[ \boldmath\begin{eqnarray*} &&\mathrm{InverseZTransform}[zCos,n,z]//\mathrm{ExpToTrig}//\mathrm{Simplify} \end{eqnarray*}

MSPSessionVariable[zCos,Null]; MSPFormat[Simplify[ExpToTrig[InverseZTransform[zCos ,z,n]]],StandardForm]
]]>
20070825125138
80 submitas Vykdyti

Kadangi vadovėliuose \twcal Z transformacija kol kas sutinkama retai, išvardinsime svarbiausias šios transformacijos savybes [Doetsch67]. Jos yra labai panašios į Laplace'o transformacijos savybes.

Pažymėkime signalų $x_1 (n)$ ir $x_2 (n)$ vaizdus atitinkamai $X_1 (z)$ ir $X_2 (z)$: $x_1 (n)\Leftrightarrow X_1 (z)$, $x_2 (n)\Leftrightarrow X_2 (z)$. Tada signalus ir jų vaizdus sieja tokie ryšiai:

  • Tiesiškumas. Jei $a$ ir $b$ yra kompleksinės pastoviosios, tada $a x_1(n)+b x_2(n)\Leftrightarrow a X_1(z)+ b X_2(z)$.
  • Vėlavimas. Jei $n_0$ atitinka tam tikrą fiksuotą laiko momentą $t_0$, t.y. $n_0=t_0\,/T$, tada \twcal Z transformacijos $x(n)\Leftrightarrow X(z)$ pirminę ir vaizdo funkcijas pradiniu ir vėlesniu laiko momentu sieja vėlinimo teorema $x(n-n_0)\Leftrightarrow z^{-n_0}X(z)$.
  • Sąsukos vaizdas. Dviejų diskrečiųjų signalų $x$ ir $w$ sąsuka (convolution) vadinama $y(n)=\sum\limits_{m=-\infty }^\infty x(n) w(n-m)$ išraiška. Sąsukos \twcal Z vaizdas yra atskirų jos sandų \twcal Z vaizdų sandauga (sąsukos teorema): $Y(z)=X(z)W(z)$.

Jei sąsukos formulėje vietoje $x(n)$ įstatysime vienetinę funkciją, $x(n)=\delta(n)$ (kur $\delta(n)=1$, kai $n=0$ ir $\delta(n)=0$ visais kitais atvejais), tada $w(n-m)$ ir $W(z)$ yra vadinamos diskrečiojo signalo perdavimo charakteristikomis.

Dažninė charakteristika

Susipažinę su \twcal Z transformacija jau galime aptarti rekursinius ir nerekursinius skaitinius filtrus bei jų savybes. Svarbiausias filtro parametras yra jo dažninė charakteristika, kuri nusako praėjusio pro filtrą harmoninio signalo kompleksinės amplitudės priklausomybę nuo dažnio. Todėl naudinga žinoti, kaip signalo \twcal Z vaizdas yra susijęs su dažnine charakteristika. Kadangi griežto pagrindimo čia pateikti negalime (tam tektų pasitelkti diskretinės matematikos metodus), apsiribosime svarbiausių žingsnių aprašymu.

Pradžioje diskretizavimo impulsų seką $S(t)$ aproksimuojame Diraco $\delta $ funkcijomis: S(t)=\sum_{n=-\infty }^\infty \delta (t-n T).\tag{3} Čia $T$ yra diskretizavimo žingsnis (periodas). Pasinaudoję šia imties funkcija, tolydų signalą $x(t)$ galime pakeisti skaitmeniniu: x_s (t)=\sum _{n=-\infty }^\infty x(t) \delta (t-n T).\tag{4} Dabar apskaičiuokime diskretinio signalo $x_s (t)$ Fourier vaizdą: X_s(\mathrm{i} \omega )=\int_{-\infty }^\infty \!x_s(t) \mathrm{e}^{-\mathrm{i} \omega t}\, \mathrm{d} x=\sum _{n=-\infty }^\infty x(n T)\, \mathrm{e}^{-\mathrm{i} \omega n T}.\tag{5}

Jei signalą diskretizuoti pradedame ne nuo be galo nutolusio laiko momento $n=-\infty$, o nuo $n=0$, tada gautą $X_s(\mathrm{i} \omega)$ išraišką galime palyginti su signalo $x(n)$ \twcal Z vaizdo $X(z)$ apibrėžimu $X(z)=\sum _{n=0}^{n=\infty}x(n) z^{-n}$, įvestu skyriaus pradžioje. Palyginę matome, kad nuo diskretinto signalo \twcal Z vaizdo galima pereiti prie signalo Fourier vaizdo ir atvirkščiai, formulėse padarius pakeitimą $z\Leftrightarrow \mathrm{e}^{\mathrm{i} \omega T}$. Tokiu būdu formules, gautas po \twcal Z transformacijos, galima interpretuoti kaip signalo spektrą, jei jose atliksime pakeitimą z\rightarrow \mathrm{e}^{\mathrm{i} \omega T}=\mathrm{e}^{2\pi\mathrm{i} f/f_d}.\tag{6} Čia $f_d$ yra signalo diskretizavimo dažnis. Spektre pasirodantis didžiausias dažnis $f$ turi tenkinti nelygybę $ff_d/2)$ spektro dalimi, kurios virpesių periodas yra mažesnis už diskretizavimo periodą. Todėl signalą diskretizavus, skaitmeninis signalas jau neatspindės tikrųjų signalo savybių. Su šiuo reikalavimu jau buvome susidūrę aptardami Fourier transformaciją. Ribinis dydis $f_N=f_d/2=1/(2 T)$ yra vadinamas Nyquisto dažniu, o pats apribojimas — Kotelnikovo teorema. Iš čia seka svarbi praktinė išvada: prieš tolydinio signalo diskretizavimą, pavyzdžiui, AD konvertoriumi, jo spektre būtina pašalinti sandus, kurių dažnis aukštesnis už $f_N.$ Jei prieš diskretizavimą dažnių $f>f_N$ neeliminuosime, jie susimaišys su žemais dažniais ir spektras, o tuo pačiu ir sudiskretintas signalas bus iškraipytas.

Nerekursiniai skaitiniai filtrai

Toks keistas pavadinimas pasidarys aiškesnis, aptarus sudėtingesnius, taip vadinamus rekursinius, filtrus. Filtro sąvoka atkeliavo iš radioelektronikos. Ten filtru vadinamas įrenginys, kuris vienus priimto signalo spektrinius sandus slopina, o kitus praleidžia. Mes skaitiniu filtru vadinsime transformacijų seką, kurią atlikus skaitmeninio signalo spektras pasikeičia. Į filtrą patenkantį signalą žymėsime $x(n)$, o iš jo išeinantį — $y(n)$. Bendru atveju nerekursiniu filtru yra vadinama sąsuka y(n)=\sum_{k=-\infty}^\infty c(k) x(n-k).\tag{7} Joje $c(k)$ koeficientai nusako konkretaus filtro savybes ir yra parenkami iš anksto. Konkrečios $c(k)$ koeficientų vertės vadinamos svoriais (weights), o jų visuma — svorio funkcija (weight function). Sritis, kurioje koeficientai nėra lygūs nuliui, vadinama langu (window). Praktikoje langą gali sudaryti kelios dešimtys ar net šimtai $c(k)$ koeficientų. Panagrinėkime paprastą pavyzdį, paėmę penkis vienodus svorius: $c(1)=c(2)=c(3)=c(4)=c(5)=1/5$. Tada nerekursinis filtras atliks tokią signalo transformaciją: y(n)=\bigl( x(n-1)+x(n-2)+x(n-3)+x(n-4)+x(n-5)\bigr)\big/ 5 .\tag{8} Literatūroje šis algoritmas vadinamas ,,vidurkinimu iš penkių''. Pasinaudoję \twcal Z transformacija ir vėlinimo teorema, gauname tokį signalo atvaizdą: Y(z)=\frac15 \bigl( X(z)z^{-1}+X(z)z^{-2}+X(z)z^{-3}+X(z)z^{-4}+X(z)z^{-5}\bigr)= \frac15 \smash{\sum_{r=1}^5} z^{-r} X(z) .\tag{9} Iš čia seka, kad filtro perdavimo charakteristika $H(z)=Y(z)/X(z)$ yra $H(z)=\frac15 \sum_{r=1}^5 z^{-r}$. Dabar tą patį rezultatą apskaičiuosime, panaudoję Mathematica sistemos \twcal Z transformaciją. Pradžioje apibrėžiame filtro transformacijos funkciją: SF09<\![CDATA[ \boldmath\begin{eqnarray*} &&y[n\_]:=(x[n-1]+x[n-2]+x[n-3]+x[n-4]+x[n-5])/5 \end{eqnarray*}
]]>
20070825130504
90 submitas Vykdyti

Po \twcal Z transformacijos turime: SF10<\![CDATA[ \boldmath\begin{eqnarray*} &&\mathrm{Factor}[(\mathrm{ZTransform}[y[n],n,z])/.\,x[n\_/;n<0]\to 0] \end{eqnarray*}

MSPSessionVariable[Yz,Null]; y[n_]:=(x[n-1]+x[n-2]+x[n-3]+x[n-4]+x[n-5])/5; Yz = Factor[ZTransform[y[n], n, z] /. x[n_ /; n < 0] -> 0]; MSPFormat[Yz,StandardForm]
]]>
20070825201458
100 submitas Vykdyti

Keitimo taisykle \boldmath$x[n\_{}/;\,n<0]\to 0$ visas pradines sąlygas prilyginome nuliui. Išraiška \boldmath$\mathrm{ZTransform}[x[n],\,n,\,z]$ žymi \twcal Z vaizdą $X(z)$. Dažninę perdavimo charakteristiką gausime kintamąjį $z$ išraiškoje pakeitę eksponente $\mathrm{e}^{2\pi\mathrm{i} f/f_d}$. SF11<\![CDATA[ \boldmath\begin{eqnarray*} &&daznineCharakteristika = Yz\big/ \mathrm{ZTransform}[x[n], n, z] /. \,z \to \mathrm{Exp}[2*\pi\ *\mathrm{i}*f/f_d] // Expand \end{eqnarray*}

MSPSessionVariable[Yz,Null]; MSPSessionVariable[daznineCharakteristika,Null]; daznineCharakteristika = Expand[Yz/ZTransform[x[n], n, z] /. z -> Exp[2*Pi*I*(f/Subscript[f, d])]]; MSPFormat[daznineCharakteristika,StandardForm]
]]>
20070825202743
110 submitas Vykdyti

Ją vizualizuosime, brėžinyje atidėję amplitudės priklausomybę nuo signalo dažnio. Signalo diskretizavimo dažnį imsime lygų vienai atvirkštinei sekundei: $f_d=1\ s{}^{-1}$. SF12<\![CDATA[ \boldmath\begin{eqnarray*} &&\mathrm{Plot}[\mathrm{Evaluate}[\mathrm{Abs}[ daznineCharakteristika /. f_d \to \fcolorbox[rgb]{1,0,0}{1,1,1}{$\vphantom{\vbox to 4pt{}}\smash{f_{d0}}$} ]], \{f, \fcolorbox[rgb]{1,0,0}{1,1,1}{$\vphantom{\vbox to 4pt{}}\smash{f_{start}}$}, \fcolorbox[rgb]{1,0,0}{1,1,1}{$\vphantom{\vbox to 4pt{}}\smash{f_{end}}$}\}, \\ &&\mathrm{AxesLabel}\to\{"\omega","Ampl."\},\fcolorbox[rgb]{1,0,0}{1,1,1}{$\vphantom{\vbox to 4pt{}}\smash{parinktys}$}]; \end{eqnarray*}

Signalo diskretizavimo dažnis \boldmath$f_{d0}$=" >
Dažnis nuo pradinio \boldmath$f_{start}$=" > iki \boldmath$f_{end}$=" >
Piešimo parinktys: All, ImageSize -> 360,PlotStyle->Thickness[0.01]}"]">

MSPSessionVariable[daznineCharakteristika,Null]; MSPBlock[{$$r1,$$r2,$$fd,$$parinktys}, MSPShow[Plot[Evaluate[Abs[daznineCharakteristika /. Subscript[f, d] -> $$fd]], {f, $$r1, $$r2},AxesLabel->{"\[Omega]","Ampl."}, Evaluate[Sequence@@$$parinktys]]] ]
]]>
20070827150531
121 submitas Vykdyti
Matome, kad toks filtras praleidžia žemus, t.y. gerokai žemesnio nei diskretizavimo dažnis $f_d$ periodinius signalus. Antra vertus, toks filtras visiškai nepraleidžia signalų, kurių dažnis yra $5$ ir $2{,}5$ karto mažesnis už $f_d$. Taigi, ,,vidurkinimo iš penkių'' filtras turi gana keistą dažninę charakteristiką, kurią intuityviai vargu ar buvo galima nuspėti. Deja, šis filtras nėra geras, nes nepakankamai slopina aukšto dažnio sandus.

Realus žemų dažnių filtras

Išnagrinėtame pavyzdyje paėmėme tik penkis nelygius nuliui svorius. Realiame filtre svorio funkciją sudaro kur kas daugiau koeficientų. Kaip reiktų šiuos koeficientus parinkti, kad skaitinis filtras turėtų pageidaujamą dažninę charakteristiką, aprašyta straipsnyje [Kaiser77]. Mes panagrinėsime žemų dažnių filtrą, t.y. tokį filtrą, kuris praleidžia tik žemus dažnius. Sekdami minėtu straipsniu, imsime Gausso pavidalo svorio funkciją su nelyginio dydžio $(2\, \textbf{np}+1)$ langu ($\textbf{np}=25$): SF13<\![CDATA[ \boldmath$ svorioFunkcijaNenormuota=\mathrm{Table}[\mathrm{N}[$ " > \boldmath$], \{m,\fcolorbox[rgb]{1,0,0}{1,1,1}{$\vphantom{\vbox to 6pt{}}\smash{np_{start}}$},\fcolorbox[rgb]{1,0,0}{1,1,1}{$\vphantom{\vbox to 6pt{}}\smash{np_{end}}$}\}]$

Duomenys nuo pradinio \boldmath$np_{start}$=" > iki \boldmath$np_{end}$=" >

MSPSessionVariable[svorioFunkcijaNenormuota,Null]; MSPBlock[{$$r1,$$r2,$$funcv}, svorioFunkcijaNenormuota =Table[N[$$funcv], {m, $$r1, $$r2}]; MSPFormat[Shallow[svorioFunkcijaNenormuota],StandardForm] ]
]]>
20070825210448
131 submitas Vykdyti

Svorio funkciją normuokime į vienetą. SF14<\![CDATA[ \boldmath\begin{eqnarray*} &&svorioFunkcija= svorioFunkcijaNenormuota\big/\mathrm{Apply}[\mathrm{Plus},svorioFunkcijaNenormuota]; \end{eqnarray*}
]]>
20070826155249
140 submitas Vykdyti

Filtras su normuota svorio funkcija nekeičia nuostoviosios (nuo laiko nepriklausančios) signalo dalies, nes visų jo koeficientų suma lygi vienetui. SF15<\![CDATA[ \boldmath\begin{eqnarray*} &&\mathrm{Apply}[\mathrm{Plus},svorioFunkcija] \end{eqnarray*}

MSPSessionVariable[svorioFunkcijaNenormuota,Null]; MSPSessionVariable[svorioFunkcija,Null]; If[svorioFunkcijaNenormuota=!=Null, svorioFunkcija= svorioFunkcijaNenormuota/Apply[Plus,svorioFunkcijaNenormuota]; MSPFormat[Apply[Plus,svorioFunkcija],StandardForm]]
]]>
20070828091521
150 submitas Vykdyti

Atvaizduokime svorio funkcijos koeficientų priklausomybę nuo jų numerio: SF17<\![CDATA[ \boldmath\begin{eqnarray*} &&\mathrm{ListPlot}[svorioFunkcija,\mathrm{AxesLabel}\to\{"n",""\},\mathrm{PlotLabel}\to"Svorio funkcija"]; \end{eqnarray*}

Piešimo parinktys: All, ImageSize -> 360, PlotStyle->{RGBColor[0,0,1],PointSize[0.02]}}"]">

MSPSessionVariable[svorioFunkcija,Null]; MSPBlock[{$$parinktys}, If[svorioFunkcija=!=Null, MSPShow[ListPlot[svorioFunkcija, AxesLabel->{"n",""},PlotLabel->"Svorio funkcija",Evaluate[Sequence@@$$parinktys]]] ]]
]]>
20070828091521
171 submitas Vykdyti

Kaip buvo aptarta, filtro su šia svorio funkcija dažninę charakteristiką rasime apskaičiavę filtro perdavimo charakteristikos \twcal Z vaizdą. Tam pradžioje komanda \boldmath$\mathrm{Table}[~]$ sudarykime vektorių sąrašą $\bigl\{x(n-1),\dotsc ,x(n-2 \textbf{np}-1)\bigr\}$. \boldmath$\mathrm{Dot}[~]$ komanda jį padauginę iš svorio funkcijos sąrašo ir padalinę iš $X(z)\equiv$ \boldmath$\mathrm{ZTransform}[x[n], n, x]$, gausime perdavimo charakteristiką. SF18<\![CDATA[ \boldmath\begin{eqnarray*} &&xVektorius=\mathrm{Table}[x[n-k],\{k,1,\mathrm{Length}[svorioFunkcija]\}];\\ &&zCharakteristika=\mathrm{Expand}[\frac{\mathrm{ZTransform}[svorioFunkcija\ .\, xVektorius,n,z]/.\,x[n\_/;n<0]\to 0}{\mathrm{ZTransform}[x[n],n,z]}] \end{eqnarray*}

MSPSessionVariable[svorioFunkcija,Null]; MSPSessionVariable[zCharakteristika,Null]; If[svorioFunkcija=!=Null, xVektorius = Table[x[n - k], {k, 1, Length[svorioFunkcija]}]; zCharakteristika = Expand[(ZTransform[svorioFunkcija . xVektorius, n, z] /. x[n_ /; n < 0] -> 0)/ZTransform[x[n], n, z]]; MSPFormat[zCharakteristika,StandardForm]]
]]>
20070828091521
180 submitas Vykdyti

Kintamąjį $z$ pakeitę eksponente $\mathrm{e}^{2\pi\mathrm{i} f/f_d}$, rasime dažninę filtro charakteristiką. Jos amplitudės priklausomybė nuo dažnio yra SF19<\![CDATA[ \boldmath\begin{eqnarray*} &&\mathrm{Plot}[\mathrm{Evaluate}[\mathrm{Abs}[ daznineCharakteristika /. f_d \to \fcolorbox[rgb]{1,0,0}{1,1,1}{$\vphantom{\vbox to 4pt{}}\smash{f_{d0}}$} ]], \{f, \fcolorbox[rgb]{1,0,0}{1,1,1}{$\vphantom{\vbox to 4pt{}}\smash{f_{start}}$}, \fcolorbox[rgb]{1,0,0}{1,1,1}{$\vphantom{\vbox to 4pt{}}\smash{f_{end}}$}\}, \\ &&\mathrm{FrameLabel}\to\{"\omega","Ampl."\},\fcolorbox[rgb]{1,0,0}{1,1,1}{$\vphantom{\vbox to 4pt{}}\smash{parinktys}$}]; \end{eqnarray*}

Signalo diskretizavimo dažnis \boldmath$f_{d0}$=" >
Dažnis nuo pradinio \boldmath$f_{start}$=" > iki \boldmath$f_{end}$=" >
Piešimo parinktys: All, ImageSize -> 360,Frame -> True,GridLines -> Automatic, PlotStyle->Thickness[0.01]}"]">

MSPSessionVariable[daznineCharakteristika,Null]; MSPSessionVariable[zCharakteristika,Null]; MSPBlock[{$$r1,$$r2,$$fd,$$parinktys}, daznineCharakteristika = zCharakteristika /. z -> Exp[2*Pi*I*(f/Subscript[f, d])]; MSPShow[Plot[Chop[Evaluate[Abs[daznineCharakteristika /. Subscript[f, d] -> $$fd]]], {f, $$r1, $$r2},FrameLabel->{"\[Omega]","Ampl."}, Evaluate[Sequence@@$$parinktys]]] ]
]]>
20070826193713
191 submitas Vykdyti

Šis filtras labai stipriai (daugiau kaip tūkstantį kartų) slopina aukšto dažnio sandus. Tuo įsitikinsime, amplitudę atidėję ne visame dažnių ruože $0\cdots 0{,}5$, o tik ties aukštesniaisiais dažniais: $f=0{,}1\cdots 0{,}5$. Skaitytojas tegu pats įrašo šias dažnių vertes ankstesnėje interaktyvioje ląstelėje.

Paeksperimentuokime su ką tik sukonstruotu nerekursiniu filtru. Paduokime į filtrą žemo dažnio signalą $1-\cos(n \pi /90)$, prieš tai jį ,,užteršę'' triukšmu, ir pažiūrėkime, kaip filtras šį triukšmą pašalina. Triukšmą imituosime atsitiktinių skaičių generatoriumi \boldmath$\mathrm{Random}[~]$, kuris pridės mažus atsitiktinius signalo iškraipymus. Gautą suminį signalą vizualizuojame. SF20<\![CDATA[ \boldmath\begin{eqnarray*} &&\mathrm{ListPlot}[(sumSignalas=\mathrm{Table}[$ ">\boldmath$,\{n,\fcolorbox[rgb]{1,0,0}{1,1,1}{$\vphantom{\vbox to 6pt{}}\smash{n_{start}}$},\fcolorbox[rgb]{1,0,0}{1,1,1}{$\vphantom{\vbox to 6pt{}}\smash{n_{end}}$}\}]),\\ &&\mathrm{AxesLabel}\to\{"n",""\},\mathrm{PlotLabel}\to"Pradinis\ signalas", \fcolorbox[rgb]{1,0,0}{1,1,1}{$\vphantom{\vbox to 4pt{}}\smash{parinktys}$}] \end{eqnarray*}

Duomenys nuo pradinio \boldmath$n_{start}$=" > iki \boldmath$n_{end}$=" >
Piešimo parinktys: All, ImageSize -> 360,PlotStyle->{RGBColor[0,0,1],PointSize[0.01]}}"]">

MSPSessionVariable[sumSignalas,Null]; MSPBlock[{$$r1,$$r2,$$funcv,$$parinktys}, sumSignalas = Table[$$funcv, {n, $$r1, $$r2}]; MSPShow[ListPlot[sumSignalas,AxesLabel->{"n",""},PlotLabel->"Pradinis signalas",Evaluate[Sequence@@$$parinktys]]] ]
]]>
20070826201017
201 submitas Vykdyti

Taigi, turime du skaičiais užpildytus sąrašus — filtro svorių sąrašą \boldmath$svorioFunkcija$ ir signalo su triukšmu sąrašą \boldmath$sumSignalas$, pavaizduotą paveiksle. Abiejuose sąrašuose esančių elementų skaičių sužinosime, komanda \boldmath$\mathrm{Length}[~]$ ,,išmatavę'' jų ilgį. SF21<\![CDATA[ \boldmath\begin{eqnarray*} &&\{\mathrm{Length}[svorioFunkcija],\mathrm{Length}[sumSignalas]\} \end{eqnarray*}

MSPSessionVariable[sumSignalas,Null]; MSPSessionVariable[svorioFunkcija,Null]; If[sumSignalas=!=Null&&svorioFunkcija=!=Null, MSPFormat[{Length[svorioFunkcija],Length[sumSignalas]},StandardForm]]
]]>
20070829093650
210 submitas Vykdyti

Kadangi signalo ir svorio funkcijos elementai gali būti dauginami tik lango ribose, komanda \boldmath$\mathrm{Take}[~]$ iš sąrašo \boldmath$sumSignalas$ imsime tik tuos narius, kurie duotu momentu pakliūva į lango sritį. Sąsuką skaičiuosime, langą ,,traukdami'' per sąrašo \boldmath$sumSignalas$ elementus nuo $\textbf{n}=1$ iki \boldmath$\mathrm{Length}[sumSignalas]-\mathrm{Length}[svorioFunkcija]$. SF22<\![CDATA[ \boldmath\begin{eqnarray*} &&sasuka[n\_]:=svorioFunkcija\ .\, \mathrm{Take}[sumSignalas,\{n,\mathrm{Length}[svorioFunkcija]+n-1\}]\\ &&signalasPoFiltro=\mathrm{Table}[sasuka[n],\{n,\mathrm{Length}[sumSignalas]-\mathrm{Length}[svorioFunkcija]\}]];\\ &&\mathrm{ListPlot}[signalasPoFiltro,\mathrm{PlotLabel}\to"Signalas\ po\ nerekurs.\ filtro",\\ &&\qquad \mathrm{AxesLabel}\to\{"n",""\}, \fcolorbox[rgb]{1,0,0}{1,1,1}{$\vphantom{\vbox to 4pt{}}\smash{parinktys}$}] \end{eqnarray*}

Piešimo parinktys: All, ImageSize -> 360,PlotStyle->{RGBColor[0,0,1],PointSize[0.015]}}"]">

MSPSessionVariable[sumSignalas,Null]; MSPSessionVariable[svorioFunkcija,Null]; If[sumSignalas=!=Null&&svorioFunkcija=!=Null, sasuka[n_] := svorioFunkcija . Take[sumSignalas, {n, Length[svorioFunkcija] + n - 1}]; signalasPoFiltro = Table[sasuka[n], {n, 1, Length[sumSignalas] - Length[svorioFunkcija]}]; MSPBlock[{$$parinktys}, MSPShow[ListPlot[signalasPoFiltro, AxesLabel -> {"n", ""},PlotLabel->"Signalas po nerekurs. filtro",Evaluate[Sequence@@$$parinktys]]] ]]
]]>
20070829093650
221 submitas Vykdyti

Matome, kad signalui praėjus pro nerekursinį filtrą iš pradinio signalo dingo aukšto dažnio triukšmas. Kadangi kraštuose sąsukos apskaičiuoti negalime, tai išfiltruoto signalo kraštus iš abiejų pusių nupjovėme per lango plotį. Reiktų pasakyti, kad toks paprastas sąsukos skaičiavimo algoritmas yra pernelyg lėtas. Ketvirtoje Mathematica versijoje panašioms transformacijoms atlikti įvesta greita komanda \boldmath$\mathrm{ListConvolve}[~]$. Apskaičiuosime mūsų pavyzdį dar kartą, panaudoję šią funkciją. SF23<\![CDATA[ \boldmath\begin{eqnarray*} &&\mathrm{ListPlot}[\mathrm{ListConvolve}[svorioFunkcija,sumSignalas],\mathrm{PlotLabel}\to"Signalas\ po\ nerekurs.\ filtro",\\ &&\qquad \mathrm{AxesLabel}\to\{"n",""\}, \fcolorbox[rgb]{1,0,0}{1,1,1}{$\vphantom{\vbox to 4pt{}}\smash{parinktys}$}] \end{eqnarray*}

Piešimo parinktys: All, ImageSize -> 360,PlotStyle->{RGBColor[0,0,1],PointSize[0.015]}}"]">

MSPSessionVariable[sumSignalas,Null]; MSPSessionVariable[svorioFunkcija,Null]; If[sumSignalas=!=Null&&svorioFunkcija=!=Null, MSPBlock[{$$parinktys}, MSPShow[ListPlot[ListConvolve[svorioFunkcija,sumSignalas], AxesLabel -> {"n", ""},PlotLabel->"Signalas po nerekurs. filtro",Evaluate[Sequence@@$$parinktys]]] ]]
]]>
20070828095659
231 submitas Vykdyti
Matome, kad pirmasis funkcijos \boldmath$\mathrm{ListConvolve}[~]$ argumentas yra svorio funkcijos koeficientų sąrašas (transformacijos branduolys), o antrasis — apdorojamas signalas. Galimi papildomi argumentai leidžia įvairiausiais būdais atsižvelgti į čia nupjautų kraštų transformacijas ir dar daugelį kitų dalykų, kuriuos perprasti paliksime pačiam skaitytojui.

Filtravimas Fourier transformacija

Praeitame eksperimente aptartą diskretinę Fourier transformaciją taip pat galima pritaikyti triukšmui filtruoti. Iliustracijai pasinaudosime S.Wolframo knygoje [Wolfram99] aprašytu pavyzdžiu. Filtravimas Fourier transformacija yra pagrįstas labai panašia sąsukos teorema, kuri teigia, kad dviejų funkcijų sąsukos Fourier vaizdas yra lygus pirminių funkcijų Fourier vaizdų sandaugai. Pasinaudoję šia savybe dar kartą išfiltruosime triukšmą iš jau nagrinėto signalo. Skaičiavimų eigoje mums teks panariui dauginti svorio funkcijos ir signalo Fourier vaizdų sąrašus. Tam turime suvienodinti abiejų sąrašų ilgius. Paprasčiausias sprendimas būtų svorio funkcijos langą iš kairės ir dešinės simetriškai papildyti nuliais: SF24<\![CDATA[ \boldmath\begin{eqnarray*} &&nuliai=\mathrm{Table}[ 0,\{\mathrm{Floor}[\frac{1}{2} (\fcolorbox[rgb]{1,0,0}{1,1,1}{$\vphantom{\vbox to 4pt{}}\smash{n_{end}}$}-\mathrm{Length}[svorioFunkcijaNenormuota])] \}];\\ &&kernPlot=\mathrm{Flatten}[\{0,nuliai,svorioFunkcijaNenormuota, nuliai\}];\\ &&\mathrm{ListPlot}[kernPlot, PlotLabel\to "Nenormuota\ svorio\ funkcija,\\ &&\qquad \mathrm{AxesLabel}\to\{"n",""\}, \fcolorbox[rgb]{1,0,0}{1,1,1}{$\vphantom{\vbox to 4pt{}}\smash{parinktys}$}] \end{eqnarray*}

\boldmath$n_{end}$=" >
Piešimo parinktys: All, ImageSize -> 360,PlotStyle->{RGBColor[0,0,1],PointSize[0.02]}}"]">

MSPSessionVariable[kernPlot,Null]; MSPSessionVariable[svorioFunkcijaNenormuota,Null]; If[svorioFunkcijaNenormuota=!=Null, MSPBlock[{$$parinktys,$$r2}, nuliai = Table[0, {Floor[($$r2 - Length[svorioFunkcijaNenormuota])/2]}]; kernPlot = Flatten[{0, nuliai, svorioFunkcijaNenormuota, nuliai}]; MSPShow[ListPlot[kernPlot, AxesLabel -> {"n", ""}, PlotLabel -> "Nenormuota svorio funkcija",Evaluate[Sequence@@$$parinktys]]] ]]
]]>
20070828093454
241 submitas Vykdyti

Tikrąjį Fourier filtro branduolį gausime svorio funkcijos elementus cikliškai perstatę taip, kad didžiausio svorio elementas atsidurtų pirmoje pozicijoje (to reikalauja Mathematica'os diskretinės Fourier transformacijos algoritmas). Pertvarkytą sąrašą, kaip ir anksčiau, normuojame į vienetą. SF25<\![CDATA[ \boldmath\begin{eqnarray*} &&kern=\frac{\mathrm{RotateLeft}[ kernPlot,\mathrm{Floor}[\frac{\mathrm{Length}[ kernPlot]}{2}] ]}{Plus@@kernPlot};\\ &&\mathrm{ListPlot}[kern, PlotLabel\to "Normuota\ svorio\ funkcija, \mathrm{AxesLabel}\to\{"n",""\}, \fcolorbox[rgb]{1,0,0}{1,1,1}{$\vphantom{\vbox to 4pt{}}\smash{parinktys}$}] \end{eqnarray*}

Piešimo parinktys: All, ImageSize -> 360,PlotStyle->{RGBColor[0,0,1],PointSize[0.02]}}"]">

MSPSessionVariable[kernPlot,Null]; MSPSessionVariable[kern,Null]; If[kernPlot=!=Null, MSPBlock[{$$parinktys}, kern = RotateLeft[kernPlot, Floor[Length[kernPlot]/2]]/Plus @@ kernPlot; MSPShow[ListPlot[kern, AxesLabel -> {"n", ""}, PlotLabel -> "Normuota svorio funkcija",Evaluate[Sequence@@$$parinktys]]] ]]
]]>
20070828095027
251 submitas Vykdyti

Dabar jau galima iš signalo eliminuoti triukšmą. Tam pradžioje rasime triukšmu užteršto signalo ir naujosios svorio funkcijos \boldmath$kern$ Fourier vaizdus. Vaizdus sudauginę ir atlikę atvirkštinę Fourier transformaciją (bei viską papildomai padauginę iš $\sqrt{N}$, kur $N$ yra signalo sąrašo ilgis) pamatysime, kad gavome (jei neatsižvelgsime į signalo kraštus) tą patį rezultatą. SF26<\![CDATA[ \boldmath\begin{eqnarray*} &&signalasPraejesFourierFiltra=\mathrm{InverseFourier}[\\ &&\qquad \sqrt{\mathrm{Length}[kern]} \mathrm{Fourier}[sumSignalas] \mathrm{Fourier}[kern]];\\ &&\mathrm{ListPlot}[\mathrm{Chop}[signalasPraejesFourierFiltra],PlotLabel\to "Signalas\ po\ Fourier\ filtro", \mathrm{AxesLabel}\to\{"n",""\}, \fcolorbox[rgb]{1,0,0}{1,1,1}{$\vphantom{\vbox to 4pt{}}\smash{parinktys}$}] \end{eqnarray*}

Piešimo parinktys: All, ImageSize -> 360,PlotStyle->{RGBColor[0,0,1],PointSize[0.02]}}"]">

MSPSessionVariable[kern,Null]; MSPSessionVariable[sumSignalas,Null]; If[sumSignalas=!=Null&&kern=!=Null, MSPBlock[{$$parinktys}, signalasPraejesFourierFiltra = InverseFourier[Sqrt[Length[kern]]*Fourier[sumSignalas]*Fourier[kern]]; MSPShow[ListPlot[Chop[signalasPraejesFourierFiltra], AxesLabel -> {"n", ""}, PlotLabel -> "Signalas po Fourier filtro",Evaluate[Sequence@@$$parinktys]]] ]]
]]>
20070828095027
261 submitas Vykdyti
Taigi, ir nerekursinis ir diskretinės Fourier transformacijos filtrai iš esmės duoda tą patį rezultatą. Tačiau skaičiavimo metodo požiūriu abi filtravimo procedūros labai skiriasi: kai naudojame nerekursinį filtrą, duomenų srautą galime filtruoti realiame laike, duomenims ,,bėgant'' pro filtrą. Tuo tarpu Fourier transformacijos atveju duomenis pradžioje reikėjo sukaupti ir tik po to galėjome juos išvalyti nuo triukšmo. Tokiu būdu nerekursinis filtras šiuo požiūriu daug pranašesnis, nes leidžia apdoroti signalą jo registravimo eigoje. Jo privalumas taps neginčytinas, jei bus priimamas ne šiaip sau, o pavojaus signalas.

Pereinamoji nerekursinio filtro charakteristika

Ištirsime pereinamąją mūsų sukonstruoto nerekursinio filtro charakteristiką. Ji nusako filtro atsaką laiko ašyje, kai į filtrą pasiunčiamas diskretinis laiptelio pavidalo signalas. Signalo formą prieš ir po filtravimo pavaizduosime, panaudoję \boldmath$\mathrm{ListPlot}[~]$. Abu signalų brėžinius sujungsime \boldmath$\mathrm{GraphicsArray}[~]$ komanda. FS27<\![CDATA[ \boldmath\begin{eqnarray*} &&laiptelis=\mathrm{Flatten}[\{\mathrm{Table}[0,\{\fcolorbox[rgb]{1,0,0}{1,1,1}{$\vphantom{\vbox to 6pt{}}\smash{kiekNuliu}$}\}],\mathrm{Table}[1,\{\fcolorbox[rgb]{1,0,0}{1,1,1}{$\vphantom{\vbox to 6pt{}}\smash{taskai}$}-\fcolorbox[rgb]{1,0,0}{1,1,1}{$\vphantom{\vbox to 6pt{}}\smash{kiekNuliu}$}\}]\}];\\ &&\mathrm{Show}[\mathrm{GraphicsArray}[\{\\ &&\qquad \mathrm{ListPlot}[laiptelis,\mathrm{PlotLabel}\to "Laiptelis", \mathrm{AxesLabel}\to \{"t", ""\}, \mathrm{DisplayFunction}\to \mathrm{Identity},\fcolorbox[rgb]{1,0,0}{1,1,1}{$\vphantom{\vbox to 6pt{}}\smash{parinktys1}$}],\\ &&\qquad \mathrm{ListPlot}[\mathrm{ListConvolve}[svorioFunkcija, laiptelis], \mathrm{PlotLabel}\to "Laiptelis\ po\ nerekursinio\ filtro", \mathrm{AxesLabel}\to \{n, ""\}, \\[2pt] &&\qquad\qquad\mathrm{DisplayFunction}\to \mathrm{Identity},\fcolorbox[rgb]{1,0,0}{1,1,1}{$\vphantom{\vbox to 6pt{}}\smash{parinktys2}$}]\\ &&\}],\mathrm{DisplayFunction}\to \\$\mathrm{DisplayFunction},\fcolorbox[rgb]{1,0,0}{1,1,1}{$\vphantom{\vbox to 6pt{}}\smash{parinktys3}$}]; \end{eqnarray*}

Laiptelio ilgis: \boldmath$taskai=$" >. Iš jų nulių ilgis: \boldmath$kiekNuliu=$" >
Pirmojo grafiko piešimo parinktys: All,PlotStyle->{RGBColor[0,0,1],PointSize[0.01]}}"]" >
Antrojo grafiko piešimo parinktys: All,PlotStyle->{RGBColor[0,0,1],PointSize[0.01]}}"]" >
Grafikų grupės piešimo parinktys: 720}"]" >
MSPSessionVariable[svorioFunkcija,Null]; If[svorioFunkcija=!=Null, MSPBlock[{$$taskai,$$nuliai,$$parinktys1,$$parinktys2,$$parinktys3}, laiptelis = Flatten[{Table[0, {$$nuliai}], Table[1, {$$taskai - $$nuliai}]}]; MSPShow[Show[GraphicsArray[{ListPlot[laiptelis, AxesLabel -> {"t", ""}, PlotLabel -> "Laiptelis",DisplayFunction -> Identity,Evaluate[Sequence@@$$parinktys1]], ListPlot[ListConvolve[svorioFunkcija, laiptelis], AxesLabel -> {"n", ""}, PlotLabel -> "Laiptelis po nerekursinio filtro", DisplayFunction -> Identity,Evaluate[Sequence@@$$parinktys2]]}], DisplayFunction -> $DisplayFunction,Evaluate[Sequence@@$$parinktys3]] ]] ]
]]>
20070830100859
271 submitas Vykdyti
Kaip matome, praėjęs filtrą laiptelio pavidalo frontas ištįso. Taip atsitiko todėl, kad staigų signalo šuoliuką lemia aukšto dažnio harmoniniai sandai. Kadangi filtras nuslopino signalo aukšto dažnio sandus, nufiltruoto signalo frontas pasidarė glodesnis.

Rekursiniai skaitiniai filtrai

Kaip matėme, nerekursinis filtras duotu laiko momentu apdoroja tik tuos duomenis, kurie pakliūna į filtro langą. Rekursiniame filtre, priešingai, duomenų apdorojimą įtakoja ne tik į langą pakliūnantys, bet ir anksčiau apdoroti duomenys. Iš čia ir pavadinimas ,,rekursinis'', t.y. grįžtantis į tą patį ,,kursą'' ar kelią. Bendru atveju rekursiniu filtru vadinama tokia transformacija: y(n)=\sum_k c(k) x(n-k) \sum_k d(k) y(n-k) .\tag{10} Čia $c(k)$ ir $d(k)$ yra filtrą apibūdinantys koeficientai (svoriai). Iš formulės matyti, kad ieškomasis $y(n)$ priklauso nuo jau apdorotų $y(n-k)$ signalo verčių, paimtų su tam tikrais svoriais $d(k)$. Priklausomai nuo $d(k)$ ženklo ir dydžio, grįžtamasis ryšis gali būti teigiamas arba neigiamas, t.y. įšėjimo signalas gali būti stiprinamas arba silpninamas. Iš rekursinio filtro apibrėžimo taip pat seka (sąsukos formulė), kad po \twcal Z transformacijos filtro perdavimo charakteristika bendru atveju turi $H(z)=C(z)/D(z)$ pavidalą. Čia $C(z)$ ir $D(z)$ yra kintamojo $z$ polinomai. Jei keitiniu $z\rightarrow \mathrm{e}^{\mathrm{i}\omega T}$ pereisime į dažnių atvaizdį, gausime formulę H(\mathrm{i}\omega )=C(\mathrm{i}\omega )/D(\mathrm{i}\omega ) . Kadangi daugelį polinomų $C(\mathrm{i}\omega )$ ir $D(\mathrm{i}\omega )$ savybių lemia jų šaknų skaičius, tai ir rekursinio filtro esmines savybes nusako perdavimo charakteristikos nuliai ir poliai, t.y. skaitiklio ir vardiklio šaknys. Todėl projektuojant filtrą pirmiausia reikia nustatyti šaknų skaičių, ir tik po to ieškoti tinkamų svorio verčių. Deja, čia neturime vietos šiems praktiniams klausimams nagrinėti [Krivickas84,Rabiner75]. Apsiribosime nesudėtingu pavyzdžiu, iliustruojančiu metodo esmę. Panagrinėkime paprastą rekursinį filtrą: $y(n)=x(n)+x(n-1)+y(n-1)/2$. Po \twcal Z transformacijos iš šios lygties gauname

SF28<\![CDATA[ \boldmath\begin{eqnarray*} &&\mathrm{Clear}[y];\\ &&zAtvaizdas=\mathrm{ZTransform}[y[n],n, z]=\mathrm{ZTransform}[x[n-1]+x[n]+\frac{1}{2} y[n-1],n, z]/.\\ &&\{y[n\_/;n<0]\to 0, x[n\_/;n<0]\to 0,\mathrm{ZTransform}[y[n],n,z]\to Y[z],\mathrm{ZTransform}[x[n],n,z]\to X[z]\} \end{eqnarray*}

MSPSessionVariable[zAtvaizdas,Null]; Clear[y]; zAtvaizdas = (ZTransform[y[n], n, z] == ZTransform[x[n] + x[n - 1] + y[n - 1]/2, n, z] /. {y[n_ /; n < 0] -> 0, x[n_ /; n < 0] -> 0, ZTransform[y[n], n, z] -> Y[z], ZTransform[x[n], n, z] -> X[z]}); MSPFormat[zAtvaizdas,StandardForm]
]]>
20070829173253
280 submitas Vykdyti

Išsprendę atsakymą $Y(z)$ atžvilgiu ir padalinę jį iš $X(z)$, randame perdavimo charakteristiką SF29<\![CDATA[ \boldmath\begin{eqnarray*} &&yz=\mathrm{Solve}[zAtvaizdas,Y[z]];\\ &&zPerdavimoCharakteristika=\mathrm{Simplify}[\frac{Y[ z]/.\mathrm{First}[yz]}{X[z]}] \end{eqnarray*}

MSPSessionVariable[zAtvaizdas,Null]; MSPSessionVariable[zPerdavimoCharakteristika,Null]; yz = Solve[zAtvaizdas, Y[z]]; zPerdavimoCharakteristika = Simplify[(Y[z] /. First[yz])/X[z]]; MSPFormat[zPerdavimoCharakteristika,StandardForm]
]]>
20070830090718
290 submitas Vykdyti

Šią \twcal Z perdavimo charakteristiką atitinka dažninė: SF30<\![CDATA[ \boldmath\begin{eqnarray*} &&fPerdavimoCharakteristika= zPerdavimoCharakteristika/.z\to \mathrm{Exp}[\frac{2 \pi \mathrm{i} f}{f_d}] \end{eqnarray*}

MSPSessionVariable[zPerdavimoCharakteristika,Null]; MSPSessionVariable[fPerdavimoCharakteristika,Null]; fPerdavimoCharakteristika = zPerdavimoCharakteristika /. z -> Exp[2*Pi*I*(f/Subscript[f, d])]; MSPFormat[fPerdavimoCharakteristika,StandardForm]
]]>
20070830092139
300 submitas Vykdyti

Atvaizduokime jos amplitudės priklausomybę nuo dažnio $f=\omega/2 \pi$, kai $f_d=1$. SF31<\![CDATA[ \boldmath\begin{eqnarray*} &&\mathrm{Plot}[\mathrm{Evaluate}[\mathrm{Abs}[ fPerdavimoCharakteristika /. f_d \to \fcolorbox[rgb]{1,0,0}{1,1,1}{$\vphantom{\vbox to 4pt{}}\smash{f_{d0}}$} ]], \{f, \fcolorbox[rgb]{1,0,0}{1,1,1}{$\vphantom{\vbox to 4pt{}}\smash{f_{start}}$}, \fcolorbox[rgb]{1,0,0}{1,1,1}{$\vphantom{\vbox to 4pt{}}\smash{f_{end}}$}\}, \\ &&\mathrm{AxesLabel}\to\{"\omega","Ampl."\},\fcolorbox[rgb]{1,0,0}{1,1,1}{$\vphantom{\vbox to 4pt{}}\smash{parinktys}$}]; \end{eqnarray*}

Signalo diskretizavimo dažnis \boldmath$f_{d0}$=" >
Dažnis nuo pradinio \boldmath$f_{start}$=" > iki \boldmath$f_{end}$=" >
Piešimo parinktys: All, ImageSize -> 360,PlotStyle->Thickness[0.01]}"]">

MSPSessionVariable[fPerdavimoCharakteristika,Null]; If[fPerdavimoCharakteristika=!=Null, MSPBlock[{$$r1,$$r2,$$fd,$$parinktys}, MSPShow[Plot[Evaluate[Abs[fPerdavimoCharakteristika /. Subscript[f, d] -> $$fd]], {f, $$r1, $$r2},AxesLabel->{"\[Omega]","Ampl."}, Evaluate[Sequence@@$$parinktys]]] ]]
]]>
20070830095757
311 submitas Vykdyti

Tokiu būdu nagrinėjamas rekursinis filtras stiprina harmoninį signalą, kai jo dažnis $f/f_d<0{,}2902$: SF32<\![CDATA[ \boldmath\begin{eqnarray*} &&FindRoot[Abs[fPerdavimoCharakteristika /. f_d \to \fcolorbox[rgb]{1,0,0}{1,1,1}{$\vphantom{\vbox to 4pt{}}\smash{f_{d0}}$}] == 1, \{f, \fcolorbox[rgb]{1,0,0}{1,1,1}{$\vphantom{\vbox to 4pt{}}\smash{f_{start}}$}, \fcolorbox[rgb]{1,0,0}{1,1,1}{$\vphantom{\vbox to 4pt{}}\smash{f_{end}}$}\},\fcolorbox[rgb]{1,0,0}{1,1,1}{$\vphantom{\vbox to 4pt{}}\smash{parinktys}$}]; \end{eqnarray*}

Signalo diskretizavimo dažnis \boldmath$f_{d0}$=" >
Dažnis nuo pradinio \boldmath$f_{start}$=" > iki \boldmath$f_{end}$=" >
\boldmath$\mathrm{FindRoot}$ komandos parinktys: ">

MSPSessionVariable[fPerdavimoCharakteristika,Null]; If[fPerdavimoCharakteristika=!=Null, MSPBlock[{$$r1,$$r2,$$fd,$$parinktys}, rts=FindRoot[Abs[fPerdavimoCharakteristika /. Subscript[f, d] -> $$fd ] == 1, {f, $$r1, $$r2},Evaluate[Sequence@@$$parinktys]]; MSPFormat[rts,StandardForm] ]]
]]>
20070830101200
321 submitas Vykdyti
ir silpnina, jei signalo dažnis yra aukštesnis. Praktikoje projektuojant tiek rekursinį, tiek nerekursinį skaitinį filtrą, pradžioje užduodama pageidaujama filtro dažninė charakteristika ir tik po to ieškomi tinkami filtro koeficientai, t.y. filtro svorio funkcija ir lango dydis. Rekursiniais filtrais galima suprojektuoti rezonansinius bei juostinius filtrus, kurie slopina arba praleidžia signalą tik tam tikroje dažnių juostoje. Dėl grįžtamojo ryšio rekursiniai filtrai kartais gali būti nestabilūs, todėl juos būtina projektuoti labai rūpestingai. Antra vertus, nerekursiniai filtrai visada yra stabilūs, tačiau jie reikalauja daugiau matematinių veiksmų, o atliekamų transformacijų klasė siauresnė.

Literatūra

R. Krivickas, "Skaitinis signalų apdorojimas", Mokslas, Vilnius (1984)

G. Doetsch "Anleitung zum Praktischen Gebrauch der Laplace-Transformation und der Z-Transformation", Oldenbourg, 1967. Yra rusiškas vertimas: G. Dёч, Руководство к практическому применению преобразования Лапласа и \twcal Z-преобразования, Наука, Москва, 1971

J.F. Kaiser,W.A. Reed, "Data smoothing using low-pass digital filters", Rev.Sci Instr.V.48, No.11, 1447-1455 (1977)

S. Wolfram, The Mathematica Book, 4-th ed., Wolfram Media/Cambridge University Press, 1999

L. R. Rabiner, B. Gold, "Theory and Application of Digital Signal Processing", Prentice-Hall, New Jersey (1975)