% Use matlab fft/ifft function to evaluate Fourier coefficients
%     c_k = sum_{j=0}^{2m-1} y_j exp(i*k*pi*j/m), k=0,...,2m-1.
%
% Notice that X=fft(x) computes
%     X_k = sum_{j=1}^N x_j exp((j-1)*(k-1)(-2*pi*i)/N, j=1,...,N
% and x=ifft(X) computes
%     x_k = 1/N*sum_{j=1}^N X_j exp((j-1)*(k-1)(2*pi*i)/N, j=1,...,N.
% We determine that the function ifft matches the form for c_k we need, but
% some index changes are needed to use this function to compute c_k. 
% We assume 2m=N.
%
% First we change the dummy index j to jp=j+1 in c_k
%     c_k = sum_{jp=1}^{N} y_{jp-1} exp(i*k*pi*(jp-1)/(N/2)), k=0,...,N-1.
% Since jp is just another dummy index, we can call it j
%     c_k = sum_{j=1}^{N} y_{j-1} exp(2*i*k*pi*(j-1)/N), k=0,...,N-1.
% On the other hand, we change k to k-1:
%     c_{k-1} = sum_{j=1}^{N} y_{j-1} exp(2*i*(k-1)*pi*(j-1)/N), k=1,...,N
% Now the right hand side matches
%     N*x_k, if y_{j-1} = X_j, j=1,...,N.
%
% This is fine with us since the formula in the text starts k=0, while in
% Matlab we usually start the index=1. We will just have to use c(k+1) for
% the coefficient c_k.
%
% Let's try an example:
%     f(x) = |x|, x in [-pi,pi]
%
% We have a total of 2*m or N data points in the array y and the function
% is assumed to be periodic so y(N+1)=y(1).
N=32; m=N/2;
x = linspace(-pi,pi,N+1);
y = abs(x);
ck = N*ifft(y,N);   %compute the comflex coefficient via ifft (inverse fft)
ak = real(ck.*repmat([1 -1],1,N/2))*2/N;  %the real part scaled to represent a_k
bk = imag(ck.*repmat([1 -1],1,N/2))*2/N;  %the imaginary part scaled to represent a_k

yapp = zeros(1,N+1); %approximation from Fourier series

% the following use the Fourier series to evaluate our approximations at
% x(j), j=1,...,N.
for j=1:N
    yapp(j) = 0.5*(ak(1) + ak(m+1)*cos(m*x(j)));
for k=2:m
    yapp(j) = yapp(j) + ak(k)*cos((k-1)*x(j)) + bk(k)*sin((k-1)*x(j));
end
end
yapp(N+1) = yapp(1);  % enforcing the periodic condition

plot(x,y,'+',x,yapp,'-');
xlabel('x'); ylabel('y'); title('FFT to generate Fourier series');
legend('data','approximation');
    


