www.pudn.com > Pulse_Airborne_PD_DownLook_HighPRF_B.rar > Pulse_Airborne_PD_DownLook_HighPRF_B.m, change:2011-11-18,size:11899b


clear all 
clc 
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%%%%%    载机平台仿真参数    %%%%% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
Pos_Aircraft=[0;0;8000];         % 载机三维坐标【m】 
V_Aircraft=[0;300;0];            % 载机速度矢量【m/s】 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%%%%%    雷达系统仿真参数    %%%%% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
c=3e8;                           % 光速 
k=1.38e-23;                      % 玻尔兹曼常数 
 
Pt=10e3;                         % 发射功率【W】 
F=10^(4/10);                     % 噪声系数 
 
Fc=10e9;                         % 中心频率【Hz】 
Wavelength=c/Fc;                 % 工作波长【m】 
 
Tp=2.5e-6;                       % 脉冲宽度【微秒】 
B=5e6;                           % 带宽【Hz】 
K=B/Tp;                          % 调频率【Hz】 
 
Fr=50e3;                         % 脉冲重复频率【Hz】 
Tr=1/Fr;                         % 脉冲重复周期【秒】 
 
Fs=5e6;                          % 采样率【Hz】 
Delta_t=1/Fs;                    % 时域采样点时间间隔【秒】 
 
D=0.6;                           % 天线孔径【m】 
Ae=1*pi*(D/2)^2;                 % 天线有效面积【m^2】 
G=4*pi*Ae/Wavelength^2;          % 天线增益 
BeamWidth=0.88*Wavelength/D;     % 天线3dB波束宽度【deg】 
 
Theta0=0*pi/180;                 % 波束主瓣方位指向【度】 
Phita0=40*pi/180;                % 波束主瓣方位指向【度】 
Beam_vector=[cos(Phita0)*sin(Theta0);... 
    cos(Phita0)*cos(Theta0);... 
    -sin(Phita0)];               % 波束主瓣指向矢量 
 
CPI=64*Tr;                       % CPI仿真持续时间【秒】 
Num_Tr_CPI=128;                   % CPI周期数 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%%%%%      目标仿真参数     %%%%% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
Pos_target=[ Pos_Aircraft+ 4.500e3*[0;       cos(Phita0);             -sin(Phita0)],...            % 目标1 
             Pos_Aircraft+ 9.000e3*[0;       cos(Phita0);             -sin(Phita0)],...            % 目标2 
             Pos_Aircraft+12.700e3*[0;       cos(Phita0-0.5*pi/180);  -sin(Phita0-0.5*pi/180)],... % 目标3 
             Pos_Aircraft+12.446e3*[0;       cos(Phita0);             -sin(Phita0)]];              % 目标4 
%   目标坐标【m】                   X轴坐标       Y轴坐标                 Z轴坐标          
 
RCS_target=[20;20;20;20].*exp(j*2*pi*rand(4,1));  
                 % 目标平均后向散射截面积【m^2】    
                 % 假定目标RCS慢起伏,即在仿真时间段内恒定 
 
V_target=[ [ 0;            280;         0],...   % 目标2  
           [ 0;            200;         0],...   % 目标2  
           [ 0;           -190;        80],...   % 目标3  
           [ 0;            20;         20] ];    % 目标4  
%          X轴速度分量  Y轴速度分量  Z轴速度分量 【m/s】   
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%%%%%    数据流仿真方式     %%%%% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
t_set_Tp=(0:Delta_t:Tp)';                 % 一个脉冲的时间采样点数组 
N=length(t_set_Tp);                       % 一个脉冲的时间采样点数 
 
R_start=(Tp+Delta_t)*c/2; 
R_stop=(Tr-Delta_t)*c/2; 
R_set=R_start:c/2/Fs:R_stop; 
 
s=zeros(length(R_set)+N,Num_Tr_CPI);        % 定义信号数组 
 
