%
%  fft314test.m  -- djm -- 28 jan 2008
%
%  RUN THIS TO LEARN ABOUT THE FFT ROUTINES
%

clear

%  interval parameters
a = -pi;  b = +pi;  N = 64;  dx = (b-a)/N;  Lx = (b-a)/2;

%  define x-points  (' means real-valued transpose)
xn = a + dx*(0:N-1)';

%  define real f-values  (!!! ELEMENTWISE ARITHMETIC !!!)
kk = 4;  fn = cos(kk*(pi/Lx)*xn) + (0.5)*cos(7*(pi/Lx)*xn);
%kk = 4;  fn = sin(kk*(pi/Lx)*xn);
%fn = 1./(1.1 - cos(xn-pi/8));

% % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
%
%  THE 3 FFT FUNCTION CALLS ARE HERE
%

%  complex fourier coeffs
[aj,bj,cj,kj] = fft314(xn,fn);

%  complex inverse fourier
[fnC] = ifftC314(xn,cj);

%  real inverse fourier
[fnR] = ifftR314(xn,aj,bj);

%
% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % 

%  check plot (include periodic points)
xnP = [xn ; b];  
fnP = [fn ; fn(1)];  fnCP = [fnC ; fnC(1)];  fnRP = [fnR ; fnR(1)];

figure(1);  clf
subplot(2,1,1)
plot(xnP,real(fnP)  ,'r:' ,xnP,real(fnP ),'k.');  hold on
plot(xnP,real(fnCP) ,'b--',xnP,real(fnRP),'g-.');
plot(b,fnP(end),'ko')

title(['\bf original f(x_n) - red & partial sum f_{N/2}(x_n) - blue'],'fontsize',14)
xlabel('\bf x-axis','fontsize',12)
ylabel('\bf f-axis','fontsize',12)

axis([a b min(real(fnC))-0.1 max(real(fnC))+0.1])

%  real fourier coefficients 
subplot(2,1,2)
plot(0:N/2-1,aj,'rx',0:N/2-1,bj,'bo');

title(['\bf fourier coefficients:  a_j (red x) & b_j (blue o)'],'fontsize',14)
xlabel('\bf n-axis','fontsize',12)
ylabel('\bf (a_j,b_j)-axis','fontsize',12)

axis([0 N/2-1 min([aj ; bj])-0.1 max([aj ; bj])+0.1])

%
%  because of the "nyquist" issue, fnC & fnR are not quite fn !!!
%

disp(' ')
disp('Errors:')
disp(['   max|fn -fnC| = ' num2str(max(abs(fn -fnC)))])
disp(['   max|fn -fnR| = ' num2str(max(abs(fn -fnR)))])
disp(['   max|fnC-fnR| = ' num2str(max(abs(fnC-fnR)))])