10种概率密度函数程序.docx

10种概率密度函数程序.docx

《10种概率密度函数程序.docx》由会员分享，可在线阅读，更多相关《10种概率密度函数程序.docx（21页珍藏版）》请在冰豆网上搜索。

10种概率密度函数程序

10种概率密度函数

functionzhifangtu（x,m）

%画数据的直方图,x表示要画的随机数,m表示所要画的条数

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

a=min（x）;

b=max（x）;

l=length（x）;

h=（b-a）/m;

%量化x

x=x/h;

x=ceil（x）;

w=zeros（1,m）;

fori=1:

l

forj=1:

m

if（x（i）==j）

%x（i）落在j的区间上，则w（j）加1

w（j）=w（j）+1;

else

continue

end

end

end

w=w/（h*l）;

z=a:

h:

（b-h）;

bar（z,w）;

title（'直方图'）

functiony=junyun（n）

%0-1的均匀分布，n代表数据量，一般要大于1024

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

y=ones（1,n）;

x=ones（1,n）;

m=100000;

x0=mod（ceil（m*rand（1,1））,m）;

x0=floor（x0/2）;

x0=2*x0+1;

u=11;

x

（1）=x0;

fori=1:

n-1

x（i+1）=u*x（i）+0;

x（i+1）=mod（x（i+1）,m）;

x（i）=x（i）/m;

end

%x（n）单位化

x（n）=x（n）/m;

y=x;

functiony=zhishu（m,n）

%指数分布,m表示指数分布的参数,m不能为0.n表示数据量,n一般要大于1024

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

x=junyun（n）;

fori=1;n

if（x（i）==0）

x（i）=0.0001;

else

continue;

end

end

u=log（x）;

y=-（1/m）*u;

functiony=ruili（m,n）

%瑞利分布,m是瑞利分布的参数,n代表数据量,n一般要大于1024

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

x=junyun（n）;

fori=1:

n

if（x（i）==0）

x（i）=0.0001;

else

continue;

end

end

u=（-2）*log（x）;

y=m*sqrt（u）;

functiony=weibuer（a,b,n）

%韦布尔分布,a,b表示参数,b不能为0.n表示数据量,一般要大于1024

%a=1时,是指数分布

%a=2时,是瑞利分布

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

x=junyun（n）;

fori=1:

n

if（x（i）==0）

x（i）=0.0001;

else

continue;

end

end

u=-log（x）;

y=b*u.^（1/a）;

functiony=swerling（n）

%swelingII分布

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

r=ones（1,n）;

u=junyun（n）;

v=junyun（n）;

fori=1:

n

if（u（i）==0）

u（i）=0.0001;

else

continue

end

end

fori=1:

n

if（u（i）==v（i））

u（i）=u（i）+0.0001

elsecontinue

end

end

t=-2*log（u）;

h=2*pi*v;

x=sqrt（t）.*cos（h）;

z=sqrt（t）.*sin（h）;

y=（r/2）.*（x.^2+z.^2）;

functiony=bernoulli（p,n）

%产生数据量为n的贝努利分布,其中p属于（0-1）之间。

%-----------------------

%

u=junyun（n）;

y=zeros（1,n）;

fori=1:

n

if（u（i）<=p）

y（i）=1;

else

y（i）=0;

end

end

functiony=duishuzhengtai（a,b,n）

%产生对数正态分布，a,b为随机分布的参数，n为数据量

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

x=gaussian（n）;

u=sqrt（b）*x+a;

y=exp（u）;

functiony=kaifeng（m,n）

%产生开丰分布，其中m代表开丰分布的自由度，n表示产生的点数量

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

y=zeros（1,n）;

if（floor（m/2）==m/2）

fori=1:

m/2

[x1,x2]=gaussian（n）;

forj=1:

n

y（j）=x1（j）^2+x2（j）^2+y（j）;

end

end

else

fori=1:

floor（m/2）

[x1,x2]=gaussian（n）;

forj=1:

n

y（j）=x1（j）^2+x2（j）^2+y（j）;

end

end

x=gaussian（n）;

forj=1:

n

y（j）=y（j）+x（j）^2;

end

end

functiony=dajiama（a,b,n）

%产生伽马随机分布的数据，a、b为随机分布的参数，数据量为n

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

%

k=1;

if（a<1）

while（k<=n）

x1=junyun

（1）;

x2=junyun

（1）;

y2=（exp

