随机信号分析实验报告.docx
《随机信号分析实验报告.docx》由会员分享，可在线阅读，更多相关《随机信号分析实验报告.docx（25页珍藏版）》请在冰点文库上搜索。
随机信号分析实验报告
HarbinInstituteofTechnology
实验报告
课程名称：
随机信号分析
院系：
电子与信息工程学院
班级：
姓名：
学号：
指导教师：
实验时间：
实验一、各种分布随机数的产生
（一）实验原理
1.均匀分布随机数的产生原理
产生伪随机数的一种实用方法是同余法，它利用同余运算递推产生伪随机数序列。
最简单的方法是加同余法
为了保证产生的伪随机数能在[0,1]内均匀分布，需要M为正整数，此外常数c和初值y0亦为正整数。
加同余法虽然简单，但产生的伪随机数效果不好。
另一种同余法为乘同余法，它需要两次乘法才能产生一个[0,1]上均匀分布的随机数
式中，a为正整数。
用加法和乘法完成递推运算的称为混合同余法，即
用混合同余法产生的伪随机数具有较好的特性，一些程序库中都有成熟的程序供选择。
常用的计算语言如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]上均匀分布的随机变量，那么
即是分布函数为FX（x）的随机变量。
式中
为
的反函数。
这样，欲求某个分布的随机变量，先产生在[0,1]区间上的均匀分布随机数，再经上式变换，便可求得所需分布的随机数。
3.高斯分布随机数的仿真
广泛应用的有两种产生高斯随机数的方法，一种是变换法，一种是近似法。
如果X1,X2是两个互相独立的均匀分布随机数，那么下式给出的Y1,Y2
便是数学期望为m，方差为
的高斯分布随机数，且互相独立，这就是变换法。
另外一种产生高斯随机数的方法是近似法。
在学习中心极限定理时，曾提到n个在[0,1]区间上均匀分布的互相独立随机变量Xi（i=1,2…,n），当n足够大时，其和的分布接近高斯分布。
当然，只要n不是无穷大，这个高斯分布是近似的。
由于近似法避免了开方和三角函数运算，计算量大大降低。
当精度要求不太高时，近似法还是具有很大应用价值的。
4.各种分布随机数的仿真
有了高斯随机变量的仿真方法，就可以构成与高斯变量有关的其他分布随机变量，如瑞利分布、指数分布和
分布随机变量。
（二）实验目的
在很多系统仿真的过程中，需要产生不同分布的随机变量。
利用计算机可以很方便地产生不同分布的随机变量，各种分布的随机变量的基础是均匀分布的随机变量。
有了均匀分布的随机变量，就可以用函数变换等方法得到其他分布的随机变量。
（三）实验结果
附：
源程序
subplot（2,2,1）;
x=random（'unif',2,5,1,1024）;
plot（x）;
title（'均匀分布随机数'）
subplot（2,2,2）;
G1=random（'Normal',0,1,1,20000）;
plot（G1）;
title（'高斯分布随机数'）
subplot（2,2,3）;
G2=random（'Normal',0,1,1,20000）;
R=sqrt（G1.*G1+G2.*G2）;
plot（R）;
title（'瑞利分布随机数'）
subplot（2,2,4）;
G3=random（'Normal',0,1,1,20000）;
G4=random（'Normal',0,1,1,20000）;
X=G1.*G1+G2.*G2+G3.*G3+G4.*G4;
plot（X）;
title（'x^2分布随机数'）
实验二、随机变量检验
（一）实验原理
1、均值的计算
在实际计算时，如果平稳随机序列满足各态历经性，则统计均值可用时间均值代替。
这样，在计算统计均值时，并不需要大量样本函数的集合，只需对一个样本函数求时间平均即可。
甚至有时也不需要计算
时的极限，况且也不可能。
通常的做法是取一个有限的、计算系统能够承受的N求时间均值和时间方差。
根据强调计算速度或精度的不同，可选择不同的算法。
设随机数序列{
}，一种计算均值的方法是直接计算下
式中，xn为随机数序列中的第n个随机数。
另一种方法是利用递推算法，第n次迭代的均值也亦即前n个随机数的均值为
迭代结束后，便得到随机数序列的均值
递推算法的优点是可以实时计算均值，这种方法常用在实时获取数据的场合。
当数据量较大时，为防止计算误差的积累，也可采用
式中，m1是取一小部分随机数计算的均值。
2、方差的计算
计算方差也分为直接法和递推法。
仿照均值的做法
方差的递推算法需要同时递推均值和方差
迭代结束后，得到随机数序列的方差为
其它矩函数也可用类似的方法得到。
3、统计随机数的概率密度直方图
假定被统计的序列
的最大值和最小值分别为a和b。
将
区间等分M（M应与被统计的序列
的个数N相适应，否则统计效果不好。
）份后的区间为
，
，…,
，…,
。
用
，表示序列
的值落在
区间里的个数，统计序列
的值在各个区间的个数
，
，则
就粗略地反映了随机序列的概率密度的情况。
用图形方式显示出来就是随机数的概率密度直方图。
（二）实验目的
随机数产生之后，必须对它的统计特性做严格的检验。
一般来讲，统计特性的检验包括参数检验、均匀性检验和独立性检验等。
事实上，我们如果在二阶矩范围内讨论随机信号，那么参数检验只对产生的随机数一、二阶矩进行检验。
我们可以把产生的随机数序列作为一个随机变量，也可以看成随机过程中的一个样本函数。
不论是随机变量还是随机过程的样本函数，都会遇到求其数字特征的情况，有时需要计算随机变量的概率密度直方图等。
（三）实验结果
附：
源程序
subplot（2,2,1）;
x=random（'unif',2,5,1,1024）;
hist（x,2:
0.2:
5）;
title（'均匀分布随机数直方图'）;
s1=0
forn1=1:
1024
s1=x（n1）+s1;
end
Mean1=s1/1024;
t1=0
forn1=1:
1024
t1=（x（n1）-Mean1）^2+t1;
end
Variance1=t1/1024;
subplot（2,2,2）;
G1=random（'Normal',0,1,1,20000）;
hist（G1,-4:
0.2:
4）;
title（'高斯分布随机数直方图'）;
s2=0
forn2=1:
20000
s2=G1（n2）+s2;
end
Mean2=s2/20000;
t2=0
forn2=1:
20000
t2=（G1（n2）-Mean2）^2+t2;
end
Variance2=t2/20000;
subplot（2,2,3）;
G2=random（'Normal',0,1,1,20000）;
R=sqrt（G1.*G1+G2.*G2）;
hist（R,0:
0.2:
5）;
title（'瑞利分布随机数直方图'）;
s3=0
forn3=1:
20000
s3=R（n3）+s3;
end
Mean3=s3/20000;
t3=0
forn3=1:
20000
t3=（R（n3）-Mean3）^2+t3;
end
Variance3=t3/20000;
subplot（2,2,4）;
G3=random（'Normal',0,1,1,20000）;
G4=random（'Normal',0,1,1,20000）;
X=G1.*G1+G2.*G2+G3.*G3+G4.*G4;
hist（X,0:
0.5:
30）;
title（'x^2分布随机数直方图'）
s4=0
forn4=1:
20000
s4=X（n4）+s4;
end
Mean4=s4/20000;
t4=0
forn4=1:
20000
t4=（X（n4）-Mean4）^2+t4;
end
实验三、中心极限定理的验证
（一）实验原理
如果n个独立随机变量的分布是相同的，并且具有有限的数学期望和方差，当n无穷大时，它们之和的分布趋近于高斯分布。
这就是中心极限定理中的一个定理。
我们以均匀分布为例，来解释这个定理。
若n个随机变量Xi（i=1,2,…,n）都为[0,1]区间上的均匀分布的随机变量，且互相独立，当n足够大时，其和
的分布接近高斯分布。
（二）实验目的
利用计算机产生均匀分布的随机数。
对相互独立的均匀分布的随机变量做和，可以很直观看到均匀分布的随机变量的和，随着做和次数的增加分布情况的变化，通过实验对中心极限定理的进行验证。
（三）实验结果
分析：
随n取值的增大，均匀分布随机序列求和的图形越发接近于高斯分布。
附：
源程序
X0=random（'unif',0,1,1,1024）;
X1=random（'unif',0,1,1,1024）;
X2=random（'unif',0,1,1,1024）;
X3=random（'unif',0,1,1,1024）;
X4=random（'unif',0,1,1,1024）;
X5=random（'unif',0,1,1,1024）;
X6=random（'unif',0,1,1,1024）;
X7=random（'unif',0,1,1,1024）;
X8=random（'unif',0,1,1,1024）;
X9=random（'unif',0,1,1,1024）;
G=random（'normal',0,1,1,1024）;
Y1=X0+X1+X2+X3+X4;
Y2=X0+X1+X2+X3+X4+X5+X6+X7+X8+X9;
subplot（2,2,1）;
hist（X0,0:
0.2:
2）;
title（'均匀分布随机数直方图'）
subplot（2,2,2）;
hist（Y1,0:
0.2:
6）;
title（'五个均匀分布之和随机数直方图'）
subplot（2,2,3）;
hist（Y2,0:
0.2:
8）;
title（'十个均匀分布之和随机数直方图'）
subplot（2,2,4）;
hist（G,-4:
0.2:
4）;
title（'高斯分布随机数直方图'）
实验四、中心极限定理的验证
（一）实验原理
在实际应用中，我们可以把产生的随机数序列看成随机过程中的一个样本函数。
如果平稳随机序列满足各态历经性，则统计自相关序列可用时间自相关序列代替。
当数据的样本数有限时，也只能用有限个数据来估计时间自相关序列，统计自相关序列的估值。
若各态历经序列X（n）的一个样本有N个数据
，由于实序列自相关序列是对称的，自相关函数的估值为
（二）实验目的
在随机信号理论中，自相关函数是非常重要的概念。
在实际系统仿真中也会经常计算自相关函数。
通过本试验学生可以亲自动手计算自相关函数，加深对概念的理解，并增强实际动手能力。
（三）实验结果
分析：
分别生成均值为0和1，方差为1的高斯随机数，由图形可以明显看出两者自相关函数的差异。
附：
源程序
N=256;
xn=random（'norm',0,1,1,N）;
Rx=xcorr（xn,'biased'）;
m=-N+1:
N-1;
subplot（2,1,1）;
plot（m,Rx）;
title（'均值为0，方差为1的高斯分布的自相关函数'）;
axis（[-NN-1-0.51.5]）;
N=256;
xn=random（'norm',1,1,1,N）;
Xk=fft（xn,2*N）;
Rx=ifft（（abs（Xk）.^2）/N）;
m=-N:
N-1;
subplot（2,1,2）;
plot（m,fftshift（Rx））;
title（'均值为1，方差为1的高斯分布的自相关函数'）;
axis（[-NN-1-0.51.5]）;
实验五、功率谱密度
（一）实验原理
一般把平稳随机序列的功率谱定义为自相关序列的傅里叶变换。
如果自相关序列是周期序列，可仿照随机过程的情况，引人适当的函数。
平稳序列X（n）的功率谱与自相关序列的关系为
与实平稳过程一样，实平稳序列的功率谱也是非负偶函数，即
可以证明，功率谱还可表示为
当X（n）为各态历经序列时，可去掉上式中的统计均值计算，将随机序列X（n）用它的一个样本序列x（n）代替。
在实际应用中，由于一个样本序列的可用数据个数N有限，功率谱密度也只能是估计
式中，X（）是x（n）的傅里叶变换。
这是比较简单的一种估计方法，这种功率谱密度的估计方法称为周期图方法。
如果直接利用数据样本做离散傅里叶变换，可得到X（）的离散值。
由于这种方法可借助FFT算法实现，所以得到了广泛的应用。
（二）实验目的
在随机信号理论中，功率谱密度和自相关函数一样都是非常重要的概念。
在实际系统仿真中也会经常计算。
通过本试验学生可以亲自动手，加深对概念的理解，并增强实际动手能力。
（三）实验结果
附：
源程序
N=256;
x1=random（'normal',0,1,1,N）;
Sx1=abs（fft（x1）.^2）/N;
subplot（2,1,1）;
plot（10*log10（Sx1））;
title（'均值为0，方差为1的高斯分布的功率谱密度'）;
xlabel（'f/Hz'）
ylabel（'Sx1/dB'）
x2=random（'normal',1,1,1,N）;
Sx2=abs（fft（x2）.^2）/N;
subplot（2,1,2）;
plot（10*log10（Sx2））;
title（'均值为1，方差为1的高斯分布的功率谱密度'）;
xlabel（'f/Hz'）
ylabel（'Sx2/dB'）
实验六、随机信号经过线性系统前后信号仿真
（一）实验原理
需要先仿真一个指定系统，再根据需要仿真输入的随机信号，然后使这个随机信号通过指定的系统。
通过对实际系统建模,计算机可以对很多系统进行仿真。
在信号处理中，一般将线性系统分解为一个全通放大器（或衰减器）和一个特定频率响应的滤波器。
由于全通放大器可以用一个常数代替，因此线性系统的仿真往往只需设计一个数字滤波器。
滤波器设计可采用MATLAB提供的函数，也可利用相应的方法自行设计。
MATLAB提供了多个设计滤波器的函数，可以很方便地设计低通、带通、高通、多带通、带阻滤波器。
（二）实验目的
系统仿真是信号仿真处理的一个重要部分，通过该实验要求学生掌握系统仿真的基本概念，并学会系统的仿真方法。
（三）实验结果
1、低通滤波器
2、带通滤波器
3、高通滤波器
4、多带通滤波器
5、带阻滤波器
附：
源程序
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=sin（2*pi*50*t+fi
（1））+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.低通滤波器
h=fir1（100,0.4）;
H=fft（h,2*N）;
HW=abs（H）.^2;
Rx=xcorr（xn,'biased'）;
Sx=abs（fftshift（fft（xn,2*N））.^2）/（2*N）;
Sy=Sx.*HW;
Ry=fftshift（ifft（Sy））;
f=（-N:
N-1）*fs/（2*N）;
m=（-N:
N-1）;
subplot（311）;plot（（-N:
N-1）/N,fftshift（abs（HW（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（[-200200-2020]）;
xlabel（'f/Hz'）
ylabel（'Sy/dB'）
title（'xn经低通滤波器的功率谱密度'）;
3.带通滤波器
h=fir1（100,[0.10.5]）;
H=fft（h,2*N）;
HW=abs（H）.^2;
Rx=xcorr（xn,'biased'）;
Sx=abs（fftshift（fft（xn,2*N））.^2）/（2*N）;
Sy=Sx.*HW;
Ry=fftshift（ifft（Sy））;
f=（-N:
N-1）*fs/（2*N）;
m=（-N:
N-1）;
subplot（311）;plot（（-N:
N-1）/N,fftshift（abs（HW（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（[-200200-2020]）;
xlabel（'f/Hz'）
ylabel（'Sy/dB'）
title（'xn经带通滤波器的功率谱密度'）;
4.高通滤波器
h=fir1（100,0.6,'high'）;
H=fft（h,2*N）;
HW=abs（H）.^2;
Rx=xcorr（xn,'biased'）;
Sx=abs（fftshift（fft（xn,2*N））.^2）/（2*N）;
Sy=Sx.*HW;
Ry=fftshift（ifft（Sy））;
f=（-N:
N-1）*fs/（2*N）;
m=（-N:
N-1）;
subplot（311）;plot（（-N:
N-1）/N,fftshift（abs（HW（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（[-200200-2020]）;
xlabel（'f/Hz'）
ylabel（'Sy/dB'）
title（'xn经高通滤波器的功率谱密度'）;
5.多带通滤波器
h=fir1（100,[0.1,0.3,0.5,0.7]）;
H=fft（h,2*N）;
HW=abs（H）.^2;
Rx=xcorr（xn,'biased'）;
Sx=abs（fftshift（fft（xn,2*N））.^2）/（2*N）;
Sy=Sx.*HW;
Ry=fftshift（ifft（Sy））;
f=（-N:
N-1）*fs/（2*N）;
m=（-N:
N-1）;
subplot（311）;plot（（-N:
N-1）/N,fftshift（abs（HW（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（[-200200-2020]）;
xlabel（'f/Hz'）
ylabel（'Sy/dB'）
title（'xn经多带通滤波器的功率谱密度'）;
6.带阻滤波器
h=fir1（100,[0.1,0.4],'stop'）;
H=fft（h,2*N）;
HW=abs（H）.^2;
Rx=xcorr（xn,'biased'）;
Sx=abs（fftshift（fft（xn,2*N））.^2）/（2*N）;
Sy=Sx.*HW;
Ry=fftshift（ifft（Sy））;
f=（-N:
N-1）*fs/（2*N）;
m=（-N:
N-1）;
subplot（311）;plot（（-N:
N-1）/N,fftshift（abs（HW（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（[-200200-2020]）;
xlabel（'f/Hz'）
ylabel（'Sy/dB'）
title（'xn经带阻滤波器的功率谱密度'）;
升级会员 