############################################
#                                          #
#  Cvičení č. 6 - Fourierovy transformace  #
#                                          #
############################################


install.packages(c("fourierin","dplyr","purrr"))
library("fourierin")
library("dplyr")
library("purrr")

#####################################
#                                   #
#  Spojitá Fourierova transformace  #
#                                   #
#####################################

# V praktickém zpracování časových řad je obvykle spojitá Fourierova transformace bez užitku, nicméně
# pro úplnost uveďme, že v R existuje knohovna "fourierin", která se umí určitým způsobem s výpočtem
# spojité Fourierovy transformace popasovat. Naneštěstí se tak neděje analyticky (jako v přpadě CAS
# systémů typu Maple) ale opět numericky. Pro uživatele je „příjemné“, že vstupem funkce fourierin
# je „spojitá“ funkce, tedy objekt třídy function, nakonec je ale výpočet stejně realizován pomocí DFT,
# což s sebou nese nežádoucí artefakty, které uvidíme níže ve výsledcích.


## Nejprve nadefinujeme funkci, jejíž Fourierovu trasformaci budeme chítí spočítat.
# V totmo případě bude mít tři harmonické složky o úhlových frekvencích 1/2, 1/4 a 1/5.
fce1<-function(t) {
  return(sin(2*t)-0.6*cos(4*t+0.5)+0.3*cos(5*t-1))  # perioda bude rovna 2*pi (protoze 2,4, a 5 jsou nesoudelne)
}

# Vykreslení funkce (použijme 400 bodů za sebou, aby byl obrázek hezky hladký)
t<-seq(0,6*pi,length.out=1000)
plot(t,fce1(t),type="l",lwd=3,col="forestgreen",main="Fce1",xlab="Čas",ylab="Hodnota")

# Nastavení „rozlišení“ pro funkci fourierin
rozliseni<-1024

# Výpocet Fourierova integrálu
spojita_four_transf<-function(rozliseni){
  out <- fourierin(f=fce1,lower_int=0,upper_int=20*pi,lower_eval=0,upper_eval=7,const_adj=-1,freq_adj=-1,resolution=rozliseni)
  out %>%
    as.data.frame() %>%
    transmute(
      x=w,
      amplituda=Mod(values),     # Mod() vraci modul komplexniho cisla
      faze=Arg(values))
}
hodnoty<-map_df(rozliseni,spojita_four_transf)

# Vykreslení frekvenčního spektra
par(mfrow=c(2,1)) # vykreslení dvou grafů nad sebe do jednoho okna
plot(hodnoty$x,hodnoty$amplituda,type="l",lwd=3,col="forestgreen",main="Amp1",xlab="Frekvence",ylab="Hodnota")
plot(hodnoty$x,hodnoty$faze,type="l",lwd=3,col="brown",main="Faz1",xlab="Frekvence",ylab="Hodnota")
par(mfrow=c(1,1))

# S určitými nepřenostmi danými použitím numerického výpočtu a tedy FFT vidíme, že v amplitudovém spektru nám vysokčily jasně
# tři frekvence na hodnotách 2, 4 a 5. Hůře čitelné je fázové spektrum, kde dochází k oscilacím u méně významných složek, nicméně
# pro tři hlavní složky s největšími amplitudami zde můžeme odečíst hodnoty zhruba -pi/2 pro sinus a 1/2 (tj. zhruba 0.16*pi) a
# -1 (tj. zhruba 0.32*pi) pro oba cosiny.



#######################################
#                                     #
#  Diskrétní Fourierova transformace  #
#                                     #
#######################################

# Běžná a v praxi velmi často užívaná je diskrétní Fourierova transformace (DFT), resp. její výpočetní implementace rychlá
# Fourierova transformace (příkaz fft()).

## Využijeme funkce z předchozího příkladu a navzorkujeme časovou řadu v délce tří period.
t<-t<-seq(0,6*pi,length.out=61) # délka bude 61, tj. dvacet vzorků na periodu a pak jeden na závěr, zopakovaný z nuly
vzorky<-fce1(t)

# Protože nejvyšší ze tří frekvenci, ktere jsou v časové řadě zastoupeny, má hodnotu 5, musíme zovlit více než dvojnásobnou
# (dle Shannonova-Nyquistova-Kotelnikovova teorému) vzorkovací frekvenci. Proto volíme 20>2*5.

# Zobrazení v grafu
plot(t,vzorky,pch=19,cex=2,col="navyblue")
lines(c(-10,20),c(0,0)) # osa x
for (i in 1:length(t)) {
  lines(c(t[i],t[i]),c(0,vzorky[i]),lwd=2,col="navyblue")
}
  
