main36.m

仿真了地球杂波的情况。考虑了地球自转带来的多普勒旋转等影响。仿真结果和理论相一致。可作为星载雷达地球杂波仿真的参考

mesh(Delta_Rst_Phi_Result);
xlabel('Phi');
ylabel('Theta');
zlabel('Delt_R/Delt_Phi');

subplot(3,2,3);
mesh(Delta_Fd_Theta_Result);
xlabel('Phi');
ylabel('Theta');
zlabel('Delt_Fd/Delt_Theta');
subplot(3,2,6);
mesh(Delta_Fd_Phi_Result);
xlabel('Phi');
ylabel('Theta');
zlabel('Delt_Fd/Delt_Phi');

%% 对相交区域进行分割
% 求得分隔区域的间隔大小
Delta_Theta_From_Rst_All = Radar_Range_Resolution ./ (Delta_Rst_Theta_Result+1e-20);
Delta_Theta_From_Fd_All  = Radar_Fd_Resolution ./ (Delta_Fd_Theta_Result+1e-20);
Delta_Phi_From_Rst_All = Radar_Range_Resolution ./ (Delta_Rst_Phi_Result+1e-20);
Delta_Phi_From_Fd_All  = Radar_Fd_Resolution ./ (Delta_Fd_Phi_Result+1e-20);
Delta_Theta_From_Rst_All = abs(Delta_Theta_From_Rst_All);
Delta_Theta_From_Fd_All = abs(Delta_Theta_From_Fd_All);
Delta_Phi_From_Rst_All = abs(Delta_Phi_From_Rst_All);
Delta_Phi_From_Fd_All =abs(Delta_Phi_From_Fd_All);
clear Delta_Rst_Theta_Result Delta_Fd_Theta_Result Delta_Rst_Phi_Result Delta_Fd_Phi_Result;

% 比较获得最终的分割间距
Min_Delta_Theta_Get = min(min(min(Delta_Theta_From_Rst_All)), min(min(Delta_Theta_From_Fd_All)));
Min_Delta_Phi_Get = min(min(min(Delta_Phi_From_Rst_All)), min(min(Delta_Phi_From_Fd_All)));
% min(min(Delta_Theta_From_Rst_All))
% min(min(Delta_Theta_From_Fd_All))
% min(min(Delta_Phi_From_Rst_All))
% min(min(Delta_Phi_From_Fd_All))

Min_Delta_Theta_Get = Min_Delta_Theta_Get*Min_Delta_Theta_Factor;
Min_Delta_Phi_Get = Min_Delta_Phi_Get*Min_Delta_Phi_Factor;
clear Delta_Theta_From_Rst_All Delta_Theta_From_Rst_All Delta_Phi_From_Rst_All Delta_Phi_From_Rst_All ...
    Delta_Phi_From_Fd_All Delta_Theta_From_Fd_All;

% 对相交区域进行分割，并得到各个分割区域的相关位置信息
% Phi应该的分割个数
Delta_Phi_Num = floor((Accurate_Intersect_Data_Phi_Radian(end) - Accurate_Intersect_Data_Phi_Radian(1)) / ...
    Min_Delta_Phi_Get);
