FFT Home     Rachunek prawdopodobieństwa i statystyka


Transformata Fouriera

Dyskretna, szybka = FFT

Potrzebny nam będzie R.
The R Project for Statistical Computing. Dlaczego R, a nie S? S już było zajęte, teraz nikt o tym już nie pamięta.
Nie jest potrzebna znajomość języka R, ale w konsoli jakoś się trzeba poruszać.
Zróbmy (w R) prosty szereg czasowy x=c(1,2,3,4,3,6,7,8) wykreślmy go plot(x) , ma trend i jakiś błąd na pozycji 5.

Zróbmy jakiś ciekawszy szereg, biały szum, gaussowski: x=rnorm(128) r - random, norm - wiadomo. Gdzie siedzi, że biały szum? W niezależności kolejnych wartości. W tym, że one nie maja tu żadnej kolejności, czyli przypadkową mają. Jeżeli mamy robić FFT, to dodajmy do naszego szeregu jednak szczyptę sinusa, o okresie ... 32. Najpierw zróbmy czas t=c(0:127) to teraz s=sin(2*pi/32*t) i teraz dodajmy sinusa do x. Ale najpierw zachowajmy sobie x w jakimś schowku x0=x teraz x=x+s na wykresie niewiele widać, ten sinus trochę się gubi w szumie, ale o to chodzi, zobaczymy do czego służy FFT. No to y=fft(x) teraz uwaga, piszemy y i ENTER.
Co to jest? Dlaczego przed i nie ma spacji? I że i co? Ok. wpiszmy 1i^2 i ENTER.

Teoria. Po co i? Zróbmy tak: z=exp(2i*pi/32*t) i plot(z)
exp(it) = cos(t) + i*sin(t) stąd kółko na wykresie. Hasło: wzór Eulera. Czy i jest tu konieczne? Nie, ale tak jest elegancko, matematycznie. Tak się robi, jak widać. W sumie chodzi o fazę sin(2*pi/T*t-faza), żeby ją mieć w jednej liczbie razem z amplitudą. Mamy transformacje sinusową i cosinusową (używaną w jpg), które obchodzą się bez liczb zespolonych.

Wracamy do naszego wyniku FFT, nazywa się y. Sprawdźmy y[2] i y[128] (y[1] to inna sprawa) otrzymujemy:
-5.490172-7.874846i
-5.490172+7.874846i Przypadek? Nie ma takich przypadków, z dokładnością do 7 cyfr. Zrobiliśmy transformatę wektora (szeregu) 128 liczb (rzeczywistych). A w wyniku otrzymaliśmy 128 zespolonych, czyli dwa razy więcej liczb (rzeczywistych, bo zesopolona składa się z dwóch rzeczywistych). Skąd one się wzięły, te dodatkowe. Nie wzięły się. Mamy tylko 128 różnych liczb rzeczywistych (w tych zespolonych). Liczby zespolone a+bi i a-bi nazywają się sprzężonymi. Różnych liczb zespolonych mamy tu 64, następne 64 są sprzężone z tymi pierwszymi, czyli nic nowego. Na wykresie plot(y) widać tę symetrię.
Czy to jest najciekawszy wykres transformaty Fouriera naszego szeregu czasowego? Nie, nie jest. Chcemy zobaczyć amplitudy, po prostu wartości liczb zebranych w wektorze y. Policzyć trzeba moduły liczb z=a+bi: √(a2+b2), ym=Mod(y) i plot(ym). Też widać symetrię, połowa wykresu niepotrzebna. Więc: ym=ym[1:65] , połowa plus 1. Zauważmy, że w y mamy dwie unikalne wartości, y[1] i y[65], obie są rzeczywiste, obu potrzebujemy. Wykreślamy ym (aktualne, przycięte), jednak musimy to zrobić porządnie. Przygotujmy zmienną f=c(0:64) (litera f nieprzypadkowa) i plot(f,ym) poniższy wykres, patrz Notatki.
Co to f? Co mamy na osi poziomej?
Ciekawa rzecz, ale nie mamy tu wykresu analizowanego szeregu czasowego. Już mamy plot(t,x,type='l'):
Tu oś pozioma nazywa się t.
Czy widać tu ślad sinusa o okresie 32? Trochę widać. Na FFT za to widać, ym dla f=4 - to ten sinus! Dlaczego na czwórce? Bo 128/32=4. Cztery okresy sinusa mieszczą się w naszym szeregu. Co to jest f=4? To nie całkiem 4 Hz, 4 na sekundę (bo co to tutaj sekunda?), tylko 4 cykle na cały odcinek czasu (128 [może i sekund, może dni]). To jest częstotliwość, frequency.

