function main_fft_speed

%DRIVER1
DRIVER2

return

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

function DRIVER1

rng(0)

nvec = 2.^(10:16);
tt0  = NaN*ones(size(nvec));
tt1  = zeros(size(nvec));
tt2  = zeros(size(nvec));
for i = 1:length(nvec)
  xx = randn(nvec(i),1);
  if (nvec(i) <=5000)
    tic
    yy = DFT_matrix(nvec(i))*xx;
    tt0(i) = toc;
  end
  tic
  yy = myFFT(xx);
  tt1(i) = toc;
  tic
  yy = myFFTrev(xx);
  tt2(i) = toc;
  fprintf(1,'n = %7d     t0 = %6.3f     t1 = %6.3f     t2 = %6.3f\n',...
          nvec(i),tt0(i),tt1(i),tt2(i))
end

figure(1)
loglog(nvec,tt0,'c',...
       nvec,tt1,'b',...
       nvec,tt2,'r',...
       'LineWidth',2)
legend('Naive - N^2','myFFT','myFFTrev')
xlabel('n')
ylabel('Execution time in seconds')

figure(2)
semilogx(nvec,tt1./nvec,'b',...
         nvec,tt2./nvec,'r',...
         nvec(1),0,...
         'LineWidth',2)
legend('myFFT','myFFTrev')
xlabel('n')
ylabel('time / n')

keyboard

return

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

function DRIVER2

rng(0)

nvec  = 2.^(10:18);
tt1 = zeros(size(nvec));
tt2 = zeros(size(nvec));
tt3 = zeros(size(nvec));
for i = 1:length(nvec)
  xx = randn(nvec(i),1);
  tic
  yy = myFFT(xx);
  tt1(i) = toc;
  tic
  yy = myFFTrev(xx);
  tt2(i) = toc;
  tic
  for j = 1:50
    yy = fft(xx);
  end
  tt3(i) = toc/50;
  fprintf(1,'n = %7d     t1 = %6.3f     t2 = %6.3f     t3 = %7.5f\n',...
          nvec(i),tt1(i),tt2(i),tt3(i))
end

figure(1)
loglog(nvec,tt1,'b',...
       nvec,tt2,'r',...
       nvec,tt3,'k',...
       'LineWidth',2)
legend('myFFT','myFFTrev','Matlab built in fft')
xlabel('n')
ylabel('Execution time in seconds')

figure(2)
semilogx(nvec,tt1./nvec,'b',...
         nvec,tt2./nvec,'r',...
         nvec,tt3./nvec,'k',...
         nvec(1),0,...
         'LineWidth',2)
legend('myFFT','myFFTrev','Matlab built in fft')
xlabel('n')
ylabel('time / n')

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