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')