%
%  FFT example
%
%- - - - - - - - - - - - - - - - - - - - - - - - - - - - 

%  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 = cos(xj) + 4*sin(2*xj);

%  FORWARD fourier transform:  compute complex fourier coeffs

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

aj = real(cj);
bj = imag(cj);

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

Yj = ifft(fftshift(cj)*N)

%  plot 

clf

subplot(2,1,1)
plot(xj,yj,'x')
hold on
plot(xj,real(Yj),'o')
title('\bf discrete function')
xlabel('\bf x-axis');  ylabel('\bf y(x) [x] & Y(x) [o]')
A1 = max(abs(yj));
axis([0 2*L -1.2*A1 1.2*A1])
text(3.5, -1,'original y(x): x')
text(3.5, -3,'reconstructed Y(x): o')

subplot(2,1,2)
plot(nj,aj,'o')
hold on
plot(nj,bj,'x')

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

text(-7,-1,'example:  y(x) = cos(x) + 4 sin(2x)')