####################################################
#                                                  #
#  Cvičení č. 8 - Vnější popis lineárních systémů  #
#                                                  #
####################################################

#install.packages("nleqslv")
library("nleqslv")


#############
# Příklad 1 #
#############

# V příkladu 1 není, co řešit, protože jde jen o slovní popis.


#############
# Příklad 2 #
#############

# Definice systému pomocí diferenční rovnice
y[n]<-3*x[n]-2*x[n-1]+2*y[n-1]-y[n-2]

# Za předpokladu nulových počátečních podmínek lze systém přímo použít pro výpočet výstupů

x<-c(0,0,rep(1,20)) # vstup x (konstantní o hodnotě 1)
y<-rep(0,3)           # počáteční podmínky pro y (hodnoty 0 pro časy -2, -1 a 0)

for(i in 0:19) {
  n<-i+3                                # posun v čase o 3, abychom zamezili negativním indexům
  y[n]<-3*x[n]-2*x[n-1]+2*y[n-1]-y[n-2]
}

plot(-2:19,y,pch=19,main="Výstup systému\r\n y[n]=3·x[n]-2·x[n-1]+2·y[n-1]-y[n-2]\r\n s konstatním vstupem x = 1, 1, 1,…",xlab="Čas")

# Odvození obrazové přenosové funkce je popsáno v textu
H<-function(z) {
  return((3-2*z^-1)/(1-2*z^-1+z^-2))
}

φ<-seq(0,2*pi,length.out=1000) # vykreslíme pro 1000 bodů na jednotkové kružnici
z<-cos(φ)+(0+1i)*sin(φ)

# Zobrazení obrazové= přenosové funkce (závislost magnitudy na normalizované frekvence děleno π)
plot(φ/pi,abs(H(z)),type="l",lwd=3,col="forestgreen",xlab="Normalizovná f [rad×π]",ylab="Magnituda [1]")                            # standardní graf pro z od 0 do 2*pi
plot((φ/pi)[1:500],abs(H(z))[1:500],type="l",lwd=3,col="forestgreen",xlab="Normalizovná f [rad×π]",ylab="Magnituda [1]")            # standardní graf pro z od 0 do pi
plot(φ/pi,20*log10(abs(H(z))),type="l",lwd=3,col="forestgreen",xlab="Normalizovná f [rad×π]",ylab="Magnituda [dB]")                 # graf v dB pro z od 0 do 2*pi
plot((φ/pi)[1:500],20*log10(abs(H(z)))[1:500],type="l",lwd=3,col="forestgreen",xlab="Normalizovná f [rad×π]",ylab="Magnituda [dB]") # graf v dB pro z od 0 do pi

# Dle očekávání konstantní vstup (tedy vstup s nulovou frekvencí, třeba náš c(1,1,1,…))
# je zesilován až do nekonečné intenzity, tak, jak je vidětr na grafu pro nulovou normalizovanou frekvenci.
# To odpovídá našemu prvotnímu výsledku s rostoucí posloupností výstupů při vstupu x = 1, 1, 1,…

#############
# Příklad 3 #
#############

# Definice systému pomocí obrazové přenosové funkce
H<-function(z) {
  return((1+z+z^2)/(2*z*(1-z)))
}

φ<-seq(0,pi,length.out=500)
z<-cos(φ)+(0+1i)*sin(φ)

plot(φ/pi,20*log10(abs(H(z))),type="l",ylim=c(-50,50),lwd=3,col="forestgreen",xlab="Normalizovná f [rad×π]",ylab="Magnituda [dB]") # graf v dB pro z od 0 do pi

# Odvození diferenční rovnice je popsáno v textu
y[n]<--0.5*x[n]-0.5*x[n-1]-0.5*x[n-2]+y[n-1]

# Vykreslení
x<-rep(1,22) # vstup x (konstantní o hodnotě 1)
y<-rep(0,2)  # počáteční podmínky pro y (hodnoty 0 pro časy -1 a 0)

for(i in 1:20) {
  n<-i+2                                # posun v čase o 2, abychom zamezili negativním indexům
  y[n]<--0.5*x[n]-0.5*x[n-1]-0.5*x[n-2]+y[n-1]
}

plot(y,pch=19,main="Výstup systému\r\n y[n]=-0.5·x[n]-0.5·x[n-1]-0.5·x[n-2]+y[n-1]\r\n s konstatním vstupem x = 1, 1, 1,…",xlab="Čas")

# Kromě počátečního úseku výstup systému při konstatním vstupu lineárně klesá.


#############
# Příklad 4 #
#############

# Odvození obrazové přenosové funkce je popsáno v textu


#############
# Příklad 5 #
#############

# Určení nulových bodů a pólů
H<-function(z) {
  return((1-z^-2)/(1-z^-1+0.5*z^-2))
}

φ<-seq(0,pi,length.out=500)
z<-cos(φ)+(0+1i)*sin(φ)

# Modulová charakteristika
plot(φ/pi,20*log10(Mod(H(z))),type="l",ylim=c(-20,20),lwd=3,col="forestgreen",xlab="Normalizovná f [rad×π]",ylab="Magnituda [dB]") # graf v dB pro z od 0 do pi

# Fázová charakteristika
plot(φ/pi,Arg(H(z)),type="l",ylim=c(-pi/2,pi/2),lwd=3,col="cornflowerblue",xlab="Normalizovná f [rad×π]",ylab="Fáze [rad]") # graf v dB pro z od 0 do pi

# Vyjádříme jako koeficienty u z seřazené dle mocnin 
citatel   <-1-z^-2          
jmenovatel<-1-z^-1+0.5*z^-2

cit<-c(1,0,1)
jme<-c(0.5,-1,1)

