function [y,nop] = myfft(x)
% Implementation of FFT
%
% Input: x, supposed to be a vector of length 2^p
%
% Output: y:   FFT(x)
%         nop: number of operations required
%
N = length(x);
l2N = log2(N);
assert(abs(l2N-round(l2N)) < 1.0E-10);
w = exp(-i*2*pi/N);
[y,nop] = domyfft(x,N,w);
% scale y
y = y/sqrt(N);


function [y,nop] = domyfft(x,N,w)
y = zeros(N,1);
if N==1,
    % in this case FFT is trivial
    nop = 0;
    y   = x;
else
    % split the sequence into even and odd parts
    x0  = x(1:2:end);
    x1  = x(2:2:end);
    % compute their FFTs
    mu  = w^2;
    [y0,nop0]  = domyfft(x0,N/2,mu);
    [y1,nop1]  = domyfft(x1,N/2,mu);
    
    % recombine them
    y(1:N/2)     = y0 + w.^(0:N/2-1).' .* y1; % 2*N/2-1 ops
    y(N/2+1:end) = y0 + w.^(N/2:N-1).' .* y1; % 2*N/2   ops
    nop = nop0 + nop1 + 2*N-1;
end