%改進處理3源10%抖動
clc;
clear;
N=1000;
taumin=0;
taumax=10;
K=201;
epsilon1=0.05;?
epsilon2=0.05*sqrt（2）;
epsilon3=0.05*sqrt（5）;
p1=1;
p2=sqrt（2）;
p3=sqrt（5）;
zeta0=0.03;
i=sqrt（-1）;
alpha=0.0375;
beta=0.1;
gamma=3;
T=465;
N1=465;
N2=328;
N3=207;
b=（taumax-taumin）/K;
for?k=1:K
????tau1（k）=（k-0.5）*b+taumin;
????b11（k）=max（[b2*epsilon1*tau1（k）]）;?
????b12（k）=max（[b2*epsilon2*tau1（k）]）;?
????b13（k）=max（[b2*epsilon3*tau1（k）]）;?
end
D1=zeros（size（tau1））;
C1=zeros（size（tau1））;???????????
O1=zeros（size（tau1））;
first1=zeros（size（tau1））;
D2=zeros（size（tau1））;
C2=zeros（size（tau1））;???????????
O2=zeros（size（tau1））;
first2=zeros（size（tau1））;
D3=zeros（size（tau1））;
C3=zeros（size（tau1））;???????????
O3=zeros（size（tau1））;
first3=zeros（size（tau1））;
threshold=zeros（size（tau1））;?????????????
for?k1=1:N1
????x1（k1）=k1*p1-epsilon1+（2*epsilon1）.*rand（1）;
end
for?k2=1:N2
????x2（k2）=k2*p2-epsilon2+（2*epsilon2）.*rand（1）;
end
for?k3=1:N3
????x3（k3）=k3*p3-epsilon3+（2*epsilon3）.*rand（1）;
end
for?n=2:N1
????for?m=1:n-1
????????delta_tau1=x1（n）-x1（m）;
????????if?（（delta_tau1>（1-epsilon1）*taumin）&&（delta_tau1<（1+epsilon1）*taumax））????
???????????for?k6=1:K???????
???????????????if?（（delta_tau1>tau1（k6）-b11（k6）/2）&&（delta_tau1??????????????????C1（k6）=C1（k6）+1;???????????????????
??????????????????if?first1（k6）==0
?????????????????????O1（k6）=x1（n）;%到達時間
?????????????????????first1（k6）=1;?????????%v????????????
??????????????????end??????
??????????????????eta01=（x1（n）-O1（k6））/tau1（k6）;?%初始相位????????
??????????????????nu1=floor（eta01+0.4999）;??????%v的迭代公式??
??????????????????if?nu1>0??????
?????????????????????if?（（（nu1==1）&&（x1（m）==O1（k6）））||（（nu1>=2）&&（abs（eta01/nu1-1）<=zeta0）））???????
????????????????????????O1（k6）=x1（n）;????
?????????????????????end?????????????????????????????????????
??????????????????end??????????
??????????????????eta1=（x1（n）-O1（k6））/tau1（k6）;????%初始相位???????
??????????????????D1（k6）=D1（k6）+exp（2*pi*i*eta1）;??????????
???????????????end?????????????
???????????end?
????????end
????end
end?
for?n=2:N2
????for?m=1:n-1
????????delta_tau2=x2（n）-x2（m）;
?????????if?（（delta_tau2>（1-epsilon2）*taumin）&&（delta_tau2<（1+epsilon2）*taumax））????
???????????for?k7=1:K???????
???????????????if?（（delta_tau2>tau1（k7）-b12（k7）/2）&&（delta_tau2??????????????????C2（k7）=C2（k7）+1;???????????????????
??????????????????if?first2（k7）==0
?