#Fourier basis  
nb=45


n=2002    #size of the equidistant one dimensional grid
h=1/n #increment size
t=matrix(seq(h,1,by=h),ncol=1) #grid

Fbase=matrix(0,ncol=n,nrow=nb) #matrix of basis values 
                                #basis functions are in rows


Sn=matrix(2*(1:floor(nb/2)),ncol=1)
Cn=matrix(2*(1:ceiling(nb/2))-1,ncol=1)

Fbase[Sn,]=sqrt(2)*sin(pi*Sn%*%t(t)) #sin in basis 
Fbase[Cn,]=sqrt(2)*cos(pi*(Cn-1)%*%t(t)) #cos in basis 
Fbase[1,]=Fbase[1,]/sqrt(2)

#Checking that the elements of the basis are orthonormal (computing inner products)
IP=Fbase%*%t(Fbase)*h
View(round(IP,2))

#Test function and its plot
f=t(8*exp(-5000*(t-1/2)^2)+sin(30*t))
plot(t,f,type="l")

snf=f%*%t(f)*h

#Computation of coefficients

Cf=f%*%t(Fbase)*h 
Cf[1]

#Spectrum of the signal
freq1=1/(2*pi*Sn)
freq2=1/(2*pi*(Cn-1))
Sp1=Cf[Sn]
Sp2=Cf[Cn]

plot(freq2,Sp2^2,type='l')
lines(freq1,Sp1^2,type='l', col='red')

#The approximation error
Err=(f-Cf%*%Fbase)

plot(t,Err,type='l')

#The squared norm of the error
sne=Err%*%t(Err)*h

#Pythagorean theorem

snf

sne

Cf%*%t(Cf)

snf-Cf%*%t(Cf)

#The squared error norm as a function of the number of the elements of the basis

EE=vector('numeric', nb)
 

for( k in 1:nb)
  {
  EE[k]=snf-sum(Cf[1,1:k]^2)
  }

plot(1:nb,EE,type='l')
abline(0,0)

#From the graph we can infer some exponential exp(-l*n) decline of the error, where n is the sample size. 
#If we know l then exp(-l*N)=exp(-l*45)/2 and from that N=log(2)/l + 45.
#To get estimate of l we draw the previous graphs for log of EE
plot(1:nb,log(EE),type='l')

#The at hoc estimator would be 0.1. Thus


#So one can try with nb=52.

nb=52

#The error is now 0.0284 which is roughly half of the previous one. 


#Implementation of Gram-Schmidt orthonormalization
gso=function(A,H)
#Input A -- the matrix of row vectors representing (in a certain basis) the vectors (functions) to be normalized
#      H -- the matrix of the inner products of the vectors represented in A 
#Output B -- the matrix of row vectors representing in the same basis as in A the orthonormalized vectors
{
  nb=dim(H)[1]
  k=dim(A)[1]
  HH=H  
  AA=A
  B=A #the k x nb, k>=nb matrix output with columns being orthonormalized columns of A
  i=1
  B[,1]=AA[,1, drop=F]/HH[1,1] #the first output vector normalized
  while(nb>1){
    H1=HH[1,,drop=F]            # row 1 x nb
    AA=AA-(H1%x%AA[,1, drop=F])/HH[1,1]          # k x nb
    HH=HH-t(H1)%*%H1/HH[1,1]
    HH=HH[2:nb,2:nb, drop=F]
    AA=AA[,2:nb,drop=F]
    i=i+1
    B[,i]=AA[,1,drop=F]/sqrt(HH[1,1]) #subsequent column to the output
    nb=nb-1
  }
  
  #G-S orthonormalization
  nb=10
  A=diag(nb)
  
  I=t(rep(1,nb))
  S=0:(nb-1)
  H=1/(I%x%S+1/2+t(I%x%S+1/2))
  S=(sqrt(2*S+1))%*%t(sqrt(2*S+1))
  H=H*S
  
  B=gso(A,H)
  
  #Plotting orthonormal polynomial basis
  
  PBase=matrix(1,ncol=nb,nrow=n)
  PBase=Mono[,1:nb]%*%B
  Max=max(PBase[,1:nb])
  Min=min(PBase[,1:nb])
  plot(t,PBase[,1],type='l',ylim=c(Min,Max),ylab='')
  for(i in 2:10)
  {
    lines(t,PBase[,i],type='l',col=i)
  }
  
  B
}