# DFT (vezmeme vzorky za jednu periodu)
transformace<-fft(vzorky[1:20])/20  # získáváme komplexní hodnoty spektra - aby amplitudy odpovídaly, je v R třeba ještě vydělit spketrum N (=10)
  
par(mfrow=c(2,1)) # vykreslení dvou grafů nad sebe do jednoho okna
plot(0:19,Mod(transformace),pch=19,cex=2,col="blue",main="Amp1",xlab="Frekvence",ylab="Hodnota")
lines(c(-10,20),c(0,0))
for (i in 0:19) {
  lines(c(i,i),c(0,Mod(transformace)[i+1]),lwd=2,col="blue")
}
plot(0:19,Arg(transformace),pch=19,cex=2,col="red3",main="Faz1",xlab="Frekvence",ylab="Hodnota")
lines(c(-10,20),c(0,0))
for (i in 0:19) {
  lines(c(i,i),c(0,Arg(transformace)[i+1]),lwd=2,col="red3")
}
par(mfrow=c(1,1)) 
  
# V diskrétním spektru jsou nyní krásně vidět nenulové harmonické složky pro hodnoty 2, 4 a 5 - tedy dvojnásobek, čtyřnásobek a pětinásobek
# základní frekvence. Čára v 0 je nulová => původní časová řada má nulový průměr. Ampitudové spektrum je symetrické okolo poloviny vzorkovací
# frekvence (10), dle očekávání. Výšky čar jsou tedy v navzájem symetrických hodnotách poloviční (0.5 + 0.5 pro obě složky odpovídající
# sinu v 2 a v 8 apod.). Fázové spektrum obsahuje fázi -pi/2 pro sinovou složku a -2.64 (ekvivalentně 0.5-pi) a -1 pro obě cosinové složky.



#####################################
#                                   #
#  Několik užitečných příkladů DFT  #
#                                   #
#####################################

## Spektrum konstantní časové řady
fce2<-c(1,1,1,1,1,1,1,1)
spektrum2<-fft(fce2)/8

par(mfrow=c(2,1)) # vykreslení dvou grafů nad sebe do jednoho okna
plot(0:7,Mod(spektrum2),pch=19,cex=2,col="blue",main="Amplitudové spektrum",xlab="Frekvence",ylab="Hodnota",ylim=c(0,max(Mod(spektrum2))))
lines(c(-10,20),c(0,0))
for (i in 0:7) {
  lines(c(i,i),c(0,Mod(spektrum2)[i+1]),lwd=2,col="blue")
}
plot(0:7,Arg(spektrum2),pch=19,cex=2,col="red3",main="Frekvenční spektrum",xlab="Frekvence",ylab="Hodnota",ylim=c(-pi,pi))
lines(c(-10,20),c(0,0))
for (i in 0:7) {
  lines(c(i,i),c(0,Arg(spektrum2)[i+1]),lwd=2,col="red3")
}
par(mfrow=c(1,1)) 

# Jediná nenulová složka je v nule, tj. konstantní harmonická složka o nulové frekvenci (na fázi pak nezáleží)



## Spektrum harmonické časové řady
fce3<-c(1,sqrt(1/2),0,-sqrt(1/2),-1,-sqrt(1/2),0,sqrt(1/2))
spektrum3<-fft(fce3)/8

par(mfrow=c(2,1)) # vykreslení dvou grafů nad sebe do jednoho okna
plot(0:7,Mod(spektrum3),pch=19,cex=2,col="blue",main="Amplitudové spektrum",xlab="Frekvence",ylab="Hodnota",ylim=c(0,max(Mod(spektrum3))))
lines(c(-10,20),c(0,0))
for (i in 0:7) {
  lines(c(i,i),c(0,Mod(spektrum3)[i+1]),lwd=2,col="blue")
}
plot(0:7,Arg(spektrum3),pch=19,cex=2,col="red3",main="Frekvenční spektrum",xlab="Frekvence",ylab="Hodnota",ylim=c(-pi,pi))
lines(c(-10,20),c(0,0))
for (i in 0:7) {
  lines(c(i,i),c(0,Arg(spektrum3)[i+1]),lwd=2,col="red3")
}
par(mfrow=c(1,1)) 

# Jediná nenulová složka je v jedničce, protože její frekvence odpovídá právě jedné periodě (vzorkovanému úseku)




## Spektrum jednotkového impulzu
fce4<-c(1,0,0,0,0,0,0,0)
spektrum4<-fft(fce4)/8

