graphics_toolkit gnuplot
pkg load signal
%=======================================================
function Y = DFT(y,t,f)
  W = exp(-j*2*pi * t' * f);                    % Nx1 × 1x8N = Nx8N
  Y = abs(y * W);                               % 1xN × Nx8N = 1x8N
% Y(1)  = SUM(n=1,2,...,N): { e^(-B × t(n)^2) × e^(-j2p ×-4096/8N × t(n)) }
% Y(2)  = SUM(n=1,2,...,N): { e^(-B × t(n)^2) × e^(-j2p ×-4095/8N × t(n)) }
% Y(8N) = SUM(n=1,2,...,N): { e^(-B × t(n)^2) × e^(-j2p × 4095/8N × t(n)) }
  Y = Y/max(Y);
endfunction  

  T = 1;                             % time resolution (arbitrary)
  Nyquist = 1/T;                     % Nyquist bandwidth
  N = 1024;                          % sample size
  I  = 8;                            % freq interpolation factor
  NI = N*I;                          % number of frequencies in Nyquist bandwidth
  freq_resolution = Nyquist/NI;
  X  = (-NI/2 : NI/2 -1);            % center the frequencies at the origin
  freqs = X * freq_resolution;       % actual frequencies to be sampled and plotted

% (https://octave.org/doc/v4.2.1/Graphics-Object-Properties.html#Graphics-Object-Properties)
  set(0, "DefaultAxesXlim",[min(freqs) max(freqs)])
  set(0, "DefaultAxesYlim",[0 1.05])
  set(0, "DefaultAxesXtick",[0])
  set(0, "DefaultAxesYtick",[])
% set(0, "DefaultAxesXlabel","frequency")
  set(0, "DefaultAxesYlabel","amplitude")

#{
Sample a funtion at intervals of T, and display only the Nyquist bandwidth [-0.5/T 0.5/T].  
Technically this is just one cycle of a periodic DTFT, but since we can't see the periodicity,
it looks the same as a continuous Fourier transform, provided that the actual bandwidth is
significantly less than the Nyquist bandwidth; i.e. no aliasing.
#}
% We choose the Gaussian function e^{-B (nT)^2}, where B is proportional to bandwidth.
  B = 0.1*Nyquist;
  x = (-N/2 : N/2 -1);              % center the samples at the origin
  t = x*T;                          % actual sample times
  y = exp(-B*t.^2);                 % 1xN  matrix
  Y = DFT(y, t, freqs);             % 1x8N matrix

% Re-sample to reduce the periodicity of the DTFT.  But plot the same frequency range.
  T = 8/3;
  t = x*T;                         % 1xN
  z = exp(-B*t.^2);                % 1xN
  Z = DFT(z, t, freqs);            % 1x8N
%=======================================================
  hfig = figure("position", [1 1 1200 900]);

  x1 = .08;                   % left margin for annotation
  x2 = .02;                   % right margin
  dx = .05;                   % whitespace between plots
  y1 = .08;                   % bottom margin
  y2 = .08;                   % top margin
  dy = .12;                   % vertical space between rows
  height = (1-y1-y2-dy)/2;    % space allocated for each of 2 rows
  width  = (1-x1-dx-x2)/2;    % space allocated for each of 2 columns
  x_origin1 = x1;
  y_origin1 = 1 -y2 -height;  % position of top row
  y_origin2 = y_origin1 -dy -height;
  x_origin2 = x_origin1 +dx +width;
%=======================================================
% Plot the Fourier transform, S(f)

  subplot("position",[x_origin1 y_origin1 width height])
  area(freqs, Y, "FaceColor", [0 .4 .6])
% xlabel("frequency")            % leave blank for LibreOffice input
%=======================================================
% Plot the DTFT

  subplot("position",[x_origin1 y_origin2 width height])
  area(freqs, Z, "FaceColor", [0 .4 .6])
  xlabel("frequency")
%=======================================================
% Sample S(f) to portray Fourier series coefficients

  subplot("position",[x_origin2 y_origin1 width height])
  stem(freqs(1:128:end), Y(1:128:end), "-", "Color",[0 .4 .6]);
  set(findobj("Type","line"),"Marker","none")
% xlabel("frequency")            % leave blank for LibreOffice input
  box on
%=======================================================
% Sample the DTFT to portray a DFT

  FFT_indices = [32:55]*128+1;
  DFT_indices = [0:31 56:63]*128+1;
  subplot("position",[x_origin2 y_origin2 width height])
  stem(freqs(DFT_indices), Z(DFT_indices), "-", "Color",[0 .4 .6]);
  hold on
  stem(freqs(FFT_indices), Z(FFT_indices), "-", "Color","red");
  set(findobj("Type","line"),"Marker","none")
  xlabel("frequency")
  box on
%=======================================================
% Output (or use the export function on the GNUPlot figure toolbar).
  print(hfig,"-dsvg", "-S1200,800","-color", 'C:\Users\BobK\Fourier transform, Fourier series, DTFT, DFT.svg')