% function c = myfft(f)
%
% Recursive implementation of the FFT
%
% Inputs:
%    f    2^m by 1 vector (signal)
%
% Outputs
%    c    2^m by 1 vector (fft of signal)
%
%    The entries of the vector f are
%
%     f(j) = F(x_{j-1}) 
%
%    where
%
%     x_j = 2*pi*j/N, j = 0..N-1, N = 2^m
%
%    The entries of the vector c are the coefficients in the polynomial
%
%      p(x) = \sum_{k=0}^{N-1} c(k+1) exp(i k x)
%
%    such that p(x_j) = F(x_j), j = 0..N-1
%
% Note: This code is not efficient at all. The overhead of
%       recursion will take over for length(x) large.
%
% Note: Matlab's fft is such that fft(f) == N*myfft(f)
%
% Fernando Guevara Vasquez 2011
function c = myfft(f)
N = length(f);
if (N==1)
 c = f;
else
 ce = myfft(f(1:2:N-1)); % even interpolant
 co = myfft(f(2:2:N));   % odd  interpolant
 ph = exp(-2i*(0:(N/2-1))'*pi/N);
 c = 0.5*[ ce + ph.*co 
           ce - ph.*co ];
end;
