%信號(hào)相乘后的正弦信號(hào)，進(jìn)而用互相關(guān)
for nn=1:1
clear

fs=30E+6;%采樣頻率
ts=1/fs;%采樣時(shí)間間隔
k=10E+10;  %調(diào)頻斜率
fo=5E+6;%載頻
T=40E-6;%脈沖寬度
B=k*T;
n=3500;%抽樣點(diǎn)數(shù)
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=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=-9;
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)
%--------------------------------------互相關(guān)四階累積量 
mix1=mix(1:2900);
mix2=mix(50:2950);

m=200;
for i=0:m
    mix_cum(:,i+1)=cum4x(mix1,mix2,mix2,mix2,m,1000,50,'unbiased',i,i);
end
for i=1:m+1
    mix_cum_chirp(i)=mix_cum(i+m,i);
end

l2=1024;
mix_cum_chirp_spect=abs(fft(mix_cum_chirp,l2));
figure
subplot(2,1,1)
plot(mix_cum_chirp)
subplot(2,1,2)
plot(mix_cum_chirp_spect)

value_max=max(mix_cum_chirp_spect);
position=find(mix_cum_chirp_spect==value_max);
k_detect=(position(1)-1)*fs/l2/(delay*ts)

%-----------直接求對(duì)角切片的的頻譜
mix_cum=cum4est(mix,m,1000,50,'unbiased',0,0);
mix_cum(m)=0;
mix_cum_spect=abs(fft(mix_cum,1024*2));
figure
subplot(2,1,1)
plot(mix_cum)
subplot(2,1,2)
plot(mix_cum_spect)

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

end