graphics_toolkit gnuplot
pkg load signal

% Characteristics common to both plots
  set(0, "DefaultAxesFontName", "Microsoft Sans Serif")
  set(0, "DefaultTextFontName", "Microsoft Sans Serif") 
  set(0, "DefaultAxesTitleFontWeight", "bold")
  set(0, "DefaultAxesFontWeight",      "bold")
  set(0, "DefaultAxesFontSize", 20)
  set(0, "DefaultAxesLineWidth", 3)
  set(0, "DefaultAxesBox", "on")
  set(0, "DefaultAxesGridLineStyle", "-")
  set(0, "DefaultAxesGridColor", [0 0 0])  % black
  set(0, "DefaultAxesGridAlpha", 0.25)     % opaqueness of grid
  set(0, "DefaultAxesLayer", "bottom")     % grid not visible where overlapped by graph
%========================================================================
function plotWindow (w, wname, wfilename = "", wspecifier = "", wfilespecifier = "")
 
  close     % If there is a previous screen image, remove it.
  M = 32;   % Fourier transform size as multiple of window length
  Q = 512;  % Number of samples in time domain plot
  P = 40;   % Maximum bin index drawn
  dr = 130; % (dynamic range) Maximum attenuation (dB) drawn in frequency domain plot
 
  L = length(w);
  B = L*sum(w.^2)/sum(w)^2;              % noise bandwidth (bins)
  
  n = [0 : 1/Q : 1];
  w2 = interp1 ([0 : 1/(L-1) : 1], w, n);
 
  if (M/L < Q)
    Q = M/L;
  endif
 
  figure("position", [1 1 1200 600])  % width = 2×height, because there are 2 plots
% Plot the window function
  subplot(1,2,1)
  area(n,w2,"FaceColor", [0 0.4 0.6], "edgecolor", [0 0 0], "linewidth", 1)
  
  g_x = [0 : 1/8 : 1];    % user defined grid X [start:spaces:end]
  g_y = [0 : 0.1 : 1];
  set(gca,"XTick", g_x)
  set(gca,"YTick", g_y)

% Special y-scale if filename includes "flat top"
  if(index(wname, "flat top"))
    ylimits = [-0.1 1.05];
  else
    ylimits = [0 1.05];
  endif
  ylim(ylimits)
  ylabel("amplitude","FontSize",28)  
  set(gca,"XTickLabel",[" 0"; " "; " "; " "; " "; " "; " "; " "; "  N"])
  grid("on")
  xlabel("samples","FontSize",28)
#{
% This is a disabled work-around for an Octave bug, if you don't want to run the perl post-processor.
  text(-.18, .4,"amplitude","rotation",90, "Fontsize", 28);
  text(1.15, .4,"decibels", "rotation",90, "Fontsize", 28);
#}

%Construct a title from input arguments.
%The default interpreter is "tex", which can render subscripts and the following Greek character codes:
%  \alpha \beta \gamma \delta \epsilon \zeta \eta \theta \vartheta \iota \kappa \lambda \mu \nu \xi \o
%  \pi \varpi \rho \sigma \varsigma \tau \upsilon \phi \chi \psi \omega.
%
  if (strcmp (wspecifier, ""))
    title(cstrcat(wname," window"), "FontSize", 28)    
  elseif (length(strfind (wspecifier, "&#")) == 0 )
    title(cstrcat(wname,' window (', wspecifier, ')'), "FontSize", 28)
  else
% The specifiers '\sigma_t' and '\mu' work correctly in the output file, but not in subsequent thumbnails.
% So UNICODE substitutes are used.  The tex interpreter would remove the & character, needed by the Perl script.
    title(cstrcat(wname,' window (', wspecifier, ')'), "interpreter", "none", "FontSize", 28)
  endif
  ax1 = gca;
 
