k_distribution.m

SIRP法相干相关K分布雷达杂波的建模与仿真

%	罗军辉，et al. Matlab 7.0在数字信号处理中的应用. 北京：机械工业出版社，2005:180-184.
%%%%K分布杂波的产生过程
clear all;
close all
%%%%%%%%%%%%雷达和海浪基本参数%%%%%%%%%%%%%%%%%%%%%%%%%
azi_num=2000;   %雷达回波帧数,一帧表示一个重复周期
 fr=1000;       %脉冲重复频率，即对回波进行离散化时的抽样频率
 lamda0=0.05;   %波长
 sigmav=1.0;    %sigmav为海浪速度散布的均方根值，其值只和杂波内部起伏的程度有关，而和雷达工作波长无关
 sigmaf=2*sigmav/lamda0;%杂波功率谱方差，可认为是杂波功率谱密度的3dB带宽。
 %%%%%%%%%%%%高斯白噪声的产生%%%%%%%%%%%%%%%%%%%%%%%%%这部分也可用函数randn()直接产生
 rand('state',sum(100*clock));%不同时间，不同的rand状态，产生不同的随机数序列
 d1=rand(1,azi_num);          %产生均匀分布序列
 rand('state',7*sum(100*clock)+3);
 d2=rand(1,azi_num);
 xi=sqrt(-2*log(d1)).*cos(2*pi*d2);%产生复高斯白噪声的实部，I通道
 xq=sqrt(-2*log(d1)).*sin(2*pi*d2);%产生复高斯白噪声的虚部，Q通道
 %关于高斯分布与均匀分布和指数分布的关系可参考
 %《雷达系统模拟》( (美)米切尔(R.L. Mitchell)著，科学出版社，1982) 第九章：随机数及库函数的产生
 %%%%%%%%%%%%傅立叶级数展开法求滤波器1的系数%%%%%%%%%%%%%%%%%%%%%%%%%
 coe_num=12;      %N=12        
 for n=0:coe_num; 
     coeff(n+1)=2*sigmaf*sqrt(pi)*exp(-4*sigmaf^2*pi^2*n^2/fr^2); %傅立叶系数展开
 end
 %%%%%%书上源程序部分%%%%%%%%%%%%%%%%%
 for n=1:2*coe_num+1;          %一般序列可用共轭对称序列和共轭反对称序列之和表示。这里求的是                 
     if n<=coe_num+1;          %coeff 序列的共轭对称部分b(n),FT[b(n)]=Re(FT[coeff(n)]);具体可参考
         b(n)=1/2*coeff(coe_num+2-n);%丁玉美《数字信号处理》（第二版）P32.
     elseif n==coe_num+2;            %这行与书上源程序略有不同
         b(n)=coeff(n-coe_num);      %求出b(n)后将其向右平移[（2*coe_num+1）-1]/2
     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:azi_num+coe_num*2);
 xxq=xxq(coe_num*2+1:azi_num+coe_num*2);
 %%%%%%%%%修改部分%%%%%%%%%%%%%%%%%%%%
%  b=coeff; 
%  xxi=conv(b,xi);
%  xxq=conv(b,xq);
%  xxi=xxi(coe_num+1:azi_num+coe_num);
%  xxq=xxq(coe_num+1:azi_num+coe_num);
 %%%%%%%%%修改部分%%%%%%%%%%%%%%%%%%%%
 vmuc=2;%niu
 xisigmac=std(xxi);
 ximuc=mean(xxi);
 xxi=(xxi-ximuc)/xisigmac; %均值为零
 xqsigmac=std(xxq);
 xqmuc=mean(xxq);
 xxq=(xxq-xqmuc)/xqsigmac;%均值为零
 xdata=xxi+j*xxq;
 tmpdat=randn(1,azi_num);
 [b,a]=butter(5,0.01);%窄带滤波器
 sk_dat=filter(b,a,tmpdat);
 sk_dat=sk_dat/std(sk_dat);
%下面的程序解非线性方程
max_z=6;
step=0.005;
table_z=0:step:max_z;
table_s=nonline_eq_sirp(table_z,vmuc);%创建非线性表
for n=1:azi_num;
    index=floor(abs(sk_dat(n))/max_z*length(table_z)+1);%计算index
    sk_dat(n)=table_s(index);%引用非线性表
end
ydata=xdata.*sk_dat;
figure, subplot(2,1,1),plot(real(ydata));
 title('K分布杂波时域波形--实部');
 subplot(2,1,2),plot(imag(ydata));
 title('K分布杂波时域波形--虚部');
%求概率密度函数的参数
 num=100;
 maxdat=max(abs(ydata));
 mindat=min(abs(ydata)); 
 NN=hist(abs(ydata),num);
 xpdf1=num*NN/((sum(NN))*(maxdat-mindat));%由直方图求pdf，具体可参考“参考文献”文件夹中
 xaxis1=mindat:(maxdat-mindat)/num:maxdat-(maxdat-mindat)/num;%“频数图与概率密度曲线的问题”这个文件。
 alpha=sqrt(std(xdata).^2./(2*vmuc));
 th_va1=2*((xaxis1/(2*alpha)).^vmuc).*BESSELK((vmuc-1),xaxis1/alpha)./(alpha*gamma(vmuc));%theory_value
 figure,
 subplot(211);
 plot(xaxis1,xpdf1);
 hold ;plot(xaxis1,th_va1,':red');
 title('杂波的幅度分布');
 xlabel('杂波的幅度')
 ylabel('概率密度')
 LEGEND('估计值','理论值');
 signal=ydata;
 signal=signal-mean(signal);
 %%%%用burg法来估计功率谱密度
%  figure(4),
 M=256;
 psd_dat=pburg(real(signal),16,M,fr);
 psd_dat=psd_dat/(max(psd_dat));
 freqx=0:0.5*M;
 freqx=freqx*fr/M;
  subplot(212);
 plot(freqx,psd_dat);
 title('杂波频谱');
 xlabel('频率/Hz')
 ylabel('功率谱密度')
 %%%作出理想的功率谱曲线
 powerf=exp(-freqx.^2/(2*sigmaf.^2));
 hold
 plot(freqx,powerf,':black');
 LEGEND('估计值','理论值')