%%%%%%%%%%%%%%%%%%%%%%%%% 
%%%%%  生成目标回波 %%%%% 
%%%%%%%%%%%%%%%%%%%%%%%%% 
for No_target=1:4 
    for No_Tr=1:Num_Tr_CPI 
        Pos_Aircraft_Tr=Pos_Aircraft+V_Aircraft*(No_Tr-1)*Tr; 
        % 雷达坐标 
        Pos_target_Tr=Pos_target(:,No_target)+V_target(:,No_target)*(No_Tr-1)*Tr; 
        % 目标坐标 
        Delta_vector=Pos_target_Tr-Pos_Aircraft_Tr; 
        % 目标与雷达之间的连线矢量 
        R=sqrt(Delta_vector.'*Delta_vector); 
        % 目标与雷达之间的距离 
        R_ambiguity=mod(R,c*Tr/2); 
        % 距离模糊后目标与雷达之间的距离 
        Angle_Different=acos((Delta_vector.'*Beam_vector)/R); 
        % 目标偏离雷达波束中心指向的角度 
        Gt=G*(sinc(Angle_Different/BeamWidth)); 
        Gr=G*(sinc(Angle_Different/BeamWidth)); 
        % 天线在目标方向上的增益 
        Lp=1./((4*pi)^2*R.^4); 
        % 目标处的电磁波双程传播损耗 
        Magnitude_echo=sqrt(Pt*RCS_target(No_target)*Gt*Gr*Wavelength^2/(4*pi)*Lp); 
        % 目标回波幅度 
        RVP=-2*pi*Fc*2*R/c; 
        % 目标回波视频检波剩余相位 
        if (R_ambiguity>=R_start)&&(R_ambiguity<=(R_stop-Tp*c/2)) 
            Echo_start=ceil((R_ambiguity-R_start)*2/c*Fs); 
            % 目标回波在数据数组中的起始采样点 
            Echo_stop=Echo_start+N-1; 
            % 目标回波结束采样点 
            Delta_t_sample=(Echo_start/Fs*c-(R_ambiguity-R_start)*2)/c; 
            % 目标回波起始采样点对应时延与目标回波包络前沿对应时延间的时间差 
            s(Echo_start:Echo_stop,No_Tr)=s(Echo_start:Echo_stop,No_Tr)+... 
                Magnitude_echo*exp(j*(pi*K*(t_set_Tp+Delta_t_sample-Tp/2).^2+RVP)); 
            % 将第n个目标的复回波的添加到其在信号数组中的对应位置; 
        elseif (R_ambiguity<R_start) 
            Echo_start=1; 
            % 目标回波在数据数组中的起始采样点 
            N_mis=floor((R_start-R_ambiguity)*2/c*Fs); 
            % 由于发射遮挡,目标回波丢失的采样点数 
            Echo_stop=N-N_mis; 
            % 目标回波结束采样点 
            Delta_t_sample=((R_ambiguity-R_start)*2-Echo_start/Fs*c)/c; 
            % 目标回波起始采样点对应时延与目标回波包络前沿对应时延间的时间差 
            s(Echo_start:Echo_stop,No_Tr)=s(Echo_start:Echo_stop,No_Tr)+... 
                Magnitude_echo*exp(j*(pi*K*(t_set_Tp(N_mis+1:N)+Delta_t_sample-Tp/2).^2+RVP)); 
            % 将第n个目标的复回波的添加到其在信号数组中的对应位置; 
        elseif (R_ambiguity>(R_stop-Tp*c/2)) 
            Echo_start=ceil((R_ambiguity-R_start)*2/c*Fs); 
            % 目标回波在数据数组中的起始采样点 
            N_mis=ceil((R_ambiguity*2/c-(Tr-Tp))*Fs); 
            % 由于发射遮挡,目标回波丢失的采样点数 
            Echo_stop=Echo_start+(N-N_mis)-1; 
            % 目标回波结束采样点 
            Delta_t_sample=(Echo_start/Fs*c-(R_ambiguity-R_start)*2)/c; 
            % 目标回波起始采样点对应时延与目标回波包络前沿对应时延间的时间差 
            s(Echo_start:Echo_stop,No_Tr)=s(Echo_start:Echo_stop,No_Tr)+... 
                Magnitude_echo*exp(j*(pi*K*(t_set_Tp(1:N-N_mis)+Delta_t_sample-Tp/2).^2+RVP)); 
            % 将第n个目标的复回波的添加到其在信号数组中的对应位置; 
        end 
    end 
end 
%%%%%%%%%%%%%%%%%%%%%%%%% 
 
 
%%%%%%%%%%%%%%%%%%%%%%%%% 
%%%%% 生成地面杂波  %%%%% 
%%%%%%%%%%%%%%%%%%%%%%%%%  
tic 
 
%  假定地面是平面  % 
%  采用等间隔同心圆划分方式划分地面  % 
R_ground_step=c/2/Fs/2; 
% 各等间隔同心圆半径斜距间隔 
R_ground_set=Pos_Aircraft(3):R_ground_step:18e3; 
% 各等间隔同心圆对应的地面到雷达斜距 
theta_delta=BeamWidth/10; 
% 各等间隔射线间的方位夹角 
 
for R_Cell=R_ground_set 
    %  按等间隔同心圆进行外循环 
    phita_ground=acos(Pos_Aircraft(3)/R_Cell); 
    %  该等距离环对应的入射角 
    if R_Cell==Pos_Aircraft(3) 
        theta_set=Theta0+(-pi:theta_delta:pi-theta_delta); 
        % 高度线杂波取-pi~pi整个同心圆 
    else 
        theta_set=Theta0+(-6*BeamWidth:theta_delta:6*BeamWidth); 
        % 各等间隔射线对应的方位角度,取波束中心方位角附近一部分角度范围 
    end 
    R_Cell_ground=sqrt(R_Cell^2-Pos_Aircraft(3)^2); 
    %  该等距离环对应的地面距离 
    Area_Cell=(((R_Cell+R_ground_step)^2-Pos_Aircraft(3)^2)-... 
        (R_Cell_ground)^2)*theta_delta/2; 
    %  该等距离环上每个网格的面积 
    RCS0=0.0133*cos(phita_ground)./(sin(phita_ground)+0.1*cos(phita_ground)).^3; 
    %  该等距离环上的单位面积平均后向散射截面积 【m^2】 
    RCS_set=Area_Cell*sqrt(((RCS0)*randn(length(theta_set),1)).^2+((RCS0)*... 
        randn(length(theta_set),1)).^2)/sqrt(pi/2); 
    %  该等距离环上的每个网格的后向散射截面积 【m^2】 
    No_Cell=1; 
    for theta=theta_set 
        %  按等间隔射线划分的网格进行内循环 
        Pos_Ground_Cell=R_Cell_ground*[sin(theta);cos(theta);0]; 
        %  该网格的三维坐标 
        Delta_vector=Pos_Ground_Cell-Pos_Aircraft; 
        % 该网格与雷达之间的连线矢量 
        R=sqrt(Delta_vector.'*Delta_vector); 
        % 该网格与雷达之间的距离 
        R_ambiguity=mod(R,c*Tr/2); 
        % 距离模糊后该网格与雷达之间的距离 
        Angle_Different=acos((Delta_vector.'*Beam_vector)/R); 
        % 该网格偏离雷达波束中心指向的角度 
        Gt=G*(sinc(Angle_Different/BeamWidth)); 
        Gr=G*(sinc(Angle_Different/BeamWidth)); 
        % 天线在该网格方向上的增益 
        Lp=1./((4*pi)^2*R.^4); 
        % 该网格处的电磁波双程传播损耗 
        Magnitude_echo=sqrt(Pt*RCS_set(No_Cell)*Gt*Gr*Wavelength^2/(4*pi)*Lp)*exp(1i*2*pi*rand(1,1)); 
        % 该网格回波幅度 
        Fd=2*(V_Aircraft.'*Delta_vector)/R/Wavelength; 
        % 该网格回波多普勒频率 
        if (R_ambiguity>=R_start)&&(R_ambiguity<=(R_stop-Tp*c/2)) 
            Echo_start=ceil((R_ambiguity-R_start)*2/c*Fs); 
            % 该网格回波在数据数组中的起始采样点 
            Echo_stop=Echo_start+N-1; 
            % 该网格回波结束采样点 
            Delta_t_sample=(Echo_start/Fs*c-(R_ambiguity-R_start)*2)/c; 
            % 该网格回波起始采样点对应时延与该网格回波包络前沿对应时延间的时间差 
            for No_Tr=1:Num_Tr_CPI 
                RVP=2*pi*Fd*(No_Tr-1)*Tr; 
                % 该网格回波视频检波剩余相位 
                s(Echo_start:Echo_stop,No_Tr)=s(Echo_start:Echo_stop,No_Tr)+... 
                    Magnitude_echo*exp(j*(pi*K*(t_set_Tp+Delta_t_sample-Tp/2).^2+RVP)); 
                % 将该网格的复回波的添加到其在信号数组中的对应位置; 
            end 
        elseif (R_ambiguity<R_start) 
            Echo_start=1; 
            % 该网格回波在数据数组中的起始采样点 
            N_mis=floor((R_start-R_ambiguity)*2/c*Fs); 
            % 由于发射遮挡,该网格回波丢失的采样点数 
            Echo_stop=N-N_mis; 
            % 该网格回波结束采样点 
            Delta_t_sample=((R_ambiguity-R_start)*2-Echo_start/Fs*c)/c; 
            % 该网格回波起始采样点对应时延与该网格回波包络前沿对应时延间的时间差 
           for No_Tr=1:Num_Tr_CPI 
                RVP=2*pi*Fd*(No_Tr-1)*Tr; 
                % 该网格回波视频检波剩余相位 
                s(Echo_start:Echo_stop,No_Tr)=s(Echo_start:Echo_stop,No_Tr)+... 
                    Magnitude_echo*exp(j*(pi*K*(t_set_Tp(N_mis+1:N)+Delta_t_sample-Tp/2).^2+RVP)); 
                % 将该网格的复回波的添加到其在信号数组中的对应位置; 
            end 
        elseif (R_ambiguity>(R_stop-Tp*c/2)) 
            Echo_start=ceil((R_ambiguity-R_start)*2/c*Fs); 
            % 该网格回波在数据数组中的起始采样点 
            N_mis=ceil((R_ambiguity*2/c-(Tr-Tp))*Fs); 
            % 由于发射遮挡,该网格回波丢失的采样点数 
            Echo_stop=Echo_start+(N-N_mis)-1; 
            % 该网格回波结束采样点 
            Delta_t_sample=(Echo_start/Fs*c-(R_ambiguity-R_start)*2)/c; 
            % 该网格回波起始采样点对应时延与该网格回波包络前沿对应时延间的时间差 
            for No_Tr=1:Num_Tr_CPI 
                RVP=2*pi*Fd*(No_Tr-1)*Tr; 
                % 该网格回波视频检波剩余相位 
                s(Echo_start:Echo_stop,No_Tr)=s(Echo_start:Echo_stop,No_Tr)+... 
                    Magnitude_echo*exp(j*(pi*K*(t_set_Tp(1:N-N_mis)+Delta_t_sample-Tp/2).^2+RVP)); 
                % 将该网格的复回波的添加到其在信号数组中的对应位置; 
            end 
        end 
        No_Cell=No_Cell+1; 
    end 
end 
toc 
%%%%%%%%%%%%%%%%%%%%%%%%% 
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%%%%%    仿真热噪声信号     %%%%% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
n=sqrt(k*B*F*290)*randn(length(R_set)+N,Num_Tr_CPI); 
s=s+n; 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
 
%%%%%%%%%%%%%%%%%%%%%%%%% 
Num_sample=length(R_set)+N; 
frange=((-Fs/2+Fs/(2*Num_sample)):Fs/Num_sample:(Fs/2-Fs/(2*Num_sample))).'; 
Win=hamming(Num_sample)*ones(1,Num_Tr_CPI); 
S=fftshift(fft(s,[],1),1); 
S=S.*(exp(1i*pi*frange.^2/K)*ones(1,Num_Tr_CPI)); 
s1=ifft(ifftshift(S.*Win,1),[],1); 
figure 
plot(abs(s1(:,1))); 
xlabel('一个脉冲内的快时间') 
ylabel('回波幅度') 
S=fftshift(fft(s1,[],2),2); 
figure 
plot(sum(abs(S),1)) 
xlabel('多普勒频率') 
ylabel('频谱幅度')