nakagamipowerf1216.m

该程序是对扩展目标所使用的模型进行了matlab仿真

clear all;clc;close all;
%定义所需的参数值
R=0:0.1:200;
%注意此处均采用归一化之后的标准正态分布，其方差相应的均为1
delta=1;
m=4;




% 以下公式参见庄铭杰《计算机仿真无线RICE衰落信道的实现方法》对Nakagami-m分布的定义
PDF3=((2/gamma(m))*(m/(2*delta^2))^m)*R.^(2*m-1).*exp(-(R.^2)*m/(2*delta^2));

% figure(2)
% % subplot(3,2,[1 2])
% plot(R,PDF3,'c+-')
% axis([0 20 0 max(PDF3)])

K=sqrt(abs(m^2-m))/(m-sqrt(abs(m^2-m)));%根据m值与莱斯分布中K值的关系，可得当m值变化时，对应的得到莱斯分布中相应的LOS通道幅度与其它散射通道幅度的比值

% rice=rice_glints_powerf(2000,Kcount);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fr=1000;
lamda0=0.05;
sigmav=1.0;
sigmaf=2*sigmav/lamda0;
azi_num=25000;
t = clock;
rand('state',sum(100*clock));
d1=rand(1,azi_num);
rand('state',7*sum(100*clock)+3);%保证了两个随机序列的独立性
d2=rand(1,azi_num);
xi=2*sqrt(-2*log(d1)).*cos(2*pi*d2);%产生服从高斯分布的随机数
xq=2*sqrt(-2*log(d1)).*sin(2*pi*d2);%产生服从高斯分布的随机数
coe_num=12;

Km = 0.2;%常用比例系数，取值范围为0.15-0.25，典型值为0.20
Dimension = 20; %目标的最大尺寸，单位为（米）
Distance =1000;
sita = Km*Dimension/Distance;%严格根据论文所提出的定义式写出如右等式
sita2 = sita^2;
belda=1;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
M=256; %均从图四中的设置部分获得相应的数据
freqx=1:1:M;
% powerfxcorr=exp(-freqx.^2/(2*sigmaf.^2));%其高斯功率谱，傅立叶变换形式不变

f_s =(2*belda*sita^2)./(belda.^2.+freqx.^2);
powerfxcorr = abs(ifft(f_s,M));
[filtercoeff,delta2] = levinson(powerfxcorr,25); 
xxi=filter(1,filtercoeff,xi);
xxq=filter(1,filtercoeff,xq);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

xisigmac=std(xxi);%方差
ximuc=mean(xxi);%均值
yyi=(xxi-ximuc)/xisigmac;
xqsigmac=std(xxq);
xqmuc=mean(xxq);
yyq=(xxq-xqmuc)/xqsigmac;
sigmac=1.2;
yyi=sigmac*yyi;
yyq=sigmac*yyq;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 将两路信号合成为一个莱斯信号
% K=100;
rice=rice_fading(K,yyi,yyq);%无虚部，仅有实部

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


rayleigh=rayleigh_glints_powerf(azi_num);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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=2*sqrt(-2*log(d1)).*cos(2*pi*d2);%产生服从高斯分布的随机数
xq=2*sqrt(-2*log(d1)).*sin(2*pi*d2);%产生服从高斯分布的随机数
coe_num=12;


Km = 0.2;%常用比例系数，取值范围为0.15-0.25，典型值为0.20
Dimension = 20; %目标的最大尺寸，单位为（米）
Distance =1000;
sita = Km*Dimension/Distance;%严格根据论文所提出的定义式写出如右等式
sita2 = sita^2;
belda=1;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
M=256; %均从图四中的设置部分获得相应的数据
freqx=1:1:M;
% powerfxcorr=exp(-freqx.^2/(2*sigmaf.^2));%其高斯功率谱，傅立叶变换形式不变

f_s =(2*belda*sita2)./(belda.^2.+freqx.^2);
powerfxcorr = abs(ifft(f_s,M));
[filtercoeff,delta2] = levinson(powerfxcorr,25); 
xxi=filter(1,filtercoeff,xi);
xxq=filter(1,filtercoeff,xq);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

xisigmac=std(xxi);%方差
ximuc=mean(xxi);%均值
yyi=(xxi-ximuc)/xisigmac;
xqsigmac=std(xxq);
xqmuc=mean(xxq);
yyq=(xxq-xqmuc)/xqsigmac;
sigmac=1.2;
yyi=sigmac*yyi;
yyq=sigmac*yyq;
ydata=yyi+j*yyq;
RayAmp=ydata

coeffient=1.4-0.01^m;
Nakagami=(RayAmp*exp(1-m)+rice*(1-exp(1-m)))*coeffient;

NakagamiTime=(abs(xi+j*xq)*exp(1-m)+rice_fading(K,xi,xq)*(1-exp(1-m)))*coeffient;
subplot(2,1,1);
hist(NakagamiTime,100);
title('Nakagami-m分布的时域序列')
ylabel('直方图');
subplot(2,1,2);
a=rand(1,25000)*2*pi;
hist(a,100)
title('相位分布')
ylabel('直方图');
Nakagamicosttime=etime(clock,t);%产生数据所需要的时间 
ydata=Nakagami;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

num=100;
maxdat=max(abs(ydata));
mindat=min(abs(ydata));
NN=hist(abs(ydata),num);
xpdf1=num*NN/((sum(NN))*(maxdat-mindat));
xaxis1=mindat:(maxdat-mindat)/num:maxdat-(maxdat-mindat)/num;
% th_val=(xaxis1./sigmac.^2).*exp(-xaxis1.^2./(2*sigmac.^2));
figure(3);plot(xaxis1,xpdf1);
hold on,plot(R,PDF3,'k');xlim([0 4])
title('随机序列幅度统计结果');xlabel('幅度');ylabel('概率密度');
 legend('输出随机信号统计结果','标准nakagami分布')
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));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fw = fft(ydata,M);
pw = fw.*conj(fw)/M;
N=M;
pxt(1:N/2-2) = pw(N/2+1:N-2);
pxt(N/2-1:N-4) = pw(3:N/2)/max(pw(3:N/2));
outputdisplay=pxt(N/2-1:N-4);

plot(outputdisplay,'b');
title('输出随机信号功率谱特性');xlabel('频率(Hz)');ylabel('功率');
hold on;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
M=2000
freqx=0:0.5*M;
freqx=freqx*fr/M;
standardf_s =(2*belda*sita2)./(belda.^2.+freqx.^2);
% plot(freqx(1:length(outputdisplay)),standardf_s(1:length(outputdisplay)),'k+-');hold off;
plot(freqx(1:length(outputdisplay)),standardf_s(1:length(outputdisplay)),'k+-');hold off;
axis([0 200 0 1])
legend('输出随机信号功率谱特性','理想功率谱特性')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(5);
plot(freqx,standardf_s);
axis([0 100 0 3e-5])

xlabel('频率');ylabel('功率');title('功率谱密度');