% simulation data is the received data in practice
% 0:direct cross-correlation
% 1: direct interpolation
% 2: direct parabolic interpolation
% 3: parabolic interpolation with bias-compensation
% 4: parabolic-fit interploation combined with linear filter interpolation
% 5: parabolic interpolation to the complex correlation function envelope
% 6: matched filtering for interpolation
% 7: undersuample-linear interpolation
% 8: cosine interpolation
% 9: reconstructive interpolation
% 10: lagrange interpolation
% 11: sinc interpolation
% 12: spine interpolation

clc; clear all;
data=load('data_2_3');
s=data.No1_0;
clear data;

Fs=1e6; L=500; D=10.345e-6;
M=100;
s=s(1:L+50);


count=0;
for loop=1:1;
    count=count+1;  
%========interpolation==========
interp_s=interp(s,M);
t1=(1/Fs)*(1/M);
delay=round(D/t1);

%======two receive signals========
snr=0.05;
x1=s(1:L)+snr*randn(1,L);
   tempt=interp_s(delay+1:length(interp_s));
   tx2=tempt(1:M:length(tempt));
x2=tx2(1:L)+snr*randn(1,L);

clear tempt tx2;

%========= 0: direct cross-correaltion ============
 cx12=xcorr(x1,x2);
 [cm,ci]=max(cx12);
 ext0(count)=(ci-length(x1))/Fs;
 
%=========1: direct interpolation=========
px1=interp(x1,M);
px2=interp(x2,M);
px12=xcorr(px1,px2);
 [pm,pi]=max(px12);
 ext1(count)=(pi-length(x1)*M)*(1/Fs)*(1/M);

 %========2: direct parabolic interpolation=======
fx1=fft(x1);
fx2=fft(x2);
h=1./sqrt(abs(fx1).*abs(fx2));
gr=ifft(fx1.*conj(fx2).*h);
gr=abs(gr)/max(abs(gr));
[gm,gi]=max((gr));

es_D=(gr(gi+1)-gr(gi-1))/(2*(-gr(gi+1)+2*gr(gi)-gr(gi-1)));
est2(count)=(gi+es_D-1)/(Fs);
clear es_D;

 %==========4:parabolic-fit interploation combined with linear filter interpolation
 P=2;
w12(1:length(cx12)*P)=0;
for k=1:length(cx12);
    w12((k-1)*P+1)=cx12(k);
end
clear k;
[b,a]=butter(10,1/P);
fw12=filter(b,a,w12);
fw12=fw12/max(fw12);
[wm,wi]=max((fw12));
wi=wi-(L*P-1);
es_D=(fw12(wi+1)-fw12(wi-1))/(2*(-fw12(wi+1)+2*fw12(wi)-fw12(wi-1)));
est4(count)=(wi+es_D-1)/(P*Fs);
%est0(count)=i0/(Fs*P);
clear es_D;
 
 %==========5: parabolic interpolation to the complex correlation function envelope
 Fc=0.2*Fs;
 dx1=ademod(x1,Fc,Fs,'fm');
 dx2=ademod(x2,Fc,Fs,'fm');
 cdx12=xcorr(dx1,dx2);
 [cdm,cdi]=max(cdx12);
 cdi=cdi-length(dx1);
 est5(count)=(cdi)/(Fs);
 
 %==========6: matched filtering for interpolation
 
 
 %==========7: undersuample-linear interpolation
 N=10;
 ux1=x1(1:N:length(x1));
 ux2=x1(1:N:length(x2));
 
 
end