Transformata Fouriera przenosi nas z domeny czasu do domeny częstotliwości time domain frequency domain. (Ale i z powrotem potrafi.)
Δf * Δt = 1 Kojarzy nam się zasada nieoznaczoności Heisenberga? U nas f=4, a t=1/4, przy jednostce czasu określonej przez długość szeregu. Do zasady nieoznaczoności jeszcze wrócimy.
Przyjrzyjmy się wykresowi ym=f(f) - moduł amplitudy (czego?) składowych harmonicznych w funkcji częstotliwości. f=4 - dodany sinus, f=0? Częstotliwość zero oznacza okres T=∞, czyli coś co nazywa się składowa stała. U nas niemal zero. Nic dziwnego szum rnorm ma μ=0, a sinus, jak to sinus, kręci się wokół zera (tym bardziej, że mamy tu cztery kompletne okresy, bez żadnych ścinków - 128/32=4 dokładnie).
Popatrzmy dokładnie: y[1]=-1.069716+0i, średnia szeregu mean(x)=-0.008357157. Podzielmy te liczby, otrzymamy 128. Tak Oni to liczą. Nie podzielili przez liczbę punktów szeregu. No to jaka jest wartość transformaty dla f=4, czyli y[5]/128=0.0326416-0.5742255i. Czynnik przy cosinusie mały, ok, cosinusa nie było w modelu, ale czynnik przy sinusie nie 1 jak w modelu, tylko 0.57. Jakby połowa? No tak. Przecież jest jeszcze druga składowa, sprzężona, w y[125] (w ym obcięta dla celów wykresu). OK! (Składowa stała występuje tylko raz.)
Wróćmy do częstotliwości f. Weźmy ostatni punkt wykresu ym=f(f=64). Popatrzmy na y[65]/128 = 0.1093646+0i. Liczba rzeczywista, nieduża, losowa, bez przyczyny. Ona też występuje tylko raz, jak f=0.
f=0 tak wygląda: -----------, a f=64 tak: /\/\/\/\/\/\/\/\/\. Co w nich podobnego, oba przebiegi są ekstremalne i nie ma w nich miejsca na fazę. Wykresu 1111111 nie da się przesunąć w lewo, ani w prawo, podobnie 0101010101, też się nie da, tzn. da się, o pi, ale to tyle, co zmiana znaku /\ na \/. Zróbmy tak: fft(c(0,1,0,1,0,1,0,1)) wynik:   4+0i   0+0i   0+0i   0+0i -4+0i   0+0i   0+0i   0+0i.
Jeszcze raz. Po co i? Zróbmy tak, przesuwajmy (niby) sinusa w lewo: Arg(fft(c(0,1,2,1,0,-1,-2,-1))[2]) = -1.57 Arg(fft(c(-1,0,1,2,1,0,-1,-2))[2]) = -2.35 Arg(fft(c(-2,-1,0,1,2,1,0,-1))[2]) = 3.141593 = π Arg(fft(c(-1,-2-1,0,1,2,1, 0))[2]) = 2.56 Arg( fft( c(0,-1,-2-1,0,1,2, 1) )[2]) = 1.66
Arg zwraca argument liczby zespolonej, czyli kąt. Kręci się. Popatrzmy na wynik takiej komendy:
plot(c(0+0i, -1+1i, -1-1i, 1-1i, 1+1i))
Arg(1+1i)/pi = 0.25 = 45°

