clc %emol
clearvars %emol
close all %emol

fs=1000; T=5; 
t=0:1/fs:T-1/fs;

  % generation of broadband white noise (emol)
randn('state',0);
so=randn(size(t));
  % generation of band-limited white noise (emol)
fc=100; 
[b,a] = butter(9,fc/(fs/2));
so=filtfilt(b,a,so); 
  % signal with zero mean value and unitary standard deviation (emol)
s=so-mean(so); % mean(s)=0
s=s/std(s); % std(s)=1;

  % generation of signals with time delay equal to 0.2s (emol)
delta=0.2;
xo=s(delta*fs+1:end); % xo(t) = s(t + delta) (emol)
yo=s(1:end-delta*fs); % yo(t) = s(t) (emol)
tr=t(1:end-delta*fs); %emol
  % generation of broadband white noise to corrupt the sequences (emol)
randn('state',1); nx=1*std(s)*randn(size(xo));
randn('state',2); ny=1*std(s)*randn(size(yo));
  % generation of noisy signals
x=xo+nx; y=yo+ny; % SNR is 0 dB for both signals

  % computation of autocorrelation and cross-correlation (emol)
maxlag1=0.25*fs; maxlag2=0.5*fs;
[Rss,tau1]=xcorr(s,s,maxlag1,'unbiased');
[Rxy,tau2]=xcorr(y,x,maxlag2,'unbiased');
tau1=tau1/fs; tau2=tau2/fs;

  % visualization of signals (emol)
figure(1) 
subplot(3,2,1) 
plot(t,so) % original profile s(t)
xlabel('\itt') 
ylabel('\itso(t)')
title('Original profile')
axis([0 5 -5 5])
grid on
subplot(3,2,2)
plot(t,s) % profile with zero mean value and unit standard deviation
xlabel('\itt') 
ylabel('\its(t)')
title('Profile of zero mean and unit standard deviation')
axis([0 5 -5 5])
grid on
subplot(3,2,3)
plot(tr,xo) % uncorrupted x(t)
xlabel('\itt') 
ylabel('\itxo(t)')
title('Uncorrupted front signal')
axis([0 5 -5 5])
grid on
subplot(3,2,4)
plot(tr,yo) % uncorrupted y(t)
xlabel('\itt') 
ylabel('\ityo(t)')
title('Uncorrupted rear signal')
axis([0 5 -5 5])
grid on
subplot(3,2,5)
plot(tr,x) % noisy x(t)
xlabel('\itt') 
ylabel('\itx(t)')
title('Noisy front signal')
axis([0 5 -5 5])
grid on
subplot(3,2,6)
plot(tr,y) % noisy y(t)
xlabel('\itt') 
ylabel('\ity(t)')
title('Noisy rear signal')
axis([0 5 -5 5])
grid on

% figure(2) 
% plot(tau1,Rss)
% axis([-0.25 0.25 -0.4 1.2])
% xlabel('Lag (\it\tau)') 
% ylabel('Autocorrelation (\itR_s_s\rm(\it\tau\rm))')
% grid on %emol

% figure(3)
% plot(tau2(maxlag2+1:end),Rxy(maxlag2+1:end))
% axis([0 0.5 -0.4 1.2])
% xlabel('Lag (\it\tau)')
% ylabel('Cross-correlation (\itR_x_y\rm(\it\tau\rm))')
% grid on %emol

% 1st stage: signals;
% 2nd stage: signals plus correlations.