function main_fft

rng(0)

N  = 8;
W  = exp(-1i*2*pi/N);

F8 = DFT_matrix(8);
F4 = DFT_matrix(4);

fprintf(1,'|DFT_matrix(8) - fft(eye(8))| = %10.3e\n',max(max(abs(F8 - fft(eye(8))))))

xx = randn(8,1);

%%% Try out the recursive formula.
I_even  = 1:2:7;  % Note how Matlab starts indexing on "1", but we use "0".
I_odd   = 2:2:8;
yy_even = F4*xx(I_even);
yy_odd  = F4*xx(I_odd );
yy      = [yy_even(1) + W^0*yy_odd(1);...
           yy_even(2) + W^1*yy_odd(2);...
           yy_even(3) + W^2*yy_odd(3);...
           yy_even(4) + W^3*yy_odd(4);...
           yy_even(1) - W^0*yy_odd(1);...
           yy_even(2) - W^1*yy_odd(2);...
           yy_even(3) - W^2*yy_odd(3);...
           yy_even(4) - W^3*yy_odd(4)];
fprintf(1,'|yy - F8*xx| = %10.3e\n',max(abs(yy - F8*xx)))

%%% Write the recursive formula more compactly.
dd = exp(-1i*pi*linspace(0,1-2/N,N/2)');
yy      = [yy_even + dd.*yy_odd;...
           yy_even - dd.*yy_odd];
fprintf(1,'|yy - F8*xx| = %10.3e\n',max(abs(yy - F8*xx)))

%%% Try the other recursive formula.
yy = zeros(8,1);
yy(I_even) = F4*(    xx(1:4) +     xx(5:8));
yy(I_odd ) = F4*(dd.*xx(1:4) - dd.*xx(5:8));
fprintf(1,'|yy - F8*xx| = %10.3e\n',max(abs(yy - F8*xx)))

%%% Check the block factorizations
dd     = exp(-1i*pi*linspace(0,1-2/N,N/2)');
Fblock = [F4, diag(dd)*F4;...
          F4,-diag(dd)*F4];
fprintf(1,'|Fblock - F8(:,[I_even,I_odd])| = %10.3e\n',...
        max(max(abs(Fblock - F8(:,[I_even,I_odd])))))
Fblock = [F4,          F4;...
          F4*diag(dd),-F4*diag(dd)];
fprintf(1,'|Fblock - F8([I_even,I_odd],:)| = %10.3e\n',...
        max(max(abs(Fblock - F8([I_even,I_odd],:)))))

xx = randn(8,1);
fprintf(1,'|myFFT(xx)    - fft(xx)| = %10.3e\n',max(max(abs(myFFT(xx) - fft(xx)))))
fprintf(1,'|myFFTrev(xx) - fft(xx)| = %10.3e\n',max(max(abs(myFFTrev(xx) - fft(xx)))))

    keyboard

return

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

function F = DFT_matrix(N)

W  = exp(-1i*2*pi/N);
ee = (0:(N-1))'*(0:(N-1));
F  = W.^ee;

return

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

function yy = myFFT(xx)

N       = numel(xx);
xx_even = xx(1:2:end);
xx_odd  = xx(2:2:end);

if (N>=8)
  yy_even = myFFT(xx_even);
  yy_odd  = myFFT(xx_odd);
  dd      = exp(-1i*pi*linspace(0,1-2/N,N/2)');
  tmp     = dd .* yy_odd;
  yy      = [yy_even + tmp;...
             yy_even - tmp];
elseif (N==2)
  yy = [1 1;1 -1]*xx;
elseif (N==4)
  yy = [1 0 1 0; 0 1 0 -1i; 1 0 -1 0;0 1 0 1i]*[1 0 1 0;1 0 -1 0;0 1 0 1;0 1 0 -1]*xx;
else
  fprintf(1,'N is not a power of two.\n');
  keyboard
end

return

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

function yy = myFFTrev(xx)

N = numel(xx);

if (N>=8)
  yy         = zeros(N,1);
  dd         = exp(-1i*pi*linspace(0,1-2/N,N/2)');
  yy(1:2:(N-1)) = myFFTrev(xx(1:(N/2)) + xx((N/2) + (1:(N/2))));
  yy(2:2:N)     = myFFTrev(xx(1:(N/2)).*dd - xx((N/2) + (1:(N/2))).*dd);
elseif (N==2)
  yy = [1 1;1 -1]*xx;
elseif (N==4)
  yy = [1 0 1 0; 0 1 0 -1i; 1 0 -1 0;0 1 0 1i]*[1 0 1 0;1 0 -1 0;0 1 0 1;0 1 0 -1]*xx;
else
  fprintf(1,'N is not a power of two.\n');
  keyboard
end

return