graphics_toolkit gnuplot
  pkg load signal
#{
The gnuplot graphics toolkit is not actively maintained and has a number
of limitations that are unlikely to be fixed.  Communication with gnuplot
uses a one-directional pipe and limited information is passed back to the
Octave interpreter so most changes made interactively in the plot window
will not be reflected in the graphics properties managed by Octave.  For
example, if the plot window is closed with a mouse click, Octave will not
be notified and will not update its internal list of open figure windows.
The qt toolkit is recommended instead.
#}
  set(0, "DefaultAxesTitleFontWeight", "bold")
  set(0, "DefaultAxesFontWeight",      "bold")
  set(0, "DefaultAxesFontSize",        14)

  P = 2000;                     % signal duration
  x = 0:P-1;			% sampling times

% Any function will do for the signal.  It needn't be periodic, but a square wave is convenient.
  period=233;			% square wave period
  s = square((x-period*.67)*2*pi/period);
  n = 4;			% a harmonic (any harmonic will do)

  num_angles=1000;		% eval the correlation function at 1000 phases
  correl = [];
  phi = (0:1:num_angles-1)/num_angles;
  for m=phi
% The right-hand side is equivalent to sum(s.*cos(2*pi*(n/P*x - m))
  correl(end+1) = real(sum(s.*exp(-i*2*pi*(n/P*x - m))));
  endfor

  figure("position",[1 1 700 700])
  plot(phi*360, correl, "color","b", "linewidth",3)
  xlim([0 360]);
  set(gca, 'xtick', [0 1 2 3 4]*90);

  [maxval,idx] = max(correl);
  cossin = sum(s.*exp(i*2*pi*n*x/P))
% The next two things are equal
  mod(atan2(imag(cossin), real(cossin))*360/(2*pi), 360)
  (idx-1)*360/num_angles
% The next two things are equal
  maxval
  abs(cossin)

  hold on
  stem((idx-1)*360/num_angles, correl(idx), "filled", "linewidth",4, "color",[1 0 1], "markerfacecolor",[1 0 1])
  stem( 0,correl( 0*num_angles/360+1),"filled", "linewidth",4, "color",[1 0 1], "markerfacecolor", [1 0 1])
  stem(90,correl(90*num_angles/360+1),"filled", "linewidth",4, "color",[1 0 1], "markerfacecolor", [1 0 1])

  title("Cross-correlation function", "fontsize",14);
  xlabel("phase lag (degrees)",       "fontsize",14)
  ylabel("correlation amplitude",     "fontsize",14)