% E02_4.m             09mar18	  06jun19   23mar23
% Computation of cross-correlation sequence; examples
clear all;  clf;

x = [1 3 -2 1 2 -1 4 4 2];
y = [2 -1 4 1 -2 3];
% y = y + randn(size(y))*std(y)/2; %%%%%%%% snr=2
Lx = length(x); Ly = length(y);  L = Lx+Ly-1;

r1 = conv(x,fliplr(y));      % crosscorr. in data domain
k = -Ly+1:Lx-1;
figure(1); stem(k,r1);
xlabel('Lag index'); ylabel('Amplitude');
v = axis;   axis([-Ly+1 Lx-1 v(3:end)]); 

[r2,lags] = xcorr(x,y);      % off-the-shelf Matlab command
figure(2); stem(lags,r2)     % (performs zero-padding)

% r3 = ifft(fft(x,L).*conj(fft(y,L)));      % crosscorr. in Fourier domain
% figure(3); stem(r3);
% xlabel('Lag index'); ylabel('Amplitude');
% 
%                           phase correlation [note +1 at den. to avoid /0]
% rph = ifft( (fft(x,L).*conj(fft(y,L))) ./ (abs((fft(x,L).*conj(fft(y,L)))) +1) );
% figure(4); stem(rph);
% xlabel('Lag index'); ylabel('Amplitude');

return;


%% --- autocorrelation
r = conv(x,fliplr(x));
Lx1 = length(x)-1;
k = -Lx1:Lx1;
stem(k,r);
xlabel('Lag index'); ylabel('Amplitude');
v = axis;   axis([-Lx1 Lx1 v(3:end)]);

%% ---  autocorr. of a noisy periodic signal
clf;
L = 1000;
n = 1:L;

x1 = sin(.01*pi*n) + randn(1,length(n))*.04;  % test with noise std in [0...5]
% x1 = sin(.01*pi*n) + cos(.07*pi*n) + randn(1,length(n))*.04;  % two freqs.
    % rem:  omega0 N = 2 pi r  --> N = 2 pi r / .01 pi = 200 r = 200
    % corr. values cycle 5 times in 1000 samples --> freq. is 5*2*pi/1000 = .01*pi
% x1 = 2*(round(sin(.01*pi*n)/2+.5)-.5) + randn(1,length(n))*0; % square wave
r = conv(x1,fliplr(x1));
subplot(211); plot(x1); grid on; 
subplot(212); plot(-L+1:L-1,r); grid on;

%% --- radar (sort of)
Lp = 25;  % TX duration 
x = linspace(3,0,Lp);  % e.g., amplitude of TX pulses has sawtooth shape
delay = 4;   % delay due to target's distance
y = [zeros(1,delay)  x  zeros(1,10)];  % ideal return, w. some final padding
Ly = Lp + delay + 10; % RX duration
yn = y + randn(1,Ly) * .5;  % noisy return, test with std in [0...5] 
x = [x  zeros(1,Ly-Lp)];  % padding of TX sequence
r = conv(yn,fliplr(x));   % crosscorr. with noisy version
k = -(Ly-1):(Ly-1);
time = 0:Ly-1;
subplot(311); stem(time,x);
subplot(312); stem(time,yn);
subplot(313); stem(k,r); hold on; stem(delay,r(Ly+delay),'filled'); hold off;
xlabel('Lag index'); ylabel('Amplitude');
v = axis;   axis([-Ly Ly v(3:end)]);

return
