%修改PRI變換法分選程序
%作用：對已有脈沖序列進行修改PRI分選，繪制PRI譜圖，以及門限
%需要先生成雷達信號
%作者：巫書航
%時間：2007年11月13日


tmin=0.5;tmax=10;              %PRI取值范圍設定 
K=200;                        %PRI箱個數，K越大，PRI變換圖越精確
doum=0.05;                    %脈沖抖動上限
b=(tmax-tmin)/K;              %每個PRi箱寬度，此處也代表求得PRI的精度
k=[1:K];
tk=(k-0.5)*b+tmin;            %每個PRI箱的中心坐標
for tk1=1:length(tk)
bk(tk1)=max([b,2*doum*tk(tk1)]);                 %每個PRI箱寬度
end
NN=length(C);
Zi0=0.03;



Dk=zeros(size(tk));           %PRI變換函數Dk
Ck=zeros(size(tk));           %標記脈沖數函數
Ok=zeros(size(tk));           %時間起點函數
JI=zeros(size(tk));           %標記函數

    for nn=2:NN
            for mm=1:nn-1
                    tao=C(nn)-C(mm);                     %對脈沖到達時間進行重復相減計算差值 
                                     
                    if tao>(1-doum)*tmin&tao<(1+doum)*tmax,               %防止超限制的處理    
                        
                        for l=1:K
                            
                        if tao>tk(l)-bk(l)/2&tao<tk(l)+bk(l)/2,
                            Ck(l)=Ck(l)+1;
                            if JI(l)==0,Ok(l)=C(nn);JI(l)=1;
                            end                                  %若PRI箱第一次使用，則選時間起點為C(nn),并標記PRI箱已使用
                            Yi=(C(nn)-Ok(l))/tk(l);
                            vv=floor(Yi+0.4999);
                           
                            if vv>0;
                            if (vv==1&C(mm)==Ok(l))||(vv>=2&abs(Yi/vv-1)<=Zi0),
                                Ok(l)=C(nn);                                    %轉換時間起點
                            end
                            end
                            Yii=(C(nn)-Ok(l))/tk(l);                            %計算相位 
                            Dk(l)=Dk(l)+exp(2*pi*i*Yii);                        %更新PRI變換值    
                        end
                        end
                    end
    
            end
    end
   

for ll=1:K
 Dk(ll)=abs(Dk(ll));                %對Dk取模
end

plot(tk,Dk,'k');                      %做出PRI變換圖 



    
    

        
Dk1=zeros(size(tk));
%設定門限
aa=0.3;
bb=0.15;
cc=3;
for v=1:length(tk) 
    
       Men1(v)=aa*TT/tk(v);                     %代表以tt為PRI的脈沖個數  150是觀察時間
       Men2(v)=bb*Ck(v);                     %為了消除子諧波
       Men3=cc*(TT*(length(C)/TT)^2*bk(v))^0.5;    %為消除噪聲
       Dk1(v)=max([Men1(v),Men2(v),Men3]);
   end

plot(tk,Dk,'k',tk,Dk1,'k');                      %做出PRI變換圖