% Problem 3

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

% 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.
y = x; % In problem 4, this should be y=x.^2, of course.

% Normalize the function to have a unit L^2 norm.
ey = norm(y);
y = y/ey;

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

c = zeros(1,N);
% c(1) = 1/12.0; % for problem 4.
for k=1:N-1
  c(k+1)=i*(-1)^k/(2*pi*k); % c(k+1)=(-1)^k/(2*(pi*k)^2); % for problem 4.
end


% Normalize the coefficients
c = c/ey;

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

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

subplot(1,2,2);
plot(imag(c(1:N/2)));
grid
hold on
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(real(c(1:N/16)),'o');
grid
hold on
plot(real(fy(1:N/16)),'r.');
title('Real Part')
hold off

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