par(mfrow=c(2,1)) # vykreslení dvou grafů nad sebe do jednoho okna
plot(0:7,Mod(spektrum4),pch=19,cex=2,col="blue",main="Amplitudové spektrum",xlab="Frekvence",ylab="Hodnota",ylim=c(0,max(Mod(spektrum4))))
lines(c(-10,20),c(0,0))
for (i in 0:7) {
  lines(c(i,i),c(0,Mod(spektrum4)[i+1]),lwd=2,col="blue")
}
plot(0:7,Arg(spektrum1),pch=19,cex=2,col="red3",main="Frekvenční spektrum",xlab="Frekvence",ylab="Hodnota",ylim=c(-pi,pi))
lines(c(-10,20),c(0,0))
for (i in 0:7) {
  lines(c(i,i),c(0,Arg(spektrum4)[i+1]),lwd=2,col="red3")
}
par(mfrow=c(1,1)) 

# Všechny frekvenční složky jsou zastoupeny stejnou měrou a mají stejnou fázi



## Spektrum posunutého jednotkového impulzu
fce5<-c(0,1,0,0,0,0,0,0)
spektrum5<-fft(fce5)/8

par(mfrow=c(2,1)) # vykreslení dvou grafů nad sebe do jednoho okna
plot(0:7,Mod(spektrum5),pch=19,cex=2,col="blue",main="Amplitudové spektrum",xlab="Frekvence",ylab="Hodnota",ylim=c(0,max(Mod(spektrum5))))
lines(c(-10,20),c(0,0))
for (i in 0:7) {
  lines(c(i,i),c(0,Mod(spektrum5)[i+1]),lwd=2,col="blue")
}
plot(0:7,Arg(spektrum5),pch=19,cex=2,col="red3",main="Frekvenční spektrum",xlab="Frekvence",ylab="Hodnota",ylim=c(-pi,pi))
lines(c(-10,20),c(0,0))
for (i in 0:7) {
  lines(c(i,i),c(0,Arg(spektrum5)[i+1]),lwd=2,col="red3")
}
par(mfrow=c(1,1)) 

# Všechny frekvenční složky jsou zastoupeny stejnou měrou a příslušně posunuty do jedničky



## Spektrum bílého šumu
fce6<-rnorm(8,0,1)
spektrum6<-fft(fce6)/8

par(mfrow=c(2,1)) # vykreslení dvou grafů nad sebe do jednoho okna
plot(0:7,Mod(spektrum6),pch=19,cex=2,col="blue",main="Amplitudové spektrum",xlab="Frekvence",ylim=c(0,max(Mod(spektrum6))))
lines(c(-10,20),c(0,0))
for (i in 0:7) {
  lines(c(i,i),c(0,Mod(spektrum6)[i+1]),lwd=2,col="blue")
}
plot(0:7,Arg(spektrum6),pch=19,cex=2,col="red3",main="Frekvenční spektrum",xlab="Frekvence",ylab="Hodnota",ylim=c(-pi,pi))
lines(c(-10,20),c(0,0))
for (i in 0:7) {
  lines(c(i,i),c(0,Arg(spektrum6)[i+1]),lwd=2,col="red3")
}
par(mfrow=c(1,1)) 

# Náhodně zastoupené složky, povšimněte si vždy sudého amplitudového a vždy lichého fázového spektra



## Spektrum jednotkového skoku
fce7<-c(0,0,0,0,1,1,1,1)
spektrum7<-fft(fce7)/8

par(mfrow=c(2,1)) # vykreslení dvou grafů nad sebe do jednoho okna
plot(0:7,Mod(spektrum7),pch=19,cex=2,col="blue",main="Amplitudové spektrum",xlab="Frekvence",ylim=c(0,max(Mod(spektrum7))))
lines(c(-10,20),c(0,0))
for (i in 0:7) {
  lines(c(i,i),c(0,Mod(spektrum7)[i+1]),lwd=2,col="blue")
}
plot(0:7,Arg(spektrum7),pch=19,cex=2,col="red3",main="Frekvenční spektrum",xlab="Frekvence",ylab="Hodnota",ylim=c(-pi,pi))
lines(c(-10,20),c(0,0))
for (i in 0:7) {
  lines(c(i,i),c(0,Arg(spektrum7)[i+1]),lwd=2,col="red3")
}
par(mfrow=c(1,1)) 

# Hyperbolický pokles v lichých frekvencích, nulová hodnota v sudých frekvencích a 1/2 v nule (odpovídající průměru časové řady).










  
  
  
  
  
  
  