% 每个Phi对应的Theta分割个数
Num_Delta_Theta_Of_EachPhi = zeros(1, Delta_Phi_Num);
Num_Tmp = 500;
% Theta角度，中心位置
Theta_Of_EachPart_Set = zeros(Delta_Phi_Num, Num_Tmp);
% Phi角度，中心位置
Phi_Of_EachPart_Set = zeros(Delta_Phi_Num, Num_Tmp);
% 擦地角，中心位置
Beta_S_Of_EachPart_Set = zeros(Delta_Phi_Num, Num_Tmp);
% 面积
Area_Of_EachPart_Set = zeros(Delta_Phi_Num, Num_Tmp);
% 径向距离
R_Of_EachPart_Set = zeros(Delta_Phi_Num, Num_Tmp);
% 最小径向距离
Rst_Min = Const_Earth_Re;
% 多普勒频移
Fd_Of_EachPart_Set = zeros(Delta_Phi_Num, Num_Tmp);
for w1 = 1 : 1 : Delta_Phi_Num
    Phi_Central_Tmp = Accurate_Intersect_Data_Phi_Radian(1) + (w1-1/2)*Min_Delta_Phi_Get;
    % 计算对于当前Phi的Theta可以划分的Theta区域个数
    Upside_Theta_Radian_Tmp = interp1(Accurate_Intersect_Data_Phi_Radian, Accurate_Intersect_Data_Upside_Theta_Radian,...
        Phi_Central_Tmp);
    Downside_Theta_Radian_Tmp = interp1(Accurate_Intersect_Data_Phi_Radian, Accurate_Intersect_Data_Downside_Theta_Radian,...
        Phi_Central_Tmp);
    Tmp = (Upside_Theta_Radian_Tmp - Downside_Theta_Radian_Tmp) / Min_Delta_Theta_Get;
       
    Tmp = round(Tmp);
    Num_Delta_Theta_Of_EachPhi(w1) = Tmp;
    for w2 = 1 : 1 : Num_Delta_Theta_Of_EachPhi(w1)
        % Theta角度
        Theta_Central_Tmp = Downside_Theta_Radian_Tmp + (w2-1/2)*Min_Delta_Theta_Get;
        Theta_Of_EachPart_Set(w1, w2) = Theta_Central_Tmp;
        % Phi角度
        Phi_Of_EachPart_Set(w1, w2) = Phi_Central_Tmp;
        % 面积
        Area_Of_EachPart_Set(w1, w2) = Const_Earth_Re^2 * Min_Delta_Phi_Get * Min_Delta_Theta_Get;
        % 径向距离和径向速度
        [Delta_Rst_Theta, Delta_Rst_Phi, Delta_Vst_Theta, Delta_Vst_Phi, Rst, Vst] = ...
            Calculate_R_V_DeltaR_DeltaV(Satellite_r, Const_Earth_Re, Satellite_Theta_s_Radian, ...
            Satellite_Phi_s_Radian, Satellite_Theta_i_Radian, Earth_We, Satellite_Wp, ...
            Theta_Central_Tmp, Phi_Central_Tmp);

        R_Of_EachPart_Set(w1, w2) = Rst;
        if(Rst_Min - Rst > 0)
            Rst_Min = Rst;
        end;
        % 多普勒频移
        Fd_Of_EachPart_Set(w1, w2) = -2*Vst / Const_Wave_Length;
        % 擦地角
        Beta_S_Of_EachPart_Set(w1, w2) = acos((Rst^2 + Const_Earth_Re^2 - Satellite_r^2) / (2*Rst*Const_Earth_Re)) - pi/2;
    end;
    
%     Phi_Central_Tmp
%     Upside_Theta_Radian_Tmp
%     Downside_Theta_Radian_Tmp
%     figure;
%     Num_Delta_Theta_Of_EachPhi(w1)
%     plot(Theta_Of_EachPart_Set(w1, 1:Num_Delta_Theta_Of_EachPhi(w1)), ...
%         Fd_Of_EachPart_Set(w1, 1:Num_Delta_Theta_Of_EachPhi(w1)));
%     grid on;
%     pause;
%     close all;
%     clc;
end;
fprintf('按任意键继续\n');
%  figure;
%  plot(Num_Delta_Theta_Of_EachPhi);
%  grid on;
pause;
% figure;
% mesh(R_Of_EachPart_Set);
% xlabel('Phi(经度)');
% ylabel('Theta(纬度)');
% grid on;
% title('相对距离');
% 
% figure;
% mesh(Fd_Of_EachPart_Set);
% xlabel('Phi(经度)');
% ylabel('Theta(纬度)');
% grid on;
% title('多普勒频移');

fprintf('地面散射点到卫星间距范围: %5.9f ——%5.9f\n', Rst_Min, max(max(R_Of_EachPart_Set)));
fprintf('地面散射点多普勒频率范围: %5.9f——%5.9f\n', min(min(Fd_Of_EachPart_Set)), max(max(Fd_Of_EachPart_Set)));
% pause;
% figure;
% plot(Num_Delta_Theta_Of_EachPhi);
% figure;
% plot(Theta_Of_EachPart_Set );
% figure;
% plot(Phi_Of_EachPart_Set );
% figure;
% plot(R_Of_EachPart_Set);
% figure;
% plot(Fd_Of_EachPart_Set);

