graphics_toolkitgnuplot
pkgloadsignal
#{
Thegnuplotgraphicstoolkitisnotactivelymaintainedandhasanumber
oflimitationsthatareunlikelytobefixed.Communicationwithgnuplot
usesaone-directionalpipeandlimitedinformationispassedbacktothe
Octaveinterpretersomostchangesmadeinteractivelyintheplotwindow
willnotbereflectedinthegraphicspropertiesmanagedbyOctave.For
example,iftheplotwindowisclosedwithamouseclick,Octavewillnot
benotifiedandwillnotupdateitsinternallistofopenfigurewindows.
Theqttoolkitisrecommendedinstead.
#}
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;
form=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",[11700700])
plot(phi*360,correl,"color","b","linewidth",3)
xlim([0360]);
set(gca,'xtick',[01234]*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)

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

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