matlab例程

开发平台：
Matlab

ca3.m：源码内容
							%%%%%%%%%%相干韦布尔杂波分布%%%%%%%%%%
azi_num=1000;
fr=1000; %雷达发射信号重频
lamda0=0.05;
sigmav=1.0;
sigmaf=2*sigmav/lamda0;
rand('state',sum(100*clock));
d1=rand(1,azi_num);
rand('state',7*sum(100*clock)+3);
d2=rand(1,azi_num);
xi=1*(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);
xsigmac=std(xxi);
xmuc=mean(xxi);
yyi=(xxi-xmuc)/xsigmac;
xq=2*(sqrt(-2*log(d1)).*cos(2*pi*d2));
xxq=conv(b,xq);
xxi=xxi(coe_num*2+1:azi_num+coe_num*2);
xxq=xxq(coe_num*2+1:azi_num+coe_num*2);
xisigmac=std(xxi);
ximuc=mean(xxi);
yyi1=(xxi-ximuc)/xisigmac;
xqsigmac=std(xxq);
xqmuc=mean(xxq);
yyq=(xxq-xqmuc)/xqsigmac;
p=1.5;
q=2.2;
sigmac=sqrt((q.^p)/2);
yyi=sigmac*yyi;
yyq=sigmac*yyq;
xdata_1=(yyi1.*yyi1+yyq.*yyq).^(1/p);
%%%%%%%%%%参数设定%%%%%%%%%%
fr=1000; %雷达发射信号重频
a=1;
fs=4e5; %采样频率
ts=1/fs;
RCS=140;
k=8.6*10E-5 ;
sigma0=10*log10(RCS); %杂波后散射系数
thetaE=2; %天线3db海拔带宽
thetaA=1; %天线3db方位角带宽
SL=-20; %天线旁瓣水平
range_n=200;
range_obj1=50;
range_obj2=75;
range_obj3=100;
range=linspace(1,range_n,1000);
hr=3; %雷达高度
ht=3; %目标高度
pt=75; %雷达峰值功率
rang_obj1=range_obj1*1000; %目标距离
rang_obj2=range_obj2*1000;
rang_obj3=range_obj3*1000;
t0=290; %有效噪声温度
l=10; %雷达损耗
ant_id=1; 
clight=3.e8; %光速
re=6371000; %地球半径
t=1/fr;
F=1;
B=1e6;
%%%%%%%%%%产生杂波回波%%%%%%%%%%
lamda=clight/fr;
thetaA_deg=thetaA;
thetaE_deg=thetaE;
thetaA=thetaA_deg*pi/180;
thetaE=thetaE_deg*pi/180;
rh=sqrt(8.0*hr*re/3.);
SLv=10.0^(SL/10);
sigma0v=10.0^(sigma0/10);
tau=1/B;
deltar=clight*tau/2.;
range_m=1000.*range;
thetar=asin(hr./range_m);
thetae=asin((ht-hr)./range_m);
propag_atten=1.+((range_m./rh).^4);
Rg=range_m.*cos(thetar);
deltaRg=deltar.*cos(thetar);
theta_sum=thetae+thetar;
if(ant_id==0)
ant_arg=(2.78*theta_sum)./(pi*thetaE);
gain=(sinc(ant_arg)).^2;
else
gain=exp(-2.776.*(theta_sum./thetaE).^2);
end
sigmac=(sigma0v.*Rg.*deltaRg).*(pi*SLv*SLv+thetaA.*gain.^2)./propag_atten;
sigmac=10*log10(sigmac);
%%%%%%%%%%回波幅度计算%%%%%%
G=150;
A1=sqrt(pt/(4*pi)^3*l)*G^2*lamda*sqrt(RCS)/rang_obj1^2;
A2=sqrt(pt/(4*pi)^3*l)*G^2*lamda*sqrt(RCS)/rang_obj2^2;
A3=sqrt(pt/(4*pi)^3*l)*G^2*lamda*sqrt(RCS)/rang_obj3^2;
%%%%%%%%%%计算目标%%%%%%%%%%
s_pc_1=[zeros(1,rang_obj1/range_n),A1,A1,A1,A1,A1,A1,A1,A1,A1,A1,zeros(1,((rang_obj2-rang_obj1)/range_n)-10),A2,A2,A2,A2,A2,A2,A2,A2,A2,A2,zeros(1,((rang_obj3-rang_obj2)/range_n)-10),A3,A3,A3,A3,A3,A3,A3,A3,A3,A3,zeros(1,(1000-rang_obj3/range_n)-10)]; 
         s_pc=xdata_1+s_pc_1+sigmac;
         Pn=sum(sum((xdata_1).^2));
         snr=pt/Pn;
        % SNR=pt*G^2*lamda^2*RCS/((4*pi)^2*(rang_obj)^4*k*t0*B*F*l);
y=log(range_m);
c=0.3;
logG=c*y+log( s_pc);
figure(1)
plot(s_pc)
grid
%figure(2)
%plot(logG)
%grid
N=1000;
%第1点恒虚警处理时噪声均值由其后面的n点的噪声决定
cfar_ca_result=zeros(1,N);
cfar_ca_result(1,1)=s_pc(1,1)/(sqrt(2)/pi*mean(s_pc(1,2:17)));%估计值要乘一个门限加权系数T？
%第2点到第n点恒虚警处理时噪声均值由其前面和后面的n点的噪声共同决定
for i=2:16
    noise_mean=sqrt(2)/pi*(mean(s_pc(1,1:i-1))+mean(s_pc(1,i+1:i+16)))/2;
    cfar_ca_result(1,i)=s_pc(1,i)/noise_mean;
end
%正常数据点恒虚警处理时噪声均值由其前面和后面各n点的噪声共同决定
    for i=17:N-17
         noise_mean=sqrt(2)/pi*(mean(s_pc(1,i-16:i-1))+mean(s_pc(1,i+1:i+16)))/2;
          cfar_ca_result(1,i)=s_pc(1,i)/noise_mean;
    end
%倒数n点到倒数第2点恒虚警处理时噪声均值由其前面n点和后面点的噪声共同决定
for i=N-16:N-1
    noise_mean=sqrt(2)/pi*(mean(s_pc(1,i-16:i-1))+mean(s_pc(1,i+1:N)))/2;
   cfar_ca_result(1,i)=s_pc(1,i)/noise_mean;
end
%最后一点的恒虚警处理的噪声均值由其前面的n点的噪声决定
noise_mean=sqrt(2)/pi*mean(s_pc(1,N-16:N-1));
cfar_ca_result(1,N)=s_pc(1,N)/noise_mean;
    
 figure      
plot(cfar_ca_result,'g');grid
xlabel('N'),ylabel('信号幅度'),title('采用CA处理的结果');

						