% Compute spectal leakage distribution
  H = abs(fft([w zeros(1,(M-1)*L)]));
  H = fftshift(H);
  H = H/max(H);
  H = 20*log10(H);
  H = max(-dr,H);
  n = ([1:M*L]-1-M*L/2)/M;
  k2 = [-P : 1/M : P];
  H2 = interp1 (n, H, k2);
 
% Plot the leakage distribution
  subplot(1,2,2)
  h = stem(k2,H2,"-");
  set(h,"BaseValue",-dr)
  xlim([-P P])
  ylim([-dr 6])
  set(gca,"YTick", [0 : -10 : -dr])
  set(findobj("Type","line"), "Marker", "none", "Color", [0.8710 0.49 0])
  grid("on")
  set(findobj("Type","gridline"), "Color", [.871 .49 0])
  ylabel("decibels","FontSize",28)
  xlabel("bins","FontSize",28)
  title("Fourier transform","FontSize",28)
  text(-5, -126, ['B = ' num2str(B,'%5.3f')],"FontWeight","bold","FontSize",14)
  ax2 = gca;

% Configure the plots so that they look right after the Perl post-processor.
% These are empirical values (trial & error).
% Note: Would move labels and title closer to axes, if I could figure out how to do it.
  x1 = .08;         %   left margin for y-axis labels
  x2 = .02;         %  right margin
  y1 = .14;         % bottom margin for x-axis labels
  y2 = .14;         %    top margin for title
  ws = .13;         % whitespace between plots
  width  = (1-x1-x2-ws)/2;
  height = 1-y1-y2;
  set(ax1,"Position", [x1         y1 width height])      % [left bottom width height]
  set(ax2,"Position", [1-width-x2 y1 width height])
  
%Construct a filename from input arguments.
  if (strcmp (wfilename, ""))
    wfilename = wname;
  endif
  if (strcmp (wfilespecifier, ""))
    wfilespecifier = wspecifier;
  endif
  if (strcmp (wfilespecifier, ""))
    savetoname = cstrcat("Window function and frequency response - ", wfilename, ".svg");
  else
    savetoname = cstrcat("Window function and frequency response - ", wfilename, " (", wfilespecifier, ").svg");
  endif
  print(savetoname, "-dsvg", "-S1200,600")
% close    % Relocated to the top of the function
endfunction

%========================================================================
global N L 
% Generate odd-length, symmetric windows
N = 2^16;          % Large value ensures most accurate value of B
n = 0:N;
L = length(n);     % Window length

%========================================================================
w = ones(1,L);
plotWindow(w, "Rectangular")

%========================================================================
w = 1 - abs(n-N/2)/(L/2);
plotWindow(w, "Triangular")

% Indistinguishable from Triangular for large N
% w = 1 - abs(n-N/2)/(N/2);
% plotWindow(w, "Bartlett")
 
%========================================================================
w = parzenwin(L).';
plotWindow(w, "Parzen");

%========================================================================
w = 1-((n-N/2)/(N/2)).^2;
plotWindow(w, "Welch");

%========================================================================
w = sin(pi*n/N);
plotWindow(w, "Sine")
 
%========================================================================
w = 0.5 - 0.5*cos(2*pi*n/N);
plotWindow(w, "Hann")
 
%========================================================================
w = 0.53836 - 0.46164*cos(2*pi*n/N);
plotWindow(w, "Hamming", "Hamming", 'a_0 = 0.53836', "alpha = 0.53836")
 
%========================================================================
w = 0.42 - 0.5*cos(2*pi*n/N) + 0.08*cos(4*pi*n/N);
plotWindow(w, "Blackman")
 
%========================================================================
w = 0.355768 - 0.487396*cos(2*pi*n/N) + 0.144232*cos(4*pi*n/N) -0.012604*cos(6*pi*n/N);
plotWindow(w, "Nuttall", "Nuttall", "continuous first derivative")

%========================================================================
w = 0.3635819 - 0.4891775*cos(2*pi*n/N) + 0.1365995*cos(4*pi*n/N) -0.0106411*cos(6*pi*n/N);
plotWindow(w, "Blackman-Nuttall", "Blackman-Nuttall")
 
