#--------------------------------#
#      Anwendungsbeispiel:       #
#    Fourierreihenentwicklung    #
#--------------------------------#


n = 1000
p = 20

x = seq(from=-pi,to=pi,length=n)
x


# Wir definieren zwei Funktionen, die dann jeweils  
# in eine Fourrierreihe entwickelt werden sollen:

# eine ungerade Funktion => nur Sinus-Terme
y = ifelse(x>=0,1,-1)
y

# eine gerade Funktion => nur Cosinus-Terme
y = ifelse(abs(x)<=pi/2,1,0)
y

plot(x,y,ylim=c(-1.5,1.5))


# Matrix der Regressoren:

X = cos(x)
for(j in 2:p)
{
  xcosj = cos(j*x)
  X = cbind( X , xcosj )
}

for(j in 1:p)
{
  xsinj = sin(j*x)
  X = cbind( X , xsinj )
}


# Damit haben wir schon alles upgesetted und koennen jetzt 
# sehr einfach die Fourierkoeffizienten bestimmen:

res = lm( y ~ X )

beta = res$coeff
beta

beta0 = beta[1]
a = beta[2:(p+1)]
b = beta[(p+2):(2*p+1)]


# Jetzt machen wir noch eine Visualisierung:

x0 = rep(1,n)
sumFR = beta0*x0 
for(j in 1:p)
{
  sumFR = sumFR + a[j]*cos(j*x) 
  info = paste("adding term a[j]*cos(j*x) with j =",j)
  plot(x,y,ylim=c(-1.5,1.5),main=info)
  points(x,sumFR,col="red")
  readline("press enter to continue..")

  sumFR = sumFR + b[j]*sin(j*x) 
  info = paste("adding term b[j]*sin(j*x) with j =",j)
  plot(x,y,ylim=c(-1.5,1.5),main=info)
  points(x,sumFR,col="red")
  readline("press enter to continue..")
}