%相采用論文中介紹的方法增加信噪比
clear

fs=60E+6;%采樣頻率
ts=1/fs;%采樣時間間隔
k=10E+10;  %調頻斜率
fo=5E+6;%載頻
T=40E-6;%脈沖寬度
B=k*T;
n=4500;%抽樣點數
t=ts:ts:n*ts;

p=zeros(1,n);
p(500+1:500+T*fs)=ones(1,T*fs);
signal=1/sqrt(T)*cos(2*pi*(fo*t+k*t.^2/2)).*p;
signal_spect=abs(fft(signal));

signal_power=sum(signal.^2);
snr=-5;
noise_power=signal_power*10^(-snr/10);
a=sqrt(noise_power/n);
%a=0;
noise=a*randn(1,n);
s=signal+noise;
s_spect=abs(fft(s));

x=[zeros(1,n/2),s,zeros(1,n/2)];
r=0;
for tao=0:n/2
    r1=zeros(1,2*n);
    r1(n/2+1:n/2+n)=x(n/2+1+tao:n/2+n+tao).*x(n/2+1-tao:n/2+n-tao);
    r=r+r1;
end
r=r(n/2+1:n/2+n);
r_spect=abs(fft(r));

figure
subplot(2,1,1)
plot(s)
subplot(2,1,2)
plot(r)
figure
subplot(2,1,1)
plot(s_spect)
subplot(2,1,2)
plot(r_spect)

delay=500;
mix=r(1:n-delay).*r(delay+1:n);
mix_spect=abs(fft(mix));

figure
subplot(2,1,1)
plot(mix)
subplot(2,1,2)
plot(mix_spect)

mix_cum=cum4est(mix,1024,1600,60,'unbiased',0,0);
%mix_cum(1025)=mix_cum(1025)/10;
mix_cum_spect=abs(fft(mix_cum(1:1024),length(mix)));
figure
subplot(2,1,1)
plot(mix_cum)
subplot(2,1,2)
plot(mix_cum_spect)

value_max=max(mix_cum_spect(20:length(mix_cum_spect)/2-20));
position=find(mix_cum_spect==value_max);
k_detect=(position(1)-1)*fs/length(mix)/(delay*ts)/2;

ratio=abs((k_detect-k)/k)