clc %emol
clearvars %emol
close all %emol

A=1; B=1; C=2; D=2; 
thetax=0; thetay=-pi/4; phi=pi/2; %emol
w=2*pi*1; n=2; %emol
fs=200; T=100; t=0:1/fs:T-1/fs; %emol

rel_time_delay=(thetax-thetay)/w; % relative time delay (emol)

x=A*sin(w*t+thetax)+B; % signal x(t) %emol
y=C*sin(w*t+thetay)+D*sin(n*w*t+phi); % signal y(t) %emol

maxlag=4*fs; % maximum time lag = 4 s
[Rxy2,tau]=xcorr(y,x,maxlag,'unbiased'); %emol
tau=tau/fs;

Rxy=0.5*A*C*cos(w*tau-(thetax-thetay)); % theoretical cross-correlation %emol

figure(1) 
plot(t(1:maxlag),x(1:maxlag),'b',t(1:maxlag),y(1:maxlag),'r')
xlabel('Time (s)'); ylabel('\itx\rm(\itt\rm) and \ity\rm(\itt\rm)')
legend('x(t)','y(t)')
grid on %emol
pause(3) %emol

figure(2)
  % time average
plot(tau(maxlag+1:end),Rxy2(maxlag+1:end),'b') % for lags >= 0 %emol
xlabel('Lag (\it\tau)'); ylabel('Cross-correlation')
grid on %emol
pause(3) %emol
hold on
  % theoretical
plot(tau(maxlag+1:end),Rxy(maxlag+1:end),'r'); % for lags >= 0 %emol
pause(2) %emol
  % indication of relative time delay (emol)
plot([rel_time_delay rel_time_delay],[-1.5 1.5],'k:') %emol
hold off
legend('time average','theoretical') %emol