%
%  lect10.m -- fourier series
%
%- - - - - - - - - - - - - - - - - - - - - - - - - - - - 

%  number of points & interval

N = 14;  L = pi;  h = 2*L/N;  K = 1;

%  define real periodic function on x=0->2L at N points

xj = 0:h:2*L-h;
yj = exp(K*cos(xj));

%  FORWARD fourier transform:  compute complex fourier coeffs

nj = -N/2:1:N/2-1;
cj = fftshift(fft(yj)/N);

%  fourier cosine/sine coefficients (NOTE ERROR IN MY NOTES!!)

aj = 2*real(cj);  aj(N/2+1) = real(cj(N/2+1));
bj = 2*imag(cj);  bj(N/2+1) = 0;

%  INVERSE fourier transform:  reconstruct N points of function 
%         from fourier coeffs

Yj = ifft(fftshift(cj)*N);

%  EXACT fourier coefficients (very mysterious)

cx = besseli(abs(nj),K);
error = real(cj)-cx;

%  plot 

figure(1);clf

subplot(2,1,1)
plot(xj,yj,'x')
hold on
plot(xj,real(Yj),'o')
title('\bf a periodic function')
xlabel('\bf x-axis');  ylabel('\bf y(x) [x] & Y(x) [o]')
A1 = max(abs(yj));
axis([0 2*L -0.2 1.2*A1])
text(2, 2.0,'[x] original y(x) = e^{cos x}')
text(2, 1.5,'[o] reconstructed Y(x)')

subplot(2,1,2)
plot(nj,real(cj),'o')
hold on
plot(nj,cx,'x')

title('\bf fourier series coefficients')
xlabel('\bf n-axis');  ylabel('\bf a_n:  fft [o] &  exact [x]')
A2 = max(abs(cj));
axis([-N/2 N/2 -0.2 1.2*A2])

figure(2);

plot(nj,error,'rx')
hold on
title('\bf coefficient errors')
xlabel('\bf n-axis');  ylabel('\bf error:  c^{fft}_n-c^{exact}_n [x]')

disp(' ')
disp(['error for c_0 = ' num2str(error(N/2+1))])
disp(' ')