（1）+a）/exp

（1）*x2;

if（y2>1）

p=-log（（（exp

（1）+a）/exp

（1）-y2）/a）;

if（x1

y（k）=p;

k=k+1;

else

continue;

end

else

p=y2^（1/a）;

if（x1

y（k）=p;

k=k+1;

else

continue;

end

end

end

elseif（a>=1）

while（k<=n）

x1=junyun

（1）;

x2=junyun

（1）;

v=（2*a-1）^（-0.5）*log（x1/（1-x2））;

x=a*exp（v）;

z=x1^2*x2;

w=a-log（4）+（a+sqrt（2*a-1））*v-x;

if（w>=log（z））

y（k）=x;

k=k+1;

else

continue;

end

end

end

y=b*y;

functiony=beitafenbu（a1,a2,n）

%产生贝他分布的随机数，其中a1、a2是贝他分布的参数，n代表数据量

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

x1=dajiama（a1,1,n）;

x2=dajiama（a2,1,n）;

y=x1./（x1+x2）;

function[y1,y2]=gaussian（n）

%产生数据量为n的两个相互独立高斯分布y1、y2

%---------------------------------------

%

k=1;

y1=zeros（1,n）;

y2=zeros（1,n）;

while（k<=n）

u1=junyun

（1）;

u2=junyun

（1）;

v1=2*u1-1;

v2=2*u2-1;

s=v1^2+v2^2;

if（s>=1）

continue;

elseif（s==0）

k=k+1;

else

y1（k）=v1*sqrt（-2*log（s）/s）;

y2（k）=v2*sqrt（-2*log（s）/s）;

k=k+1;

end

end

functiony=canshu（x）;

y=ones（1,2）;

n=length（x）;

y

（1）=sum（x）/n;

z=x-y

（1）;

z=z.^2;

y

（2）=sum（z）/（n-1）;

functiony=correlation（x）

%计算x的自相关函数

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

n=length（x）;

fori=1:

n

x1（i）=x（n+1-i）;

end

y=conv（x,x1）;

二．三种相关杂波

functiony=gaussianpu（x）

%由数据量为n的高斯白噪声产生向量为n，功率谱为高斯型的高斯随机向量

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

N=0:

20;

f=20;

T=1/256;

c=2*f*T*sqrt（pi）*exp（-4*f^2*pi^2*T^2*N.^2）;

n=length（x）;

y=zeros（1,n）;

fork=1:

n

fori=20:

-1:

0

if（（k-i）<=0）

continue;

else

y（k）=y（k）+c（21-i）*x（k-i）;

end

end

fori=20:

40

if（（k-i）<=0）

continue;

else

y（k）=y（k）+c（i-19）*x（k-i）;

end

end

end

y=0.5*y;

functiony=weibuerpu（a,b,n）

%由数据量为n的高斯白噪声产生向量为n，功率谱为高斯型的韦布尔分布的随机向量

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%[z1,z2]=gaussian（n）;

z1=5*z1;

z2=5*z2;

y1=sqrt（b^a/2）*z1;

y2=sqrt（b^a/2）*z2;

x1=gaussianpu（y1）;

x2=gaussianpu（y2）;

x1=sqrt（b^a/2）*x1;

x2=sqrt（b^a/2）*x2;

y=x1.^2+x2.^2;

b=canshu（y）;

y=y-b

（1）;

functiony=duishuzhengtaipu（a,b,n）

%由数据量为n的高斯白噪声产生向量为n，功率谱为高斯型的对数正态随机向量

%a表示标准方差,b表示均值

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

z1=gaussian（n）;

x=gaussianpu（z1）;

y=a*x;

y=exp（y）;

y=b*y;

b=canshu（y）;

y=y-b

（1）;%去掉直流分量

functiony=swerling2pu（n）

%由数据量为n的高斯白噪声产生向量为n，功率谱为高斯型的斯维凌II型随机向量

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

r=6;

[z1,z2]=gaussian（n）;

x1=gaussianpu（z1）;

x2=gaussianpu（z2）;

y=x1.^2+x2.^2;

y=r*y;

b=canshu（y）;

y=y-b

（1）;%去掉直流分量

functiony=kexipu（m,n）

%由数据量为n的高斯白噪声产生向量为n，功率谱为柯西谱的高斯随机向量

wc=2*pi*256;

T0=1/（256*m）;

x=gaussian（n）;