Transformata odwrotna. Nie dlatego, że musimy ją poznać (dlatego, też), ale żeby jeszcze lepiej zrozumieć wynik transformaty prostej.
Zróbmy nowy szereg, krótszy, prostszy: p=c(0,1,2,3,4,5,6,7,0,1,2,3,4,5,6,7) (p od piła) transformatę nazwijmy pt=fft(p) teraz z powrotem p1=fft(pt) i patrzymy, ale na p1/16 ENTER no zespolone są, ale wartości te same. Gdyby kogoś kusiło uprościć sobie życie i policzyć p-p1/16, to owszem można, ale odczytanie takiej liczby wymaga uwagi: 0-3.330669e-16i za dużo liter. x1=fft(y, inverse = TRUE) też jest ok, równe x. Tylko po co nam to?

Filtrowanie. Chcemy z danych wydłubać znalezionego tam sinusa (zapominamy, że to nasz model). Gdzie w transformacie siedzi sinus? W y[5] i y[125]. Resztę trzeba wyzerować. Tak to robimy.
filtr=rep(0,128); tr[5]=1; filtr[125]=1; y1=y*filtr; x1=fft(y1, inverse = TRUE); x1=Re(x1)/128; plot(x1)
Dlaczego Re zamiast dotychczas używanego Mod, do powrotu do realu? Po prostu, sinus bywa ujemny, a przy i i tak są zera. Powyższa, jednoliniowa komenda, jeżeli ktoś chciałby coś zapamiętać, powinna być realizowana krok po kroku, z wypisywaniem sobie wyników i ewentualnie ich plot'owaniem.
plot(x,type='l'); par(new = TRUE); plot(x1,type='l', col='red')    wynik:
Z osią pionową coś się komuś nie udało. Nie jest to ładne, ale błędne nie jest. To NIE jest narysowana funkcja sinus. Owszem, to jest dokładnie sinus (bo wzięliśmy tylko jeden punkt, jedną częstotliwość), ale niedokładnie ten z naszego modelu. Niedokładnie, bo to tylko statystyka. Transformacja Fouriera, to nie statystyka, tylko przekształcenie bez straty informacji. Natomiast filtrowanie, gdzie wyrzuciliśmy cały szum (cały?), to JEST statystyka.
Szersze okno filtru: fil1=c(rep(1, 10), rep(0, 108), rep(1, 10)) i dalej jak poprzednio y1=y*fil1; x1=fft(y1, inverse = TRUE); x1=Re(x1)/128; plot(x,type='l'); par(new = TRUE); plot(x1,type='l', col='red')

Jeszcze szersze: fil1=c(rep(1, 30), rep(0, 68), rep(1, 30))

Górnoprzepustowy (zera na jedynki i odwrotnie): fil1=c(rep(0, 30), rep(1, 68), rep(0, 30))

Aliasing. Włączmy ten programik Uwaga, przy najwyższych częstotliwościach czerwonego sinusa zaczyna się jeszcze jedno próbkowanie, ujawniają się piksele ekranu. Ten programik pokazuje sprawę szprych albo helikoptera. Koło ma tu jedną szprychę, a szybkość filmowania (czy flesza - szprycha robi się biała) jest taka, że jest, mniej więcej jedna klataka (frame) na obrót koła. Można to zmieniać.
Wykonajmy tę komendę:
t=c(0:31); x=sin(2*pi*4/32*t)+sin(2*pi*8/32*t); y=fft(x); ym=Mod(y); plot(ym[2:17])
Mamy tu krótszy szereg, żadnego szumu, a w sinusach zamiast okresu podajemy czestotliwość (zamiast 1/8 (T=8) piszemy 4/32 (f=4 cykle na szereg (na 32 jednostki czasu))). Nie rysujemy też ym[1] - niepotrzebnej składowej stałej (której tu i tak nie ma). Wykres jest klarowny. "Index" to częstotliwość, mamy 4 i 8.
Dodajmy jeszcze jeden, bardzo szybki sinus (najszybszy możliwy, Nyquista)
y=fft(x+sin(2*pi*16/32*(t-0.5))); ym=Mod(y); plot(ym[2:17])
Dlaczego tu akurat (t-0.5)? Wypróbujmy bez fazy 0.5 i pomyślmy nad wynikiem (np. z pomocą pierwszego programiku). Dlaczego amplituda przy f=16 dwa razy większa niż dla f=4 i 8 (i każdej innej częstotliwości)? Bo ta nie ma fazy i jest tylko w t[17], a tamte, np. w t[16] i sprzężone, w t[18] (było o tym).
Skorygujmy powyższą formułę w konsoli z 16 na 17. Amplituda przeniosła się z pozycji 16 na pozycję 15! Oraz 17, dlatego amplituda spadła dwukrotnie. Badamy f=18. Sygnał mamy w f=14. Itd. Za wysokie czestotliwości nie giną, odbijają się od granicy nazywanej częstotliwością Nyquista. To może wprowadzac w błąd. Psuć nagranie. ... the anti-aliasing filter blocks frequencies above the Nyquist frequency from being converted. This is going to prevent any signals from changing frequency during the conversion process. źródło

