Edit the function in fun, bound, and the value of n, and plot the fourier series for fun on the interval [-bound,bound] with n terms. Used with assignment 1 in A First Course in Wavelets with Fourier Analysis, http://www.math.tamu.edu/~francis.narcowich/m414/s16/m414s16_hwc.html.

fourier_series.R

#Given a function, calculates the
#Fourier series up to n terms

#Clear previous variables
rm(list=ls())

#For complex integration
library(elliptic)

#The function to approximate
fun = function(x){
    
    #9
    #return(cos(x))
    
    #10a
    #Use the fact that $f_e is even, and f(-x)=f(x)
    #https://stat.ethz.ch/pipermail/r-help/2005-July/076128.html
    #x = ifelse(x<0,-x,x)
    #return((pi-x)%%(2*pi))
    
    #10b
    #Use the fact that $f_o is odd, and -f(-x)=f(x)
    #m = ifelse(x<0,-1,1)
    #x = m*x
    #return(m*((pi-x)%%(2*pi)))
    
    #10c
    #x = ifelse(x<0,-x,x)
    #return((pi-x)%%(pi))
    
    #11a
    #gig = (x<=pi)&(x>=-pi)
    #x = ifelse(gig,x,(x+pi)%%(2*pi)-pi)
    #return(exp(-x/3))
    
    #11b
    #gig = (x<=2*pi)&(x>=0)
    #x = ifelse(gig,x,x%%(2*pi))
    #return(exp(-x/3))   
    
    #8
    #gig = (x > -1/2) & (x <= 1/2)
    #em =  ((x > -1) & (x <= -1/2)) | (( x > 1/2) & (x <= 1))
    #y = ifelse(gig,1,ifelse(em,0,NaN))
    #return(y)
    
    #7
    #return(abs(sin(x)))
    
    #21
    #gig = (x<=pi)&(x>=-pi)
    #x = ifelse(gig,x,x%%(2*pi))
    #return(x)
    
    #23a
    #gig = (x<=1)&(x>=-1)
    #x = ifelse(gig,x,(x+1)%%2-1)
    #ifelse(x!=0,return(exp(x)),return(exp(1)/2))
    
    #23e
    #return(cos(x)+abs(cos(x)))
    
    #23f
    gig = (x<=pi)&(x>=-pi)
    x = ifelse(gig,x,(x+pi)%%(2*pi)-pi)
    return(ifelse(x!=0,sin(x)/x,1))
    
    #return((1/12)*(x^3-pi^2*x))
}

#Upper bound (lower bound is -bound)
bound = pi

#a_0
a_0 = 1/(2*bound)*integrate(Vectorize(fun),-bound,bound)$value
#Half-interval
#a_0 = 1/(bound)*integrate(Vectorize(fun),0,bound)$value

#The number of repitions
n = 20 

#Vectors for the coefficients
a = rep(0,n)
b = rep(0,n)
c = rep(0,n)

#Calculate the coefficients
for(k in 1:n){
    a[k]=1/bound * integrate(function(x){return(fun(x)*cos(k*pi*x/bound))},-bound,bound)$value
    b[k]=1/bound * integrate(function(x){return(fun(x)*sin(k*pi*x/bound))},-bound,bound)$value
    
    #For half interval
    #a[k]=2/bound * integrate(function(x){return(fun(x)*cos(k*pi*x/bound))},0,bound)$value
    #b[k]=2/bound * integrate(function(x){return(fun(x)*sin(k*pi*x/bound))},0,bound)$value
    
    #Complex coefficients
    c[k] = 1/(2*bound) * myintegrate(function(x){return(fun(x)*exp(complex(imaginary=-k*x)))},-bound,bound)
}

#Fourier expansion
s_n = function(x){
    s = a_0
    for(i in 1:n)
        s=s+a[i]*cos(i*pi*x/bound)+b[i]*sin(i*pi*x/bound)
    return(s)
}

#http://stackoverflow.com/questions/7144118/how-to-save-a-plot-as-image-on-the-disk

#x-coordinates
mult = 3 #How wide should the plot be?
xc = seq(-mult*bound,mult*bound,by=0.01)

#plot fun
plot(xc,lapply(xc,fun),type="l",xlab="x",ylab="y",main="f(x)")

#Plot s_n
plot(xc,lapply(xc,s_n),type="l",xlab="x",ylab="s_n",main=paste("n=",n))