#---------------------------------------------------------------------- 
# Risposta dell'RC (ai capi di R) per onda quadra
# di frequenza arbitraria rispetto a quella di taglio
# (e non solo il caso elementare di 'carica completa')
# mediante analisi di Fourier
#
# GdA, 15/4/2015
#----------------------------------------------------------------------

# onda quadra 'pari'
n = 101
A.oq.pari <- rep(0,n)
phi.oq.pari <- rep(0,n)    # sono nulli! (messa per generatiltà)
for (k in seq(1,n,by=2)) {
  A.oq.pari[k] <- 2/(pi*k)*sin(k*pi/2)
}

# metodo 'alternativo' => STESSO RISULTATO  (-> settare 'TRUE' per provare)
metodo.alternativo <- TRUE # vedi nota 19, p. 174 dispensa 
if (metodo.alternativo) {
  A.oq.pari <- rep(0,n)
  phi.oq.pari <- rep(0,n)   
  for (k in seq(1,n,by=2)) {
    A.oq.pari[k]   <- 2/(pi*k)
    phi.oq.pari[k] <- (k-1)/2*pi
  }
}

T.oq <- 1 # secondo (arbitrario...  - conta la scala nuT !! vedo dopo)
nu.oq <- 1/T.oq

t <- seq(-1*T.oq, +1*T.oq, by=T.oq/10000)    # molto fitti !!
nt <- length(t)

fk <- matrix(rep(0, n*nt), c(nt,n))
ft <- rep(0,nt) 

for(k in 1:n) {
  fk[,k] <- A.oq.pari[k]*cos(k*2*pi*nu.oq*t + phi.oq.pari[k])
  ft     <-  ft + fk[,k]
}

plot(t,ft, ty='l')

# uscita ai capi di un RC di frequenza di taglio nuT
# cambiare nuT a piacere (in unita di nu.oq!!)
nuT = 0.2*nu.oq  #        
fk.C <- matrix(rep(0, n*nt), c(nt,n))
ft.C <- rep(0,nt) 
for(k in 1:n) {
  fk.C[,k] <- A.oq.pari[k]*1/sqrt(1+(k*nu.oq/nuT)^2)*cos(k*2*pi*nu.oq*t + phi.oq.pari[k] - atan(k*nu.oq/nuT))
  ft.C     <- ft.C + fk.C[,k]
}

points(t,ft.C, ty='l', col='red')

# Fare altri plot di controllo!
#
# Si noti come, anche se con n=101 l'onda quadra non è ben riprodotta
# il segnale in uscita somiglia a quello visto all'oscilloscopio
# in quanto le alte frequenze, cruciali per riprodurre l'onda quadra,
# sono filtrate dall'RC e quindi il loro contributo all'output è irrilevante
# (già n=11 dà qualcosa di accettabile, mentre l'onda quadra è riprodotta malamente)

# Come esercizio usare lo stesso metodo per l'onda triangolare.