MOKSLASplius.lt
Į pradžią » virtualūs eksperimentai » Optika ir bangos

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.

Nustatymai

Pirmiausia nustatome atsakymų šrifto dydį šiame puslapyje taip, kad gerai matytume atsakymus
Šrifto dydis:

  

\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,\quadn=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.

IKI CIA

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 $.

\boldmath$\mathrm{ListPlot}[\mathrm{Table}[\{n,$\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:

  

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:

\boldmath\begin{eqnarray*}&&\mathrm{Sum}[z^{-n},\{n,0,\infty\}]\end{eqnarray*}

Jei diskretinė funkcija $ x(n) $ turi gęstančios eksponentės pavidalą, pavyzdžiui, $ \mathrm{e}^{-n/5} $,

\boldmath$\mathrm{ListPlot}[\mathrm{Table}[\{n,$\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:

  

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) $:

\boldmath\begin{eqnarray*}&&\mathrm{Sum}[\ee^{-\frac15 n} z^{-n},\{n,0,\infty\}]\end{eqnarray*}

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:

\boldmath\begin{eqnarray*}&&\mathrm{ZTransform}[\mathrm{Exp}[-b n] ,n,z]\end{eqnarray*}

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) $:

\boldmath$\mathrm{ListPlot}[\mathrm{Table}[\{n,$\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:

  

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

\boldmath\begin{eqnarray*}&&zCos=\mathrm{ZTransform}[\mathrm{Cos}[-n*\omega*T] ,n,z]//\mathrm{Simplify}\end{eqnarray*}

Atvirkštinę \twcal Z transformaciją atlieka komanda \boldmath$\mathrm{InverseZTransform}[X(z),z,n], kuria iš gauto \twcal Z vaizdo atstatysime pirminę diskretinę funkciją:

\boldmath\begin{eqnarray*}&&\mathrm{InverseZTransform}[zCos,n,z]//\mathrm{ExpToTrig}//\mathrm{Simplify}\end{eqnarray*}

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ę $ f<f_d/2 $. Priešingu atveju spektras persiklos su aukštadažnumine $ (f>f_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ą:
\boldmath\begin{eqnarray*}&&y[n\_]:=(x[n-1]+x[n-2]+x[n-3]+x[n-4]+x[n-5])/5\end{eqnarray*}

Po \twcal Z transformacijos turime:

\boldmath\begin{eqnarray*}&&\mathrm{Factor}[(\mathrm{ZTransform}[y[n],n,z])/.\,x[n\_/;n<0]\to 0]\end{eqnarray*}

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} $.

\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*}

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} $.

\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:

  

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 $):

\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}$=

  

Svorio funkciją normuokime į vienetą.

\boldmath\begin{eqnarray*}&&svorioFunkcija= svorioFunkcijaNenormuota\big/\mathrm{Apply}[\mathrm{Plus},svorioFunkcijaNenormuota];\end{eqnarray*}

Filtras su normuota svorio funkcija nekeičia nuostoviosios (nuo laiko nepriklausančios) signalo dalies, nes visų jo koeficientų suma lygi vienetui.

\boldmath\begin{eqnarray*}&&\mathrm{Apply}[\mathrm{Plus},svorioFunkcija]\end{eqnarray*}

Atvaizduokime svorio funkcijos koeficientų priklausomybę nuo jų numerio:

\boldmath\begin{eqnarray*}&&\mathrm{ListPlot}[svorioFunkcija,\mathrm{AxesLabel}\to\{"n",""\},\mathrm{PlotLabel}\to"Svorio funkcija"];\end{eqnarray*}

Piešimo parinktys:

  

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ą.

\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*}

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

\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:

  

Š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.

\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:

  

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į.

\boldmath\begin{eqnarray*}&&\{\mathrm{Length}[svorioFunkcija],\mathrm{Length}[sumSignalas]\}\end{eqnarray*}

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]$.

\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:

  

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ą.

\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:

  

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:

\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:

  

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ą.

\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:

  

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ą.

\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:

  

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.

\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:
Antrojo grafiko piešimo parinktys:
Grafikų grupės piešimo parinktys:

  

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

\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*}

Išsprendę atsakymą $ Y(z) $ atžvilgiu ir padalinę jį iš $ X(z) $, randame perdavimo charakteristiką

\boldmath\begin{eqnarray*}&&yz=\mathrm{Solve}[zAtvaizdas,Y[z]];\\&&zPerdavimoCharakteristika=\mathrm{Simplify}[\frac{Y[ z]/.\mathrm{First}[yz]}{X[z]}]\end{eqnarray*}

Šią \twcal Z perdavimo charakteristiką atitinka dažninė:

\boldmath\begin{eqnarray*}&&fPerdavimoCharakteristika= zPerdavimoCharakteristika/.z\to \mathrm{Exp}[\frac{2 \pi \mathrm{i} f}{f_d}]\end{eqnarray*}

Atvaizduokime jos amplitudės priklausomybę nuo dažnio $ f=\omega/2 \pi $, kai $ f_d=1 $.

\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:

  

Tokiu būdu nagrinėjamas rekursinis filtras stiprina harmoninį signalą, kai jo dažnis $ f/f_d<0{,}2902 $:

\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:

  

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)

‹ Fourier transformacija ir Fourier spektroskopijaBangų interferencija ›viršun
spausdinti