% Script for plotting the Fourier power spectrum DFT power spectrum, and
% DFT approximation with N = 16 modes to three [0, 1]-periodic functions, 
% given in the m-files stp.m, saw.m, swell.m.

  N = 16;


  x = (0:(N-1))/N;  % Define the grid points
  M = [0:(N/2-1) (-N/2):(-1)];
  k = 2*pi*M;

  Nhi = 8*N;  % Define hi-res grid to plot functions and DFT approximations 
  xhi = (0:(Nhi-1))/Nhi;
  Mhi = [0:(Nhi/2-1) (-Nhi/2):(-1)];

%-----------------------------------------------------------
%  Step function

  [y,yhat] = stp(x,M);
  Y = fft(y);

%  Use DFT approximation to interpolate to grid points xhi.
%  Could replace the two lines below with the MATLAB function
%  yintp = interpft(y,Nhi)

  Yintp = (Nhi/N)*[Y(1:(N/2)) zeros(1,Nhi-N) Y((N/2 + 1):N)];
  yintp = real(ifft(Yintp));

%  Plot function and DFT approximation on hi-res grid

  subplot(2,2,1)
  ip = [1:N 1];
  iphi = [1:Nhi 1];
  [yhi,yhihat] = stp(xhi,Mhi);
  plot(xhi,yhi,'-',x,y,'.',xhi,yintp,'--')
  xlabel('x')
  title(['Step function, N = ' num2str(N)])
  legend('y(x)','y_i','y_N(x)',0)

%  Plot exact and DFT power spectrum

  subplot(2,2,2)
  semilogy(M, abs(yhat).^2,'+', M, (abs(Y)/N).^2,'o')
  xlabel('M')
  ylabel('Squared amplitude')
  title('Fourier power spectrum')
  legend('yhat', 'DFT',0)

  pause
%-----------------------------------------------------------
%  Sawtooth function

  [y,yhat] = saw(x,M);
  Y = fft(y);

%  Interpolate to grid points xhi 

  Yintp = (Nhi/N)*[Y(1:(N/2)) zeros(1,Nhi-N) Y((N/2 + 1):N)];
  yintp = real(ifft(Yintp));
  
%  Plot function and DFT approximation on hi-res grid

  subplot(2,2,1)
  ip = [1:N 1];
  iphi = [1:Nhi 1];
  [yhi,yhihat] = saw(xhi,Mhi);
  plot(xhi,yhi,'-',x,y,'.',xhi,yintp,'--')
  xlabel('x')
  title(['Sawtooth, N = ' num2str(N)])
  legend('y(x)','y_i','y_N(x)',0)

%  Plot exact and DFT power spectrum

  subplot(2,2,2)
  semilogy(M, abs(yhat).^2,'+', M, (abs(Y)/N).^2,'o')
  xlabel('M')
  ylabel('Squared amplitude')
  title('Fourier power spectrum')
  legend('yhat', 'DFT',0)

  pause

%-----------------------------------------------------------
%  Swell function 1/(1-0.6*cos(x))

  [y,yhat] = swell(x,M);
  Y = fft(y);

%  Interpolate to grid points xhi 

  Yintp = (Nhi/N)*[Y(1:(N/2)) zeros(1,Nhi-N) Y((N/2 + 1):N)];
  yintp = real(ifft(Yintp));
  
%  Plot function and DFT approximation on hi-res grid

  subplot(2,2,1)
  ip = [1:N 1];
  iphi = [1:Nhi 1];
  [yhi,yhihat] = swell(xhi,Mhi);
  plot(xhi,yhi,'-',x,y,'.',xhi,yintp,'--')
  xlabel('x')
  title(['1/(1 - 0.6(cos(2\pix)), N = ' num2str(N)])
  legend('y(x)','y_i','y_N(x)',0)

%  Plot exact and DFT power spectrum

  subplot(2,2,2)
  semilogy(M, abs(yhat).^2,'+', M, (abs(Y)/N).^2,'o')
  xlabel('M')
  ylabel('Squared amplitude')
  title('Fourier power spectrum')
  legend('yhat', 'DFT',0)