% Problem 5

% define the basic parameters.
N=1024;
a=-0.5;
b=0.5;
sigma=1;

% Create an array of N equidistant points over [a,b].
x=linspace(a,b,N);

% Create a function.
y = exp(-x.^2/(2*sigma))/sqrt(2*pi)/sigma;

% 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.

fy = fft(y)/N;
fy(2:2:N)=-fy(2:2:N);

% Now, prepare the analytical Fourier coefficients you derived by
% hand.

xi = 0:(N-1);
c = exp(-2*(pi*sigma*xi).^2);

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

figure(1)
clf;
subplot(1,2,1);
plot(real(c(1:N/2)));
hold on
%plot(real(fy(1:N/2))*c(1)/real(fy(1)),'r.');
plot(real(fy(1:N/2)),'r.');
title('Real Part')
hold off

subplot(1,2,2);
plot(imag(c(1:N/2)));
hold on
%plot(imag(fy(1:N/2))*c(1)/real(fy(1)),'r.');
plot(imag(fy(1:N/2)),'r.');
title('Imaginary Part')
hold off

% Let's look at the more details around the origin.
figure(2)
clf;
subplot(1,2,1);
plot(c(1:N/16),'o')
hold on
%plot(real(fy(1:N/16))*c(1)/real(fy(1)),'r.');
plot(real(fy(1:N/16)),'r.');
title('Real Part')
hold off

subplot(1,2,2);
plot(imag(c(1:N/16)),'o');
hold on
%plot(imag(fy(1:N/16))*c(1)/real(fy(1)),'r.');
plot(imag(fy(1:N/16)),'r.');
title('Imaginary Part')
hold off

