function c=myfft(f)
% primitive fft implementation for powers of two
n=length(f);
if n==1
    c=f;
    return
end
m=n/2;
fplus=(f(1:m)+f(m+1:n))/2;
fminus=exp(-2*pi*(0:m-1)'/n).*(f(1:m)-f(m+1:n))/2;
ceven=myfft(fplus);
codd =myfft(fminus);
c(1+2*(0:m-1),1)=ceven;
c(2+2*(0:m-1),1)=codd;
return
