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

if rem(N,2) == 0  % Recursive divide and conquer
   n = (0:N/2-1)';
   yhat = recfft(v(1:2:N-1));
   zhat = recfft(v(2:2:N));
   fac = omega.^n;
   vhat = [yhat+fac.*zhat;yhat-fac.*zhat];
else              % direct Fourier matrix-vector product
   j = 0:N-1; k = j';
   F = omega.^(k*j);
   vhat = F*v;
end