%
%  hw03a.m -- fourier cosine series
%
%- - - - - - - - - - - - - - - - - - - - - - - - - - - - 

%  number of points & interval

N = 32;  L = pi;  h = 2*L/N;

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

xj = 0:h:2*L-h;
yj = pi - abs(pi-xj);

%  FORWARD fourier transform:  compute complex fourier coeffs

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

%  fourier cosine/sine coefficients

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 cosine coefficients (from graphs)

ax = zeros(size(aj));
for k = 1:2:N
  ax(k)   = 0;
  ax(k+1) = -(4/pi)/(nj(k+1)^2);
end
ax(N/2+1) = pi/2;

%  plot 

figure(1);clf

subplot(2,1,1)
plot(xj,yj,'x')
hold on
plot(xj,real(Yj),'o')
title('\bf zigzag 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, 1.0,'original y(x): x')
text(2, 0.5,'reconstructed Y(x): o')

subplot(2,1,2)
plot(nj,aj,'o')
hold on
plot(nj,ax,'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 -1.2*A2 1.2*A2])

figure(2);

error = aj-ax;

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

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