%========================================================================
w = 0.35875 - 0.48829*cos(2*pi*n/N) + 0.14128*cos(4*pi*n/N) -0.01168*cos(6*pi*n/N);
plotWindow(w, "Blackman-Harris", "Blackman-Harris")
 
%========================================================================
% Matlab coefficients
a = [0.21557895 0.41663158 0.277263158 0.083578947 0.006947368];
% Stanford Research Systems (SRS) coefficients
% a = [1  1.93  1.29  0.388  0.028];
% a = a / sum(a);
w = a(1) - a(2)*cos(2*pi*n/N) + a(3)*cos(4*pi*n/N) -a(4)*cos(6*pi*n/N) +a(5)*cos(8*pi*n/N);
plotWindow(w, "flat top")
 
%========================================================================
% The version using \sigma no longer renders correct thumbnail previews.
% Ollie's older version using &#963; seems to solve that problem.
sigma = 0.4;
w = exp(-0.5*( (n-N/2)/(sigma*N/2) ).^2);
% plotWindow(w, "Gaussian", "Gaussian", '\sigma = 0.4', "sigma = 0.4")
  plotWindow(w, "Gaussian", "Gaussian", "&#963; = 0.4", "sigma = 0.4")
 
%========================================================================
% Confined Gaussian
global T P abar target_stnorm
N = 512;                % Reduce N to avoid excessive computation time
n = 0:N;
L = length(n);          % Window length

target_stnorm = 0.1;
function [g,sigma_w,sigma_t] = CGWn(alpha, M)
  % determine eigenvectors of M(alpha)
  global L P T
  opts.maxit = 10000;
  if(M ~= L)
    [g,lambda] = eigs(P + alpha*T, M, 'sa', opts);
  else
    [g,lambda] = eig(P + alpha*T);
  end
  sigma_t = sqrt(diag((g'*T*g) / (g'*g)));
  sigma_w = sqrt(diag((g'*P*g) / (g'*g)));
end
function [h1] = helperCGW(anorm)
  global L abar target_stnorm
  [~,~,sigma_t] = CGWn(anorm*abar,1);
  h1 = sigma_t - target_stnorm * L;
end
% define alphabar, and matrices T and P
T = zeros(L,L);
P = zeros(L,L);
for m=1:L
  T(m,m) = (m - (L+1)/2)^2;
  for l=1:L
    if m ~= l
      P(m,l) = 2*(-1)^(m-l)/(m-l)^2;
    else
      P(m,l) = pi^2/3;
    end
  end
end
abar = (10/L)^4/4;
[anorm, aval] = fzero(@helperCGW, 0.1/target_stnorm);
[CGWg, CGWsigma_w, CGWsigma_t] = CGWn(anorm*abar,1);
sigma_t = CGWsigma_t/L          % Confirm sigma_t
w = CGWg * sign(mean(CGWg));
w = w'/max(w);
% \sigma_t works correctly in actual file, but not in thumbnail versions.
% plotWindow(w, "Confined Gaussian", "Confined Gaussian", '\sigma_t = 0.1', "sigma_t = 0.1");
  plotWindow(w, "Confined Gaussian", "Confined Gaussian", "&#963;&#8348; = 0.1", "sigma_t = 0.1");

  N = 2^16;          % restore original N
  n = 0:N;
  L = length(n);     % Window length

%========================================================================
global denominator;
sigma = 0.1;
denominator = (2*L*sigma).^2;
function [gaussout] = gauss(x)
  global N denominator
  gaussout = exp(- (x-N/2).^2 ./ denominator);
