%% Fourier Series

%%
% You can use the following commands to calculate the $n\mbox{th}$ partial sum of 
% the Fourier series of the expression $f$ on the interval $[-L,L]$.


syms x k L n
evalin(symengine,'assume(k,Type::Integer)');
a = @(f,x,k,L) int(f*cos(k*sym('pi')*x/L)/L,x,-L,L);
b = @(f,x,k,L) int(f*sin(k*sym('pi')*x/L)/L,x,-L,L);
fs = @(f,x,n,L) a(f,x,0,L)/2 + ... 
    symsum(a(f,x,k,L)*cos(k*pi*x/L) + b(f,x,k,L)*sin(k*pi*x/L),k,1,n);
pfs = @(f,x,n,L) pretty(simple(fs(f,x,n,L)));


%%
% I'll do it for the function $f(x)=|x|$ on the interval $[-1,1]$.
f = abs(x)

%%
% The $n\mbox{th}$ partial sum is 

fs(f,x,n,1)

%%
% Here's a more readable form.

pfs(f,x,n,1)


%%
% The Fourier series is
pfs(f,x,inf,1)


%%
% Here are the plots of the partial sums for $n=2,5,10$.  The plot also shows
% the function $f$.

ezplot(fs(f,x,2,1),-1,1)
hold on
ezplot(f,-1,1)
hold off
title('Partial sum with n=2')

%%

ezplot(fs(f,x,5,1),-1,1)
hold on
ezplot(f,-1,1)
hold off
title('Partial sum with n=5')

%%

ezplot(fs(f,x,10,1),-1,1)
hold on
ezplot(f,-1,1)
hold off
title('Partial sum with n=10')

%%
% Now I do it with the function $f(x) = x$ on $[-1,1]$.

f = x

fs(f,x,n,1)
%%
pfs(f,x,n,1)

%%
% Here are plots of the partial sums for $n=2,5,10,20,50$.

ezplot(fs(f,x,2,1),-1,1)
hold on
ezplot(f,-1,1)
hold off
title('Partial sum with n=2')

%%

ezplot(fs(f,x,5,1),-1,1)
hold on
ezplot(f,-1,1)
hold off
title('Partial sum with n=5')

%%

ezplot(fs(f,x,10,1),-1,1)
hold on
ezplot(f,-1,1)
hold off
title('Partial sum with n=10')


%%

ezplot(fs(f,x,20,1),-1,1)
hold on
ezplot(f,-1,1)
hold off
title('Partial sum with n=20')

%%

ezplot(fs(f,x,50,1),-1,1)
hold on
ezplot(f,-1,1)
hold off
title('Partial sum with n=50')

%%
% If you zoom in here, you notice that the graph of the partial sum of the 
% Fourier series looks very jagged.  That is because *ezplot* does not plot 
% enough points compared to the number of oscillations in the functions in 
% the partial sum.  We can fix this problem using plot.  This also allows 
% us to use a different colors for the plot of the original function and 
% the plot of the partial sum.  To use plot, we need to turn the partial 
% sum into an inline vectorized function and specify the points where it 
% will be evaluated.

g = inline(vectorize(fs(f,x,50,1)));
X = -1:.001:1;
plot(X,g(X),'r')
hold on
ezplot(f,-1,1)
hold off
title('Partial sum with n=50')

%%

g = inline(vectorize(fs(f,x,100,1)));
X = -1:.0005:1;
plot(X,g(X),'r')
hold on
ezplot(f,-1,1)
hold off
title('Partial sum with n=100')

%%
% Notice that no matter how many terms we take, the partial sum always
% overshoots the function $x$ at the end points of the interval.  If we do
% our two plots on the interval $[-2,2]$, we see that outside $[-1,1]$, the
% partial sum doesn't have much to do with the function.  

ezplot(fs(f,x,20,1),-2,2)
hold on
ezplot(f,-2,2)
hold off

%%
% A Fourier series on $[-L,L]$ is $2L$ periodic, and so are all its partial
% sums.  So, what we are really doing when we compute the Fourier series of
% a function $f$ on the interval $(-L,L)$ is computing the Fourier series of
% a $2L$ periodic extension of f.  To do that in MATLAB, we have to make
% use of the unit step function u(x), which is 0 if $x < 0$ and 1 if $x \ge 0$.
% It is known as the Heaviside function, and the MATLAB command for it is
% *heaviside*.  The function $h(x) =  u(x-a)u(b-x)$ is 1 on the interval 
% [a,b] and 0 outside the interval. So the formula

 
f = heaviside(x+3)*heaviside(-1-x)*(x+2) + heaviside(x+1)*heaviside(1-x)*x ...
    + heaviside(x-1)*heaviside(3-x)*(x-2);

%%
% extends $f(x) = x$ to be periodic on $[-3,3]$, with period 2.  To check that 
% we've extended it correctly, we plot it.

ezplot(f,-3,3)

%%
% Here is a plot of the function and its Fourier series.

ezplot(fs(f,x,20,1),-3,3)
hold on
ezplot(f,-3,3)
hold off

%%
% It is somewhat clearer if the plots aren't the same color, so I'll use
% *plot* for the partial sum.

X = -3:0.001:3;
g = inline(vectorize(fs(f,x,20,1)));

ezplot(f,-3,3)
hold on
plot(X,g(X),'r')
hold off

%%
% The overshoots are at the discontinuities of the 2 periodic extension of
% f.  This is called the *Gibbs Phenomenon.*