ruili.m

这是雷达幅度符合瑞利分布的杂波仿真程序

clear all;
close all
num=4000; 
fr=1000;
lamda0=0.05;
sigmav=1.0;
sigmaf=2*sigmav/lamda0;
randn('state',sum(100*clock));
d1=randn(1,num);
rand('state',7*sum(100*clock)+3);
d2=randn(1,num);
xi=2*sqrt(-2*log(d1)).*cos(2*pi*d2);
xq=2*sqrt(-2*log(d1)).*cos(2*pi*d2);
coe_num=12; %用傅里叶级数展开法求滤波器系数
for n=0:coe_num; 
coeff(n+1)=2*sigmaf*sqrt(pi)*exp(-4*sigmaf^2*pi^2*n^2/fr^2)/fr;
end
for n=1:2*coe_num+1
if n<=coe_num+1
b(n)=1/2*coeff(coe_num+2-n);
else
b(n)=1/2*coeff(n-coe_num);
end
end
xxi=conv(b,xi);
xxq=conv(b,xq);
xxi=xxi(coe_num*2+1:num+coe_num*2);
xxq=xxq(coe_num*2+1:num+coe_num*2);
xisigmac=std(xxi);
ximuc=mean(xxi);
xxi=(xxi-ximuc)/xisigmac;
xqsigmac=std(xxq);
xqmuc=mean(xxq);
xxq=(xxq-xqmuc)/xqsigmac;
sigmac=1.2;
xxi=sigmac*xxi;
xxq=sigmac*xxq;
xdata=xxi+j*xxq;
figure(1),
subplot(2,1,1),plot(real(xdata));
title('瑞利杂波时域波形--实部')
subplot(2,1,2),plot(imag(xdata));
title('瑞利杂波时域波形--虚部')

Nc=1024;         %   杂波分辨单元
Nm=4096;         %   相干相关随机序列个数
v=0.5;
sita=2*(v+1);
v1=zeros(sita,Nc);
for i=1:sita
    v1(i,1:Nc)=randn(1,Nc)*(v/2);    %  N(0,v/2)
end
a=0.054; %空间相关函数参数
b=0.52;
c=20.69;
t=100;
d=0:t/Nc:t-t/Nc;
Rd=a*exp(-d.^2/b^2)+(1-a)*exp(-d/c);   %空间相关函数
rn=sqrt(Rd);
pr=fft(rn);
Hw=sqrt(pr);
hw=Hw.*exp(j*2*pi*rand(1,Nc));

w=zeros(sita,Nc);
s2=zeros(1,Nc);
for i=1:sita
    w(i,1:Nc)=ifft(fft(v1(i,1:Nc)).*hw);
    s2=s2+w(i,1:Nc).^2;
end
ss=sqrt(s2);
figure(2);
hist(abs(ss),30);%产生gamma分布随机序列
figure(3);
plot(abs(ss));
ps=psd(ss,Nc,'hamming(16)',8);
figure(4);
plot(abs(Rd));%实际相关函数曲线

n=500;%任取一杂波分辨单元
ydata=xdata*ss(n);
figure(5);
hist(abs(ydata),30);%－－－取第500个分辨单元的杂波序列
figure(6);
plot(abs(ydata));
figure(7);
hist(real(ydata),30);
%figure(12); hist(imag(yadta),30);




Tp=80e-6; %---脉冲宽度
B=1e6;%--调频带宽
N=B*Tp;%--脉压比
tao=1/B;%---子脉冲宽度
lamda=0.23;
h=5e3;%--雷达高度
c=3e8;

t=0:tao/num:tao-tao/num;
for i=0:N-1
    K(1,i+1)=sqrt(lamda^2/(64*pi^3*(h+0.5*i*c*tao).^4));%距离衰减因子
end
mq=ss(1:N).*K;%---乘以相应的距离衰减因子作为Nc个雷达分辨单元回波幅度时间序列均值mq
for i=0:N-1
    %H(i+1,1:num)=x* K(1,i+1)*ss(i+1).*exp(j*pi*B*t.^2/Tp);  %---常规脉冲雷达每个距离单元的幅度
    H(i+1,1:num)=(xdata-mean(xdata)*ones(1,num)+mq(i+1)*ones(1,num))* K(1,i+1);
end


%----------线性调频脉压雷达每个距离单元的幅度---------%
n=N;
for k=1:n
    Y(k,1:num)=zeros(1,num);
    for i=1:k
         Y(k,1:num)=Y(k,1:num)+[H(i,k-(i-1):end) H(i,1:k-i)];    %线性调频脉压雷达每个距离单元的幅度
     end
 end
 

figure(8);                  % 各距离波门的回波幅度
hist(abs(Y(1,:)),30);

figure(9);
hist(abs(Y(15,:)),30);

figure(10);      
hist(abs(Y(40,:)),30);

figure(11);
hist(abs(Y(60,:)),30);
%------对每个距离波门杂波信号进行时间域平均------%
for k=1:n
    mqY(k)=mean(abs(Y(k,:)));%----线性调频脉压雷达
end
ycor=xcorr(mqY);
figure(12);
plot(real(ycor(80:159)/max(ycor(80:159))),'b--');
hold on




num=100;
maxdat=max(abs(xdata));
mindat=min(abs(xdata)); 
NN=hist(abs(xdata),num);
xpdf1=num*NN/((sum(NN))*(maxdat-mindat));
xaxis1=mindat:(maxdat-mindat)/num:maxdat-(maxdat-mindat)/num;
th_va1=(xaxis1./sigmac^2).*exp(-xaxis1.^2./(2*sigmac^2));
figure(3),
plot(xaxis1,xpdf1);
hold ;plot(xaxis1,th_va1,':r');
title('杂波的幅度分布');
xlabel('杂波的幅度')
ylabel('概率密度')
signal=ydata;
signal=signal-mean(signal);
figure(4),
M=256;
psd_dat=pburg(real(signal),32,M,fr);
psd_dat=psd_dat/(max(psd_dat));
freqx=0:0.5*M;
freqx=freqx*fr/M;
plot(freqx,psd_dat);
title('杂波功率谱密度');
xlabel('频率/Hz')
ylabel('功率谱密度')

powerf=exp(-freqx.^2/(2*sigmaf.^2));
hold
plot(freqx,powerf,':r');