system.file('scripts', package='fda')
# [1] "/Users/mats-ksp/Library/R/3.4/library/fda/scripts"
dir( "/Users/mats-ksp/Library/R/3.4/library/fda/scripts")
# [1] "afda-ch01.R"  "afda-ch02.R"  "afda-ch03.R"  "afda-ch04.R"  "afda-ch05.R"  "afda-ch06.R" 
# [7] "afda-ch07.R"  "fda-ch01.R"   "fda-ch02.R"   "fda-ch03.R"   "fda-ch07.R"   "fda-ch08.R"  
#[13] "fda-ch13.R"   "fda-ch17.R"   "fdarm-ch01.R" "fdarm-ch02.R" "fdarm-ch03.R" "fdarm-ch04.R"
#[19] "fdarm-ch05.R" "fdarm-ch06.R" "fdarm-ch07.R" "fdarm-ch08.R" "fdarm-ch09.R" "fdarm-ch10.R"
#[25] "fdarm-ch11.R"


daybasis65 <- create.fourier.basis(rangeval=c(0, 365), nbasis=65)

#  -----------  set up the harmonic acceleration operator  ----------

harmaccelLfd365 <- vec2Lfd(c(0,(2*pi/365)^2,0), c(0, 365))

#  ---------  create fd objects for temp. and prec. ---------------

# First check the distribution
qqnorm(CanadianWeather$dailyAv[,,"Temperature.C"], datax=TRUE)
# Consistent with a strong annual cycle
# plus weaker normal noise

daytempfd <- with(CanadianWeather, smooth.basis(day.5,
                                                dailyAv[,,"Temperature.C"],
                                                daybasis65, fdnames=list("Day", "Station", "Deg C"))$fd )
plot(daytempfd, axes=FALSE)
axisIntervals(1)
axis(2)

a=CanadianWeather$dailyAv[,,"Temperature.C"]
plot(a[,1],type="l")



#Fourier basis  
nb=45


n=2002    #size of the equidistant one dimensional grid
h=1/n #increment size
t=matrix(seq(h,1,by=h),ncol=1) #grid

Fbase=matrix(0,ncol=n,nrow=nb) #matrix of basis values

Sn=matrix(2*(1:floor(nb/2)),ncol=1)
Cn=matrix(2*(1:ceiling(nb/2))-1,ncol=1)

Fbase[Sn,]=sqrt(2)*sin(pi*Sn%*%t(t)) #sin in basis 
Fbase[Cn,]=sqrt(2)*cos(pi*(Cn-1)%*%t(t)) #cos in basis 
Fbase[1,]=Fbase[1,]/sqrt(2)



plot(t,Fbase[2,],type="l", col=2,ylab="") 
lines(t,Fbase[1,],type="l", col="red") 
for(i in 3:9)
{
lines(t,Fbase[i,],type="l", col=i) 
}



#Test function and its plot
t=seq(0,1,by=0.01)
f=t(8*exp(-5000*(t-1/2)^2)+sin(30*t))
plot(t,f,type="l")

#Computation of coefficients

Cf=f%*%t(Fbase)*h 

#Computation of expansion (projections)
Pf=Cf%*%Fbase

#Plot of the projection
lines(t,Pf,type="l",col='red')


#Plotting monomial basis
nb=45
n=2002
h=1/n
t=matrix(seq(h,1,by=h),ncol=1) #grid 

n=dim(t)[1]

i=0
Mono=matrix(1,ncol=nb,nrow=n)
Mono[,1]=sqrt(2*i+1)*matrix(1,ncol=1,nrow=n)