end
w = gauss(n) - gauss(-1/2).*(gauss(n+L) + gauss(n-L))./(gauss(-1/2 + L) + gauss(-1/2 - L));
% \sigma_t works correctly in actual file, but not in thumbnail versions
% plotWindow(w, "App. conf. Gaussian", "Approximate confined Gaussian", '\sigma_t = 0.1', "sigma_t = 0.1");
  plotWindow(w, "App. conf. Gaussian", "Approximate confined Gaussian", "&#963;&#8348; = 0.1", "sigma_t = 0.1");

%========================================================================
alpha = 0.5;
a = alpha*N/2;
w = ones(1,L);
m = 0 : a;
if( max(m) == a )
  m = m(1:end-1);
endif
M = length(m);
w(1:M) = 0.5*(1-cos(pi*m/a));
w(L:-1:L-M+1) = w(1:M);
% plotWindow(w, "Tukey", "Tukey", '\alpha = 0.5', "alpha = 0.5")
  plotWindow(w, "Tukey", "Tukey", "&#945; = 0.5", "alpha = 0.5")

%========================================================================
epsilon = 0.1;
a = N*epsilon;

w = ones(1,L);
m = 0 : a;
if( max(m) == a )
  m = m(1:end-1);
endif
% Divide by 0 is handled by Octave.  Results in w(1) = 0.
  z_exp = a./m - a./(a-m);
  M = length(m);
  w(1:M) = 1 ./ (exp(z_exp) + 1);
  w(L:-1:L-M+1) = w(1:M);

#{
% The original method is harder to understand:
t_cut = N/2 - a;
T_in = abs(n - N/2);
z_exp = (t_cut - N/2) ./ (T_in - t_cut)...
      + (t_cut - N/2) ./ (T_in - N/2);
% The numerator forces sigma = 0 at n = 0:
  sigma =  (T_in < N/2) ./ (exp(z_exp) + 1);        
% Either the 1st term or the 2nd term is 0, depending on n:
  w = 1 * (T_in <= t_cut)  +  sigma .* (T_in > t_cut);
#}

% plotWindow(w, "Planck-taper", "Planck-taper", '\epsilon = 0.1', "epsilon = 0.1")
  plotWindow(w, "Planck-taper", "Planck-taper", "&#949; = 0.1", "epsilon = 0.1")

%========================================================================
N = 2^12;               % Reduce N to avoid excess memory requirement
n = 0:N;
L = length(n);          % Window length

alpha = 2;
s = sin(alpha*2*pi/L*[1:N])./[1:N];
c0 = [alpha*2*pi/L,s];
A = toeplitz(c0);
[V,evals] = eigs(A, 1);
[emax,imax] = max(abs(diag(evals)));
w = abs(V(:,imax));
w = w.';
w = w / max(w);
% plotWindow(w, "DPSS", "DPSS", '\alpha = 2', "alpha = 2")
  plotWindow(w, "DPSS", "DPSS", "&#945; = 2", "alpha = 2")

%========================================================================
alpha = 3;
s = sin(alpha*2*pi/L*[1:N])./[1:N];
c0 = [alpha*2*pi/L,s];
A = toeplitz(c0);
[V,evals] = eigs(A, 1);
[emax,imax] = max(abs(diag(evals)));
w = abs(V(:,imax));
w = w.';
w = w / max(w);
% plotWindow(w, "DPSS", "DPSS", '\alpha = 3', "alpha = 3")
  plotWindow(w, "DPSS", "DPSS", "&#945; = 3", "alpha = 3")

N = 2^16;          % Restore original N
n = 0:N;
L = length(n);     % Window length
%========================================================================
alpha = 2;
w = besseli(0,pi*alpha*sqrt(1-(2*n/N -1).^2))/besseli(0,pi*alpha);
% plotWindow(w, "Kaiser", "Kaiser", '\alpha = 2', "alpha = 2")
  plotWindow(w, "Kaiser", "Kaiser", "&#945; = 2", "alpha = 2")
 
