哈工大随机信号实验报告

H a r b i n I n s t i t u t e o f T e c h n o l o g y
实验报告
课程名称：随机信号分析
院系：电信学院
班级： 1205201 姓名：
学号：
指导教师：郑薇
实验时间： 2014年 11月
哈尔滨工业大学
实验一 各种分布随机数的产生
一、 实验目的
在很多系统仿真的过程中，需要产生不同分布的随机变量。

利用计算机可以很方便地产生不同分布的随机变量，各种分布的随机变量的基础是均匀分布的随机变量。

有了均匀分布的随机变量，就可以用函数变换等方法得到其他分布的随机变量。

二、 实验内容
产生均匀分布的随机数、高斯分布的随机数和其它分布的随机数。

三、 实验原理
1. 均匀分布随机数的产生原理
产生伪随机数的一种实用方法是同余法，它利用同余运算递推产生伪随机数序列。

最简单的方法是加同余法
)(m od 1M c y y n n +=+
M
y x n n 1
1++=
为了保证产生的伪随机数能在[0,1]内均匀分布，需要M 为正整数，此外常数c 和初值y0亦为正整数。

加同余法虽然简单，但产生的伪随机数效果不好。

另一种同余法为乘同余法，它需要两次乘法才能产生一个[0,1]上均匀分布的随机数
)(m od 1M ay y n
n =+
M y x n n 1
1++=
式中，a 为正整数。

用加法和乘法完成递推运算的称为混合同余法，即
)(m od 1M c ay y n n +=+
M
y x n n 1
1++=
用混合同余法产生的伪随机数具有较好的特性，一些程序库中都有成熟的程序供选择。

常用的计算语言如Basic 、C 和Matlab 都有产生均匀分布随机数的函数可以调用，只是用各种编程语言对应的函数产生的均匀分布随机数的范围不同，有的函数可能还需要提供种子或初始化。

Matlab 提供的函数rand()可以产生一个在[0,1]区间分布的随机数，rand(2,4)则可以产生一个在[0,1]区间分布的随机数矩阵，矩阵为2行4列。

Matlab 提供的另一个产生随机数的函数是random('unif',a,b,N,M)，unif 表示均匀分布，a 和b 是均匀分布区间的上下界，N 和M 分别是矩阵的行和列。

2. 随机变量的仿真
根据随机变量函数变换的原理，如果能将两个分布之间的函数关系用显式表达，那么就可以利用一种分布的随机变量通过变换得到另一种分布的随机变量。

若X 是分布函数为F(x)的随机变量，且分布函数F(x)为严格单调升函数，令Y=F(X)，则Y 必为在[0,1]上均匀分布的随机变量。

反之，若Y 是在[0,1]上均匀分布的随机变量，那么
)(1
Y F X X -=
即是分布函数为FX(x)的随机变量。

式中F X -⋅1
()为F X ()⋅的反函数。

这样，欲求
某个分布的随机变量，先产生在[0,1]区间上的均匀分布随机数，再经上式变换，便可求得所需分布的随机数。

3. 高斯分布随机数的仿真
广泛应用的有两种产生高斯随机数的方法，一种是变换法，一种是近似法。

如果X1,X2是两个互相独立的均匀分布随机数，那么下式给出的Y1,Y2
⎪⎩⎪⎨
⎧+-=+-=m X X Y m
X X Y )π2sin(ln 2)π2cos(ln 2212
211σσ 便是数学期望为m 2
σ
变换法。

另外一种产生高斯随机数的方法是近似法。

在学习中心极限定理时，曾提到n 个在[0,1]区间上均匀分布的互相独立随机变量Xi (i=1,2…,n)，当n 足够大时，其和的分布接近高斯分布。

当然，只要n 不是无穷大，这个高斯分布是近似的。

由于近似法避免了开方和三角函数运算，计算量大大降低。

当精度要求不太高时，近似法还是具有很大应用价值的。

4. 各种分布随机数的仿真
有了高斯随机变量的仿真方法，就可以构成与高斯变量有关的其他分布随机
2χ分布随机变量。

四、 实验过程和结果分析
1、均匀分布、高斯分布随机数的产生与仿真
① 思路：利用已知matlab 函数直接产生随机数。

② 程序
x=random('unif',5,10,1,1000);%产生1000个服从于U （5,10）的随机数
y=random('normal',0,1,1,3000);%产生3000个服从于N （0,1）的随机数 subplot(211),plot(x);title('均匀分布随机数') subplot(212),plot(y);title('高斯分布随机数')
③仿真图形
④分析：产生的随机数呈现中间多，两头少的趋势，普遍集中于期望附近。

2、瑞利分布、指数分布及2 分布随机数的产生与仿真
①思路：利用已知matlab函数的变换加和产生随机数。

②程序
N=5000;
G1=random('Normal',0,1,1,N);
G2=random('Normal',0,1,1,N);
G3=random('Normal',0,1,1,N);
G4=random('Normal',0,1,1,N);
R=sqrt(G1.*G1+G2.*G2);
E=G1.*G1+G2.*G2;
X=G1.*G1+G2.*G2+G3.*G3+G4.*G4;
subplot(311);plot(R);title('瑞利分布随机数')
subplot(312);plot(E);title('指数分布随机数')
subplot(313);plot(X);title('4自由度x^2分布随机数')
③仿真图形
④分析：经变换后的随机数生成规律满足所需要的随机数如指数、瑞利分布。

实验二随机变量检验
一、实验目的
随机数产生之后，必须对它的统计特性做严格的检验。

一般来讲，统计特性的检验包括参数检验、均匀性检验和独立性检验等。

事实上，我们如果在二阶矩范围内讨论随机信号，那么参数检验只对产生的随机数一、二阶矩进行检验。

我们可以把产生的随机数序列作为一个随机变量，也可以看成随机过程中的一个样本函数。

不论是随机变量还是随机过程的样本函数，都会遇到求其数字特征的情况，有时需要计算随机变量的概率密度直方图等。

二、 实验内容
1. 对实验一产生的各种分布的随机数进行均值和方差的检验。

2. 对实验一产生的各种分布的随机数概率分布进行统计，并在计算机屏幕上显示实际统计的概率密度直方图。

三、 实验原理
1. 均值的计算
在实际计算时，如果平稳随机序列满足各态历经性，则统计均值可用时间均值代替。

这样，在计算统计均值时，并不需要大量样本函数的集合，只需对一个样本函数求时间平均即可。

甚至有时也不需要计算N →∞时的极限，况且也不可能。

通常的做法是取一个有限的、计算系统能够承受的N 求时间均值和时间方差。

根据强调计算速度或精度的不同，可选择不同的算法。

设随机数序列{N x x x ,,,21 }，一种计算均值的方法是直接计算下式
∑==
N
n n x N m 1
1 式中，xn 为随机数序列中的第n 个随机数。

另一种方法是利用递推算法，第n 次迭代的均值也亦即前n 个随机数的均值为
)(1
11111----+=+-=
n n n n n n m x n
m x n m n n m 迭代结束后，便得到随机数序列的均值
m m N =
递推算法的优点是可以实时计算均值，这种方法常用在实时获取数据的场合。

当数据量较大时，为防止计算误差的积累，也可采用
)(111
1m x N m m n N
n -+=∑=
式中，m1是取一小部分随机数计算的均值。

2. 方差的计算
计算方差也分为直接法和递推法。

仿照均值的做法
21
2)(1m x N N
n n -=
∑=σ 21
2
2
1m x N N n n -=∑=σ
方差的递推算法需要同时递推均值和方差
m m n
x m n n n n =+
---111
() ])(1
[121212---+-=
n n n n m x n
n n σσ 迭代结束后，得到随机数序列的方差为
22N
σσ= 其它矩函数也可用类似的方法得到。

3. 统计随机数的概率密度直方图
假定被统计的序列)(n x 的最大值和最小值分别为a 和b 。

将),(b a 区间等分M （M 应与被统计的序列)(n x 的个数N 相适应，否则统计效果不好。

）份后的区间为
))(,(M a b a a -+
，))
(*2,)((M
a b a M a b a -+-+，… , )*)(*2,)1)(((M
i
a b a M i a b a -+--+，… , ),)1)(((b M M a b a --+。

用)(i f ，表示序列)(n x 的值落在)*)(*2,)1)(((M i
a b a M i a b a -+--+
区间里的个数，统计序列)(n x 的值在各个区间的个数)(i f ，1,,2,0-=M i ，则)(i f 就粗略地反映了随机序列的概率密度的情况。

用图形方式显示出来就是随机数的概率密度直方图。

四、 实验过程和结果分析
1、均匀分布、高斯分布随机数均值、方差的检验及概率密度直方图
① 思路：随机产生一组数算出均值、方差，与理论值比较。

② 程序
x=random('unif',5,10,1,20000);%产生20000个服从于U （5,10）的随机数 y=random('normal',0,1,1,3000);%产生3000个服从于N （0,1）的随机数 subplot(211),hist(x,5:0.1:10);title('均匀分布随机数') subplot(212),hist(y,-3:0.1:3);title('高斯分布随机数') m1=mean(x)v1=var(x) m2=mean(y) v2=var(y)
③ 仿真图形
2、瑞利、指数、2 分布随机数均值、方差的检验及概率密度直方图
①思路：随机产生一组数算出均值、方差，与理论值比较。

②程序
N=5000;
G1=random('Normal',0,1,1,N);
G2=random('Normal',0,1,1,N);
G3=random('Normal',0,1,1,N);
G4=random('Normal',0,1,1,N);
R=sqrt(G1.*G1+G2.*G2);
E=G1.*G1+G2.*G2;
X=G1.*G1+G2.*G2+G3.*G3+G4.*G4;
subplot(311);hist(R,0:0.05:4);title('瑞利分布随机数')
subplot(312);hist(E,0:0.1:15);title('指数分布随机数')
subplot(313);hist(X,0:0.2:21);title('4自由度x^2分布随机数')
m1=mean(R)
v1=var(R)
m2=mean(E)
v2=var(E)
m3=mean(X)
v3=var(X)
③仿真图形
实验三中心极限定理的验证
一、实验目的
利用计算机产生均匀分布的随机数。

对相互独立的均匀分布的随机变量做和，可以很直观看到均匀分布的随机变量的和，随着做和次数的增加分布情况的变化，通过实验对中心极限定理的进行验证。

二、 实验内容
产生多组[0,1]区间上的均匀分布的随机数序列，各序列的对应元素做和，够成的和序列再进行随机数的概率密度直方图的统计，并作图显示。

三、 实验原理
如果n 个独立随机变量的分布是相同的，并且具有有限的数学期望和方差，当n 无穷大时，它们之和的分布趋近于高斯分布。

这就是中心极限定理中的一个定理。

我们以均匀分布为例，来解释这个定理。

若n 个随机变量Xi (i=1,2,…,n)都为[0,1]区间上的均匀分布的随机变量，且互相独立，当n 足够大时，其和∑==n
i i
X Y 1的分布接近高斯分布。

四、 实验过程和结果分析
① 思路：产生n 个[0,1]区间上的均匀分布的随机数序列并作和，n 取三组值，此外再产生一个高斯分布随机数，对四组随机数进行比较。

② 程序
X1=random('unif',0,1,1,2000); X2=random('unif',0,1,1,2000); X3=random('unif',0,1,1,2000); X4=random('unif',0,1,1,2000); X5=random('unif',0,1,1,2000); X6=random('unif',0,1,1,2000); G=random('normal',0,1,1,2000);
Y1=X1+X2+X3; Y2=X1+X2+X3+X4+X5+X6; subplot(411);hist(X1,0:0.05:2); subplot(412);hist(Y1,0:0.05:4); subplot(413);hist(Y2,0:0.05:6); subplot(414);hist(G,-3:0.05:3);
③ 仿真图形
④分析：
随n取值的增大，均匀分布随机序列求和的图形越发接近于高斯分布。

实验四自相关函数的计算
一、实验目的
在随机信号理论中，自相关函数是非常重要的概念。

在实际系统仿真中也会经常计算自相关函数。

通过本试验学生可以亲自动手计算自相关函数，加深对概念的理解，并增强实际动手能力。

二、实验内容
用一个数学期望为零和非零，方差为某值的高斯分布随机数，作为样本序列求自相关函数的估值，并用图形显示。

三、 实验原理
在实际应用中，我们可以把产生的随机数序列看成随机过程中的一个样本函数。

如果平稳随机序列满足各态历经性，则统计自相关序列可用时间自相关序列代替。

当数据的样本数有限时，也只能用有限个数据来估计时间自相关序列，统计自相关序列的估值。

若各态历经序列X(n)的一个样本有N 个数据
)}1(,),1(),0({-N x x x ，由于实序列自相关序列是对称的，自相关函数的估值为
|)|()(1)(ˆ1||0
m n x n x N
m R
m N n +=∑--=
四、 实验过程和结果分析
① 思路：利用matlab 函数直接产生所需自相关函数。

② 程序
N=500;
x1=random('normal',0,1,1,N); Rx1=xcorr(x1,'biased'); m1=-N+1:N-1;
subplot(211),plot(m1,Rx1); xlabel('m1')
ylabel('Rx1(m1)')
title('均值为0，方差为1的高斯分布的自相关函数'); axis([-N N -0.5 1.5]);
x2=random('normal',1,1,1,N); Rx2=xcorr(x2,'biased'); m2=-N+1:N-1;
subplot(212),plot(m2,Rx2); xlabel('m2')
ylabel('Rx2(m2)')
title('均值为1，方差为1的高斯分布的自相关函数'); axis([-N N -0.5 2]);
③ 仿真图形
④分析：
分别生成均值为0和1，方差为1的高斯随机数，由图形可以明显看出两者自相关函数的差异。

实验五功率谱密度
一、实验目的
在随机信号理论中，功率谱密度和自相关函数一样都是非常重要的概念。

在实际系统仿真中也会经常计算。

通过本试验学生可以亲自动手，加深对概念的理解，并增强实际动手能力。

二、实验内容
用实验四计算出的自相关函数的估值，作为样本序列求功率谱密度的估值，并用图形显示。

三、 实验原理
一般把平稳随机序列的功率谱定义为自相关序列的傅里叶变换。

如果自相关
X(n)的
功率谱与自相关序列的关系为
m X m X m R S ωωj e )()(-∞
-∞
=∑
=
ωωωd e )(π21)(j π
π
m
X X S m R ⎰-= 与实平稳过程一样，实平稳序列的功率谱也是非负偶函数，即
)
()(0)(ωωω-=≥X X X S S S
可以证明，功率谱还可表示为
}e )(1
21
{
)(2
j lim
n
N
N
n N X n X N E S ωω--=∞→∑
+=
当X(n)为各态历经序列时，可去掉上式中的统计均值计算，将随机序列X(n)用它的一个样本序列x(n)代替。

在实际应用中，由于一个样本序列的可用数据个数N 有限，功率谱密度也只能是估计值
22
j 1
|)(|1
e )(1)(ˆωωωX N
n x N
S n
N n X
=
=--=∑
式中，X()是x(n)的傅里叶变换。

这是比较简单的一种估计方法，这种功率谱密度的估计方法称为周期图方法。

如果直接利用数据样本做离散傅里叶变换，可得到X()的离散值。

由于这种方法可借助FFT 算法实现，所以得到了广泛的应用。

四、 实验过程和结果分析
① 思路：利用实验四中的自相关函数与功率谱密度的关系产生或用matlab 函 数直接产生所需功率谱密度。

② 程序
N=500;
x1=random('normal',0,1,1,N); Sx1=abs(fft(x1).^2)/N;
subplot(211),plot(10*log10(Sx1)); axis([0 300 -40 10]) xlabel('f/Hz')
ylabel('Sx1/dB')
title('均值为0，方差为1的高斯分布的功率谱密度'); x2=random('normal',1,1,1,N);
Sx2=periodogram(x2);
subplot(212),plot(10*log10(Sx2));
xlabel('f/Hz')
ylabel('Sx2/dB')
title('均值为1，方差为1的高斯分布的功率谱密度');
③仿真图形
④分析：由波形知，两种方法均可产生功率谱密度。

实验六 随机信号经过线性系统前后信号仿真
一、 实验目的
系统仿真是信号仿真处理的一个重要部分，通过该实验要求学生掌握系统仿真的基本概念，并学会系统的仿真方法。

二、 实验内容
仿真信号和加性噪声经过各种系统前后的自相关函数和功率谱密度并图示。

三、 实验原理
需要先仿真一个指定系统，再根据需要仿真输入的随机信号，然后使这个随机信号通过指定的系统。

通过对实际系统建模, 计算机可以对很多系统进行仿真。

在信号处理中，一般将线性系统分解为一个全通放大器（或衰减器）和一个特定频率响应的滤波器。

由于全通放大器可以用一个常数代替，因此线性系统的仿真往往只需设计一个数字滤波器。

滤波器设计可采用MATLAB 提供的函数，也可利用相应的方法自行设计。

MATLAB 提供了多个设计滤波器的函数，可以很方便地设计低通、带通、高通、多带通、带阻滤波器。

四、 实验过程和结果分析
① 思路：实验产生的随机信号
1122()
cos(2)3cos(2)()X n f t f t N n πφπφ=++++
，其中1φ、2φ为[0,π2]内均匀分布的随机变量，()n N 是数学期望为0、方差为1的高斯白噪声，通过各种系统得到所需仿真图形。

② 程序
1、X （n ）信号的自相关函数及功率谱密度
N=2000;fs=400;
Nn=random('normal',0,1,1,N); t=(0:N-1)/fs;
fi=random('unif',0,1,1,2)*2*pi;
xn=cos(2*pi*50*t+fi(1))+3*cos(2*pi*150*t+fi(2))+Nn; Rx=xcorr(xn,'biased'); m=-N+1:N-1;
Sx=abs(fft(xn).^2)/N; f=(-N/2:N/2-1)*fs/N;
subplot(211),plot(m,Rx); xlabel('m')
ylabel('Rx(m)')
title('xn 的自相关函数');
subplot(212),plot(f,fftshift(10*log10(Sx(1:N)))); xlabel('f/Hz') ylabel('Sx/dB')
title('xn 的功率谱密度');
2、X （n ）通过低通滤波器
N=2000;fs=400;
Nn=random('normal',0,1,1,N); t=(0:N-1)/fs;
fi=random('unif',0,1,1,2)*2*pi;
xn=cos(2*pi*50*t+fi(1))+3*cos(2*pi*150*t+fi(2))+Nn; h=fir1(100,0.4); H=fft(h,2*N); H2=abs(H).^2;
Rx=xcorr(xn,'biased');
Sx=abs(fftshift(fft(xn,2*N)).^2)/(2*N); Sy=Sx.*H2;
Ry=fftshift(ifft(Sy)); f=(-N:N-1)*fs/(2*N); m=(-N:N-1)/N*(N/2000);
subplot(311);plot((-N:N-1)/N,fftshift(abs(H2(1:2*N)))); title('低通滤波器');
subplot(312),plot(m,Ry); xlabel('m')
ylabel('Ry(m)')
title('xn经低通滤波器的自相关函数');
subplot(313),plot(f,fftshift(10*log10(Sy(1:2*N))));
axis([-200 200 -20 50]);
xlabel('f/Hz')
ylabel('Sy/dB')
title('xn经低通滤波器的功率谱密度');
3、带通
N=2000;fs=400;
Nn=random('normal',0,1,1,N);
t=(0:N-1)/fs;
fi=random('unif',0,1,1,2)*2*pi;
xn=cos(2*pi*50*t+fi(1))+3*cos(2*pi*150*t+fi(2))+Nn;
h=fir1(100,[0.1 0.5]);
H=fft(h,2*N);
H2=abs(H).^2;
Rx=xcorr(xn,'biased');
Sx=abs(fftshift(fft(xn,2*N)).^2)/(2*N);
Sy=Sx.*H2;
Ry=fftshift(ifft(Sy));
f=(-N:N-1)*fs/(2*N);
m=(-N:N-1)/N*(N/2000);
subplot(311);plot((-N:N-1)/N,fftshift(abs(H2(1:2*N))));
title('带通滤波器');
subplot(312),plot(m,Ry);
xlabel('m')
ylabel('Ry(m)')
title('xn经带通通滤波器的自相关函数');
subplot(313),plot(f,fftshift(10*log10(Sy(1:2*N))));
axis([-200 200 -20 50]);
xlabel('f/Hz')
ylabel('Sy/dB')
title('xn经带通滤波器的功率谱密度');
4、高通
N=2000;fs=400;
Nn=random('normal',0,1,1,N);
t=(0:N-1)/fs;
fi=random('unif',0,1,1,2)*2*pi;
xn=cos(2*pi*50*t+fi(1))+3*cos(2*pi*150*t+fi(2))+Nn;
h=fir1(100,0.6,'high');
H=fft(h,2*N);
H2=abs(H).^2;
Rx=xcorr(xn,'biased');
Sx=abs(fftshift(fft(xn,2*N)).^2)/(2*N);
Sy=Sx.*H2;
Ry=fftshift(ifft(Sy));
f=(-N:N-1)*fs/(2*N);
m=(-N:N-1)/N*(N/2000);
subplot(311);plot((-N:N-1)/N,fftshift(abs(H2(1:2*N))));
title('高通滤波器');
subplot(312),plot(m,Ry);
xlabel('m')
ylabel('Ry(m)')
title('xn经高通通滤波器的自相关函数');
subplot(313),plot(f,fftshift(10*log10(Sy(1:2*N))));
axis([-200 200 -20 50]);
xlabel('f/Hz')
ylabel('Sy/dB')
title('xn经高通滤波器的功率谱密度');
5、多带通
N=2000;fs=400;
Nn=random('normal',0,1,1,N);
t=(0:N-1)/fs;
fi=random('unif',0,1,1,2)*2*pi;
xn=cos(2*pi*50*t+fi(1))+3*cos(2*pi*150*t+fi(2))+Nn;
h=fir1(100,[0.2,0.4,0.6,0.8]);
H=fft(h,2*N);
H2=abs(H).^2;
Rx=xcorr(xn,'biased');
Sx=abs(fftshift(fft(xn,2*N)).^2)/(2*N);
Sy=Sx.*H2;
Ry=fftshift(ifft(Sy));
f=(-N:N-1)*fs/(2*N);
m=(-N:N-1)/N*(N/2000);
subplot(311);plot((-N:N-1)/N,fftshift(abs(H2(1:2*N))));
title('多带通滤波器');
subplot(312),plot(m,Ry);
xlabel('m')
ylabel('Ry(m)')
title('xn经多带通通滤波器的自相关函数');
subplot(313),plot(f,fftshift(10*log10(Sy(1:2*N))));
axis([-200 200 -20 50]);
xlabel('f/Hz')
ylabel('Sy/dB')
title('xn经多带通滤波器的功率谱密度');
6、带阻
N=2000;fs=400;
Nn=random('normal',0,1,1,N);
t=(0:N-1)/fs;
fi=random('unif',0,1,1,2)*2*pi;
xn=cos(2*pi*50*t+fi(1))+3*cos(2*pi*150*t+fi(2))+Nn;
h=fir1(100,[0.2,0.5],'stop');
H=fft(h,2*N);
H2=abs(H).^2;
Rx=xcorr(xn,'biased');
Sx=abs(fftshift(fft(xn,2*N)).^2)/(2*N);
Sy=Sx.*H2;
Ry=fftshift(ifft(Sy));
f=(-N:N-1)*fs/(2*N);
m=(-N:N-1)/N*(N/2000);
subplot(311);plot((-N:N-1)/N,fftshift(abs(H2(1:2*N))));
title('带阻滤波器');
subplot(312),plot(m,Ry);
xlabel('m')
ylabel('Ry(m)')
title('xn经带阻滤波器的自相关函数');
subplot(313),plot(f,fftshift(10*log10(Sy(1:2*N))));
axis([-200 200 -20 50]);
xlabel('f/Hz')
ylabel('Sy/dB')
title('xn经带阻滤波器的功率谱密度');
③仿真图形
1、X（n）信号的自相关函数及功率谱密度
2、低通
3、带通
5、多带通
④分析：由图形可得，给定的信号频率在频谱图上明显存在峰值，当信号经过各系统时只要频率在通带内即可输出峰值。