graphics_toolkit gnuplotpkg 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 issignificantly 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')