% Problem 4.1 in Ifeachor and Jervis DSP text

% Rich Kozick, March 1997

x1 = [6.02 -5.98 7.92 -7.96 -0.78 -8.34 9.22 -2.65 -3.7 9.51 5.53 3.5 -3.18 -8.85 8.21 1.69 -0.06 6.65 -8 -9.21 -0.78 7.27 -5.98 -3.97 9.11 4.23 2.99 -1.85 -5.27 3.81 6.62 -2.64 2.08 -5.91 -3.58 -1.65 3.64 -8.19 -3.5 4.84 7.25 2.93 -4.42 -8.21 3.61];

x2 = [8.93 -7.2 -0.82 3.23 1.44 5.43 -9.88 -1.13 0.79 9.83 -8.73 4.64 -8.49 -4.66 -8.84 5.55 -8.24 -0.37 2.71 4.63 1.88 -0.92 -5.33 9.01 9.23 -3.7 5.08 -0.72 -5.08 -2.6 9.67 -8.55 -3.08 4.18 8.11 0.74 -3.87 -4.09 8.03 6.91 -9.87 -3.62 -8.29 -5.8 -7.04];

N = length(x1);
n=1:N;

figure(1)
subplot(211)
plot(n, x1)
ax = axis;
axis([1 N ax(3) ax(4)])
title('Signal 1')
subplot(212)
plot(n, x2)
ax = axis;
axis([1 N ax(3) ax(4)])
title('Signal 2')

r11b = xcorr(x1,'biased');
r22b = xcorr(x2,'biased');
r12b = xcorr(x1, x2, 'biased');

k = -(N-1):1:(N-1);

figure(2)
subplot(211)
plot(k, r11b);
ax = axis;
axis([-(N-1) (N-1) ax(3) ax(4)])
title('Autocorrelation r11: biased')
subplot(212)
plot(k, r22b);
ax = axis;
axis([-(N-1) (N-1) ax(3) ax(4)])
title('Autocorrelation r22: biased')

figure(3)
plot(k, r12b);
ax = axis;
axis([-(N-1) (N-1) ax(3) ax(4)])
title('Cross-correlation r12: biased')

% Problem 4.2:  Correct for the "end effect" (previous is without correction)

r11u = xcorr(x1,'unbiased');
r22u = xcorr(x2,'unbiased');
r12u = xcorr(x1, x2, 'unbiased');

k = -(N-1):1:(N-1);

figure(4)
subplot(211)
plot(k, r11u);
ax = axis;
axis([-(N-1) (N-1) ax(3) ax(4)])
title('Autocorrelation r11: unbiased')
subplot(212)
plot(k, r22u);
ax = axis;
axis([-(N-1) (N-1) ax(3) ax(4)])
title('Autocorrelation r22: unbiased')

figure(5)
plot(k, r12u);
ax = axis;
axis([-(N-1) (N-1) ax(3) ax(4)])
title('Cross-correlation r12: unbiased')

% Problem 4.3

rho12 = (x1*x2') / sqrt( (x1*x1')*(x2*x2') );
fprintf('Percent correlation between x1 and x2 at zero lag = %f %%\n', rho12*100);

% How can you answer part (2) of problem 4.1 ?