%% 根据所分的小区间仿真得到回波信号
% 得到每个散射区间的回波信号幅度和相位信息
[Amp_Of_EachPart_Set, PhaseOrg_Of_EachPart_Set] = Get_Amp_Phase_Of_EachPart(...
    Num_Delta_Theta_Of_EachPhi, Theta_Of_EachPart_Set, Phi_Of_EachPart_Set ,...
    Beta_S_Of_EachPart_Set , Area_Of_EachPart_Set, R_Of_EachPart_Set, ...
    Const_Earth_Re, FA2D_Factor_Matrix, FA2D_Adder_Matrix, ...
    Antenna_Delta_Phi_w_Radian, Antenna_Delta_Theta_w_Radian, Factor_0D, IsRealSimulate);
% 为了避免产生数据量过大，获取感兴趣的距离范围，注意，这里没有考虑距离模糊问题
% 每个PRT回波的基准点
R_Max = max(max(R_Of_EachPart_Set));
R_Min = Rst_Min;
R_Max_Point = ceil(2 * R_Max/3e8 * Radar_Fs);
R_Min_Point = floor(2 * R_Min/3e8 * Radar_Fs);
R_Based_Point = R_Min_Point - 1;
% 每个PRT对应的有效采样点数
PRT_Valid_Point_Num = R_Max_Point - R_Min_Point;
% 仿真得到的回波信号，1列对应一个PRT回波
Simulate_Echo_Signal_Set = zeros(PRT_Valid_Point_Num, Radar_Coherent_PRT_Num);
% 用于数据处理时的缓存
Simulate_Echo_Signal_Tmp = zeros(PRT_Valid_Point_Num, Radar_Coherent_PRT_Num);
% 分别求得每个散射区间反射的信号，并且相加后得到总体的回波信号
for w1 = 1 : 1 : Delta_Phi_Num
    Num_Delta_Theta_Of_EachPhi_Tmp = Num_Delta_Theta_Of_EachPhi(w1);
    fprintf('(%d) / %d\n',w1,Delta_Phi_Num);
    for w2 = 1 : 1 : Num_Delta_Theta_Of_EachPhi_Tmp
        Theta_Of_EachPart_Tmp = Theta_Of_EachPart_Set(w1, w2);
        Phi_Of_EachPart_Tmp = Phi_Of_EachPart_Set(w1, w2);
        Area_Of_EachPart_Tmp = Area_Of_EachPart_Set(w1, w2);
        R_Of_EachPart_Tmp = R_Of_EachPart_Set(w1, w2);
        Fd_Of_EachPart_Tmp = Fd_Of_EachPart_Set(w1, w2);
        Amp_Of_EachPart_Tmp = Amp_Of_EachPart_Set(w1, w2);
        PhaseOrg_Of_EachPart_Tmp = PhaseOrg_Of_EachPart_Set(w1, w2);   
        
        % 清空缓存中的数据
        Simulate_Echo_Signal_Tmp = Simulate_Echo_Signal_Tmp * 0;
        % 求得每个散射区的间距对应的延迟，这里求得的延迟是以R_Based_Point作为参考的
        Delay_Point = round(2*R_Of_EachPart_Tmp/3e8*Radar_Fs) - R_Based_Point;
        % 这里将其转换为索引
        Delay_Point_Index = Delay_Point;
        Simulate_Echo_Signal_Tmp(Delay_Point_Index : 1 : Delay_Point_Index+Radar_Sample_Point_Of_Pluse_Width-1,...
            :) = Amp_Of_EachPart_Tmp;

        % 产生多普勒频移对应的信号
        Time_Tmp = (0 : 1 : Radar_Coherent_PRT_Num-1) * Radar_PRT;
        Fd_Signal_Tmp = exp(j*(2*pi*Fd_Of_EachPart_Tmp*Time_Tmp + PhaseOrg_Of_EachPart_Tmp));
        for w3 = 1 : 1 : Radar_Coherent_PRT_Num
            Simulate_Echo_Signal_Tmp(:, w3) = Simulate_Echo_Signal_Tmp(:, w3) * Fd_Signal_Tmp(w3);
        end;
        % 将每个散射单元的反射信号合成相加
        Simulate_Echo_Signal_Set = Simulate_Echo_Signal_Set + Simulate_Echo_Signal_Tmp;