Potęga dwójki. O co chodzi z tą długością szeregu, że ma być 2n? A jak nie jest? FFT w R działać będzie, dopóki szereg nie będzie za długi. W Excelu (tam jest FFT) długość szeregu musi być potegą dwójki. Dawniej komputery nie były tak szybkie. Chodzi o liczbę mnożeń. Jeżeli mamy szereg o dowolnej długości, a algorytm przełyka tylko potęgę 2, to robimy zero padding, uzupełnianie zerami końca szeregu, aż uzyskamy długość 1024, czy coś takiego. Oczywiście, najpierw odejmujemy średnią od danych, jak to zwykle przy transformacji Fouriera dobrze jest robić.

Garść transformat

Fala prostokątna (tu nie potrzebujemy fft, sami sobie złożymy funkcję z sinusów (1/1+1/3+1/5+1/7+...)
t=c(0:999)
j=1
x=rep(0,1000); for(i in 1:j){f=i*2-1; x=x+sin(2*pi/500*f*t)/f}; plot(x,type='l')
Dla j=1 oczywiście mamy sinus. Dajmy j=2, 3, ..., 100. Jeżeli chcemy dokładnie zobaczyć ten skok, to wydłużmy okres z 500 do 10000, itd. To w rzeczywistości, w elektronice, na prawdę tak wygląda (Googlujemy square wave oscillogram).

Gauss
czyli zasada nieoznaczoności Heisenberga
t=c(-100:100)
par(mfrow=c(4,1))
x=exp(-0.5*(t/20)^2); plot(x,type='l'); plot(Mod(fft(x))[2:100],type='l')
x=exp(-0.5*(t/3)^2); plot(x,type='l'); plot(Mod(fft(x))[2:100],type='l')
Szeroki Gauss i wąska transformata (też Gauss). I odwrotnie.
Idźmy dalej. Najwęższy możliwy "Gauss" (nazywa się δ Diraca), tak się go robi:
x=0*t; x[100]=1    a tak transformuje i wszystko wykreśla
par(mfrow=c(2,1)); plot(x); plot(Mod(fft(x)))     "Nieskończenie" wąski pik wymaga wszystkich możliwych częstotliwości, z tą samą wagą. Jak wyglądałaby transformata funkcji stałej? Wiadomo, tylko f=0 by było, raszta zero.
Do kompletu jeszcze (jak już mowa o mechanice kwantowej) paczka falowa:
t=c(-500:500); x=rep(0,1001); for(f in -50:50){x=x+exp(-0.5*(f/20)^2)*sin(2*pi/30*(1+f/300)*t)}; plot(x,type='l')
czyli garść (101) sinusów o podobnych częstotliwościach, ważonych Gaussem, żeby ładnie wyszło. Zauważmy, że exp(-x2) nie jest tu funkcją t, tylko f, ma ona bardzo niebezpośredni związek z osią poziomą wykresu (czyli osią czasu, t). A jednak paczka ma kształt gaussowski. Transformata Gaussa jest Gaussem (wąskiego, szerokim i vice versa), a mechanika kwantowa (cząstek elementarnych) nazywana jest też mechaniką falową. O dualiźmie korpuskularno-falowym uczy się w szkole.
W szkole (muzycznej) jest też o dudnieniach
t=c(0:999); plot(sin(2*pi*0.03*t)+sin(2*pi*0.033*t),type='l')
Strunę stroimy poprzez porównanie wysokości dźwięku, który ona wydaje, ze struną wzorcową, np. koncertmistrza. Dopóki obie struny nie grają z dokładnie tą samą częstotliwością, słyszymy dudnienia. (Dokładnie ta sama częstotliwość oznacza, że różnica jest tak mała, że dudnienie są tak wolne fdudnień=(f1-f2)/2, że dźwięk szybciej zanika niż zdążymy dudnienia usłyszeć.)

Szum
Biały szum:    par(mfrow=c(2,1)); x=rnorm(100); plot(x,type='l'); plot(Mod(fft(x))[1:50],type='l')     Drugi wykres ma mniej punktów, bo wykreślamy moduł, gubiąc fazę (symetryczny wykres źle wygląda). Tak, czy inaczej z białego szumu mamy biały szum.
Szum różowy. Jak się robi różowy szum, czy czerwony? Obcina się, albo osłabia wysokie częstotliwości. Czyli tak, bierzemy biały szum, transformujemy go, wycinamy wysokie częstotliwości i retransformujemy. Czy ma sens tę retransformatę z powrotem transformować, żeby zobaczyć, jak wygląda transformata różowego szumu? Nie ma sensu.
Robimy filtr:     fil=c(-50:49); fil=abs(fil)/50     od 1 do 0 i z powrotem,
filtrujemy transformatę     y=fil*fft(x)     i ją wykreślamy    plot(Mod(y),type='l')
retransformujemy y     xf=fft(y,inverse=TRUE)    i wykreślamy    plot(Re(xf),type='l')
Na końcu warto mieć dwa wykresy:     plot(x,type='l'); plot(Re(xf),type='l')     białego szumu i tego samego ale już czerwonawego, gładszego. Re przerabia liczby formalnie zespolone (na prawdę rzeczywiste) na rzeczywiste.

Faza
Nazwijmy fazę p (phase, a f już zajęte na frequency), na początek     p=0     robimy czas     t=c(0:99)     i sinusa o okresie 50     s=sin(2*pi/50*(t-p))     Fazę można by wstawić inaczej, ale tak mamy ją w tych samych jednostkach, co czas. Jeżeli okres sinusa jest 50, to p=50 nic nie zmieni, ale p=25, albo p=12.5 ...
Idźmy od początku, ale zamieńmy sinus na cosinus:
p=0; s=cos(2*pi/50*(t-p)); y=fft(s); Arg(y[3])     wychodzi 0 (5e-16?), zgadza się z p! Dajmy p=12.5, wynik: -π/2. Liczmy Arg(y[3])*50/(2*pi) od razu będzie p. Dlaczego cosinus? Tak jest zdefiniowana transformata.
FFT szeregu czasowego pozwala nam odczytać amplitudę i fazę, czyli kompletną funkcję harmoniczną. Wnioskowanie statystyczne jest następujące. Widzimy w transformacie wysoką wartość, sprawdzamy jej istotność statystyczną. Być może uwzględniamy jakąś wiedzę a priori istosujemy twierdzenie Bayesa. Jeżeli uznamy, że dany punkt transformaty jest istotny, to odczytujemy jego amplitudę i fazę i mamy wszystko. Ale z jaką dokładnością?

Statystyka
Widzimy w transformacie jakieś maksimum, takie sobie, ale jednak. Każdy szum i jego transformata mają jakieś maksimum. Kiedy, przy jakiej amplitudzie możemy uznać, że jednak stoi za tym jakiś sinus?
t=c(0:99); x0=rnorm(100);
A=2
s=A*cos(2*pi/10*(t-2.5)); x=x0+s; plot(x,type='l')     raczej to sinus z dodatkiem szumu, niż odwrotnie. A jest dość duże.
y=fft(x)
Mod(y[11])/100*2     =1.88, nie najgorzej
Arg(y[11])*10/(2*pi)     =-1.128, j.w.
Przy amplitudzie cosinusa A=2 i szumie o σ=1, składowa harmoniczna odtworzona jest dobrze, błąd amplitudy 6%, błąd fazy (-2.5+2.43)/10*100%=0.7%, może przypadkowo tak mały. No właśnie szum, to szum. Musimy to wszystko powtórzyć z nowymi szumami. I, oczywiście z innymi Amplitudami.
x=rnorm(100)+A*cos(2*pi/10*(t-2.5)); y=fft(x); Mod(y[11])/100*2; Arg(y[11])*10/(2*pi)     to możemy parę razy powtórzyć. Jednak trzeba przejśc do większel liczby powtórzeń, może 100. Zdefiniujmy sobie dwa schowki:
am=c(1:100); ph=c(1:100)
I zróbmy tak:
for(i in 1:100){
x=rnorm(100)+A*cos(2*pi/10*(t-2.5)); y=fft(x)
am[i]=Mod(y[11])/100*2; ph[i]=Arg(y[11])*10/(2*pi)}
Następnie
mean(am); sd(am)
to jest średnia z otrzymanych amplitud i ich odchylenie standardowe. Wyników było 100, więc odchylenie standardowe średniej=sd/√100=sd/10. A więc np. tak: 2.002±0.015, trochę jakby za dobrze, ale tak wyszło. Faza?
mean(ph); sd(ph) Wyszło (w pewnym, konkretnym przypadku) -2.489±0.010 tu już nie jest za dobrze, jest dobrze, różnica od oczekiwanej wartości (-2.5) niemal dokładnie równa odchyleniu standardowemu wyniku. Belive me.
Wszystko zbieramy razem, żeby obniżać A i przyglądać się zestandaryzowanym błędom amplitudy i fazy (powinny być 0±1):
A=2
for(i in 1:100){x=rnorm(100)+A*cos(2*pi/10*(t-2.5)); y=fft(x); am[i]=Mod(y[11])/100*2; ph[i]=Arg(y[11])*10/(2*pi)}; (mean(am)-A)/(sd(am)/10); (mean(ph)+2.5)/(sd(ph)/10)
ta linia powinna się wkopiować do konsoli jako jedna linia, wtedy łatwo zmieniać A i ponownie ją wykonywać. Wyniki sa takie:
A, błąd A, błąd fazy (błędy standaryzowane)
2   0.97 -1.4
1  -0.77 -2.0
0.5 3.24 -1.6
0.2 5.0   2.7
0.1 9.8   4.9

Wszystko jasne, amplituda składowej harmonicznej 1, czyli równa odchyleniu standardowemu szumu, zapewnia jeszcze dobre parametry. Przy A=0.5 rzeczy zaczynają się psuć, przy A=0.2 jest źle, a przy A=0.1 beznadziejnie (konkretnie, średnia amplituda wyszła 0.19, dwa razy za duża, ale trzeba pamietać, że jest tam też szum, którego amplituda ma już tu duże znaczenie).
To oczywiście był model. Wszystko pod kontrolą (niemal wszysto, w końcu tam był rand). W życiu jest inaczej. Mamy szereg, robimy jego transformatę i widzimy jakiś wystający punkt. Zawsze będzie on miał jakąś amplitudę i fazę, ale czy te liczby coś znaczą? Miejmy świadomość problemu, który ma pewne rozwiązania, wymagające pewnych założeń. Najczęściej oczekujemy jakiejś składowej periodyczne, mamy jakąś wiedzę a priori, której pewność trudno skwantyfikować.

 

 

 

Notatki x=rnorm(128) t=c(0:127) rep(0,128) s=sin(2*pi/32*t) x0=x; x=x+s y=fft(x) y[2]; y[128] z=exp(2i*pi/32*t) ym=Mod(y); ym=ym[1:65] f=c(0:64) plot(f,ym,type='b',col='blue',axes=FALSE); axis(side=1, at=c(0:64)) plot(t,x,type='l') plot(x); par(new = TRUE); plot(x1) - dwa na jednym Arg(y[2]) length(x) x=rep(0,1000); for(i in 1:j){f=i*2-1; x=x+1/f}; plot(x) - cztery wykresy. save(x,file="C:\\Users\\Laurem\\Desktop\\zajecia\\x.RData") load("C:\\Users\\Laurem\\Desktop\\zajecia\\x.RData") ls() - listowanie nazw zmiennych Ctrl+L - czyści konsolę (nie czyści zmiennych)