
n=128;
x=zeros(n,1);
y=zeros(n,1);

%generate shifted cosines

for(i=1:n)
    x(i)=1+cos(i*2*pi/n);
    y(i)=cos(((i+20)*2*pi/n));

end


%generate two gaussian peaks
%{
for(i=1:n)
    c=50;
    sigma=5;
    x(i)=1*exp(-(i-c)^2/(2*sigma*sigma));
    c=90;
    sigma=5;    
    y(i)=1*exp(-(i-c)^2/(2*sigma*sigma));

end
%}
%plot them 
figure
p1=plot((1:n)',x)
hold on
p2=plot(((1+20):(n+20))',y)
%p2=plot(((1):(n))',y)
legend([p1 p2],'f', 'g')
%xlim([1 n])
%subtract mean value
dx=x-mean(x);
dy=y-mean(y);
%do not subtract mean value
%dx=x;
%dy=y;
[r1, lags]=xcorr(dx,dy);
%normalize cross-correlation (1)
normcoeff= sqrt(sum(dx.^2)*sum(dy.^2));

%normalize cross-correlation (2)
normarr=zeros(2*n-1,1);
for i=1:(n-1)
    xm=sqrt(sum(dx(i:end).^2));
    ym=sqrt(sum(dy(1:(n-i+1)).^2));
    normarr(round(n)+i-1)=xm*ym;
    ym=sqrt(sum(dy(i:end).^2));
    xm=sqrt(sum(dx(1:(n-i+1)).^2));    
    normarr(round(n)-i+1)=xm*ym;

end


%plot central part
figure
corrx=r1(round(n/2):round(3*n/2));
corrx=corrx/normcoeff;
plotx=(1:(numel(corrx)))';
plotx=plotx-0.5*n;

h1=plot(plotx, corrx);
hold on
corrxren=r1./normarr;
corrxren=corrxren(round(n/2):round(3*n/2));

h2=plot(plotx, corrxren);

legend([h1 h2],'normalized cross-correlation (1)','normalized cross-correlation (2)')