Uvod U Spektralnu Analizu PDF

Uvod u spektralnu analizu Predavanja, III. godina preddiplomskog studija Geofizike prof. dr. sc. Marijan Herak, prof. dr. sc. Dragutin Skoko Geofizički odsjek PMF-a, 28. studenoga 2022.  M. Herak i D. Skoko: Uvod u spektralnu analizu, skripta Uvod Ovaj kratki kurs osno...

Uvod u spektralnu analizu Predavanja, III. godina preddiplomskog studija Geofizike prof. dr. sc. Marijan Herak, prof. dr. sc. Dragutin Skoko Geofizički odsjek PMF-a, 28. studenoga 2022.  M. Herak i D. Skoko: Uvod u spektralnu analizu, skripta Uvod Ovaj kratki kurs osnova spektralne analize namijenjen je studentima geofizike. Nastao je na temelju predavanja prof. dr. sc. Dragutina Skoke koja je naknadno tijekom vremena dopunio prof. dr. sc. Marijan Herak. Predavanja su uglavnom bazirana na knjigama: Bendat, J. S., Piersol, A. G. (2010): Random data. Analysis and Measurement procedures (Fourth edition), John Wiley & Sons, Inc., Hoboken, New Jersey, pp. 604. Bracewell, R. N. (2000): The Fourier Transform and its Applications (3rd edition), McGraw-Hill, pp. 616. Brigham, E. O. (1974): The Fast Fourier Transform, Prentice-Hall, pp. 252. Havskov, J., Alguacil, G. (2004): Instrumentation in Earthquake Seismology, Springer, pp. 358. Osgood, B. G. (2019): Lectures on the Fourier Transform and its Applications, American Mathematical Society, Providence, Rhode Island, pp. 693. Papoulis, A. (1961): The Fourier integral and its applications, McGraw-Hill, New York, pp. 318. Scherbaum, F. (1996): Of Poles and Zeros, Springer, pp. 268. Literatura o spektralnoj analizi je izuzetno bogata, pa ovaj izbor nije nikako reprezentativan. Ohrabrujem studente da konzultiraju i druge knjige, kao i bogate resurse na internetu. Potreba za takav uvodni kolegij proizlazi iz same prirode većine geofizičkih podataka – to su najčešće vremenski nizovi, ponekad i vrlo dugoga trajanja, pa je njihov prikaz u prostoru vremena često nepregledan, a mnoge važne informacije koje oni nose nije lagano uočiti. Zbog toga su alati harmoničke i spektralne analize geofizičarima osobito važni – oni omogućavaju analizirane pojave promotriti i analizirati u području frekvencija, čime se postiže kompaktnost prikaza pojave bez gubitka informacija, te se postiže komplementaran uvid u njihova svojstva. Iako ćemo ovdje podatke shvaćati isključivo kao funkcije vremena, ista teorija i formalizam vrijedit će i ako vrijeme zamijenimo bilo kojom drugom varijablom (npr. prostorom). 2  M. Herak i D. Skoko: Uvod u spektralnu analizu, skripta Osnovni opis empirijskih podataka Empirijski podaci koji predstavljaju određenu fizikalnu pojavu mogu se općenito klasificirati u dvije skupine – determinističke i nedeterminističke. Deterministički podaci su oni koje možemo opisati eksplicitnim matematičkim izrazom. Zamislimo, na primjer, kruto tijelo mase m obješeno na zavojnici kao na slici. Neka je k konstanta elastičnosti zavojnice. Pretpostavimo da je tijelo izmaknuto iz položaja ravnoteže u udaljenost X, a zatim ga pustimo, te od tog trenutka počnemo odbrojavati vrijeme. Na temelju osnovnih zakona mehanike ili ponavljanja mjerenja može se utvrditi da postoji ovaj odnos položaja tijela x(t) i vremenskog trenutka t: 𝑘 𝑥(𝑡) = 𝑋 cos 𝑡, 𝑡≥0 𝑚 Tom jednadžbom određen je egzaktno položaj toga tijela u bilo koje vrijeme 𝑡 ≥ 0. Zato kažemo da su podaci kojima je opisano gibanje te mase deterministički. Gornji izraz, naravno, vrijedi ako uklonimo sve ostale sile osim sile elastičnosti zavojnice. Ima fizikalnih pojava u prirodi koje su opisane podacima za koje možemo smatrati da se mogu sa zadovoljavajućom točnosti izraziti matematičkom jednadžbom. Primjerice to su gibanje satelita po svojoj stazi, nastup pomrčine Sunca ili Mjeseca, i sl. Međutim, postoji niz prirodnih pojava koje se opisuju podacima koji nisu deterministički. Na primjer, nije moguće egzaktno odrediti visinu valova uzburkanog mora na određenom mjestu za zadani trenutak u budućnosti. Ti su podaci u svojoj osnovi slučajni, te se zato opisuju statističkim veličinama i parametrima vjerojatnosti, a ne eksplicitnim izrazima. Ovakva klasifikacija fizikalnih podataka na determinističke i nedetermi- nističke (stohastičke) može biti diskutabilna s više aspekata. Može se, npr., reći da ne postoje podaci koji bi u svim uvjetima mjerenja bili striktno 3  M. Herak i D. Skoko: Uvod u spektralnu analizu, skripta deterministički jer je uvijek moguće da događaj koji nismo uzeli u obzir promijeni tok pojave. Neka, npr. kod spomenute zavojnice ona pukne, ili se masa odvoji od zavojnice, ili sile koje nismo uzeli u obzir (npr. trenje) igraju preveliku ulogu... S druge, pak, strane, možemo ustvrditi da u prirodi ne mogu postojati nedeterminističke pojave jer bi ih se sve moglo opisati egzaktnim izrazima ako bismo dovoljno znali o njihovom osnovnom mehanizmu i rubnim uvjetima. U praksi ćemo zato lučiti determinističke podatke od nedeterminističkih prema tome možemo li ih reproducirati određenim eksperimentom unutar granice eksperimentalne pogreške (odnosno, kod prirodnih pojava, unutar granica koje zadovoljavaju naš interes). Ako takvo što nije moguće, podatke ćemo smatrati slučajnim. Deterministički podaci Podaci koji predstavljaju determinističku pojavu mogu se klasificirati prema ovoj shemi: Naravno, svaka kombinacija tih oblika determinističkih (pa i nedetermini- stičkih) podataka može doći u praksi. Radi jasnoće, opisat će se primjerom svaki od tih oblika. 4  M. Herak i D. Skoko: Uvod u spektralnu analizu, skripta Sinusoidalni deterministički podaci Sinusoidalni (ili kosinusoidalni) podaci su takvi tip periodičkih podataka koji se matematički mogu izraziti funkcijom oblika: 𝑥(𝑡) = 𝑋 cos(2𝜋𝑓 𝑡 − 𝜃), gdje je X – amplituda, f0 – frekvencija (broj oscilacija u jedinici vremena), – početna faza, x(t) – vrijednost funkcije x u vrijeme t. Kako predočiti takve podatke? Možemo to napraviti crtanjem vrijednosti x(t) u ovisnosti o vremenu t za odabrani interval vremena: Takva slika zove se vremenska slika pojave. Ukoliko ona traje vrlo dugo, ta slika postaje nepregledna. Istu informaciju možemo dobiti i tzv. spektralnom slikom signala, kod koje je nezavisna varijabla jednaka frekvenciji (s dimenzijom s–1 = Hz): 5  M. Herak i D. Skoko: Uvod u spektralnu analizu, skripta To je slika koja se sastoji od dva dijela koji pokazuju ovisnost amplitude X i faznog kuta (faze)  o frekvenciji f0. Prvi dio slike zove se spektar amplitude, a drugi dio spektar faze. Primijetimo da ovakav prikaz pojave pruža o njoj punu informaciju. Umjesto amplitude i faze često se spektralna slika predočuje amplitudama sinusnog i kosinusnog člana s nultim pomakom u fazi. U našem slučaju imamo: 𝑥(𝑡) = 𝑋 cos(2𝜋𝑓 𝑡 − 𝜃) = 𝑋 cos 𝜃 cos 2𝜋𝑓 𝑡 + 𝑋 sin 𝜃 sin 2𝜋𝑓 𝑡 𝑥(𝑡) = 𝑎 cos 2𝜋𝑓 𝑡 + 𝑏 sin 2𝜋𝑓 𝑡 U odnosu na spektralnu sliku podataka uočavamo: Spektralna slika sinusoidalnih periodičkih podataka sastoji se, dakle, od dva dijela i to:  od spektara amplitude i faznog kuta, ili  od spektara kosinusove amplitude i sinusove amplitude. 6  M. Herak i D. Skoko: Uvod u spektralnu analizu, skripta Ta dva načina prikaza spektralne slike (spektra) posve su ekvivalentna jer znajući jednu možemo izračunati drugu. Naime: 𝑎 = 𝑋 cos 𝜃, 𝑏 = 𝑋 sin 𝜃, pa je: tg 𝜃 = , 𝑋 = √𝑎 + 𝑏. Primijetimo, ipak, da konverzija [a, b] [𝑋, 𝜃] nije jednoznačna jer je = arctg(a/b)  n. Kako računalni programi u pravilu kao rezultat daju kosinus- i sinus-spektar, a najčešće je potrebno poznavati spektre amplitude i faze, posebnim algoritmima se dodaje n kada je potrebno da bi se fazni kut s frekvencijom mijenjao kontinuirano, bez skokova. U ovom, najjednostavnijem slučaju, svaka se spektralna slika (amplitude, faze, kosinus- ili sinus-spektar) sastoji od jedne komponente pridružene frekvenciji. To je najjednostavniji slučaj diskretnog ili linijskog spektra. Za razliku od vremenske slike podataka, koju nazivamo i dinamičkom slikom, spektralnu sliku nazivamo i statičnom slikom podataka. 7  M. Herak i D. Skoko: Uvod u spektralnu analizu, skripta Kompleksno-periodične pojave (ili samo periodičke pojave) Ovdje ćemo bez izvoda ustvrditi (a na ispitu ću pozitivno valorizirati dokaz, pronađite u knjigama i na webu!): Svaka periodična funkcija x(t) s osnovnim periodom Tp može se razviti u Fourierov red: 𝑥(𝑡) = 𝑎 + 𝑎 cos 𝜔 𝑡 + 𝑎 cos 2𝜔 𝑡 + ⋯ + 𝑎 cos 𝑛𝜔 𝑡 + +𝑏 sin 𝜔 𝑡 + 𝑏 sin 2𝜔 𝑡 + ⋯ + 𝑏 sin 𝑛𝜔 𝑡 + ⋯ Ovdje je 𝜔 = 2𝜋𝑓 = osnovna kružna frekvencija, f1 = 1/Tp je osnovna frekvencija, a an i bn su Fourierovi koeficijenti: 2 𝑎 = 𝑥(𝑡) cos 2𝜋𝑛𝑓 𝑡 𝑑𝑡, 𝑛 = 0, 1, 2, … 𝑇 2 𝑏 = 𝑥(𝑡) sin 2𝜋𝑛𝑓 𝑡 𝑑𝑡, 𝑛 = 1, 2, … 𝑇 Integrira se po jednom osnovnom periodu pa granice mogu biti i npr. . Prema tome, spektralna slika takvih podataka ima beskonačno mnogo komponenata od kojih je svaka pridružena frekvenciji koja je cjelobrojni višekratnik osnovne frekvencije f1, tj. 𝑓 = 𝑛𝑓. 8  M. Herak i D. Skoko: Uvod u spektralnu analizu, skripta Koeficijent an za n = 0 (f = 0 Hz) je 𝑎 = ∫ 𝑥(𝑡) 𝑑𝑡. Dakle on je definiran srednjom vrijednošću od 𝑥(𝑡). Po analogiji sa sinusoidalnim podacima, i ovdje možemo spektralnu sliku predočiti u obliku spektra amplitude i spektra faze, jer prethodni izraz za razvoj u Fourierov red možemo kompaktno napisati kao 𝑥(𝑡) = 𝑋 + 𝑋 cos(2𝜋𝑛𝑓 𝑡 − 𝜃 ). Ovdje je: 𝑎 𝑋 = , 2 𝑋 = (𝑎 + 𝑏 ) 𝑎 = 𝑋 cos 𝜃 𝑏 𝑛 = 1, 2, 3, …. 𝑡𝑔 𝜃 = 𝑏 = 𝑋 sin 𝜃 𝑎 Članovi 𝑋 cos(2𝜋𝑛𝑓 𝑡 − 𝜃 ) nazivaju se harmonici ili harmonijske komponente, a analiza razvojem periodičke funkcije u Fourierov red zove se harmonijska analiza. Kao i kod sinusoidalnih periodičkih podataka i ovi su spektri diskretni, linijski. Ovo je karakteristika spektara svih u vremenu periodičnih funkcija. Kasnije ćemo vidjeti da vrijedi i obrat – diskretni podaci u vremenu imat će periodične spektre. Dakle periodičnost u jednoj domeni implicira diskretnost u drugoj. 9  M. Herak i D. Skoko: Uvod u spektralnu analizu, skripta Zadatak Razvij u Fourierov red funkciju zadanu grafom (puna crvena linija), uz pretpostavku da se radi o periodičnoj, neparnoj funkciji! Kako je funkcija neparna [𝑓(−𝑡) = −𝑓(𝑡)], jasno je da dio koji odgovara negativnom vremenu mora biti kao što je naznačeno crtkanom linijom. Također s grafa vidimo da je period funkcije jednak 2𝜋, te se ona ponavlja na lijevo i na desno. Funkciju možemo algebarski izraziti kao: 𝜋 𝑓(𝑡) = 𝑡, 0≤𝑡≤ 7 1 𝜋 𝜋. 𝑓(𝑡) = − 𝑡 + , ≤𝑡≤𝜋 6 6 7 Najprije pokažimo da za neparne funkcije vrijedi: 4 𝑏 = 𝑓(𝑡) sin 𝑛𝜔 𝑡 𝑑𝑡 𝑇 Dokaz: 𝑓(𝑡) sin 𝑛𝜔 𝑡 𝑑𝑡 = {𝑡 → −𝑡} = 𝑓(−𝑡) sin 𝑛𝜔 (−𝑡) (−𝑑𝑡) = 𝑓(−𝑡) sin 𝑛𝜔 (𝑡) 𝑑𝑡 = − 𝑓(𝑡) sin 𝑛𝜔 (𝑡) 𝑑𝑡 = 𝑓(𝑡) sin 𝑛 𝜔 (𝑡) 𝑑𝑡. 10  M. Herak i D. Skoko: Uvod u spektralnu analizu, skripta Dakle: 2 𝑏 = 𝑓(𝑡) sin 𝑛𝜔 𝑡 𝑑𝑡 𝑇 2 2 4 𝑏 = 𝑓(𝑡) sin 𝑛𝜔 𝑡 𝑑𝑡 + 𝑓(𝑡) sin 𝑛𝜔 𝑡 𝑑𝑡 = 𝑓(𝑡) sin 𝑛𝜔 𝑡 𝑑𝑡. 𝑇 𝑇 𝑇 Lako je također dokazati da za neparne funkcije vrijedi i 𝑎 = 0, ∀𝑛. Dakle, parne funkcije sastoje se samo od parnih harmonika (dokažite!), a neparne od neparnih harmonika. U našem zadatku imamo da je 𝑇 = 2𝜋, pa je 𝑓 = : 4 2 𝑏 = 𝑓(𝑡) sin 2𝜋𝑛𝑓 𝑡 𝑑𝑡 = 𝑓(𝑡) sin 𝑛𝑡 𝑑𝑡. 2𝜋 𝜋 Da bismo riješili zadatak valja dakle izračunati samo 𝑏 prema gornjoj formuli, jer su svi 𝑎 = 0. Integral ćemo izračunati za svaki od dva dijela funkcije 𝑓(𝑡): 2 2 𝜋 1 𝑏 = 𝑡 sin 𝑛𝑡 𝑑𝑡 + − 𝑡 sin 𝑛𝑡 𝑑𝑡 = 𝐼 + 𝐼 𝜋 𝜋 6 6 𝜋 𝜋 𝜋 2 2 sin 𝑛𝑡 𝑡 cos 𝑛𝑡 2 sin 𝑛 7 7 cos 𝑛 7 𝐼 = 𝑡 sin 𝑛𝑡 𝑑𝑡 = − = − 𝜋 𝜋 𝑛 𝑛 𝜋 𝑛 𝑛 2 𝜋 1 1 1 𝐼 = − 𝑡 sin 𝑛𝑡 𝑑𝑡 = − 𝑡 sin 𝑛𝑡 𝑑𝑡 + sin 𝑛𝑡 𝑑𝑡 = 𝜋 6 6 3𝜋 3 11  M. Herak i D. Skoko: Uvod u spektralnu analizu, skripta 𝐼 =⋯ 𝑛𝜋 𝜋 𝑛𝜋 1 sin 𝑛𝜋 𝜋 cos 𝑛𝜋 sin 7 cos 1 1 1 𝑛𝜋 =− − − +7 7 + − cos 𝑛𝜋 + cos 3𝜋 𝑛 𝑛 𝑛 𝑛 3 𝑛 𝑛 7 𝟕 𝒏𝝅 𝐼 + 𝐼 = ⋯ = 𝒃𝒏 = 𝐬𝐢𝐧 𝟑𝒏𝟐 𝝅 𝟕 (Provedite sami račun tamo gdje je = ⋯ = !) Uočimo da amplituda koeficijenata 𝑏 opada s kvadratom broja n, a 𝑏 = 0 za n = 7N, N = 0, 1, 2, 3... Dakle i b7 = 0. Naš zadatak zapravo opisuje spektar zvuka koji proizvodi žica (npr. na gitari) kada ju se trzne na 1/7 njezine duljine. Glazbenici kažu da suzvuk osnovnog harmonika (n = 0) i sedmog harmonika (n = 7) zvuči disonantno, neugodno. Kako trzanje na 1/7 duljine žice ne proizvodi sedmi harmonik (b7 = 0), harmonička analiza kaže da se tako dobije najugodniji zvuk. S druge strane, amplituda 7. harmonika je 49 puta manja od amplitude osnovnog harmonika, pa treba zbilja dobro uho da ga se čuje! Ako imate gitaru (i izvrstan sluh!), pokušajte trzati žicu na raznim mjestima i procijenite na kojem mjestu će zvuk biti najčišći. Izaberite jednu od debljih žica da bi sedmi harmonik (7 puta veća frekvencija od osnovne frekvencije te žice!) uopće bio u Vašem slušnom području! Kompleksni oblik Fourierovog reda U spektralnoj analizi uobičajeno je korištenje kompleksnih brojeva kako bi se formule napisale u kompaktnoj formi te se tako skratili dugački izrazi. To je donekle kontraintuitivno jer se realni problemi opisuju kompleksnim brojevima, pa imaginarni dio moramo na neki način 'izmisliti', a formalno 12  M. Herak i D. Skoko: Uvod u spektralnu analizu, skripta ćemo uvesti i negativne frekvencije... To će se na kraju kompenzirati time da će se gledati samo realni dio rješenja. I Fourierov red možemo zapisati u kompleksnom obliku. Krenimo od definicije: 𝑎 𝑦(𝑡) = + [𝑎 cos 𝑛𝜔 𝑡 + 𝑏 sin 𝑛𝜔 𝑡] (1) 2 2 𝑎 = 𝑦(𝑡) cos 𝑛𝜔 𝑡 𝑑𝑡 (2) 𝑇 2 𝑏 = 𝑦(𝑡) sin 𝑛𝜔 𝑡 𝑑𝑡 (3) 𝑇 Kako je cos 𝑥 = 𝑒 +𝑒 i sin 𝑥 = 𝑒 −𝑒 = 𝑒 −𝑒 , to (1) postaje: 𝑎 1 𝑖 𝑦(𝑡) = + 𝑎 (𝑒 +𝑒 )+ 𝑏 (𝑒 −𝑒 ) 2 2 2 𝑎 1 1 = + (𝑎 − 𝑖𝑏 )𝑒 + (𝑎 + 𝑖𝑏 )𝑒 (4) 2 2 2 Kako bismo ovo pojednostavili, uvedimo negativne indekse, n < 0, u (2) i (3). Lako je vidjeti da vrijedi: 𝑎 =𝑎 𝑛 = 0, 1, 2, … i 𝑏 = −𝑏 𝑛 = 1, 2, … Provedimo zamjenu n –n kod posljednjeg člana u (4), pa dobijemo: 𝑎 1 1 𝑦(𝑡) = + (𝑎 − 𝑖𝑏 )𝑒 + (𝑎 − 𝑖𝑏 )𝑒 2 2 2 1 𝑦(𝑡) = (𝑎 − 𝑖𝑏 )𝑒 = 𝛼 𝑒. 2 𝛼 su kompleksni Fourierovi koeficijenti, 𝛼 = (𝑎 − 𝑖𝑏 ), 𝑛 = 0, ±1, ±2, … 13  M. Herak i D. Skoko: Uvod u spektralnu analizu, skripta Gotovo-periodične pojave Vidjeli smo da se periodični podaci mogu izraziti nizom sinusoidalnih valova čije su frekvencije cjelobrojni višekratnici osnovne frekvencije. Ponekad se podaci mogu izraziti sumom (ko)sinusoida, ali oni ipak neće biti periodični. Npr: 𝑥(𝑡) = 𝑋 cos(2𝜋2𝑡 + 1.3) + 𝑋 cos(2𝜋3𝑡 − 0.5) + 𝑋 cos 2𝜋√50 𝑡 + 3.0. Pripadni spektar amplituda izgleda ovako: Podaci u ovom slučaju nisu periodični, jer ne postoji frekvencija konačne vrijednosti koja bi pomnožena adekvatnim cijelim brojevima dala frekvencije 2, 3, √50. Vremensku sliku takvih podataka nazivamo 'gotovo-periodičkom'. Prema tome, ta klasa neperiodičkih podataka može se prikazati sumom koja podsjeća na Fourierov red 𝑥(𝑡) = 𝑋 + 𝑋 cos(2𝜋𝑓 𝑡 − 𝜃 ), gdje nije racionalno za sve slučajeve. Tranzijentne (neperiodične) pojave Tranzijentni su podaci oni koji imaju definiran početak i kraj. Često je početak jasan, a kraj nije, ali se može naći takvo vrijeme da se smije reći kako je pojava utrnula, te je nakon toga vremena više nema. To su, dakle, signali 14  M. Herak i D. Skoko: Uvod u spektralnu analizu, skripta ograničenoga trajanja, pa analiza razvojem u Fourierov red nije moguća. Kako ne postoji konačan osnovni period, formalno možemo uzeti da u ovom slučaju 𝑇 → ∞. Tipičan primjer iz geofizike je seizmogram potresa koji ima jasan početak te trne u vremenu. Npr., ovo je vertikalna komponenta seizmograma potresa koji je 22.03.2020. pogodio Zagreb (M = 5.5), zapisana na seizmološkoj postaji Dubrovnik: Vidi se jasan početak, te postupno trnjenje energije kako vrijeme protječe (minute su označene sivim vertikalnim linijama). Evo još tri primjera takvih signala: Fizikalne pojave predočene takvim podacima su brojne. Na primjer, druga slika prikazuje slobodno gibanje seizmografa uz mirenje. Važno svojstvo takvih podataka, nasuprot periodičkih i „gotovo- periodičkih“ jest činjenica da u njih diskretan spektar nije moguć. U većini, pak slučajeva spektralna slika može se dobiti Fourierovim integralom: 15  M. Herak i D. Skoko: Uvod u spektralnu analizu, skripta 𝑋(𝑓) = 𝑥(𝑡)𝑒 𝑑𝑡. Spektar 𝑋(𝑓) tako defniran općenito je kompleksan broj koji se, kao i svaki drugi kompleksan broj, može izraziti s: 𝑋(𝑓) = |𝑋(𝑓)|𝑒 ( ) , gdje je |𝑋(𝑓)| – apsolutna vrijednost (magnituda) od 𝑥(𝑡) 𝜃(𝑓) – argument. Fourierov spektar magnitude (amplitude) |𝑋(𝑓)| za tri prikazana primjera izgleda ovako: Izvedimo, dakle, Fourierov integral! Fourierov integral Do sada je interval u kojem smo trebali predočiti podatke 𝑥(𝑡) Fourierovim redom bio konačne duljine, npr. 𝑇. Razumljivo je stoga, ukoliko odaberemo taj interval i izračunamo Fourierov red, tada će on prikazivati naše podatke unutar toga intervala, dok će se izvan intervala slika periodički ponavljati u beskonačnost, bez obzira kakvi su podaci 𝑥(𝑡) izvan izvornog intervala duljine 𝑇. Izrazi 𝑎 = ∫ 𝑥(𝑡) cos 2𝜋𝑛𝑓 𝑡 𝑑𝑡, 𝑛 = 0, 1, 2, … 𝑏 = ∫ 𝑥(𝑡) sin 2𝜋𝑛𝑓 𝑡 𝑑𝑡, 𝑛 = 1, 2, … ukazuju da članovi 𝑎 i 𝑏 reda imaju u nazivniku 𝑇 , te ukoliko bismo htjeli te izraze primijeniti i na funkcije čiji se interval zadanosti proteže u 16  M. Herak i D. Skoko: Uvod u spektralnu analizu, skripta beskonačnost (neperiodičke funkcije), bilo bi potrebno tražiti graničnu vrijednost nepravog integrala duljine intervala integracije (tj. 𝑇 ), kada taj interval teži u beskonačnost. Fourier je to riješio općenito i našao jedan od načina prikaza i onih vrijednosti 𝑥(𝑡) koje su zadane za svaku realnu vrijednost od t, dakle za −∞ < 𝑡 < +∞. Pretpostavimo da je 𝑥(𝑡) zadano u intervalu [−𝑑, +𝑑], gdje je 𝑇 = 2𝑑 = , pa se podatak 𝑥(𝑡 ) u trenutku 𝑡 = 𝑡 unutar [−𝑑, +𝑑] može izraziti razvojem u Fourierov red: 𝑎 𝑥(𝑡 ) = + [𝑎 cos 2𝜋𝑛𝑓 𝑡 + 𝑏 sin 2𝜋𝑛𝑓 𝑡 ] 2 Uvrstimo 𝑎 i 𝑏 : 1 2 𝑥(𝑡 ) = 𝑥(𝑡)𝑑𝑡 + 2 2𝑑 2 + 𝑥(𝑡) cos 2𝜋𝑛𝑓 𝑡 𝑑𝑡 cos 2𝜋𝑛𝑓 𝑡 2𝑑 2 + 𝑥(𝑡) sin 2𝜋𝑛𝑓 𝑡 𝑑𝑡 sin 2𝜋𝑛𝑓 𝑡 2𝑑 (Izrazi u okruglim zagradama su 𝑎 i 𝑏.) 1 1 𝑥(𝑡 ) = 𝑥(𝑡)𝑑𝑡 + 𝑥(𝑡) cos 2𝜋𝑛𝑓 (𝑡 − 𝑡 )𝑑𝑡. 2𝑑 𝑑 a) Uzmimo sada da duljina intervala 2d raste i teži beskonačnosti. Tada veličina 𝑓 = 1/2𝑑 opada i može se učiniti po volji malenom, ako je d dovoljno velik. b) Neka postoji konačan, određen limes lim ∫ |𝑥(𝑡)|𝑑𝑡, tj. neka integral → konvergira apsolutno (nepravi integral). c) Kako prvi integral na desnoj strani konvergira apsolutno za 𝑑 → ∞ (ima konačnu vrijednost kada d teži u beskonačnost) prvi član na desnoj strani teži nuli. 17  M. Herak i D. Skoko: Uvod u spektralnu analizu, skripta d) Tada zbroj preostalih članova na desnoj strani možemo pisati u novom obliku (1/𝑑 = 2𝑓 ): 𝑥(𝑡 ) = 2𝑓 𝑥(𝑡) cos 2𝜋𝑛𝑓 (𝑡 − 𝑡 )𝑑𝑡. Budući da je 𝑛𝑓 = 𝑓 , 2𝑓 , 3𝑓 , … frekvencija koja odgovara pojedinom broju n (harmoniku), to ćemo je općenito označiti s f, tj. 𝑛𝑓 = 𝑓. Pustit ćemo nadalje da granice integracije rastu preko po volji velikog broja. Tada vrijednost osnovne frekvencije 𝑓 = 1/2𝑑 pada ispod po volji male veličine. Kako 𝑓 ujedno znači i razliku frekvencija (𝑛 + 1)-og i n-tog člana, tj. (𝑛 + 1)𝑓 − 𝑛𝑓 = 𝑓 , to ćemo tu malu veličinu označiti s ∆𝑓. Rekli smo ranije da ∫ 𝑥(𝑡)𝑑𝑡 konvergira apsolutno, pa će i ∫ 𝑥(𝑡) cos 2𝜋𝑓(𝑡 − 𝑡) 𝑑𝑡 konvergirati, te izraz za 𝑥(𝑡 ) pri graničnom prijelazu (d → ∞, ∆𝑓 → 𝑑𝑓) postaje: 𝑥(𝑡 ) = 2 𝑑𝑓 𝑥(𝑡) cos 2𝜋𝑓(𝑡 − 𝑡 )𝑑𝑡. Izraz za 𝑥(𝑡 ) možemo pisati i u obliku: 𝑥 (𝑡 ) = 2 𝑥 (𝑡 ) cos 2𝜋𝑓𝑡 cos 2𝜋𝑓𝑡 𝑑𝑡 + 𝑥 (𝑡 ) sin 2𝜋𝑓𝑡 sin 2𝜋𝑓𝑡 𝑑𝑡 𝑑𝑓 = =2 cos 2𝜋𝑓𝑡 𝑥 (𝑡 ) cos 2𝜋𝑓𝑡 𝑑𝑡 + sin 2𝜋𝑓𝑡 𝑥(𝑡 ) sin 2𝜋𝑓𝑡 𝑑𝑡 𝑑𝑓 ( )/ ( )/ Dakle, vrijednost podataka u trenutku 𝑡 može se izraziti Fourierovim integralom: 𝑥(𝑡 ) = [𝑎(𝑓) cos 2𝜋𝑓𝑡 + 𝑏(𝑓) sin 2𝜋𝑓𝑡 ] 𝑑𝑓, gdje su 𝑎(𝑓) = 2 ∫ 𝑥(𝑡) cos 2𝜋𝑓𝑡 𝑑𝑡 i 𝑏(𝑓) = 2 ∫ 𝑥(𝑡) sin 2𝜋𝑓𝑡 𝑑𝑡. 18  M. Herak i D. Skoko: Uvod u spektralnu analizu, skripta 𝑎(𝑓) i 𝑏(𝑓) zovu se kosinusova, odnosno sinusova Fourierova transformacija. Kao i kod Fourierovog reda, i ovdje definiramo spektar amplitude i spektar faze: 𝑏(𝑓) 𝑋(𝑓) = 𝑎 (𝑓) + 𝑏 (𝑓) tg 𝜃(𝑓) =. 𝑎(𝑓) Fourierov integral podsjeća na Fourierov red. Međutim, sada su koeficijenti 𝑎(𝑓) i 𝑏(𝑓), odnosno amplitudni i fazni spektri 𝑋(𝑓) i 𝜃(𝑓) kontinuirane funkcije frekvencije, a zbroj harmonika u Fourierovom redu zamijenjen je integralom protegnutim od 0 do +. Ovdje se na radi o izražavanju zadanih podataka zbrojem diskretnih harmonijskih oscilacija, nego kontinuiranim prijelazom neprekinutih harmonijskih funkcija kojima frekvencije f poprimaju realne vrijednosti od 0 do +. Analiza podataka pomoću računanja Fourierovih integrala zove se spektralna analiza (za razliku od harmonijske analize koju smo prije upoznali). Jasno je da Fourierov integral ne može biti striktno upotrijebljen, jer zahtijeva integriranje po cijeloj vremenskoj osi (0 do +!), što, pogotovo kod eksperimentalno mjerenih podataka, ne može biti zadovoljeno. Toga valja biti svjestan, pogotovo ako želimo na osnovi mjerenih podataka ekstrapolirati pojavu u budućnost. Kod Fourierovog reda i periodičkih podataka, te ako su mjerenja zahvatila barem jedan potpuni osnovni period, tako nešto može imati smisla. Kod računa Fourierovog integrala, međutim, on u realnim primjenama može biti izračunat samo uz konačnu donju (početak) i gornju granicu (kraj našega mjerenja), što implicira da pojave koju promatramo nema niti prije niti poslije mjerenja, pa će svaka ekstrapolacija izvan intervala mjerenja rezultirati nulom. Dakle, rezultati harmoničke i spektralne analize podudarat će se unutar intervala mjerenja, ali će se razlikovati izvan njega – kod harmoničke analize signal će se periodički ponavljati na jednu i na drugu stranu, dok će Fourierov integral po definiciji tamo biti jednak nuli. Drugim riječima, Fourierov red primjenjuje se samo na (periodičke) funkcije beskonačne energije (infinite energy), dok se Fourierova transformacija odnosi samo na pojave čija je energija konačna (sjetimo se: lim ∫ |𝑥(𝑡)|𝑑𝑡 ≠ ∞, pa će u praktički svim nama interesantnim slučajevima → biti i lim ∫ 𝑥 𝑑𝑡 ≠ ∞!). → 19  M. Herak i D. Skoko: Uvod u spektralnu analizu, skripta Zbog toga valja primjenjivati onu vrst analize koja je primjerena za podatke kojima raspolažemo, a spektre izračunate harmonijskom i spektralnom analizom ne možemo apriori uspoređivati. Zadatak 2: Izračunajte spektar funkcije 𝑥(𝑡) = A sin 𝑡 koja je zadana u intervalu [− ,+ ] prvo harmoničkom pa zatim i spektralnom analizom! Kako je to obični sinus, a harmonička analiza pretpostavlja periodičnost, funkcija koju ćemo razviti u red je prikazana crtkanom linijom na slici gore. To je neparna funkcija, pa – kao što smo dokazali – svi koeficijenti 𝑎 = 0. Postoji samo jedna neparna harmonička komponenta, ona s frekvencijom 𝑓 = 1/𝑇. Zato je 𝑏 = 0, ∀ 𝑛 ≠ 1; 𝑏 = 𝐴. Spektri amplitude i faze su: 𝑏 𝜋 𝑋 = 𝑏 = 𝑏, 𝜃 = atan = 𝑎 2 Spektri imaju samo jednu liniju: 20  M. Herak i D. Skoko: Uvod u spektralnu analizu, skripta Ako želimo koristiti Fourierov integral, kosinusova transformacija 𝑎(𝑓) = 0 (baš kao i kod Fourierovih koeficijenata pri harmonijskoj analizi – dokažite!), a gra- nice integracije pri računu koeficijenta 𝑏(𝑓) moramo postaviti na interval zada- nosti funkcije, ±𝑇 /2, umjesto na ±∞: / 2𝜋 𝑏(𝑓) = 2𝐴 sin 𝑡 𝑑𝑡 / 𝑇... (provedite potreban račun!)... 𝑓 2𝐴 sin 𝜋 1 − 𝑓 𝑏(𝑓) = , 𝑓 𝜋𝑓 1− 𝑓 pa spektri izgledaju ovako: Prilična razlika o odnosu na harmoničku analizu! 21  M. Herak i D. Skoko: Uvod u spektralnu analizu, skripta Kompleksan oblik Fourierovog integrala, Fourierova transformacija Fourierov integral funkcije 𝑥(𝑡 ) prikazivali smo u obliku: 𝑥(𝑡 ) = ∫ [𝑎(𝑓) cos 2𝜋𝑓𝑡 + 𝑏(𝑓) sin 2𝜋𝑓𝑡 ] 𝑑𝑓, ili 𝑥(𝑡 ) = 𝑋(𝑓) cos[2𝜋𝑓𝑡 − 𝜃(𝑓)] 𝑑𝑓, gdje su: 𝑎(𝑓) = 2 ∫ 𝑥(𝑡) cos 2𝜋𝑓𝑡 𝑑𝑡 i 𝑏(𝑓) = 2 ∫ 𝑥(𝑡) sin 2𝜋𝑓𝑡 𝑑𝑡, 𝑎(𝑓) = 𝑋(𝑓) cos 𝜃, 𝑏(𝑓) = 𝑋(𝑓) sin 𝜃 (*) 𝑏(𝑓) 𝑋(𝑓) = 𝑎 (𝑓) + 𝑏 (𝑓) , tg 𝜃(𝑓) =. 𝑎(𝑓) Odnos 𝑎(𝑓), 𝑏(𝑓), 𝑋(𝑓) i 𝜃(𝑓) može se predočiti grafički kao na slici: Za matematičko izvođenje često je, međutim, zgodno uvesti kompleksan oblik spektra, odnosno Fourierovog integrala. Kompleksan spektar 𝐶(𝑓) definira se kosinus i sinus spektrom tako da je 1 𝐶(𝑓) = [𝑎(𝑓) − 𝑖 𝑏(𝑓)] (∗∗) 2 kompleksan broj, tj. spektar prikazujemo točkom kompleksne ravnine. Vidimo da je apsolutna vrijednost, tj. modul kompleksnog spektra jednak spektru amplitude 𝑋(𝑓), a argument mu je jednak negativnoj vrijednosti spektra faze. Prema tome, uvođenje kompleksnog spektra ne uvodi ništa bitno novo u razmatranje. U praksi se primjenjuje onaj oblik spektra koji je jednostavniji za izvođenje. Uvrstimo (*) u (**), i dobijemo: 22  M. Herak i D. Skoko: Uvod u spektralnu analizu, skripta ( ) 2𝐶(𝑓) = 𝑋(𝑓)[cos 𝜃(𝑓) − 𝑖 sin 𝜃(𝑓)] = 𝑋(𝑓)𝑒 Uvrstimo li izraze a 𝑎(𝑓) i 𝑏(𝑓) u 𝐶(𝑓) (**) izlazi: 2 𝐶(𝑓) = 𝑎(𝑓) − 𝑖 𝑏(𝑓) = 2 𝑥(𝑡) cos 2𝜋𝑓𝑡 − 2𝑖 𝑥(𝑡) sin 2𝜋𝑓𝑡 𝐶(𝑓) = 𝑥(𝑡)[cos 2𝜋𝑓𝑡 − 𝑖 sin 2𝜋𝑓𝑡]𝑑𝑡 𝐶(𝑓) = 𝑥(𝑡) 𝑒 𝑑𝑡 = 𝑥(𝑡) 𝑒 𝑑𝑡. Ovo je Fourierova transformacija, što je osnovna formula za određivanje spektra u kompleksnom obliku. Fourierov integral 𝑥(𝑡 ) = [𝑎(𝑓) cos 2𝜋𝑓𝑡 + 𝑏(𝑓) sin 2𝜋𝑓𝑡 ] 𝑑𝑓 prilagodit ćemo kompleksnom obliku prikaza spektra amplitude i faze tako da kosinuse i sinuse izrazimo Moivreovim formulama: 1 −𝑖 𝑥(𝑡 ) = 𝑎(𝑓) 𝑒 + 𝑒 + 𝑏(𝑓) 𝑒 − 𝑒 𝑑𝑓 = 2 2 1 = 𝑒 [𝑎(𝑓) − 𝑖𝑏(𝑓)] + 𝑒 [𝑎(𝑓) + 𝑖𝑏(𝑓)] 𝑑𝑓 = 2 1 1 = 2𝐶(𝑓) 𝑒 𝑑𝑓 + [𝑎(𝑓) + 𝑖𝑏(𝑓)]𝑒 𝑑𝑓 2 2 𝑓 → −𝑓 = DRUGI INTEGRAL: 𝑑𝑓 → −𝑑𝑓 = granice → −granice 1 1 = 2𝐶(𝑓) 𝑒 𝑑𝑓 + [𝑎(−𝑓) + 𝑖𝑏(−𝑓)]𝑒 (−𝑑𝑓) = 2 2 1 1 = 2𝐶(𝑓) 𝑒 𝑑𝑓 + [𝑎(𝑓) − 𝑖𝑏(𝑓)] 𝑒 𝑑𝑓 = 2 2 ( ) 1 1 = 2𝐶(𝑓) 𝑒 𝑑𝑓 + 2𝐶(𝑓) 𝑒 𝑑𝑓 2 2 23  M. Herak i D. Skoko: Uvod u spektralnu analizu, skripta 𝑥(𝑡 ) = 𝐶(𝑓) 𝑒 𝑑𝑓, pri čemu za negativne frekvencije treba uzeti konjugirano kompleksnu veliči- nu od 𝐶(𝑓), tj. 𝐶 ∗ (𝑓), jer je 𝐶(−𝑓) = [𝑎(−𝑓) − 𝑖𝑏(−𝑓)] = [𝑎(𝑓) + 𝑖𝑏(𝑓)] = 𝐶 ∗ (𝑓). [a je parna, b je neparna funkcija] Ovo je inverzna Fourierova transformacija, koju možemo izraziti i koristeći kružnu frekvenciju, 𝜔, umjesto f (primijetimo da je 𝑑𝑓 = 𝑑𝜔/2𝜋): 1 𝑥(𝑡 ) = 𝐶(𝜔) 𝑒 𝑑𝜔. 2𝜋 24  M. Herak i D. Skoko: Uvod u spektralnu analizu, skripta O karakteristikama fizikalnih sustava u odnosu na frekvenciju Osvrnut ćemo se na neke karakteristike fizikalnih sustava u ovisnosti o frekvenciji. To, primjerice, može biti mjerni instrument (seizmograf, mareograf...), ili sredstvo kroz koje se valovi rasprostiru (npr. cijela Zemlja, površinski slojevi tla, morska voda, zrak,...) – sve što može promijeniti spektralni sastav pojave koju promatramo. Sustav ćemo predočiti 'crnom kutijom' – ne zanima nas što je u njoj niti kako npr. instrument radi, ali nas zanima kako se prolaskom kroz taj sustav promijenio ulazni signal 𝑥 (𝑡) čiji spektar označimo s 𝐶 (𝑓). Na slici dvostruke strelice označavaju Fourierov par, tj. proces (inverzne) Fourierove transformacije. Dakle 𝐶 (𝑓) je spektar izlaznog signala 𝑥 (𝑡). Razmotrit ćemo samo idealne sustave. Idealni je sustav onaj koji:  ima konstantne parametre  linearan je u odnosu izlaz/ulaz. Sustav ima konstantne parametre ako mu se bitna svojstva ne mijenjaju s vremenom. Do promjene svojstava može doći npr. starenjem elektronskih komponenata, promjenama elastičkih svojstava opruge, popravcima, i sl. Linearan sustav ima frekventnu karakteristiku koja je aditivna i homogena. Aditivnost znači: izlaz zbroja dvaju različitih pojava jednak je zbroju izlaza svake od njih (svejedno je zbraja li se prije ili poslije sustava). Homogenost znači: izlaz n-strukog ulaza jednak je 𝑛 × 1 izlaz (svejedno je množi li se prije ili poslije sustava). Linearnost je obično razumna pretpostavka, ali ipak instrumenti ponekad nisu linearni – njihov utjecaj na signal tada ovisi o amplitudi ulaznog signala. Primjer 25  M. Herak i D. Skoko: Uvod u spektralnu analizu, skripta je mehanički seizmograf kod kojega trenje igra veću ulogu kod male nego kod velike pobude. Sustav koji promatramo, npr. seizmograf, će – ovisno o tipu seizmografa – više ili manje modificirati ulazni signal (trešnju tla), pa njegov izlaz (seizmogram) neće vjerno predstavljati ulaznu pojavu. Tako ćemo definirati dinamičku karakteristiku sustava kao omjer spektara izlaznog i ulaznog signala: 𝐶 (𝑓) 𝑍(𝑓) =. 𝐶 (𝑓) Kod idealnih sustava Z ovisi jedino o frekvenciji i parametrima sustava, a ne ovisi niti o ulaznoj pojavi niti o vremenu. Dinamička karakteristika sustava općenito je kompleksan broj, te se može napisati kao ( ) 𝑍(𝑓) = |𝑍(𝑓)|𝑒. Apsolutnu vrijednost |𝑍(𝑓)| nazivamo dinamičko povećanje, a 𝛿(𝑓) faza [eng. gain factor; phase factor]. Poznavajući dinamičku karakteristiku sustava možemo na osnovi poznatog izlaza odrediti ulaz i obrnuto. Kako se određuje karakteristika nekog instrumenta? 1. Većina modernih instrumenata kalibrirana je u tvornici, te se njihova karakteristika tijekom uporabnog vijeka vrlo malo (zanemarivo) mijenja. U tom slučaju njihova je karakteristika navedena u dokumentaciji instrumenta, te ju je potrebno (vrlo rijetko) provjeravati kalibracijom. 2. Mjerenjem pojave na ulazu i izlazu. Ovo je ponekad lako izvedivo, a ponekad jako teško. 3. Teorijski, na temelju poznatih parametara sustava. Neka je npr. sustav seizmograf s mehaničkom registracijom. U seizmometriji ćete naučiti da je njegova dinamička karakteristika zadana parametrima: 𝑉 |𝑍(𝑓)| = 𝑓 𝑓 −1 + 4ℎ 𝑓 𝑓 26  M. Herak i D. Skoko: Uvod u spektralnu analizu, skripta 𝜖𝑓 𝛿(𝑓) =. 𝜋(𝑓 − 𝑓 ) Parametri sustava su: 𝑉 − 𝑠𝑡𝑎𝑡𝑖č𝑘𝑜 𝑝𝑜𝑣𝑒ć𝑎𝑛𝑗𝑒 ℎ − 𝑘𝑜𝑛𝑠𝑡𝑎𝑛𝑡𝑎 𝑚𝑖𝑟𝑒𝑛𝑗𝑎 𝑓 − 𝑣𝑙𝑎𝑠𝑡𝑖𝑡𝑎 𝑓𝑟𝑒𝑘𝑣𝑒𝑛𝑐𝑖𝑗𝑎 𝑠𝑒𝑖𝑧𝑚𝑜𝑔𝑟𝑎𝑓𝑎 𝜖 − 𝑓𝑎𝑘𝑡𝑜𝑟 𝑚𝑖𝑟𝑒𝑛𝑗𝑎. Dakle, želimo li iz zapisa seizmografa s poznatom karakteristikom 𝑍(𝑓) rekonstruirati pravo gibanje tla valja: 1. Izračunati spektar seizmograma: 𝐶 (𝑓) = 𝑥 (𝑓) 𝑒 𝑑𝑡; 2. Izračunati spektar ulaza (gibanje tla): 𝐶 (𝑓) 𝐶 (𝑓) = ; 𝑍(𝑓) 3. Izračunati 𝑥 (𝑡) kao Fourierov par od 𝐶 (𝑓): 𝑥 (𝑡) = 𝐶 (𝑓) 𝑒 𝑑𝑓. Ovdje valja biti oprezan u dijelu frekvencija koje instrument vrlo slabo bilježi, jer će tamo karakteristika 𝑍(𝑓) biti jako mala, pa će nazivnik u 2. koraku biti jako velik. Zato valja primjenu 2. koraka ograničiti samo na onaj frekventni pojas u kojemu je unutrašnji šum instrumenta znatno manji od razine ulaznog signala (tada je instrument linearan). 27  M. Herak i D. Skoko: Uvod u spektralnu analizu, skripta Slučajni podaci Već je rečeno da se podaci koji su slučajne vrijednosti ne mogu opisati eksplicitnim matematičkim izrazom, jer svako opažanje predstavlja jedinstven događaj između više mogućih. Isječak iz vremenske slike slučajnih podataka je uzorak. Slučajne se pojave mogu svrstati u stacionarne i nestacionarne, koje se dalje dijele kako je prikazano slikom. Ako se određene vrijednosti pridružene nekoj pojavi mogu smatrati slučajnima, tada se pojava može opisati npr. srednjacima uzoraka. Neka je na slikama predočeno nekoliko (N) uzoraka (to se naziva ensemble). Srednjak slučajnih vrijednosti u nekom trenutku 𝑡 možemo odrediti kao srednjak trenutnih vrijednosti svih uzoraka koji su nam na raspolaganju. Slično vrijedi i za korelaciju slučajnih vrijednosti u dva različita trenutka (autokorelacija). Možemo pisati: 1 𝜇 (𝑡 ) = 𝑥 (𝑡 ) 𝑁 1 𝑅 (𝑡 , 𝑡 + 𝜏) = 𝑥 (𝑡 )𝑥 (𝑡 + 𝜏), 𝑁 gdje smo pretpostavili da su težine uzoraka jednake. 28  M. Herak i D. Skoko: Uvod u spektralnu analizu, skripta Općenito 𝜇(𝑡 ) i 𝑅(𝑡 , 𝑡 + 𝜏) ovisit će o trenutku 𝑡 i razlici trenutaka 𝜏. Za takvu pojavu kažemo da je nestacionarna. U posebnom slučaju kad 𝜇(𝑡 ) i 𝑅(𝑡 , 𝑡 + 𝜏) ne ovise o trenutku 𝑡 , već samo o 𝜏, za pojavu kažemo da je slabo stacionarna ili stacionarna u širem smislu. Prema tome, tada je srednjak konstantan, a funkcija autokorelacije ovisi samo o vremenskoj razlici 𝜏, te je 𝜇 (𝑡 ) = 𝜇 , 𝑅 (𝑡 , 𝜏) = 𝑅 (𝜏). U slučaju da su svi momenti i pridruženi momenti neovisni o vremenu 𝑡 , za pojavu kažemo da je strogo stacionarna. Ergodička slučajna pojava Ranije smo rekli da se određena svojstva slučajne pojave mogu odrediti kao srednjak vrijednosti uzoraka u određenom trenutku 𝑡. Međutim, često je moguće opisati svojstva stacionarne slučajne pojave jednostavnim određiva- njem srednjaka jednog od uzoraka. Zamislimo npr. k-ti uzorak. Tada su srednjak i autokorelacija: 1 𝜇 (𝑡 ) = lim 𝑥 (𝑡) 𝑑𝑡, → 𝑇 1 𝑅 (𝑡 , 𝑡 + 𝜏) = lim 𝑥 (𝑡)𝑥 (𝑡 + 𝜏) 𝑑𝑡. → 𝑇 Ukoliko je slučajna pojava stacionarna, a 𝜇 (𝑥) i 𝑅 (𝜏, 𝑘) ne variraju kada ih računamo za različite uzorke i jednaki su 𝜇 i 𝑅 (𝑡) za stacionarnu slučajnu pojavo kažemo da je ergodička. Ergodički slučajni procesi su važna klasa slučajnih pojava jer sva njihova bitna obilježja možemo dobiti na osnovi jednog, odnosno, radi kontrole, nekoliko uzoraka. U praksi se nešto odstupa od stroge definicije stacionarnosti. Da bismo uzorak tretirali kao uzorak stacionarnih podataka, parametri (srednjak i autokorela- cija) određeni iz podataka kraćih vremenskih intervala ne bi smjeli „značajno“ odstupati od jednog intervala do drugog. U daljem se nećemo detaljno baviti svojstvima slučajnih pojava (to ćete detaljnije naučiti u drugim kolegijima), ali ćemo ipak posvetiti pažnju svojstvima koja se važna za spektralnu analizu. To su funkcije (auto)korelacije i funkcije spektra gustoće snage. 29  M. Herak i D. Skoko: Uvod u spektralnu analizu, skripta Funkcija autokorelacije Funkcija autokorelacije slučajnih vrijednosti opisuje opću ovisnost vrijednosti podataka u određenom trenutku, o vrijednostima podataka u neko drugo vrijeme. Neka je dana vremenska slika uzorka kao na slici, a vrijeme kada smo počeli prikupljati podatke uzmimo kao 𝑡 = 0. Ocjenu autokorelacije između vrijednosti 𝑥(𝑡) u vrijeme 𝑡 i 𝑥(𝑡 + 𝜏) u trenutku 𝑡 + 𝜏 možemo dobiti kao srednju vrijednost produkta 𝑥(𝑡)𝑥(𝑡 + 𝜏) u intervalu vremena T. Ta srednja vrijednost produkta teži egzaktnoj vrijed- nosti funkcije autokorelacije kada T teži u beskonačnost: 𝑅 (𝜏) = 𝑥(𝑡)𝑥(𝑡 + 𝜏) 𝑑𝑡. To je, kako vidimo, realna funkcija. Lako je vidjeti da je to i parna funkcija: 𝑢 =𝑡+𝜏 𝑥(𝑡)𝑥(𝑡 + 𝜏)𝑑𝑡 = = 𝑥(𝑢 − 𝜏)𝑥(𝑢)𝑑𝑢 𝑡 =𝑢−𝜏 𝑅 (𝜏) = 𝑅 (−𝜏), a može poprimiti pozitivne i negativne vrijednosti. Budući da se granica integracija proteže u beskonačnost može se postaviti pitanje o uvjetima konvergencije toga integrala. Ne ulazeći duboko u taj problem možemo reći da postoje ova dva pravila: 30  M. Herak i D. Skoko: Uvod u spektralnu analizu, skripta 1. Za teorijski definirane funkcije integral konvergira ako postoji konačna granična vrijednost, tj. ispunjeni su Dirichlet-ovi uvjeti. 2. Za empirijske vrijednosti treba navesti da postoje alternativne definicije funkcije autokorelacije, npr.: 1 𝑅 (𝜏) = 𝑥(𝑡)𝑥(𝑡 + 𝜏) 𝑑𝑡. 𝑇−𝜏 Razumljivo je da je u praktičnom postupku interval integracije ograničen, u našem slučaju T, te prethodni izraz znači prosječnu vrijednost produkta 𝑥(𝑡)𝑥(𝑡 − 𝜏) unutar intervala T. Uključivanje faktora 1/T služi kao neka vrst normiranja, što je posebno potrebno kada uspoređujemo podatke različitih intervala T. Konačno, nas i zanimaju relativni odnosi, ali određeni na način da veličine možemo uspoređivati. Primjer: Odredi autokorelogram funkcije 𝑥(𝑡) = 𝑋 sin 2𝜋𝑓 𝑡 ! Pretpostavimo 𝑇 = 𝑘𝑇 = 𝑘 , dakle gledamo funkciju u intervalu koji odgovara cjelobrojnom višekratniku perioda. 𝑥(𝑡)𝑥(𝑡 + 𝜏) 𝑑𝑡 = = 𝑋 sin 2𝜋𝑓 𝑡 𝑋 sin 2𝜋𝑓 (𝑡 + 𝜏) 𝑑𝑡 = 𝑋 sin 2𝜋𝑓 𝑡 sin 2𝜋𝑓 (𝑡 + 𝜏) 𝑑𝑡 1 1 ( ) ( ) =𝑋 𝑒 −𝑒 𝑒 −𝑒 𝑑𝑡 = 2𝑖 2𝑖 𝑋 ( ) ( ) =− 𝑒 +𝑒 −𝑒 −𝑒 𝑑𝑡 = 4 𝑋 = [2 cos 2𝜋𝑓 (2𝑡 + 𝜏) − 2 cos 2𝜋𝑓 𝜏]𝑑𝑡 = 4 [Prvi član jednak je nuli (integral kosinusa po punom periodu)] 𝑋 𝑋 = cos 2𝜋𝑓 𝜏 𝑑𝑡 = 𝑇 cos 2𝜋𝑓 𝜏, 2 2 pa je funkcija autokorelacije: 1 1 𝑋 𝑋 𝑅 (𝜏) = lim 𝑥(𝑡)𝑥(𝑡 − 𝜏) 𝑑𝑡 = lim 𝑇 cos 2𝜋𝑓 𝜏 = cos 2𝜋𝑓 𝜏. → 𝑇 → 𝑇 2 2 31  M. Herak i D. Skoko: Uvod u spektralnu analizu, skripta Dakle, korelogram sinusa je kosinus iste frekvencije, a amplituda je jednaka polovici kvadrata amplitude sinusa. Kako smo i očekivali, to je realna i parna funkcija. Informacija o fazi je izgubljena. Autokorelogram slučajnih podataka ima oštar vrh za 𝜏 = 0, a dalje brzo opada prema nuli. Dakle, slučajni podaci dobro su korelirani samo 'sami sa sobom', što je i očekivano. Autokorelogram zbroja sinusoidalnih podataka i slučajnog šuma jednak je zbroju pripadnih autokorelograma, kao na slici: Određivanje funkcije autokorelacije ponekad omogućava otkrivanje determinističke pojave koja inače može biti maskirana jakim šumom. To je jedna od mogućih primjena autokorelacije. Druge, manje očite, povezane su sa spektrima gustoće snage (vidi dolje!). 32  M. Herak i D. Skoko: Uvod u spektralnu analizu, skripta Funkcije međukorelacije Tek radi potpunosti, navedimo da je funkcija međukorelacije (cross-correlation) definirana slično, samo se gleda međusobna ovisnost podataka dvije funkcije, 𝑥(𝑡) i 𝑦(𝑡): 1 𝑅 (𝜏) = lim 𝑥(𝑡)𝑦(𝑡 + 𝜏) 𝑑𝑡 → 𝑇 Ona u općem slučaju nije parna funkcija. U geofizici, posebno u seizmologiji, funkcije međukorelacije služe i za precizno određivanje vremenskih razlika nastupa iste grupe valova na dvije seizmološke postaje u istoj mreži (vidi sliku dolje), ili npr. za precizno mjerenje dvoloma S-valova koji prolaze kroz anizotropno sredstvo. 33  M. Herak i D. Skoko: Uvod u spektralnu analizu, skripta Funkcije spektra gustoće energije (snage) (Power spectral density functions) Funkcija spektra gustoće energije za slučajne podatke opisuje općenito sastav podataka u odnosu na spektralnu gustoću vrijednosti njihovih kvadrata. Budući da je spektralna analiza kao matematička disciplina bila u prvom redu primjenjivana na električke pojave, to su se i mnogi izrazi u spektralnoj analizi zadržali u analogiji s elektrotehničkim nazivima. Tako je definirana i funkcija spektra gustoće energije, ili kako se često također zove – spektra gustoće snage. Uzmimo da se u krugu električne struje jakosti 𝐼(𝑡) nalazi otpor R od 1 . Tada je gubitak energije kruga u vremenu t: ∆𝑊 = 𝐼 (𝑡)(𝑅 = 1)∆𝑡, a ukupni gubitak energije tijekom vremena 0 𝑇 𝑊 = ∫ 𝐼 (𝑡) 𝑑𝑡. U jedinici vremena prosječni gubitak energije tj. prosječna snaga je 1 𝑊= 𝐼 (𝑡) 𝑑𝑡. 𝑇 To nije ništa drugo do prosječna vrijednost kvadrata jakosti električne struje. Tako se i u spektralnoj analizi zadržao izraz spektar gustoće snage (energije), bez obzira o kakvoj se funkciji 𝐼(𝑡) radi. Označimo funkciju vremena s 𝑥(𝑡) i neka je definirana u intervalu [0, 𝑇 ]. Analogno nazivu u elektrotehnici i ovdje možemo reći da je prosječna snaga 1 𝑥 (𝑡) = 𝑥 (𝑡) 𝑑𝑡, 𝑇 no odmah uočavamo da je to statistička prosječna vrijednost kvadrata naše funkcije. Označimo je 𝑥 (𝑡). Funkciju 𝑥(𝑡) možemo izraziti Fourierovim integralom: 𝑥(𝑡) = 𝑋(𝑓) cos(2𝜋𝑓𝑡 − 𝜃) 𝑑𝑓. 34  M. Herak i D. Skoko: Uvod u spektralnu analizu, skripta Ograničimo se sada samo na onaj dio spektra koji se nalazi unutar jediničnog intervala oko neke frekvencije f0, dakle unutar 𝑓 − , 𝑓 +. Sada možemo izračunati onaj dio funkcije 𝑥(𝑡) koji sadrži frekvencije samo unutar tog intervala, 𝑓 ∈ 𝑓 − , 𝑓 + , tako da primjereno promijenimo granice integracije Fourierovog integrala: / 𝑥(𝑡; 𝑓 , 1) = 𝑋(𝑓) cos(2𝜋𝑓𝑡 − 𝜃) 𝑑𝑓. / Ovdje je 𝑥(𝑡; 𝑓 , 1) 'dio' funkcije 𝑥(𝑡) koji sadrži frekvencije (f) samo unutar jediničnog (1) intervala oko frekvencije 𝑓. Pretpostavit ćemo da je taj interval frekvencija u odnosu na cijeli interval frekvencija dovoljno uzak da se niti [𝑋(𝑓), 𝜃(𝑓)] niti kosinus ne mijenjaju znatno, pa ih možemo zamijeniti njihovim vrijednostima u 𝑓 = 𝑓. 𝑥(𝑡; 𝑓 , 1) = 𝑋(𝑓 ) cos(2𝜋𝑓 𝑡 − 𝜃(𝑓 )). Prosječna vrijednost od 𝑥 (𝑡; 𝑓 , 1) je: 1 1 𝑥 (𝑡; 𝑓 , 1) = 𝑥 (𝑡; 𝑓 , 1) 𝑑𝑡 = 𝑋 (𝑓 ) cos (2𝜋𝑓 𝑡 − 𝜃(𝑓 )) 𝑑𝑡 = 𝑇 𝑇 𝑋 (𝑓 ) 𝑋 (𝑓 ) = cos 𝜙 𝑑𝑡 = {neka je 𝑇 = 𝑘/𝑓 } = (1 + cos 2𝜙) 𝑑𝑡 = 𝑇 2𝑇 35  M. Herak i D. Skoko: Uvod u spektralnu analizu, skripta 𝑋 (𝑓 ) 𝑋 (𝑓 ) {integral cos 2𝜙 po periodu = 0} = 𝑑𝑡 =. 2𝑇 2 Ovu veličinu zvat ćemo spektar gustoće energije (snage). spektar – jer se odnosi na određenu frekvenciju; gustoće – jer se odnosi na jedinicu frekvencije i vremena; energije – jer označuje srednju vrijednost kvadrata funkcije. Pisat ćemo: 𝑋 (𝑓) 𝐺(𝑓) =. 2 Razumljivo da bismo ukupnu prosječnu snagu dobili zbrojem doprinosa snaga pojedinih frekvencija, tj. 𝑥 (𝑡) = 𝐺(𝑓) 𝑑𝑓. Grafički prikaz funkcije 𝐺(𝑓) o frekvenciji naziva se spektralna slika snage. Spektri gustoće snage uglavnom se koriste za karakterizaciju širokopojasnih slučajnih signala (npr. mikroseizmičkog nemira), dakle pojava koje dugo traju (npr. stacionarne pojave svojstva ne mijenjaju s vremenom što implicira beskonačno dugo trajanje i beskonačnu ukupnu energiju) i nisu periodičke, dakle ne udovoljavaju kriterijima da ih se 'napadne' alatima niti harmoničke niti spektralne analize. Uočite i da je ovakvom reprezentacijom spektralnog sastava neke pojave potpuno izgubljena informacija o fazi. 36  M. Herak i D. Skoko: Uvod u spektralnu analizu, skripta Veza funkcije autokorelacije i spektra gustoće snage Funkciju autokorelacije definirali smo ovako: 𝑅(𝜏) = 𝑥(𝑡)𝑥(𝑡 + 𝜏) 𝑑𝑡 , 𝑡 ∈ [0, ∞]. Kao i svaka funkcija vremena, i ona ima svoj spektar. Izrazimo ga kosinusovom i sinusovom transformacijom. Kako je R parna funkcija, sinus-spektar joj je nula: ↗ 𝑎 (𝑓) = 2 𝑅(𝜏) cos 2𝜋𝑓𝜏 𝑑𝜏 𝑅(𝜏) ↘ 𝑏 (𝑓) = 0 pa se 𝑅(𝜏) može pisati kao 𝑅(𝜏) = 𝑎 (𝑓) cos 2𝜋𝑓𝜏 𝑑𝑓. (∗) Analogno ćemo spektar funkcije 𝑥(𝑡) prikazati također kosinus i sinus spektrom: 𝑎 (𝑓) = 2 𝑥(𝑡) cos 2𝜋𝑓𝑡 𝑑𝑡 ↗ 𝑥(𝑡). ↘ 𝑏 (𝑓) = 2 𝑥(𝑡) sin 2𝜋𝑓𝑡 𝑑𝑡 {Ovdje je donja granica integracije nula, jer smo naprijed definirali da vrijeme mjerimo od nule, 𝑡 ∈ [0, ∞].} 𝑥(𝑡) = [𝑎(𝑓) cos 2𝜋𝑓𝑡 + 𝑏(𝑓) sin 2𝜋𝑓𝑡] 𝑑𝑓. Sada imamo [𝑢𝑧 𝜔 = 2𝜋𝑓]: 𝑅(𝜏) = 𝑥(𝑡)𝑥(𝑡 + 𝜏) 𝑑𝑡 = 𝑥(𝑡) [𝑎(𝜔) cos 𝜔(𝑡 + 𝜏) + 𝑏(𝜔) sin 𝜔(𝑡 + 𝜏)]𝑑𝑓 𝑑𝑡 = {𝑧𝑎𝑚𝑖𝑗𝑒𝑛𝑖𝑚𝑜 𝑑𝑓 𝑖 𝑑𝑡} = 37  M. Herak i D. Skoko: Uvod u spektralnu analizu, skripta = 𝑥(𝑡) [𝑎(𝜔) cos 𝜔𝑡 cos 𝜔𝜏 − 𝑎(𝜔) sin 𝜔𝑡 sin 𝜔𝜏 + 𝑏(𝜔) sin 𝜔𝑡 cos 𝜔𝜏 + 𝑏(𝜔) cos 𝜔𝑡 sin 𝜔𝜏]𝑑𝑡 𝑑𝑓 = = 𝑎(𝜔) cos 𝜔𝜏 𝑥(𝑡) cos 𝜔𝑡 𝑑𝑡 − 𝑎(𝜔) sin 𝜔𝜏 𝑥(𝑡) sin 𝜔𝑡 𝑑𝑡 + 𝑏(𝜔) cos 𝜔𝜏 𝑥(𝑡) sin 𝜔𝑡 𝑑𝑡 + 𝑏(𝜔) sin 𝜔𝜏 𝑥(𝑡) cos 𝜔𝑡 𝑑𝑡 𝑑𝑓. Četiri integrala u vitičastim zagradama prepoznajemo kao 𝑎(𝜔)/2 i 𝑏(𝜔)/2, pa dalje imamo: 1 𝑅(𝜏) = [𝑎(𝜔)𝑎(𝜔) cos 𝜔𝜏 − 𝑎(𝜔)𝑏(𝜔) sin 𝜔𝜏 + 2 + 𝑏(𝜔)𝑏(𝜔) cos 𝜔𝜏 + 𝑏(𝜔)𝑎(𝜔) sin 𝜔𝜏]𝑑𝑓. 1 𝑅(𝜏) = [𝑎 (𝑓) + 𝑏 (𝑓)] cos 2𝜋𝑓𝜏 𝑑𝑓 = {(∗)} = 𝑎 (𝑓) cos 2𝜋𝑓𝜏 𝑑𝑓. 2 ( ) Slijedi da je 𝑋 (𝑓) = 𝑎 (𝑓) = 𝐺(𝑓). 2 Time smo izveli važno svojstvo spektra gustoće snage 𝐺(𝑓) – on je jednak spektru (Fourierovoj transformaciji) funkcije autokorelacije te pojave! U literaturi se zato spektri gustoće snage ponekad zovu i autospektralne funkcije gustoće (autospectral density functions). Na taj smo način povezali domenu nedeterminističkih, slučajnih, podataka (autokorelaciju i spektar gustoće snage) i alata koje smo upoznali kod determinističkih tranzijentnih pojava, konkretno Fourierove transformacije. 38  M. Herak i D. Skoko: Uvod u spektralnu analizu, skripta Diracova 𝜹-funkcija Diracova 𝛿-funkcija jedan je od osnovnih alata primijenjene matematike, a osobitu primjenu ima u spektralnoj analizi. Pripada tzv. singularnim funkcijama (singularity functions). Iako je to važna funkcija definirana je na nekoliko fundamentalno različitih načina. Spomenimo ove: 𝛿(𝑡) = 0 za 𝑡 ≠ 0 A. ∫ 𝛿(𝑡)𝑑𝑡 = 1 𝛿(𝑡) 𝑗e nedefinirana za 𝑡 = 0 'Površina' pod 𝛿-funkcijom jednaka je 1, pa je to dakle beskonačni puls u nuli. B. Kao limes: 𝛿(𝑡) = lim 𝑓 (𝑡), → gdje je 𝑓 red funkcija koje zadovoljavaju: 𝑓 (𝑡)𝑑𝑡 = 1, lim 𝑓 (𝑡) = 0 za 𝑡 ≠ 0 → C. Svojstvom: ∫ 𝛿(𝑡)𝑓(𝑡)𝑑𝑡 = 𝑓(0), gdje je f bilo kakva funkcija kontinuirana u ishodištu. D. Svojstvom da je 𝛿-funkcija jednaka inverznoj Fourierovoj transformaciji konstante: 1 𝛿(𝑡) = 1 ∙ 𝑒 𝑑𝜔. 2𝜋 Sve su ove definicije 'besmislene' ako δ-funkciju promatramo kao običnu funkciju. Zato ćemo uvesti novi koncept generalizirane funkcije (funkcionala) ili razdiobe (ne u statističkom smislu). Razdiobom 𝑔(𝑡) ćemo zvati proces koji prozvoljnoj funkciji 𝜑(𝑡) klase C (to su funkcije koje imaju bilo koju derivaciju i trnu prema nuli brže od 𝑡 , n je bilo koji cijeli broj, kad 𝑡 → ∞) pridružuje broj 𝑁 [𝜑(𝑡)]. Za 𝛿-funkciju može se pokazati da je dovoljan uvjet da je neprekinuta u ishodištu. Broj 𝑁 [𝜑(𝑡)] može biti bilo što ovisno o 𝜑(𝑡). Općenito, obično se uzima 𝑔(𝑡)𝜑(𝑡)𝑑𝑡 = 𝑁 [𝜑(𝑡)]. 39  M. Herak i D. Skoko: Uvod u spektralnu analizu, skripta 𝛿-funkciju definira prema definiciji C integral ∫ 𝛿(𝑡)𝜑(𝑡)𝑑𝑡 = 𝜑(0). Za obične funkcije ovaj integral ima smisao sam po sebi. 𝛿-funkcija je, međutim, skupa s integralom definirana brojem 𝜑(0), pa 𝛿 i integral sami po sebi nemaju nikakvo značenje. Promotrimo načas integral jediničnog pravokutnog pulsa 𝑟 (𝑡) definiranog grafom, i neke funkcije 𝜑(𝑡). Očito je da će iznos ovog integrala biti to sličniji 𝜑(0), što je 𝜖 manji, pa imamo: 1 𝜑(0) = lim 𝑟 (𝑡)𝜑(𝑡)𝑑𝑡 = lim 𝜑(𝑡)𝑑𝑡 , → → 𝜖 pa je dakle (iz C) 𝛿(𝑡) = lim 𝑟 (𝑡). → Ovdje je red funkcija 𝑓 (𝑡) iz definicije B koji ima jedan član jednak lim 𝑟. → Tako smo razdiobu 𝛿(𝑡) izrazili limesom pravokutnog pulsa. Mogli smo upotrijebiti i kakvu drugu oko nule simetričnu funkciju, npr. Gaussovu ili kardinalni sinus (vidi kasnije). Kako je površina pravokutnika na slici uvijek jednaka (P = 1), njegovim sužavanjem kako 𝜖 → 0 on mora rasti, pa je očito da 𝛿-funkciju možemo zamisliti kao puls u t = 0, dok je ona jednaka 0 za 𝑡 ≠ 0. Kakva još svojstva ima 𝛿-funkcija? Što je npr. Fourierova transformacija od 𝛿(𝑡), ∫ 𝛿(𝑡)𝑒 𝑑𝑡? Ovdje je 𝜑(𝑡) = 𝑒 , pa je 𝛿(𝑡) 𝑒 𝑑𝑡 = {definicija C} = 𝑒 = 1. Znači spektar 𝛿(𝑡) je konstanta! Ova činjenica ima veliku primjenu, pa ćemo uvijek kada imamo posla s nekom impulsnom silom očekivati bijeli šum kao spektar (bijeli šum je signal u kojem sve frekvencije imaju jednaku amplitudu). 40  M. Herak i D. Skoko: Uvod u spektralnu analizu, skripta Ovo znači da 𝛿-funkciju možete dobiti zbrajajući beskonačno mnogo kosinusa svih frekvencija. Da to vidite, pokušajte zbrajati kosinuse različitih frekvencija (ali bez faznog pomaka) – oni se svi zbroje u t = 0, a za sva druga vremena se kosinusi u sumi ponište. Npr. ovo je suma 100.000 kosinusa s frekvencijama između 0.0 Hz i 500 Hz.: Vrijedi i obrat – bijeli šum ima spektar koji je jednak 𝛿-funkciji. Pokažite to! Inverzna Fourierova transformacija konstante, dakle, daje 𝛿-funkciju, 1 1𝑒 𝑑𝜔 = 𝛿(𝑡), 2𝜋 što je i naša definicija D. Kako je 𝛿(𝑡) = 0 za 𝑡 ≠ 0, to možemo pisati: 𝛿(𝑡)𝑒 𝑑𝑡 = 1 = + + 𝛿(𝑡)𝑒 𝑑𝑡. Prvi i treći integral jednaki su nuli jer je tamo 𝛿 = 0, dok je za drugi integral 𝑒 = 1, pa je 𝛿(𝑡)𝑑𝑡 = 1. Površina ispod 𝛿-funkcije jednaka je 1, što je njezino važno svojstvo po definiciji A. Odavde izravno slijedi i svojstvo C: 𝛿(𝑡)𝑓(𝑡)𝑑𝑡 = 𝑓(0) 𝛿(𝑡)𝑑𝑡 = 𝑓(0). 41  M. Herak i D. Skoko: Uvod u spektralnu analizu, skripta Lako je pokazati da vrijede i ova svojstva: 𝑎) 𝛿(𝑡)𝑓(𝑡 − 𝜏)𝑑𝑡 = {𝑓(𝑡 − 𝜏 = 0)} = 𝑓(𝜏) 𝑏) 𝛿(𝑡 − 𝜏)𝑓(𝑡)𝑑𝑡 = {𝑡 − 𝜏 ≠ 0 → 𝛿(𝑡 − 𝜏) = 0} = = 𝛿(𝑡 − 𝜏)𝑓(𝑡)𝑑𝑡 = 𝑓(𝜏) 𝛿(𝑡 − 𝜏)𝑑𝑡 = {𝑑𝑒𝑓𝑖𝑛𝑖𝑐𝑖𝑗𝑎 𝐴) = 𝑓(𝜏) Svojstva a) i b) zovu se svojstva prosijavanja (engl: sifting property), jer iz funkcije 𝑓(𝑡) izvlače točno određenu vrijednost, onu na koju pokazuje 𝛿- funkcija. 1 𝑐) 𝛿(𝑎𝑡) = 𝛿(𝑡). Dokaz: 𝑎 ∫ 𝛿(𝑎𝑡)𝜑(𝑡) 𝑑𝑡 = {𝑢 = 𝑎𝑡, … } = ∫ 𝛿(𝑢)𝜑(𝑢/𝑎) 𝑑𝑢 = 𝜑(0), a znamo: 𝛿(𝑡)𝜑(𝑡)𝑑𝑡 = {𝑑𝑒𝑓. 𝐶} = 𝜑(0), pa slijedi: 𝑎 𝛿(𝑎𝑡)𝜑(𝑡) 𝑑𝑡 = 𝜑(0) = 𝛿(𝑡)𝜑(𝑡) 𝑑𝑡. Dakle, usporedbom podintegralnih funkcija: 1 𝑎𝛿(𝑎𝑡) = 𝛿(𝑡), 𝛿(𝑎𝑡) = 𝛿(𝑡). 𝑎 Zbog toga je npr. 𝛿 = 𝛿(𝑓) = 2𝜋𝛿(𝜔). 𝑑) 𝑓(𝑡)𝛿(𝑡) = 𝑓(0)𝛿(𝑡). Dokaz: 𝛿(𝑡)𝑓(𝑡) 𝜑(𝑡)𝑑𝑡 = {def. 𝐶} = 𝑓(0)𝜑(0) = 𝛿(𝑡)𝑓(0)𝜑(𝑡)𝑑𝑡 Potcrtani izrazi moraju dakle biti jednaki, pa je dokaz proveden! 42  M. Herak i D. Skoko: Uvod u spektralnu analizu, skripta Pomoću 𝛿-funkcije možemo izračunati i Fourierovu transformaciju funkcija kosinus i sinus u beskonačnom intervalu. Napišimo npr. prvo kosinus pomoću Moivreove formule: cos 𝜔 𝑡 = (𝑒 +𝑒 ). Trebamo izračunati Fourierovu transformaciju: 1 1 cos 𝜔 𝑡 ↔ 𝑒 𝑒 𝑑𝑡 + 𝑒 𝑒 𝑑𝑡 = 2 2 1 ( ) 1 ( ) = 𝑒 𝑑𝑡 + 𝑒 𝑑𝑡 = 2 2 1 ( ) 1 ( ) = 𝑒 𝑑𝑡 + 𝑒 𝑑𝑡 = 2 2 1 [𝛿(𝑓 − 𝑓 ) + 𝛿(𝑓 + 𝑓 )] = {𝛿(𝑓) = 2𝜋𝛿(𝜔)} = 𝜋[𝛿(𝜔 − 𝜔 ) + 𝛿(𝜔 + 𝜔 )] 2 Na isti bismo način dobili da je Fourierov par funkcije sin 𝜔 𝑡: sin 𝜔 𝑡 ↔ [𝛿(𝑓 − 𝑓 ) − 𝛿(𝑓 + 𝑓 )] = 𝑖𝜋[𝛿(𝜔 − 𝜔 ) − 𝛿(𝜔 + 𝜔 )]. 43  M. Herak i D. Skoko: Uvod u spektralnu analizu, skripta Pokažimo sada da se 𝛿-funkcija može dobiti i kao ovakav limes: sin 𝜔𝑡 𝛿(𝑡) = lim. (𝐴) → 𝜋𝑡 Dovoljno je pokazati da desna strana ima osnovno svojstvo 𝛿-funkcije: sin 𝜔𝑡 lim 𝜑(𝑡)𝑑𝑡 = 𝜑(0). → 𝜋𝑡 Rastavimo integral ovako: sin 𝜔𝑡 sin 𝜔𝑡 lim 𝜑(𝑡)𝑑𝑡 = lim + + 𝜑(𝑡)𝑑𝑡. → 𝜋𝑡 → 𝜋𝑡 Uz razumnu pretpostavku da je 𝜑(𝑡)/𝑡 integrabilno u (−∞, −𝜖) i (𝜖, ∞) slijedi prema Riemann-Lebesque-ovoj lemi da su i prvi i treći integral jednaki nuli. Ta funkcija oko nule nije integrabilna jer je t u nazivniku. Digresija – Riemann-Lebesque-ova lema: lim 𝑒 𝜑(𝑡)𝑑𝑡 = 𝜑(0), za 𝑎, 𝑏 bilo koji broj (uključujući i ± ∞), → ako je 𝜑(𝑡) apsolutno integrabilna. Sinus možemo pisati kao sin 𝛼 = (𝑒 +𝑒 ), pa je lema izravno primjenljiva. 𝜑(𝑡) je kontinuirana u nuli, pa se smije pisati 𝜑(𝑡) = 𝜑(0) za 𝑡 ∈ [−𝜖, +𝜖], ako je 𝜖 dovoljno mali: 1 𝜔 ⎧𝜔𝑡 = 𝑥, = ⎫ sin 𝜔𝑡 sin 𝜔𝑡 ⎪ 𝑡 𝑥⎪ 𝜑(𝑡)𝑑𝑡 = 𝜑(0) 𝑑𝑡 = 𝑑𝑥 = 𝜋𝑡 𝜋𝑡 ⎨ 𝑑𝑡 = ⎬ ⎪ 𝜔 ⎪ ⎩ 𝑔𝑟𝑎𝑛𝑖𝑐𝑒! ⎭ 𝜔 sin 𝑥 𝑑𝑥 sin 𝑥 𝜑(0) = 𝜑(0) 𝑑𝑥. 𝑥 𝜋 𝜔 𝜋𝑥 44  M. Herak i D. Skoko: Uvod u spektralnu analizu, skripta Ovo je tablični integral i može se pokazati da u limesu vrijedi: sin 𝑥 sin 𝑥 lim 𝑑𝑥 = 𝑑𝑥 = 1. → 𝜋𝑥 𝜋𝑥 [Izvod vidi u Papoulis, (II-57).] Tako konačno imamo: sin 𝜔𝑡 lim 𝜑(𝑡)𝑑𝑡 = 𝜑(0), → 𝜋𝑡 čime je dokaz proveden. * * * Ovo ćemo primijeniti pri traženju Fourierova para tzv. jedinične funkcije češlja, tj. niza ekvidistantnih 𝛿-funkcija. Ona izgleda ovako: U literaturi ova se funkcija često označava i simbolom III koji podsjeća i na ćirilićno veliko slovo 'Š', pa ga se ponekad 'engleski' čita 'sha' (ša)! Dokazat ćemo da je Fourierov par također funkcija češlja ovoga oblika: Znači dokazat ćemo da je: 45  M. Herak i D. Skoko: Uvod u spektralnu analizu, skripta 𝑠 (𝑡) = 𝛿(𝑡 − 𝑛𝑇 ) ↔ 𝜔 𝑠 (𝜔) = 𝛿(𝜔 − 𝑛𝜔 ). (∗) Dakle Fourierov par je također funkcija češlja s razmakom među zupcima = 𝜔. 'Visina' šiljaka kako je naznačena na slikama jednostavno ukazuje koliko puta je djelovanje svake 𝛿-funkcije (svakoga zupca češlja) jače od 'obične', jedinične. Primijetite da smo gornjom formulom funkciju češlja izrazili sumom 𝛿-funkcija kakve smo do sada upoznali. Indeks označava razmak zubaca pripadne funkcije češlja. DIGRESIJA – možemo li si na temelju onoga što smo do sada naučili i vizualno predočiti da je naša tvrdnja istinita? 𝜔 𝑠 (𝜔) izgleda ovako: Kako smo malo prije vidjeli, spektar funkcije cos 𝜔𝑡 su zapravo dvije 𝛿- funkcije (pomnožene s 𝜋) na frekvencijama ±𝜔. Tako, npr. dvije 𝛿-funkcije na poziciji 𝜔 (plave) zapravo predstavljaju spektar funkcije cos 𝜔 𝑡, a zelene su spektar od cos 2𝜔 𝑡. Srednja (crvena) je spektar od cos 0 = = 𝑘𝑜𝑛𝑠𝑡. – to smo već prije pokazali, 𝛿-funkcija je spektar konstante. Zbrojimo li sve kosinuse cos 𝑛𝜔 𝑡 za 𝑛 → ∞, oni će se uvijek zbrajati za 𝑡 = = 𝑛𝑇 , a za druge vrijednosti će se poništavati, pa ćemo dobiti 'češalj' s razmakom šiljaka jednakim 𝑇. Dokažimo to sada i egzaktno! 46  M. Herak i D. Skoko: Uvod u spektralnu analizu, skripta Da bismo dokazali tvrdnju (*), dovoljno je pokazati da je inverzna Fourierova transformacija od 𝜔 𝑠 (𝜔) jednaka 𝑠 (𝑡): 1 2𝜋 1 𝛿(𝜔 − 𝑛𝜔 ) 𝑒 𝑑𝜔 = 𝛿(𝜔 − 𝑛𝜔 ) 𝑒 𝑑𝜔 = 2𝜋 𝑇 𝑇 𝑥 = 𝜔 − 𝑛𝜔 1 1 = 𝑑𝑥 = 𝑑𝜔 = 𝑒 𝛿(𝑥) 𝑒 𝑑𝑥 = 𝑒. 𝜔 = 𝑥 + 𝑛𝜔 𝑇 𝑇 ( ) Ovo je kompleksni oblik Fourierovog reda (s je periodička funkcija!). Uvedimo sada realnu funkciju 𝑘 (𝑡) kao sumu: 1 𝑘 (𝑡) = 𝑅𝑒 𝑒. To je dakle suma kosinusa. 𝑇 Suma geometrijskog reda općenito je jednaka: 1−𝑞 𝑎 + 𝑎 𝑞 + ⋯+ 𝑎 𝑞 = 𝑎. 1−𝑞 Ovdje je 𝑞 = 𝑒 ,𝑎 = 𝑒 , a članova ima 𝑛 = 2𝑁 + 1. Zato je: ( ) ( ) 1 1−𝑒 1𝑒 −𝑒 𝑘 (𝑡) = 𝑅𝑒 𝑒 = 𝑅𝑒. 𝑇 1−𝑒 𝑇 𝑒 −1 Realni dio ovog izraza je: 1 sin 𝑁 + 𝜔 𝑡 𝑘 (𝑡) = 2 𝜔 𝑡 , 𝑇 sin 2 i zove se kernel (jezgra) Fourierovog reda. To je periodična funkcija s periodom 𝑇 = 2𝜋/𝜔 i izgleda otprilike ovako: 47  M. Herak i D. Skoko: Uvod u spektralnu analizu, skripta Ako dokažemo da u 𝑘 (𝑡) u intervalu [− , + ] teži prema 𝛿(𝑡) za 𝑁 → ∞, slijedit će: lim 𝑘 (𝑡) = 𝑠 (𝑡), što i tražimo. → , 1 𝑡 sin 𝑁 + 𝜔 𝑡 𝑡 𝑘 (𝑡) = 2 𝑡 𝑇 𝑡 𝜔 𝑡 sin 2 1 sin 𝑁 + 𝜔 𝑡 𝜋 lim 2 = 𝛿(𝑡). To smo dokazali ranije − jednadžba (𝐴). → 𝑇 𝑡 𝑇 Dakle, kako je 𝑡/ sin ograničeno u [−𝑇 /2, +𝑇 /2], za |𝑡| < , vrijedi: 𝜋 𝑡 lim 𝑘 (𝑡) = 𝛿(𝑡) 𝜔 𝑡 = {iz svojstva 𝑓(𝑡)𝛿(𝑡) = 𝑓(0)𝛿(𝑡)} = → , 𝑇 sin 2 𝜋 𝑡 = 𝛿(𝑡). 𝑇 sin 𝜔 𝑡 2 Kako je 𝜋 𝜔 𝑡 𝜋 𝑡 𝑡 𝑇 2 = 𝜔 𝑡= = 1, 𝑇 sin 𝜔 𝑡 sin sin 𝜔 𝑡 2 2 2 imamo konačno: / lim [𝑘 (𝑡)] / = 𝛿(𝑡), → pa je i inverzna Fourierova transformacija od 𝜔 𝑠 (𝜔) jednaka 𝑠 (𝑡). 48  M. Herak i D. Skoko: Uvod u spektralnu analizu, skripta Analiza digitalnih podataka Ovdje ćemo razmotriti specifičnosti analize i probleme s kojima se možemo susresti pri analizi digitalnih podataka. Pretpostavit ćemo da je vremenski diskretni slijed podataka takav da je vremenski razmak dva susjedna podatka (interval digitalizacije ili interval uzorkovanja) uvijek jednak. Prije nego se priđe spektralnoj analizi podataka, često su nužne određene prethodne obrade 'sirovih' podataka. Na početku ćemo samo ukratko razmotriti jednu od osnovnih koja se praktički uvijek provodi, eliminiranje trenda. Eliminiranje trenda Eliminiranje trenda spomenut ćemo ovdje isključivo radi potpunosti predavanja, jer se nadam da su tehnike kojima se to postiže već poznate. Trend se definira kao komponenta podataka čiji je period veći od vremenskog intervala u kojemu poznajemo podatke. Za to se koriste standardni statistički postupci. Uglavnom je to postupak najmanjih kvadrata, koji se rabi najčešće za eliminiranje linearnog trenda, ali i trendova koji se opisuju polinomima drugog ili višeg reda. Na slici je prikazan seizmogram potresa (crno) koji očito ima pozitivni linearni trend. Do toga može doći npr. zbog polaganog zagrijavanja ili hlađenja seizmometra u zoru ili sumrak. To rezultira superpozicijom seizmograma potresa na kosinusoidu vrlo dugačkoga perioda (≈ 24 sata). Lokalno, tijekom intervala vremena prikazanoga slikom (koji je mnogo manji od 24 sata) taj trend moguće je aproksimirati pravcem (plavo). Postupkom uklanjanja trenda dobije se korigirani seizmogram (crveno), koji oscilira strogo oko nule. Da je trend važno ukloniti pokazuje donja slika – spektri amplituda 'sirovog' (crno) i korigiranog seizmograma (crveno). Vidimo da se oni znatno razlikuju za niske frekvencije, jer je FFT algoritam u obzir uzeo da u signalu postoje i komponente vrlo malih frekvencija (dugačkih perioda). Ovdje je zapravo riječ o superpoziciji seizmograma na signal beskonačno velikog perioda (pravac). 49  M. Herak i D. Skoko: Uvod u spektralnu analizu, skripta Dakle, uklanjanje trenda je važan, često i neizostavan korak pripreme podataka za analizu. Najčešće se to postiže jednostavnom prilagodbom pravca regresije (u smislu najmanjih kvadrata) 'sirovim' podacima (plavo), kojega se oduzme od mjerenja. Dakle, simbolički pisano u našem slučaju: [crveno] = [crno] – [plavo]. Naravno, trend može biti dobro opisan i npr. kvadratnom (ili kojom drugom) funkcijom, pa ćemo postupkom najmanjih kvadrata naći takvu parabolu koja najbolje opisuje mjerene podatke, te je konačno oduzeti od podataka. U detalje postupka određivanja pravca regresije ovdje nećemo ulaziti. Postupak uklanjanja trenda često je u poznatim programskim paketima posve olakšan jer su dostupne gotove prikladne rutine (npr. detrend u Matlabu). Priprema podataka obično uključuje i modificiranje signala (eng. tapering) kako bi se minimizirao utjecaj činjenice da uvijek na raspolaganju imamo signal konačne duljine. O tome nešto kasnije. 50  M. Herak i D. Skoko: Uvod u spektralnu analizu, skripta Digitalno filtriranje Filtriranje podataka podrazumijeva izdvajanje samo onih spektralnih komponenti koje nas zanimaju. To se obično postiže tako da se iz signala nastoje što potpunije izbaciti sve one frekvencije koje nas ne zanimaju ili smetaju pri obradi podataka. Filtriranju se mogu podvrgnuti kako analogni, tako i digitalni signali. Iako i analogno filtriranje ima svoju važnu primjenu (vidi npr. kasnije u poglavlju o uzorkovanju), ovdje ćemo se baviti samo najosnovnijim tehnikama digitalnog filtriranja, dakle filtriranja digitalizi- ranog, diskretnog, signala. Filtriranje podataka izvodi se radi:  izglađivanja podataka,  izdvajanja interesantnog pojasa frekvencija (npr. smanjenja utjecaja šuma),  usporedbe svojstava pojave u pojedinim frekventnim pojasevima. Neka su originalni podaci označeni s 𝑥(𝑡), kao na slici: Ova slika može prikazivati, npr. mareogram (zapis mareografa, instrumenta koji bilježi razine mora). Vidi se da se sastoji od dugih oscilacija (uočava se period od oko 24 sata, a vidi se još barem jedna spektralna komponenta perioda oko 12 sati) na koje su superponirane oscilacije mnogo kraćih perioda (najizraženije su one od oko 20 minuta). Ako je ovo mareogram, tada one duge oscilacije prikazuju morska doba (plimu i oseku), dok superponirani visoko- frekventni signal možda odgovara sešima (to su stojni valovi ili vlastite oscilacije zaljeva ili bazena u kojem se mareograf nalazi; npr. Bakarski zaljev). 51  M. Herak i D. Skoko: Uvod u spektralnu analizu, skripta Ovisno o tome na čemu je trenutno Vaš fokus istraživanja, možda Vas zanimaju morska doba (pa su seši neka vrsta šuma), a možda baš želite proučiti spektralni sastav seša pa Vam smetaju dugoperiodičke oscilacije. U oba slučaja pomoći će digitalno filtriranje, kao što možete vidjeti na slici. Ovdje su na srednjem i donjem panelu izdvojene dugoperiodičke (niskofrekventne) i kratkoperiodičke (visokofrekventne) komponente opaženog signala. U prvom slučaju primijenili bismo neki od niskopropusnih (low-pass, LP) filtara, dok bi se u drugom slučaju koristili visokopropusnim (high-pass, HP) filtriranjem. Dakle filtri imena dobivaju prema dijelu frekvencija koji propuštaju. Kvantitativno, postupak filtriranja možemo opisati izrazom: 52  M. Herak i D. Skoko: Uvod u spektralnu analizu, skripta 𝑌(𝑓) = 𝑍(𝑓)𝑋(𝑓), gdje su: 𝑋(𝑓) – spektralna slika osnovnih podataka, 𝑌(𝑓) – spektralna slika podataka nakon filtriranja, 𝑍(𝑓) – karakteristika filtra. Sjetimo se da smo istu formulu već jednom napisali – u poglavlju o linearnim fizikalnim sustavima. Kako je to linearna jednadžba, i filtre koji se njome opisuju zvat ćemo linearnim filtrima. Jasno je da će i u ovom slučaju vrijediti principi aditivnosti i homogenosti. O karakteristici filtra ovisit će koje će frekvencije biti istaknute, a koje umanjene. Primjeri: 𝑍(𝑓) = 1.0 Nema filtriranja, 𝑋(𝑓) = 𝑌(𝑓). 𝑍(𝑓) = 25.5 Pojačalo 1.0 … |𝑓| ≤ 5 𝑍(𝑓) = Niskopropusni (LP, low-pass) filtar 0.0 … |𝑓| > 5 1.0 … |𝑓| ≥ 10 𝑍(𝑓) = Visokopropusni (HP, high-pass) filtar 0.0 … |𝑓| < 10 1.0 … 7 ≤ |𝑓| ≤ 15 𝑍(𝑓) = Pojasno propusni (BP, band-pass) filtar 0.0 … 𝑜𝑠𝑡𝑎𝑙𝑜 0.0 … 7 ≤ |𝑓| ≤ 15 𝑍(𝑓) = Pojasno nepropusni (BS, band-stop) filtar 1.0 … 𝑜𝑠𝑡𝑎𝑙𝑜 Filtriranje u prostoru frekvencija je dakle vrlo jednostavno – spektar pojave modificiramo tako da sadrži samo one frekvencije koje želimo, napravimo inverznu Fourierovu transformaciju – i gotovo! Postavlja se, međutim, pitanje možemo li zadani vremenski niz filtrirati u prostoru vremena, dakle bez prethodnog računanja njegovog spektra. Pogledajmo dakle, kakav je odnos vremenskih slika podataka 𝑥(𝑡), 𝑦(𝑡) i karakteristike filtra. 53  M. Herak i D. Skoko: Uvod u spektralnu analizu, skripta Primijenit ćemo i ovdje Fourierove transformacije, pa tako za vremensku sliku filtriranih podataka vrijedi: 𝑦(𝑡 ) = 𝑌(𝑓)𝑒 𝑑𝑓 = 𝑋(𝑓)𝑍(𝑓) 𝑒 𝑑𝑓. Označimo sada s 𝑧(𝜏) vremensku sliku karakteristike filtra [𝑧(𝜏) ↔ 𝑍(𝑓)], pa imamo: 𝑧(𝜏) = 𝑍(𝑓) 𝑒 𝑑𝑓 , 𝑍(𝑓) = 𝑧(𝜏) 𝑒 𝑑𝜏, a posljednji izraz uvršten u onaj za 𝑦(𝑡 ) daje 𝑦(𝑡 ) = 𝑋(𝑓) 𝑧(𝜏)𝑒 𝑑𝜏 𝑒 𝑑𝑓. Zamjenom redoslijeda integriranja imamo: ( ) 𝑦(𝑡 ) = 𝑧(𝜏) 𝑋(𝑓) 𝑒 𝑑𝑓 𝑑𝜏, odnosno 𝑦(𝑡 ) = 𝑧(𝜏)𝑥(𝑡 − 𝜏)𝑑𝜏 = 𝑧(𝑡) ∗ 𝑥(𝑡) ili, za vrijeme kao kontinuiranu varijablu: 𝑦(𝑡) = 𝑧(𝜏)𝑥(𝑡 − 𝜏)𝑑𝜏 = 𝑧(𝑡) ∗ 𝑥(𝑡) Ovo je integral konvolucije, jedan od najvažnijih integrala primijenjene matematike. Vidimo da za zadanu karakteristiku filtra 𝑧(𝜏), filtriranu vrijednost niza u trenutku 𝑡 možemo izračunati kao integral u koji ulaze sve točke originalnog niza 𝑥(𝑡)! Konvoluciju dva niza se uobičajeno označava zvjezdicom (∗), da nas podsjeti kako množenju u prostoru frekvencija odgovara konvolucija u prostoru vremena: 𝑍(𝑓)𝑋(𝑓) ↔ 𝑧(𝑡) ∗ 𝑥(𝑡), a vrijedi i obrat: 54  M. Herak i D. Skoko: Uvod u spektralnu analizu, skripta 𝑧(𝑡)𝑥(𝑡) ↔ 𝑍(𝑓) ∗ 𝑋(𝑓). Dakle, umnošku spektara dva signala odgovara konvolucija ta dva signala u vremenu, a spektar umnoška dvije vremenske funkcije jednak je konvoluciji njihovih spektara. Množenje i konvolucija su Fourierov par! Pogledajmo i neka važnija svojstva konvolucije: a) Konvolucija s 𝛿-funkcijom: 𝑓(𝑡) ∗ 𝛿(𝑡) = 𝑓(𝜏)𝛿(𝑡 − 𝜏)𝑑𝜏 = 𝑓(𝑡) pa je i 𝑓(𝑡) ∗ 𝛿(𝑡 − 𝑡 ) = 𝑓(𝜏)𝛿(𝑡 − 𝑡 − 𝜏)𝑑𝜏 = 𝑓(𝑡 − 𝑡 ). Dakle konvolucija s 𝛿-funkcijom jednostavno premješta funkciju f na mjesto 𝛿-funkcije! Primijetimo da ako za ulazni signal odaberemo 𝛿-funkciju [𝑥(𝑡) = 𝛿(𝑡)], filtrirani će signal 𝑦(𝑡) = ∫ 𝑧(𝜏)𝛿(𝑡 − 𝜏)𝑑𝜏 prema ovome svojstvu biti jednak 𝑧(𝑡). Zato 𝑧(𝑡) često nazivamo i odzivom filtra (ili sustava) na impuls (engl. impulse response). b) Svojstvo komutacije: 𝑥(𝑡) ∗ 𝑦(𝑡) = 𝑦(𝑡) ∗ 𝑥(𝑡) 𝑥(𝑡) ∗ 𝑦(𝑡) = 𝑥(𝜏)𝑦(𝑡 − 𝜏)𝑑𝜏 = {𝜆 = 𝑡 − 𝜏, 𝑑𝜏 = −𝑑𝜆, granice → –granice} =− 𝑥(𝑡 − 𝜆)𝑦(𝜆)𝑑𝜆 = 𝑦(𝜆)𝑥(𝑡 − 𝜆)𝑑𝜆 = 𝑦(𝑡) ∗ 𝑥(𝑡) c) Svojstvo asocijativnosti: [𝑥(𝑡) ∗ 𝑦(𝑡)] ∗ 𝑧(𝑡) = 𝑥(𝑡) ∗ [𝑦(𝑡) ∗ 𝑧(𝑡)]. Izvod je dosta dugačak, ostavljam ga Vama! 55  M. Herak i D. Skoko: Uvod u spektralnu analizu, skripta d) Svojstvo distributivnosti: 𝑥(𝑡) ∗ [𝑦 (𝑡) + 𝑦 (𝑡)] = 𝑥(𝑡) ∗ 𝑦 (𝑡) + 𝑥(𝑡) ∗ 𝑦 (𝑡) Dokaz: 𝑥(𝑡) ∗ [𝑦 (𝑡) + 𝑦 (𝑡)] = 𝑥(𝜏)[𝑦 (𝑡 − 𝜏) + 𝑦 (𝑡 − 𝜏)]𝑑𝜏 = 𝑥(𝜏)𝑦 (𝑡 − 𝜏)𝑑𝜏 + 𝑥(𝜏)𝑦 (𝑡 − 𝜏)𝑑𝜏 = 𝑥(𝑡) ∗ 𝑦 (𝑡) + 𝑥(𝑡) ∗ 𝑦 (𝑡) e) Derivacija konvolucije: 𝑦 (𝑡) = 𝑥 (𝑡) ∗ 𝑧(𝑡) 𝑦(𝑡) = 𝑥(𝜏) 𝑧(𝑡 − 𝜏)𝑑𝜏 Vremenska derivacija od 𝑦(𝑡) je 𝑦 (𝑡) = 𝑥(𝜏) 𝑧 (𝑡 − 𝜏)𝑑𝜏 = 𝑥(𝑡) ∗ 𝑧 (𝑡) = {komutacija!} = 𝑥 (𝑡) ∗ 𝑧(𝑡). f) Svojstvo površine: Ako je 𝑦(𝑡) = 𝑥(𝑡) ∗ 𝑧(𝑡), površina ispod y je produkt površina ispod 𝑥 i z. 𝐴(𝑦) = 𝐴(𝑥) 𝐴(𝑧) 𝐴(𝑦) = 𝑦(𝑡)𝑑𝑡 = 𝑥(𝜏) 𝑧(𝑡 − 𝜏)𝑑𝜏 𝑑𝑡 = 𝑥(𝜏)𝑑𝜏 𝑧(𝑡 − 𝜏)𝑑𝑡 = 𝐴(𝑥) 𝐴(𝑧). g) Svojstvo skaliranja: Ako je 𝑦(𝑡) = 𝑥(𝑡) ∗ 𝑧(𝑡) tada vrijedi: 𝑦(𝑎𝑡) = |𝑎|𝑥(𝑎𝑡) ∗ 𝑧(𝑎𝑡). 56  M. Herak i D. Skoko: Uvod u spektralnu analizu, skripta Kako grafički provesti konvoluciju dvije funkcije, z(𝜏) i x(𝜏)? Konvolucija u trenutku 𝑡 je 𝑦(𝑡 ) = ∫ 𝑧(𝜏)𝑥(𝑡 − 𝜏)𝑑𝜏 Valja izračunati površinu ispod produkta dvije funkcije: 𝑧(𝜏) i 𝑥(𝑡 − 𝜏). Funkcija 𝑥(𝑡 − 𝜏) je funkcija negativnog vremena 𝜏. Ako znamo kako izgleda 𝑥(𝜏), kako izgleda 𝑥(𝑡 − 𝜏)? Najprije postavimo 𝑡 = 0, i uočimo da ćemo zrcaljenjem 𝑥(𝜏) oko ordinate dobiti funkciju 𝑥′(𝜏) = 𝑥(−𝜏). Ako 𝑥′ pomaknemo u desno za 𝑡 dobit ćemo funkciju 𝑥 (𝜏) = 𝑥 (−𝑡 + 𝜏) = 𝑥(𝑡 − 𝜏). Dakle, 𝑥(𝑡 − 𝜏) dobijemo iz 𝑥(𝜏) tako da 𝑥(𝜏) zrcalimo oko ordinate i pomaknemo je za 𝑡 u desno (engl. 'flip and shift'). Primjer 1. 57  M. Herak i D. Skoko: Uvod u spektralnu analizu, skripta Primjer 2. Zadatak: Izračunajte konvoluciju funkcije 𝑧(𝜏) iz prethodnog zadatka same sa sobom! 58  M. Herak i D. Skoko: Uvod u spektralnu analizu, skripta Idealan niskopropusni (low-pass) filtar Definirajmo takav filtar njegovom karakteristikom 𝑍(𝑓): 1, |𝑓| ≤ 𝑓 𝑍(𝑓) = 0, |𝑓| > 𝑓 Dakle, ovaj će filtar propustiti samo frekvencije koje su manje ili jednake 𝑓 , dok će ostale u potpunosti potisnuti. Funkcija obli- ka 𝑍(𝑓) zove se i boxcar funkcija, i idealna je da se iz nekog niza 'izvuče' samo određeni interval. Kako izgleda funkcija 𝑧(𝜏) koja je, kako smo vidjeli, jednaka Fourierovoj transformaciji 𝑍(𝑓)? 𝑧(𝜏) = 𝑍(𝑓)𝑒 𝑑𝑓 = 𝑒 𝑑𝑓 1 1 1 = 𝑒 −𝑒 = sin 2𝜋𝑓 𝜏. 𝜋𝜏 2𝑖 𝜋𝜏 Ovdje je sinc(𝑥) = , funkcija koja se zove kardinalni sinus. 59  M. Herak i D. Skoko: Uvod u spektralnu analizu, skripta 𝑧(𝜏) je funkcija sinus frekvencije 𝑓 modificirana funkcijom 1/𝜏. Najviša je u nuli, gdje iznosi 2𝑓 (jer je lim sinc(𝑥) = 1!), poprima i pozitivne i negativne → vrijednosti, a oscilacije trnu kako 𝜏 → ±∞. Ovo je u spektralnoj analizi jako važna funkcija upravo zbog svoje veze s boxcar funkcijom: Zapamtimo: Kardinalni sinus i boxcar funkcija su Fourierov par. Uočimo i da će za 𝑓 → ∞, karakteristika filtra 𝑍(𝑓) → 1 ∀𝑓, visina funkcije 𝑧(𝜏) će rasti u beskonačnost, a njezina širina će se smanjivati prema nuli – 𝑧(𝜏) će tada težiti 𝛿-funkciji! Ranije smo pokazali da to svojstvo vrijedi i za jedinični pravokutni puls (boxcar funkciju) kada 𝑓 → 0. Iako se zovu ‘idealnim’ (jer ‘idealno’ režu neželjene frekvencije) ovakvi filtri se u praksi rijetko upotrebljavaju. Da vidimo zašto, pogledajmo kako izgleda filtrirana 𝛿-funkcija. Izračunajmo, dakle konvoluciju 𝑧(𝜏) i δ(𝜏): 60  M. Herak i D. Skoko: Uvod u spektralnu analizu, skripta Dakle, idealnim LP filtriranjem 𝛿-funkcije, tj. pulsa u 𝜏 = 𝑡, dobije se na istom mjestu kardinalni sinus. Međutim, iako je impuls nastao u 𝜏 = 𝑡, filtrirani signal postoji i prije toga vremena, što je, naravno fizikalno nemoguće, dakle signal postoji prije uzroka! To se zove akauzalnost, i ona će biti to jače izražena što je filtar uži. U praksi ovo može jako smetati – ako je npr. riječ o početku seizmograma, i ako na njega prije analize primijenimo idealni filtar da se smanji utjecaj visokofrekventnog šuma, često ćemo nastupno vrijeme potresa pročitati prerano jer će se pokrajnji ekstremi kardinalnog sinusa činiti kao početak potresa. Vidjeli smo, dakle, da ćemo filtriranje postići pomoću konvolucije mjerenog signala 𝑥(𝑡) i vremenske karakteristike filtra 𝑧(𝑡): 𝑦(𝑡) = 𝑧(𝜏)𝑥(𝑡 − 𝜏)𝑑𝜏. Primijenimo sada taj izraz na diskretne vrijednosti mjerenoga niza: 61  M. Herak i D. Skoko: Uvod u spektralnu analizu, skripta Neka je naš niz digitaliziran svakih h sekundi (kasnije ćemo naučiti kako izabrati h), i neka se sastoji od N podataka. Tada općenito, diskretno vrijeme, možemo označiti s t = nh, dok ćemo 𝜏 označiti s 𝜏 = 𝑘ℎ. Sada integral konvolucije možemo pisati u diskretnom obliku, pretvorivši ga u sumu: 𝑧(𝜏)𝑥(𝑡 − 𝜏)𝑑𝜏 → 𝑦 = 𝑧(𝑘ℎ)𝑥[(𝑛 − 𝑘)ℎ]ℎ = ℎ𝑧 𝑥 (∗) gdje je