function g = cooleytukey(f)
% recursive FFT step, as in fftw paper. expects f row vector of even length

N = numel(f);
if N==1
  g = f;
  return                        % bottom level, stop recursion
end
w = exp(-2i*pi*(0:N/2-1)/N);    % twiddle factors
e = cooleytukey(f(1:2:end));    % even elements
o = cooleytukey(f(2:2:end));    % odd elements
g = [e + w.*o, e - w.*o];