PRI变换的matlab程序

%改进处理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