% Problem 5

% define the basic parameters.
N=1024;
a=-0.5;
b=0.5;
sigma1=1;
sigma2=0.1;
sigma3=0.01;

% Create an array of N equidistant points over [a,b].
% Trying to exclude the point x=b=0.5 from the samples.
x=linspace(a,b,N+1);
x=x(1:end-1);

% Create a function.
y1 = exp(-x.^2/(2*sigma1));
y2 = exp(-x.^2/(2*sigma2));
y3 = exp(-x.^2/(2*sigma3));

% Normalize the function to have a unit L^2 norm.
a1 = 1/norm(y1);
a2 = 1/norm(y2);
a3 = 1/norm(y3);
y1 = a1*y1;
y2 = a2*y2;
y3 = a3*y3;

% Do the fft to approximate the Fourier series coefficients over this
% interval.  Note that we need to have 1/N here.  You need to go back
% to the original definition of the Fourier coefficients and its
% approximation by the trapezoidal rule.
% Note that fft essentially view the input data is defined over the
% interval on [0,1], instead of [-1/2,1/2].  So you need to do either
% of the following two:
% 1) Apply fftshift to the input vector before taking fft; or
% 2) Apply the complex exponential factor exp(pi*k)=(-1)^k to the output
%   of fft, which is equivalent to changing the signum of the fft results
%   alternatively as fy(2:2:N)=-fy(2:2:N), where fy=fft(y)/N.

fy1 = fft(y1)/N;
fy1(2:2:N)=-fy1(2:2:N);
fy2 = fft(y2)/N;
fy2(2:2:N)=-fy2(2:2:N);
fy3 = fft(y3)/N;
fy3(2:2:N)=-fy3(2:2:N);

% Now, prepare the analytical Fourier coefficients you derived by
% hand.  The length of the original spatial domain is 1, thus,
% the frequency grid spacing is also 1 (recall the reciprocity relation).

xi = 0:(N-1);
c1 = a1*sqrt(2*pi)*sigma1*exp(-2*(pi*sigma1*xi).^2);
c2 = a2*sqrt(2*pi)*sigma2*exp(-2*(pi*sigma2*xi).^2);
c3 = a3*sqrt(2*pi)*sigma3*exp(-2*(pi*sigma3*xi).^2);

% Now plot real and imaginary part separately using the semilog plot.

figure(1)
clf;
plot(real(c1(1:N/2)),'-');
grid
hold on
plot(real(c2(1:N/2)),'r-');
plot(real(c3(1:N/2)),'g-');
%plot(real(fy(1:N/2))*c(1)/real(fy(1)),'r.');
plot(real(fy1(1:N/2)),'b.-');
plot(real(fy2(1:N/2)),'r.-');
plot(real(fy3(1:N/2)),'g.-');
title('Real Part')
%hold off

figure(2)
clf;
plot(imag(c1(1:N/2)),'-');
grid
hold on
plot(imag(c2(1:N/2)),'r-');
plot(imag(c3(1:N/2)),'g-');
%plot(imag(fy(1:N/2))*c(1)/imag(fy(1)),'r.');
plot(imag(fy1(1:N/2)),'b.-');
plot(imag(fy2(1:N/2)),'r.-');
plot(imag(fy3(1:N/2)),'g.-');
title('Imaginary Part')
hold off

% Let's look at the more details around the origin.
figure(3)
clf;
plot(real(c1(1:N/16)),'-');
grid
hold on
plot(real(c2(1:N/16)),'r-');
plot(real(c3(1:N/16)),'g-');
%plot(real(fy(1:N/16))*c(1)/real(fy(1)),'r.');
plot(real(fy1(1:N/16)),'b.-');
plot(real(fy2(1:N/16)),'r.-');
plot(real(fy3(1:N/16)),'g.-');
title('Real Part')
%hold off

figure(4)
clf;
plot(imag(c1(1:N/16)),'-');
grid
hold on
plot(imag(c2(1:N/16)),'r-');
plot(imag(c3(1:N/16)),'g-');
%plot(imag(fy(1:N/16))*c(1)/imag(fy(1)),'r.');
plot(imag(fy1(1:N/16)),'b.-');
plot(imag(fy2(1:N/16)),'r.-');
plot(imag(fy3(1:N/16)),'g.-');
title('Imaginary Part')
hold off