%采用信號相乘并用2 1/2 譜求調頻斜率
for nn=1:1
clear

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

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

figure
subplot(2,1,1)
plot(signal)
subplot(2,1,2)
plot(signal_spect)

signal_power=sum(abs(signal).^2);
snr=-7.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));

%figure
%subplot(2,1,1)
%plot(s)
%subplot(2,1,2)
%plot(s_spect)

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

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

value_max=max(mix_spect(10:round(length(mix_spect)/2)));
position=find(mix_spect==value_max)
k_detect=(position(1)-1)*fs/length(mix)/(delay*ts);
ratio=abs(k-k_detect)/k

%-------------------------------
m=400;
mix_cum=cum4est(mix,m,1000,60,'unbiased',0,0);

mix_cum(m+1)=0;
mix_cum_spect=abs(fft(mix_cum,length(mix)));

figure
subplot(2,1,1)
plot(mix_cum)
subplot(2,1,2)
plot(mix_cum_spect)

value_max1=max(mix_cum_spect(10:round(length(mix_cum_spect)/2)));
position1=find(mix_cum_spect==value_max1)
k_detect1=(position1(1)-1)*fs/length(mix)/(delay*ts);
ratio1=abs(k-k_detect1)/k
%-------------------------------------
ss=signal(delay+1:n).*signal(1:n-delay);
ss_cum=cum4est(ss,m,1000,60,'unbiased',0,0);
ss_cum_spect=abs(fft(ss_cum,length(mix)));
figure
subplot(2,1,1)
plot(ss_cum)
subplot(2,1,2)
plot(ss_cum_spect)
figure
hist(ss,100)

nn1=signal(1:n-delay).*noise(delay+1:n)+signal(delay+1:n).*noise(1:n-delay);
nn1_cum=cum4est(nn1,m,1000,60,'unbiased',0,0);
nn1_cum_spect=abs(fft(nn1_cum,length(mix)));
figure
subplot(2,1,1)
plot(nn1_cum)
subplot(2,1,2)
plot(nn1_cum_spect)
figure
hist(nn1,100)

nn2=noise(delay+1:n).*noise(1:n-delay);
nn2_cum=cum4est(nn2,m,1000,60,'unbiased',0,0);
nn2_cum_spect=abs(fft(nn2_cum,length(mix)));

figure
subplot(2,1,1)
plot(nn2_cum)
subplot(2,1,2)
plot(nn2_cum_spect)
figure
hist(nn2_cum_spect,100)
%------------------------------

end