#############################################
####     SUPERPOSICION DE SINUSOIDES     ####
#############################################

# frecuencia de muestreo
fs <- 441000
# frecuencia fundamental
f0 = 440
# período
T <- 1/f0
# vector de tiempo
t <- seq(1/fs, 1, length.out=fs)
# frecuencia angular
w0 <- 2*pi/T
# coseno de frecuencia fundamental (440 Hz)
h0 <- 0.3*cos(w0*t)
# primer armónico (880 Hz)
h1 <- 2*cos(2*w0*t)
# segundo armónico (1320 Hz)
h2 <- 3*sin(3*w0*t)
# tercer armónico (1760 Hz)
h3 <- 0.5*cos(4*w0*t)
# superposición
s <- 0.5 + h0 + h1 + h2 + h3

# período T en muestras
Tm <- round(T*fs)
# muestras en un período
dt <- 1:Tm

# grafica de un periodo
# plot(h0[1:(Tm)], type='l')
plot(s[1:(Tm)], type='l')

#############################################
####         ANALISIS DE FOURIER         ####
#############################################

# cálculo de A0
A0 <- mean(s)

# cálculo de coefficientes de Fourier
# cálculo de An 
coefA <- function(n) {
  2/Tm * sum(s[dt] * cos((2*pi*n/Tm)*dt))
}
# cálculo de Bn
coefB <- function(n) {
  2/Tm * sum(s[dt] * sin((2*pi*n/Tm)*dt))
}

# sapplay para repetir el cálculo
A <- sapply(1:5, FUN=coefA)
B <- sapply(1:5, FUN=coefB)
# data frame con los coeficientes obtenidos
coefs <- data.frame(A=round(A,1), B=round(B,1))

#############################################
####         SINTESIS DE FOURIER         ####
#############################################

# reconstrucción: suma de componentes sinusoidales
res <- matrix(rep(NA, fs*5), nrow=5)
for(n in 1:5){
  res[n, ] <- A[n]*cos(2*pi*n*t/T) +
    B[n]*sin(2*pi*n*t/T)
}
s.rec <- A0 + colSums(res)

# gráfica de señal original y reconstrucción
plot(s[1:(Tm)], type='l', col=1)
par(new=T)
plot(s.rec[1:(Tm)], type='l', col=2)