y=zeros（1,n）;

y

（1）=wc*T0*x

（1）;

fori=2:

n

y（i）=wc*T0*x（i）+exp（-wc*T0）*y（i-1）;

end

b=canshu（y）;

%y=y-b

（1）;　　　　　%去掉直流分量　　　

y=conv（y,y）;

y=fft（y）;

y=abs（y）;

i=1:

2*n-1;

plot（i,y）

functionplotpu（x）

%绘出随机数的功率谱密度函数频域的图形。

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

w=fft（x）;

w=abs（w）;

v=2*pi/length（w）;

i=0:

v:

（2*pi-v）;

plot（i,w）;

三.雷达系统仿真

function[t,s,g,f0,fs,f1]=huibo

%产生目标回波信号x,系统噪声y，地物杂波z以及回波p

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

f0=3*10^7;%发射信号频率

w=0;%发射信号初始相位

c=3*10^8;%光速

l=c/f0;%雷达信号波长（载波波长）

R=40000;%目标范围

Vd=200;%雷达与目标之间的径向速度

fd=2*Vd/l;%多普勒频率

Tr=600/f0;%脉冲重复周期

N=10;%雷达脉冲串长度

f1=f0/5;%调频带宽是发射信号频率的1/5

k=10*f1/Tr;

fs=3*f0;%仿真采样频率

Ts=1/fs;

%Tt=2*R/c;

Btar=4*pi*R/l;%

M=floor（Tr*fs）;%一个脉冲重复周期内的采样点数M=600

mt1=floor（2*Tr*fs/5）;%mt1=720

mt2=floor（3*Tr*fs/10）;%mt2=540

mt3=floor（7*Tr*fs/10）;%mt3=1260

mt4=floor（3*Tr*fs/5）;%mt4=1080

mt5=mt1-mt2;%mt5=180

Vgain=6;

s=zeros（M,N）;%回波幅度起伏

form=1:

M

forn=1:

N

v（m,n）=（u（mt1-m）-u（mt2-m））*Vgain*cos（2*pi*f0*（m-mt2）*Ts+2*pi*k*Ts*（m-mt2）/2*（m-mt2）*Ts+2*pi*fd*n*Tr）;

g（m,n）=（u（mt3-m）-u（mt4-m））*Vgain*cos（2*pi*（f0*（m-mt4）*Ts+k*Ts*（m-mt4）/2*（m-mt4）*Ts））;

%u（t）是发射信号包络

end

end

fori=1:

M/10

t（i）=Vgain*cos（2*pi*（k*Ts*（M/10-i）/2*（M/10-i）*Ts））;

end

forn=1:

N

y（1:

M,n）=gaussian（M）';%系统噪声服从高斯分布

end

forn=1:

N

z（1:

M,n）=swerling2pu（M）';%地物杂波服从swerling2型分布

end

s=v+g+y+z;

%s=v+g;

[m,n]=size（s）;

x=zeros（1,m*n）;

q=zeros（1,m*n）;

x=s（:

）';

q=g（:

）';

p=x;

i=0:

length（p）-1;

subplot（2,1,1）;

plot（i,x）,title（'目标回波信号'）;%目标回波信号x

subplot（2,1,2）;

l=canshu（x）;

b=x-l

（1）;

plotpu（b）,title（'目标回波频谱'）;

functiony=gaofang（s,f0,fs,f1）

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%高频放大器

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

[M,N]=size（s）;

Vgain=3;%高放增益

w1=（f0-4*f1/3）*2*pi/fs;%w1=26pi/45

w3=（f0-f1/3）*2*pi/fs;%w3=29pi/45

w2=（f0+7*f1/3）*2*pi/fs;%w2=37pi/45

w4=（f0+4*f1/3）*2*pi/fs;%w4=34pi/45

th=min（（w3-w1）,（w2-w4））;%w3-w1=pi/15,w2-w4=pi/15

M1=ceil（6.6*pi/th）+1;%M1=99

%n=[0:

M1-1];

w5=（w1+w3）/2;w6=（w4+w2）/2;%w5=11pi/18,w6=71pi/90

h=wide（w6,M1）-wide（w5,M1）;

w=（hamming（M1））';

h=h.*w;

L=length（h）;

s=Vgain*s;

z=zeros（M+L-1,N）;

y=zeros（M,N）;

K=ceil（L/2）;

fori=1:

N

forj=1:

10