for(i in 1:(nb-1))
{
  Mono[,i+1]=(((2*i+1)^(1/(2*i)))*t)^i
}
Max=max(Mono[,1:15])
Min=min(Mono[,1:15])
plot(t,Mono[,1],type='l',ylim=c(Min,Max),ylab='')
for(i in 2:15)
{
  lines(t,Mono[,i],type='l',col=i)
}


#Implementation of Gram-Schmidt orthonormalization
gso=function(A,H)
{
  nb=dim(H)[1]
  k=dim(A)[1]
  HH=H  
  AA=A
  B=A #the k x nb, k>=nb matrix output with columns being orthonormalized columns of A
  i=1
  B[,1]=AA[,1, drop=F]/HH[1,1] #the first output vector normalized
  while(nb>1){
    H1=HH[1,,drop=F]            # row 1 x nb
    AA=AA-(H1%x%AA[,1, drop=F])/HH[1,1]          # k x nb
    HH=HH-t(H1)%*%H1/HH[1,1]
    HH=HH[2:nb,2:nb, drop=F]
    AA=AA[,2:nb,drop=F]
    i=i+1
    B[,i]=AA[,1,drop=F]/sqrt(HH[1,1]) #subsequent column to the output
    nb=nb-1
    }
B
}


#G-S orthonormalization
nb=10
A=diag(nb)

I=t(rep(1,nb))
S=0:(nb-1)
H=1/(I%x%S+1/2+t(I%x%S+1/2))
S=(sqrt(2*S+1))%*%t(sqrt(2*S+1))
H=H*S

B=gso(A,H)

#Plotting orthonormal polynomial basis

PBase=matrix(1,ncol=nb,nrow=n)
PBase=Mono[,1:nb]%*%B
Max=max(PBase[,1:nb])
Min=min(PBase[,1:nb])
plot(t,PBase[,1],type='l',ylim=c(Min,Max),ylab='')
for(i in 2:10)
{
  lines(t,PBase[,i],type='l',col=i)
}


#The test function
plot(t,f,type="l")

#Computation of coefficients

Cf=f%*%t(PBase)*h 

#Computation of expansion (projections)
Pf=Cf%*%PBase

#Plot of the projection
lines(t,Pf,type="l",col='red')

#Plotting polynomial non-orthonormal basis
nb=45
Poly=matrix(1,ncol=nb,nrow=n) 
Poly[,1]=matrix(1,ncol=1,nrow=n)
Poly[,2]=t-1/2
Poly[,2]=Poly[,2]/sqrt(sum(Poly[,2]^2)*h)
for(k in 3:nb) {
  for(i in 0:(k-2))
    {
    Poly[ ,k] =Poly[,k]*(2*sqrt(exp(1)*(k-1))*(t-i/(k-2)))
  }
  Poly[,k]=Poly[,k]*sqrt(2)/(2*exp(1))
  Poly[,k]=Poly[,k]/sqrt(sum(Poly[,k]^2)*h)
}
Max=max(Poly[ ,1:9])
Min=min(Poly[ ,1:9]) 
plot(t,Poly[,1],type='l',ylim=c(Min,Max),ylab='') 

for ( i in 2:9)
{
  lines(t,Poly[,i],type='l',col=i)
}


#G-S orthonormalization

#Computing the inner product matrix 
nb=45
A=diag(nb)

H=t(Poly)%*%Poly*h

B=gso(A,H)

#Plotting orthonormal polynomial basis

PBase=matrix(1,ncol=nb,nrow=n)
PBase=Poly%*%B
Max=max(PBase[,1:16])
Min=min(PBase[,1:16])
plot(t,PBase[,1],type='l',ylim=c(Min,Max),ylab='')
for(i in 2:16)
{
  lines(t,PBase[,i],type='l',col=i)
}


#The test function
plot(t,f,type="l")

#Computation of coefficients

Cf=f%*%t(PBase)*h 

#Computation of expansion (projections)
Pf=Cf%*%PBase

#Plot of the projection
lines(t,Pf,type="l",col='red')


install.packages("wavelets")
#Test function and its plot
f=t(8*exp(-5000*(t-1/2)^2)+sin(30*t))
plot(t,f,type="l")


WD=dwt(f,filter="haar")
plot(WD)