%         [R_Num, C_Num] = size(Simulate_Echo_Signal_Set);
%         Signal_Tmp = reshape(Simulate_Echo_Signal_Set, 1, R_Num*C_Num);
%         figure;
%         plot(Signal_Tmp);
%         pause;
%         close all;
    end;
end;

% [R_Num, C_Num] = size(Simulate_Echo_Signal_Set);
% Signal_Tmp = reshape(Simulate_Echo_Signal_Set, 1, R_Num*C_Num);
% figure;
% plot(Signal_Tmp);

%% 对回波信号进行相干累加处理
% 一列对应一个距离波门
Simulate_Echo_Signal_Set = Simulate_Echo_Signal_Set.';
% 进行FFT变换
Simulate_Echo_Signal_Set_FFT = fft(Simulate_Echo_Signal_Set);
% 求得FFT变换结果每列对应的距离
[R_Num, C_Num] = size(Simulate_Echo_Signal_Set_FFT);
R_Point_Index = R_Based_Point + (1 : 1 : C_Num);
R_Index = R_Point_Index * 1/Radar_Fs * 3e8 / 2;

% 不进行距离模糊处理的结果
Freq_Index = (-Radar_Coherent_PRT_Num/2 : 1 : Radar_Coherent_PRT_Num/2-1)*...
    1/(Radar_PRT*Radar_Coherent_PRT_Num);
[R_Index_Org, Freq_Index_Org] = meshgrid(R_Index, Freq_Index);
Simulate_Echo_Signal_Set_FFT_ABS_Org = fftshift(abs(Simulate_Echo_Signal_Set_FFT), 1);
% figure;
% mesh(R_Index_Org, Freq_Index_Org, Simulate_Echo_Signal_Set_FFT_ABS_Org);
figure;
contour(R_Index_Org, Freq_Index_Org, Simulate_Echo_Signal_Set_FFT_ABS_Org);
xlabel('距离波门(m)');
ylabel('频率(Hz)');
zlabel('幅度谱');
title('地球杂波的距离频谱二维图（原始的，没有距离模糊）');




% 最大不模糊距离
PRT_Distance = Radar_PRT * 3e8 / 2;
% 得到雷达所检测到的距离，对距离进行模糊处理
R_Index = mod(R_Index, PRT_Distance);
figure;
plot(R_Index);
grid on;
title('通过模糊处理得到的距离');
% 对原始回波数据进行距离模糊处理，即将同一个距离单元内的信号相加
Simulate_Echo_Signal_Set_FFT_Tmp = zeros(size(Simulate_Echo_Signal_Set_FFT));
R_Index_Tmp = unique(R_Index);
for w1 = 1: 1 : length(R_Index_Tmp)
    [Index_Tmp] = find(R_Index == R_Index_Tmp(w1));
    Simulate_Echo_Signal_Set_FFT_Tmp(:, w1) = sum(Simulate_Echo_Signal_Set_FFT(:, Index_Tmp), 2);
end;
Simulate_Echo_Signal_Set_FFT = Simulate_Echo_Signal_Set_FFT_Tmp(:, 1 : length(R_Index_Tmp));
R_Index = R_Index_Tmp;
clear Simulate_Echo_Signal_Set_FFT_Tmp;
% 求模
Simulate_Echo_Signal_Set_FFT_ABS = fftshift(abs(Simulate_Echo_Signal_Set_FFT), 1);

Freq_Index = (-Radar_Coherent_PRT_Num/2 : 1 : Radar_Coherent_PRT_Num/2-1)*...
    1/(Radar_PRT*Radar_Coherent_PRT_Num);
[R_Index, Freq_Index] = meshgrid(R_Index, Freq_Index);
figure;
mesh(R_Index, Freq_Index, Simulate_Echo_Signal_Set_FFT_ABS);
xlabel('距离波门(m)');
ylabel('频率(Hz)');
zlabel('幅度谱');
title('地球杂波的距离频谱二维图');
grid on;
fprintf('距离-多普勒二维图距离范围: %5.9f——%5.9f\n', min(min(R_Index)), max(max(R_Index)));
fprintf('距离-多普勒二维图多普勒频率范围: %5.9f——%5.9f\n', min(min(Freq_Index)), max(max(Freq_Index)));