%========================================================================
alpha = 3;
w = besseli(0,pi*alpha*sqrt(1-(2*n/N -1).^2))/besseli(0,pi*alpha);
% plotWindow(w, "Kaiser", "Kaiser", '\alpha = 3', "alpha = 3")
  plotWindow(w, "Kaiser", "Kaiser", "&#945; = 3", "alpha = 3")
 
%========================================================================
alpha = 5; % Attenuation in 20 dB units
w = chebwin(L, alpha * 20).';
% plotWindow(w, "Dolph-Chebyshev", "Dolph-Chebyshev", '\alpha = 5', "alpha = 5")
  plotWindow(w, "Dolph&#8211;Chebyshev", "Dolph-Chebyshev", "&#945; = 5", "alpha = 5")
 
%========================================================================
w = ultrwin(L, -.5, 100, 'a')';
% \mu works correctly in actual file, but not in thumbnail versions
% plotWindow(w, "Ultraspherical", "Ultraspherical", '\mu = -0.5', "mu = -0.5")
  plotWindow(w, "Ultraspherical", "Ultraspherical", "&#956; = -0.5", "mu = -0.5")

%========================================================================
tau = (L/2);
w = exp(-abs(n-N/2)/tau);
% plotWindow(w, "Exponential", "Exponential", '\tau = N/2', "half window decay")
  plotWindow(w, "Exponential", "Exponential", "&#964; = N/2", "half window decay")
 
%========================================================================
tau = (L/2)/(60/8.69);
w = exp(-abs(n-N/2)/tau);
% plotWindow(w, "Exponential", "Exponential", '\tau = (N/2)/(60/8.69)', "60dB decay")
  plotWindow(w, "Exponential", "Exponential", "&#964; = (N/2)/(60/8.69)", "60dB decay")

%========================================================================
w = 0.62 -0.48*abs(n/N -0.5) -0.38*cos(2*pi*n/N);
plotWindow(w, "Bartlett-Hann", "Bartlett-Hann")

%========================================================================
alpha = 4.45;
epsilon = 0.1;
t_cut = N * (0.5 - epsilon);
t_in = n - N/2;
T_in = abs(t_in);
z_exp = ((t_cut - N/2) ./ (T_in - t_cut) + (t_cut - N/2) ./ (T_in - N/2));
sigma =  (T_in < N/2) ./ (exp(z_exp) + 1);
w = (1 * (T_in <= t_cut) + sigma .* (T_in > t_cut)) .* besseli(0, pi*alpha * sqrt(1 - (2 * t_in / N).^2)) / besseli(0, pi*alpha);
% plotWindow(w, "Planck-Bessel", "Planck-Bessel", '\epsilon = 0.1, \alpha = 4.45', "epsilon = 0.1, alpha = 4.45")
  plotWindow(w, "Planck&#8211;Bessel", "Planck-Bessel", "&#949; = 0.1, &#945; = 4.45", "epsilon = 0.1, alpha = 4.45")

%========================================================================
alpha = 2;
w = 0.5*(1 - cos(2*pi*n/N)).*exp( -alpha*abs(N-2*n)/N );
% plotWindow(w, "Hann-Poisson", "Hann-Poisson", '\alpha = 2', "alpha = 2")
  plotWindow(w, "Hann&#8211;Poisson", "Hann-Poisson", "&#945; = 2", "alpha = 2")
 
%========================================================================
w = sinc(2*n/N - 1);
plotWindow(w, "Lanczos")

%========================================================================
% optimized Nutall
ak = [-1.9501232504232442 1.7516390954528638 -0.9651321809782892 0.3629219021312954 -0.0943163918335154 ...
0.0140434805881681 0.0006383045745587 -0.0009075461792061 0.0002000671118688 -0.0000161042445001];

n = -N/2:N/2;
n = n/std(n);
w = 1;
for k = 1 : length(ak)
% This is an array addition, which expands the dimension of w[] as needed, and the value "1" is replicated.
   w = w + ak(k)*(n.^(2*k));
endfor
w = w/max(w);
plotWindow(w, "GAP optimized Nuttall")