function fc=dft(N,f)

i=sqrt(-1);
h=2*pi/N; x=[0:h:2*pi*(1-1/N)]; fx=eval(f);wn=exp(-i*h);
for k=0:N-1,
  s= 0;
  for j=0:N-1
%    s = s + fx(j+1)*wn^((k-N/2)*j);
     s = s + fx(j+1)*exp( (2*pi*i*j*k)/N );
  end
 
  fc(k+1)=s/N;
end