---------------------------------
#      Fourierreihen als        #
#  lineares Regressionsproblem  #
---------------------------------


n = 1000     # devide [-pi,pi] into n or n-1 pieces
p = 20       # number of basis functions (sines and cosines)

x = seq(from=-pi,to=pi,length=n)    # Delta_x = 2*pi/(n-1)

#
# set up regressors:
#
X = outer(x,1:p,FUN="*")            # n times p matrix mit 
                                    # Eintraegen x[i]*j
Xsin = sin(X)
Xcos = cos(X)

# that's all to be done, now we can Fourier expand any y with
# lm(y ~ Xsin + Xcos),  this is super-easy!!


#
# example 1: y(x)=x, ungerade Funktion
#
y = x
plot(x,y,type="l")

res = lm(y ~ Xsin)
summary(res)
plot(x,res$fit,type="l")
lines(x,y,col="red")

plot(res$coeff)
plot(2/res$coeff[-1])
lines(1:20,col="red")
lines((-1):(-20),col="red")

res = lm(y ~ Xcos)
summary(res)                # alle betas 0 (oder e-16) 
                            # da ungerade Funktion
res = lm(y ~ Xsin + Xcos)
summary(res)

# let's see how the terms add up:
plot(x,y,type="l",col="red")
for(pp in 1:p)
{
  X = outer(1:pp,x,FUN="*")
  X = t(X)
  Xsin = sin(X)
  res = lm(y ~ Xsin)
  lines(x,res$fit)
  readline("press enter to continue..")
}

# die Bilder werden vielleicht uebersichtlicher, 
# wenn man jeweils nur eine Fourier-Summe plottet:
for(pp in 1:p)
{
  X = outer(1:pp,x,FUN="*")
  X = t(X)
  Xsin = sin(X)
  res = lm(y ~ Xsin)
  info = paste("p =",pp)
  plot(x,y,type="l",col="red",main=info)
  lines(x,res$fit)
  readline("press enter to continue..")
}



#
# example 2: ungerade Rechteck-Funktion:
#
y = ifelse(x>=0,1,-1)
plot(x,y,type="l")
for(pp in 1:p)
{
  X = outer(1:pp,x,FUN="*")
  X = t(X)
  Xsin = sin(X)
  res = lm(y ~ Xsin)
  info = paste("p =",pp)
  plot(x,y,type="l",col="red",ylim=c(-1.5,1.5),main=info)
  lines(x,res$fit)
  readline("press enter to continue..")
}


#
# example 3: y=|x|, gerade Funktion
#
y = abs(x)
plot(x,y,type="l")

res = lm(y ~ Xcos)
summary(res)
plot(x,res$fit,type="l")
lines(x,y,col="red")      # nearly perfect fit
plot(res$coeff)

# let's see how the terms add up:
for(pp in 1:p)
{
  X = outer(1:pp,x,FUN="*")
  X = t(X)
  Xcos = cos(X)
  res = lm(y ~ Xcos)
  info = paste("p =",pp)
  plot(x,y,type="l",col="red",main=info)
  lines(x,res$fit)
  readline("press enter to continue..")
}


#
# example 4: y=exp(x), weder gerade noch ungerade 
#                    -> we need sines and cosines
#
y = exp(x)
plot(x,y,type="l",col="red")
res = lm(y ~ Xsin + Xcos)
summary(res)
lines(x,res$fit,type="l")
plot(res$coeff)            # plottet erst die Konstante,
                           # dann die sin-Koeffizienten, 
                           # dann die cos-Koeffizienten
# let's see how the terms add up:
for(pp in 1:p)
{
  X = outer(1:pp,x,FUN="*")
  X = t(X)
  Xsin = sin(X)
  Xcos = cos(X)
  res = lm(y ~ Xsin + Xcos)
  info = paste("p =",pp)
  plot(x,y,type="l",col="red",ylim=c(-5,25),main=info)
  lines(x,res$fit)
  readline("press enter to continue..")
}


# alles in einem Bild:
plot(x,y,type="l",col="red",ylim=c(-5,25))
for(pp in 1:p)
{
  X = outer(1:pp,x,FUN="*")
  X = t(X)
  Xsin = sin(X)
  Xcos = cos(X)
  res = lm(y ~ Xsin + Xcos)
  lines(x,res$fit)
}
# die letzte Fourier-Summe in gruen und etwas
# fetter zum bestehenden plot hinzufuegen:
points(x,res$fit, pch=".", cex=4, col="green")

# pch = point character
# cex = character expansion, scheint wie folgt
# zu wirken: 
#  geplotteteGroesse = cex * defaultGroesse




---------------------------------------------
#   Noch ein paar Infos zu den Parametern   #
#               pch und cex                 #
---------------------------------------------

plot(1:10, pch=".")
plot(1:10, pch=".", cex=5)
plot(1:10, pch=".", cex=10)
# "." scheint offensichtlich ein Quadrat zu 
# sein, kein Kreis. Schauen wir uns mal einige 
# plot characters an:

for(i in 1:30)
{
  info = paste("pch =",i)            
  plot(1:10,pch=i,main=info)
  readline("press enter to continue..")
}
# also: pch=20 sind kleine ausgefuellte Kreise,
# pch=21 sind kleine unausgefuellte Kreise,
# kann man die mit cex beliebig vergroessern 
# und verkleinern?

# vergroessern:
for(scale in 1:20)
{
  info = paste("cex =",scale)            
  plot(1:10,pch=20,cex=scale,main=info)
  readline("press enter to continue..")
}

# verkleinern:
for(scale in 1/(1:10) )
{
  info = paste("cex =",scale)            
  plot(1:10,pch=20,cex=scale,main=info)
  readline("press enter to continue..")
}

# fuer den default plot character:
# vergroessern:
for(scale in 1:20)
{
  info = paste("cex =",scale)            
  plot(1:10,cex=scale,main=info)
  readline("press enter to continue..")
}
# verkleinern:
for(scale in 1/(1:10) )
{
  info = paste("cex =",scale)            
  plot(1:10,cex=scale,main=info)
  readline("press enter to continue..")
}