% eda09_04
% using the cross correlation to align two timeseries
% this example uses synthetic data

% time series paramaters
N=1000;
Dt=0.1;
t=Dt*[1:N]';

% synthetic time series, u and v
t1 = Dt*4*N/8;
t2 = Dt*5*N/8;
tlagtrue = t2-t1;
w0 = 2*pi*0.5;
u = cos( w0 * (t-t1) ) .* exp( -((t-t1).^2)/2 ) + random('Normal',0,0.02,N,1);
v = cos( w0 * (t-t2) ) .* exp( -((t-t2).^2)/2 ) + random('Normal',0,0.02,N,1);;

% plot original timeseries
figure(1);
clf;
subplot(2,1,1);
set(gca,'LineWidth',2);
hold on;
axis( [t(1), t(N), -1.1, 1.1] );
plot( t, u+0.1, 'k-', 'LineWidth', 2);
plot( t, v-0.1, 'k-', 'LineWidth', 2);
xlabel('time')
ylabel('u(t) and v(t)');

% cross correlation
c = xcorr(u, v);
Nc = length(c);
Ncenter = (Nc+1)/2;
[cmax, icmax] = max(c);
tlag = -Dt * (icmax-N);

% plot lagged timeseries
subplot(2,1,2);
set(gca,'LineWidth',2);
hold on;
axis( [t(1), t(N), -1.1, 1.1] );
plot( t, u+0.1, 'k-', 'LineWidth', 2);
plot( t-tlag, v-0.1, 'k-', 'LineWidth', 2);
xlabel('time')
ylabel('u(t) and v(t)');

% plot cross-correlation
figure(2);
clf;
set(gca,'LineWidth',2);
hold on;
axis( [-20, 20, 1.1*min(c), 1.1*max(c)] );
plot( Dt*[-Ncenter+1:Ncenter-1]', c, 'k-', 'LineWidth', 2);
xlabel('time')
ylabel('cross-correlation');

disp(sprintf('lag: true: %.2f estimated: %.2f', tlagtrue, tlag ));