

%                          Script Mfile  dft.m
%  This Mfile calculates the dft coefficients and plots the dft approximation
%  together with the given function f.  The function f must be given as
%  inline function or in an Mfile.  User must enter the number N of dft
%  coefficients.

N = input('enter N (N must be even)  ')

sample = linspace(0, 2*pi, N+1);
% these are the N sample points in the interval [0, 2*pi]
tt= sample(1:N);

% these are the sampled values of f
ff = f(tt);
% We use the fast fourier transform to find the dft coefficients d(k)
d = fft(ff)/N;
d0 = d(1);
dd = d(2:N);
n = [0:N-1];

disp(' these are the coefficients  d')
disp('    index      real(d)    imag(d)   ')

for k = 1:N
    fprintf('     %2d    %10.5f %10.5f \n', k-1,d(k), -i*d(k))
end


% We do not need to find the coefficients a(k) to form the interpolating
% function g.  We use the fact that a(k) = d(k) for k = 0,1,2, . . N/2, and
% a(-k) = d(N-k) for k = 1,2,. . . , N/2 .

% these are the plotting points

 t= linspace( 0, 2*pi ,201);

% this is the function g

sum = d0 + dd(N/2)*exp(i*t*N/2);
for k = 1:N/2  -1
    sum = sum + dd(k)*exp(i*k*t) + dd(N-k)*exp(-i*k*t);
end


M = max(abs(f(t))) +.5;

% plot the real parts of f and g

subplot(2,1,1)
  plot(t, real(f(t)), 'r')
  hold on
  plot(t, real(sum) )
  axis([0, 2*pi, -M,M])
  hold  off
  title(' Real parts of f (in red) and DFT approximation ')
  
  % plot the imaginary parts of f and g
  
subplot(2,1,2)
  plot(t, imag(f(t)), 'r' )
  hold on
  plot(t, imag(sum))
  axis([0 2*pi, -M,M])
  hold off
  title(' Imaginary parts of f (in red) and DFT approximation' )