polyroot(cit) # nulové body -1 a 1
polyroot(jme) # póly 0.5+0.5i a 0.5-0.5i

abs(polyroot(jme)) # pöly uvnitř jednotkové kružnice

# V souladu s grafem vidíme, že systém je stabilní.


#############
# Příklad 6 #
#############

# Určení frekvenčních charakteristik
H<-function(z) {
  return(1+z^-1)
}
z<-function(Ω) {
  return(exp(1i*Ω))
}

Ω<-seq(0,pi,length.out=500)

# Modulová charakteristika
plot(Ω/pi,20*log10(Mod(H(z(Ω)))),type="l",ylim=c(-40,10),lwd=3,col="orange",xlab="Frekvence Ω [fvz/2]",ylab="Modul [dB]")

# Fázová charakteristika
plot(Ω/pi,Arg(H(z(Ω))),type="l",ylim=c(-pi/2,0),lwd=3,col="navyblue",xlab="Frekvence Ω [fvz/2]",ylab="Fáze [rad]")


#############
# Příklad 7 #
#############

# Definice systému pomocí obrazové přenosové funkce
H<-function(z) {
  return(1-z^-1)
}

φ<-seq(0,pi,length.out=500)
z<-cos(φ)+(0+1i)*sin(φ)

plot(φ/pi,20*log10(Mod(H(z))),type="l",ylim=c(-40,10),lwd=3,col="forestgreen",xlab="Normalizovná f [rad×π]",ylab="Magnituda [dB]") # graf v dB pro z od 0 do pi

# Odvození diferenční rovnice je popsáno v textu
y[n]<-x[n]-x[n-1]

x<-c(0,1,rep(0,9)) # vstup x (jednotkový impulz v čase 0)
y<-c(0)  # počáteční podmínky pro y (hodnota 0 pro čas 0)

for(i in 0:9) {
  n<-i+2                                # posun v čase o 2, abychom zamezili negativním indexům
  y[n]<-x[n]-x[n-1]
}

plot(-1:9,y,pch=19,main="Výstup systému\r\n y[n]=x[n]-x[n-1]\r\n se vstupem x = 1, 0, 0,…",xlab="Čas")


#############
# Příklad 8 #
#############

# Definice systému pomocí obrazové přenosové funkce
H<-function(z) {
  return((1-2*z^-2+z^-3)/(1-z^-1))
}

φ<-seq(0,pi,length.out=500)
z<-cos(φ)+(0+1i)*sin(φ)

plot(φ/pi,20*log10(Mod(H(z))),type="l",ylim=c(0,10),lwd=3,col="forestgreen",xlab="Normalizovná f [rad×π]",ylab="Magnituda [dB]") # graf v dB pro z od 0 do pi

# Odvození diferenční rovnice je popsáno v textu
y[n]<-x[n]-2*x[n-2]+x[n-3]+y[n-1]

x<-c(0,0,0,1,rep(0,9)) # vstup x (jednotkový impulz v čase 0)
y<-rep(0,4)  # počáteční podmínky pro y (hodnoty 0 pro časy -3, -2, -1 a 0)

for(i in 0:9) {
  n<-i+4                                # posun v čase o 4, abychom zamezili negativním indexům
  y[n]<-x[n]-2*x[n-2]+x[n-3]+y[n-1]
}

plot(-3:9,y,pch=19,main="Výstup systému\r\n y[n]=x[n]-x[n-1]\r\n se vstupem x = 1, 0, 0,…",xlab="Čas")


#############
# Příklad 9 #
#############

# Výpočet je uveden v textu
H<-function(z,a) {
  return(z/(z-a))
}

φ<-seq(0,pi,length.out=500)
z<-cos(φ)+(0+1i)*sin(φ)

# Modulová charakteristika pro a=0,5
a<-0.5
plot(φ/pi,20*log10(Mod(H(z,a))),type="l",ylim=c(-10,10),lwd=3,col="forestgreen",main="Modulová charakteristika pro a=0,5",xlab="Normalizovná f [rad×π]",ylab="Magnituda [dB]") # graf v dB pro z od 0 do pi

# Modulová charakteristika pro a=1,5
a<-1.5
plot(φ/pi,20*log10(Mod(H(z,a))),type="l",ylim=c(-10,10),lwd=3,col="forestgreen",main="Modulová charakteristika pro a=1,5",xlab="Normalizovná f [rad×π]",ylab="Magnituda [dB]") # graf v dB pro z od 0 do pi


##############
# Příklad 10 #
##############

# Definice systému pomocí obrazové přenosové funkce
H<-function(z) {
  return((1-2*z^-2+z^-3)/(1-z^-1))
}

φ<-seq(0,pi,length.out=500)
z<-cos(φ)+(0+1i)*sin(φ)

plot(φ/pi,20*log10(Mod(H(z))),type="l",ylim=c(0,10),lwd=3,col="forestgreen",xlab="Normalizovná f [rad×π]",ylab="Magnituda [dB]") # graf v dB pro z od 0 do pi

# Odvození diferenční rovnice je popsáno v textu
y[n]<-x[n]-2*x[n-2]+x[n-3]+y[n-1]

x<-c(0,0,0,rep(1,10)) # vstup x (jednotkový přechod počínající v čase 0)
y<-rep(0,4)           # počáteční podmínky pro y (hodnoty 0 pro časy -3, -2, -1 a 0)

for(i in 0:9) {
  n<-i+4                                # posun v čase o 4, abychom zamezili negativním indexům
  y[n]<-x[n]-2*x[n-2]+x[n-3]+y[n-1]
}

plot(-3:9,y,pch=19,main="Výstup systému\r\n y[n]=x[n]-x[n-1]\r\n se vstupem x = 1, 0, 0,…",xlab="Čas")
