對機載雷達雜波的仿真，當然考慮了最簡單的情況，正側視陣，沒有考慮模糊、陣元誤差、雜波起伏等等因素，仿真出雜波數據，然后畫出功率譜圖、特征譜圖以及距離多普勒圖。
程序如下：
%碩士期間的雜波仿真，准確地說是雜波+干擾
%先產生噪聲背景下的干擾，再考慮機載雷達的情況，將雜波考慮進去.主瓣方向為theta=90°，phi=1°
% 統一時域卷積空域
clc;
clear all;
close all;
%% 參數設置
%-- --系統參數------
Pt=230e3; %峰值發射功率
lamda=0.23; %工作波長
fr=2434.8; %重頻
B=5e6; %接收機帶寬
F=3; %噪聲系數dB
tao=15e-6; %脈沖寬度
K=16; %相參脈沖數
sigma2=1; %噪聲功率 w
k1=1.38e-23; %玻爾茲曼常數
T=290; %溫度
c=3e8; %光速
Ls=10^0.4; %接收機損耗
%-----天線參數-----
%陣元采用余弦方向圖
Gel=4; %陣元增益dB
M=4; %列向陣元數，
N=16; %子陣數
d=lamda/2; %陣元間距
%------載機參數----
H=8000; %載機高度
v=140; %載機速度
%------主瓣指向----
theta=90;
phi=1; %目標俯視角，方位角為theta，錐角為psi 主瓣方向俯仰角1° 列向合成存在相位誤差，需要補償
%------雜波參數----
gamma=-3; %歸一化散射系數dB
theta_c=(0:359)*(pi/180); %雜波塊個數為360
N_RangeSample = 2000; %距離采樣數 看似是人為計算，實際上是計算出來的，一個不模糊區間的距離門個數Ru/deltaR
Nc=360;
deltaTheta = 2*pi/Nc; %方位分辨率/rad
open_phasecompensation_c=exp(-1i*2*pi*d/lamda*sin(1*pi/180)*(0:M-1)); %雜波（俯仰）相位補償算子
open_phasecompensation_c=open_phasecompensation_c/norm(open_phasecompensation_c);
Mtx_synthesisoperator_c=kron(eye(N),open_phasecompensation_c); %雜波俯仰合成算子
%% ------------目標生成------------
ksi_t=10^(3/10); %信噪比3db
a=exp(1i*2*pi*(0:N-1)'*0);
b=exp(1i*2*pi*(0:K-1)'*0.3); %2*vt/(lamda*fr)
Xt=sqrt(sigma2*ksi_t/2)*kron(b,a);
%% -----------干擾生成--------------
Sj=1e-3; %干擾機的有效輻射功率密度（W/Hz） 1000
R_j=150e3; %干擾源所在斜距,考慮干擾源在地面 150km
phi_j=asin(H/R_j); %干擾源的俯仰角
Fr_j=sum(exp(1i*2*pi*d/lamda*(0:M-1).'*(sin(phi_j)-sin(1*pi/180)))); %干擾的接收方向圖 .*Fe_j?
Fr_j=Fr_j/max(max(Fr_j));
Gr_j=M*(10^(Gel/10))*(abs(Fr_j)).^2; %干擾的接收增益
S_j=exp(1i*2*pi*(0:N-1)'*0.05); %16*1
J0=Sj*Gr_j*lamda^2/((4*pi)^2*R_j^2*Ls); %干擾功率
ksi_j=J0/(k1*T*10^(F/10)); %干噪比 分母原本乘了一個B，陳師兄
%% ------------雜波生成----------------
Ru=c/(2*fr); %最大不模糊距離約為61.6公里
d0=4.12*(sqrt(H)+0)*1000; %雷達視距 368.5km
fs = 2.5e6;
deltaR = c/(2*B); % 距離分辨率m 30m
Re = 8.49e6; %地球曲率半徑km
Ft_c= zeros(1,Nc); %發射方向圖
Fr_c= zeros(1,Nc); %接收方向圖
Xcjn=zeros(N*K,N_RangeSample); %雜波+干擾接收數據
Xcn=zeros(N*K,N_RangeSample); %雜波接收數據
S_c=zeros(N*K,Nc);
Gt_c=zeros(1,Nc); %發射增益
Gr_c=zeros(1,Nc); %接收增益
rcjn=zeros(N*K);
rcn=zeros(N*K);
Xc=zeros(N*K,N_RangeSample);
Rc_real=zeros(N*K);
Rcjn_real=zeros(N*K);
for l = 1:N_RangeSample %距離采樣遍歷
Rcell = (l-1)*deltaR;
for m=3 %對於第三個模糊區間 150km
Rc =(m-1)*Ru + Rcell; %不同模糊區間的雜波距離 (m-1)*Ru + Rcell
if Rc>=H && Rc<= d0
sin_psi = (2*H*Re+H^2-Rc^2)/2/Rc/Re; %擦地角 若不嚴格考慮擦地角而將俯仰角視為擦地角，最終無法得出結果！！！！
if sin_psi>0 && (1-sin_psi)>(1e-10)
phi_c=asin(-(2*H*Re+H^2+Rc^2)/2/Rc/(Re+H)); %雜波俯仰角
psi_c=acos(cos(theta_c)*cos(phi_c)); %計算錐角
cos_psi = sqrt(1-sin_psi^2); %擦地角余弦
Fe_c=abs(cos(theta_c+90*pi/180)); %陣元方向圖采用余弦方向圖
Ft_c=sum(exp(1i*2*pi*d/lamda*(0:N-1).'*(cos(psi_c)-cos(90*pi/180)*cos(1*pi/180))))*sum(exp(1i*2*pi*d/lamda*(0:M-1)'*(sin(phi_c)-sin(1*pi/180)))).*Fe_c; %整個陣面總的發射方向圖
Fr_c=sum(exp(1i*2*pi*d/lamda*(0:M-1).'*(sin(phi_c)-sin(1*pi/180)))).*Fe_c;
Ft_c=Ft_c/max(max(Ft_c));
Fr_c=Fr_c/max(max(Fr_c));
Gt_c=M*N*(10^(Gel/10))*(abs(Ft_c)).^2;
Gr_c=M*(10^(Gel/10))*(abs(Fr_c)).^2;
sigma0=10^(gamma/10)*sin_psi; %散射系數,等gama模型
% ClutterPatchArea = Rc*deltaTheta*deltaR/cos_psi;
% sigma = sigma0*ClutterPatchArea;
Ag=c*tao/2*1/cos_psi*Rc*sin(1/360*2*pi); %雜波塊的地面可分辨面積 長*寬
ksi_c=sqrt((Pt*B*tao*Gt_c.*Gr_c*lamda^2*Ag*sigma0)/((4*pi)^3*Rc^4*k1*T*B*10^(F/10)*Ls)); %雜噪比
fd_c=2*cos(phi_c)*cos(theta_c)*v/(lamda*fr);
St_c=exp(1i*(0:K-1)'*2*pi*fd_c); %時域導向矢量
fs_cn=cos(phi_c)*cos(theta_c)*2*pi*d/lamda;
fs_cm=sin(phi_c)*2*pi*d/lamda;
Ss_c=(ones(N,1)*ksi_c).*(Mtx_synthesisoperator_c*kron(exp(1i*(0:N-1)'*fs_cn),exp(1i*(0:M-1)'*fs_cm))); %空域導向矢量,面陣的導向矢量與線陣不同
for i=1:Nc
S_c(:,i)=kron(St_c(:,i),Ss_c(:,i)); %時域卷積空域 NK*360
end
if l == 1000
Rc_real = S_c*S_c'; %每一個距離門上真實(理論)的Rc
end
Xc(:,l)=S_c*(rand(1,Nc)+1i*randn(1,Nc))'; %每一個距離采樣點上的雜波數據
Xn=sqrt(sigma2/2)*(randn(N*K,1)+1i*randn(N*K,1)); %Xn=wgn(N*K,1,0) 0db是統計意義上的值，在非大樣本情況下實際功率不是0db
% Xjt=sqrt(sigma2*ksi_j/2)*(randn(K,1)+1i*randn(K,1)); %干擾的時域數據，全多普勒域
% if l>=1000&&l<1200
Xj=sqrt(sigma2*ksi_j/2).*kron((randn(K,1)+1i*randn(K,1)),S_j); % 時域卷積空域
% else
% Xj=zeros(N*K,1);
% end
Xcjn(:,l) = Xcjn(:,l)+Xc(:,l)+Xn+Xj; %第l個距離單元的空時快拍數據 256*1
Xcn(:,l) = Xcn(:,l)+Xc(:,l)+Xn;
end
end
end
end
%% 畫圖
Rcjn_real=Rc_real+Xj*Xj'+eye(N*K);
Rcn_real=Rc_real+eye(N*K);
for p=1:3*N*K %用3NK個樣本來進行估計
rcjn=rcjn+Xcjn(:,p)*Xcjn(:,p)'; %+eye(N*K)
rcn=rcn+Xcn(:,p)*Xcn(:,p)';
end
Rcjn=rcjn/(3*N*K); %估計值
Rcn=rcn/(3*N*K);
%計算雜波功率譜P並繪圖
fs_num=101;
fd_num=101;
fs_n=linspace(-0.5,0.5,fs_num);
fd_n=linspace(-0.5,0.5,fd_num);
Rcjn_inv=inv(Rcjn);
P=zeros(fd_num,fs_num);
for i=1:fd_num
St=exp(1i*(0:K-1)'*2*pi*fd_n(i));
for j=1:fs_num
Ss=exp(1i*(0:N-1)'*2*pi*fs_n(j)); %此處虛數原先用j，混淆了。。
S=kron(St,Ss); %先時域再空域！！！！
S=S/norm(S); %歸一化使噪聲為功率在圖中顯示為零
P(j,i)=10*log10(abs(1/(S'*Rcjn_inv*S)));
end
end
figure(1);
colormap(jet);
surf(fd_n,fs_n,P);
shading interp;lighting gouraud;colorbar;
title('雜波功率譜');xlabel('fd');ylabel('fs');zlabel('P/dB');hold on; %繪制雜波功率譜圖
%特征譜圖
EigVal = eig(Rcjn);
EigVal = db(sort(abs(EigVal),'descend'))/2;
EigVal2 = eig(Rcn);
EigVal2 = db(sort(abs(EigVal2),'descend'))/2;
figure(2);
plot(EigVal','*-');
title('特征譜');
grid on;
hold on
plot(EigVal2','+-');
legend('雜波+干擾','雜波','Location','NorthWest');
%距離多普勒圖
%先對每個陣元接收數據做DFT，后空域合成
clutter_echo_matrix = zeros(N,K,N_RangeSample);
Vec_Qs = chebwin(N,60).'; %
Vec_Qs = Vec_Qs/(sqrt(Vec_Qs*Vec_Qs')); %歸一化的空域
u_s = Vec_Qs; %1*N的空域合成算子
Vec_Qt = chebwin(K,60).'; %1*K
Vec_Qt = Vec_Qt/(sqrt(Vec_Qt*Vec_Qt'));
fd_t=-0.5:0.01:0.49; %(-K/2:K/2-1)/K;
fdata = zeros(N,100,N_RangeSample);
fdata2 = zeros(100,N_RangeSample);
for l=1:N_RangeSample
clutter_echo_matrix(:,:,l)=reshape(Xcjn(:,l),N,K); %NK*L矩陣變維為N*K*L
for n = 1:N %先對每個陣元數據進行DFT,得到每個陣元數據仍為1*K(頻域)
fdata(n,:,l)=(Vec_Qt.*clutter_echo_matrix(n,:,l)*exp(1i*2*pi*(0:K-1).'*fd_t)); %(1*K).*(1*K)*(K*K)第l個距離門中第n個陣元進行DFT處理以后的回波數據，exp(j*2*pi*(0:K-1).'*fd))為DFT權值
end
fdata2(:,l) = u_s*fdata(:,:,l); %第l個距離門中對每個陣元數據進行DFT后再用空域合成算子進行空域合成
end
x=1:length(fd_t); %1:size(fdata2,1)
y=1:N_RangeSample; %1:size(fdata2,2)
figure(2);colormap(jet);
mesh(x,y,20*log10(abs(fdata2.')));
shading interp;
lighting gouraud;colorbar;
title('處理前距離-多普勒圖');
% axis([1 length(fd_t) 1 N_RangeSample]);
xlabel('多普勒通道');ylabel('距離門');
下面是程序作的三張圖：
1.功率譜圖
