% Script for plotting the exact and DFT approximated derivative of the 
% Gaussian function exp(-x^2/2), truncated to x = [-4, 4], as well as its 
% DFT power spectrum.

  X1 = -4;
  X2 = 4;
  N = 16;

  L = X2 - X1;
  x = X1 + L*(0:(N-1))/N;
  M = [0:(N/2-1) (-N/2):(-1)];
  k = 2*pi*M/L;

  Nhi = 8*N;                    % Hi-res grid for plotting
  xhi = X1 + L*(0:(Nhi-1))/Nhi;

  y = exp(-x.^2/2);
  yhi = exp(-xhi.^2/2);
  Y = fft(y);
  dydx = real(ifft(1i*k.*Y)); % DFT approximation to dy/dx
  dydx_ex = -x.*y;            % Exact dy/dx on DFT grid
  dydx_hi = -xhi.*yhi;        % Exact dy/dx on hi-res grid

%  Plot exact and DFT dy/dx
  subplot(2,2,1)
  plot(xhi,dydx_hi,'-',x,dydx,'x')
  xlabel('x')
  ylabel('dy/dx')
  title(['Gaussian over [-4, 4], N = ' num2str(N)])
  legend('Exact','DFT',0)

%  Plot DFT power spectrum

  subplot(2,2,2)
  semilogy(M,abs(Y).^2,'x')
  xlabel('M_m')
  ylabel('|Y_m|^2')
  title('DFT Power spectrum')