function vhat = recfft(v)
%
%RECFFT - recursive radix-2 FFT.  If length(v) is
%         not divisible by 2 then a direct matrix
%         multiply is performed.
%
v = v(:);
N = length(v);
omega = exp(-2*pi*i/N);

if rem(N,2) == 0     % Recursive divide and conquer

   k = (0:N/2-1)'; twidle = omega.^k;  % twidle factors

   % transform 2 vectors that are half the length of v   
   yhat = recfft(v(1:2:N-1)); % recursive call 
   zhat = recfft(v(2:2:N));   % recursive call 

   % recombine yhat & zhat to form vhat; note only O(N) work!
   vhat = [yhat + twidle.*zhat;  % indices 0,...,N/2-1    
           yhat - twidle.*zhat]; % indices N/2,...,N-1

else  % just use direct multiplication by the Fourier matrix

   j = 0:N-1; 
   k = j';
   F = omega .^ (k*j);
   vhat = F*v;

end