z（1+（j-1）*M/10:

j*M/10+L-1,i）=conv（h,s（1+（j-1）*M/10:

j*M/10,i）'）';

y（1+（j-1）*M/10:

j*M/10,i）=z（1+（j-1）*M/10+K:

j*M/10+K,i）;

end

end

w=y（:

）';

subplot（3,1,1）;

i=0:

length（w）-1;

plot（i,w）,title（'高放'）;

subplot（3,1,2）;

plotpu（w）,title（'频谱'）;

subplot（3,1,3）;

plotpu（h）,title（'幅频特性'）;

functiony=hunpin（s,f0,fs,f1）

%%%%%%%%%%%%%%%%%%%%%%%%%%进行混频，输出为中频信号

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

[M,N]=size（s）;

i=0:

M-1;

z1=cos（2*pi*2*f0*i/3/fs）;

y1=zeros（M,N）;

fori=1:

N

y1（1:

M,i）=（z1.*s（1:

M,i）'）';

end

w1=（f0-5*f1）*2*pi/（3*fs）;%w1=pi/9

w3=（f0-2*f1）*2*pi/（3*fs）;%w3=1.6pi/9

w2=（f0+8*f1）*2*pi/（3*fs）;%w2=0.4pi

w4=（f0+5*f1）*2*pi/（3*fs）;%w4=pi/3

th=min（（w3-w1）,（w2-w4））;%w3-w1=0.2pi/3,w2-w4=0.2pi/3

M1=ceil（6.6*pi/th）+1;%M1=99

%n=[0:

M1-1];

w5=（w1+w3）/2;w6=（w4+w2）/2;%w5=1.3pi/9,w6=1.1pi/3

h=wide（w6,M1）-wide（w5,M1）;

wth=（hamming（M1））';

h=h.*wth;

L=length（h）;

z=zeros（M+L-1,N）;

y=zeros（M,N）;

K=ceil（L/2）;

fori=1:

N

forj=1:

10

z（1+（j-1）*M/10:

j*M/10+L-1,i）=conv（h,s（1+（j-1）*M/10:

j*M/10,i）'）';

y（1+（j-1）*M/10:

j*M/10,i）=z（1+（j-1）*M/10+K:

j*M/10+K,i）;

end

end

w=y（:

）';

v=y1（:

）';

subplot（4,1,1）;

i=0:

length（w）-1;

plot（i,w）,title（'混频'）;

subplot（4,1,2）;

plotpu（v）,title（'频谱'）;

subplot（4,1,3）;

plotpu（w）;

subplot（4,1,4）;

plotpu（h）,title（'幅频特性'）;

functiony=zhongfang（s,f0,fs,f1）

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%中频放大器

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

[M,N]=size（s）;

Vgain=5;%中放增益

w1=（f0-5*f1）*2*pi/（3*fs）;%w1=pi/9

w3=（f0-2*f1）*2*pi/（3*fs）;%w3=1.6pi/9

w2=（f0+8*f1）*2*pi/（3*fs）;%w2=0.4pi

w4=（f0+5*f1）*2*pi/（3*fs）;%w4=pi/3

th=min（（w3-w1）,（w2-w4））;%w3-w1=0.2pi/3,w2-w4=0.2pi/3

M1=ceil（6.6*pi/th）+1;%M1=99

%n=[0:

M1-1];

w5=（w1+w3）/2;w6=（w4+w2）/2;%w5=1.3pi/9,w6=1.1pi/3

h=wide（w6,M1）-wide（w5,M1）;

w=（hamming（M1））';

h=h.*w;

L=length（h）;

s=Vgain*s;

z=zeros（M+L-1,N）;

y=zeros（M,N）;

K=ceil（L/2）;

fori=1:

N

forj=1:

10

z（1+（j-1）*M/10:

j*M/10+L-1,i）=conv（h,s（1+（j-1）*M/10:

j*M/10,i）'）';

y（1+（j-1）*M/10:

j*M/10,i）=z（1+（j-1）*M/10+K:

j*M/10+K,i）;

end

end

w=y（:

）';

subplot（3,1,1）;

i=0:

length（w）-1;

plot（i,w）,title（'中放'）;

subplot（3,1,2）;

plotpu（w）,title（'频谱'）;

subplot（3,1,3）;

plotpu（h）,title（'幅频特性'）;

function[I,Q]=xiangganjianbo（s,fs,f0,f1）