随机信号分析.doc

Harbin Institute of Technology实验报告课程名称： 随机信号分析 院 系： 电信学院 班 级： 姓 名： 学 号： 指导教师： 郑 薇 实验时间： 2014年 11月 哈尔滨工业大学 实验一 各种分布随机数的产生实验目的在很多系统仿真的过程中，需要产生不同分布的随机变量。利用计算机可以很方便地产生不同分布的随机变量，各种分布的随机变量的基础是均匀分布的随机变量。有了均匀分布的随机变量，就可以用函数变换等方法得到其他分布的随机变量。实验内容产生均匀分布的随机数、高斯分布的随机数和其它分布的随机数。实验原理 均匀分布随机数的产生原理产生伪随机数的一种实用方法是同余法，它利用同余运算递推产生伪随机数序列。最简单的方法是加同余法为了保证产生的伪随机数能在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分别是矩阵的行和列。 随机变量的仿真根据随机变量函数变换的原理，如果能将两个分布之间的函数关系用显式表达，那么就可以利用一种分布的随机变量通过变换得到另一种分布的随机变量。若X是分布函数为F(x)的随机变量，且分布函数F(x)为严格单调升函数，令Y=F(X)，则Y必为在0,1上均匀分布的随机变量。反之，若Y是在0,1上均匀分布的随机变量，那么即是分布函数为FX(x)的随机变量。式中为的反函数。这样，欲求某个分布的随机变量，先产生在0,1区间上的均匀分布随机数，再经上式变换，便可求得所需分布的随机数。 高斯分布随机数的仿真广泛应用的有两种产生高斯随机数的方法，一种是变换法，一种是近似法。如果X1,X2是两个互相独立的均匀分布随机数，那么下式给出的Y1,Y2便是数学期望为m，方差为s2的高斯分布随机数，且互相独立，这就是变换法。另外一种产生高斯随机数的方法是近似法。在学习中心极限定理时，曾提到n个在0,1区间上均匀分布的互相独立随机变量Xi (i=1,2,n)，当n足够大时，其和的分布接近高斯分布。当然，只要n不是无穷大，这个高斯分布是近似的。由于近似法避免了开方和三角函数运算，计算量大大降低。当精度要求不太高时，近似法还是具有很大应用价值的。各种分布随机数的仿真有了高斯随机变量的仿真方法，就可以构成与高斯变量有关的其他分布随机变量，如瑞利分布、指数分布和c2分布随机变量。实验过程和结果分析1. 产生均匀分布的随机数 for n=1:1024 y=rand(); x(n)=y*(6-3)+3; end plot(x); 2.产生高斯分布的随机数 x=random(Normal,0,2,1,1024)； 3.产生瑞利分布和分布 N=30000; g=-6:0.1:6; 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); X2=G1.*G1+G2.*G2+G3.*G3+G4.*G4;实验结论使用Matlab产生均匀分布的随机数、高斯分布的随机数、瑞利分布和分布的随机数。 实验二 随机变量检验实验目的随机数产生之后，必须对它的统计特性做严格的检验。一般来讲，统计特性的检验包括参数检验、均匀性检验和独立性检验等。事实上，我们如果在二阶矩范围内讨论随机信号，那么参数检验只对产生的随机数一、二阶矩进行检验。我们可以把产生的随机数序列作为一个随机变量，也可以看成随机过程中的一个样本函数。不论是随机变量还是随机过程的样本函数，都会遇到求其数字特征的情况，有时需要计算随机变量的概率密度直方图等。实验内容1. 对实验一产生的各种分布的随机数进行均值和方差的检验。2. 对实验一产生的各种分布的随机数概率分布进行统计，并在计算机屏幕上显示实际统计的概率密度直方图。实验原理1. 均值的计算在实际计算时，如果平稳随机序列满足各态历经性，则统计均值可用时间均值代替。这样，在计算统计均值时，并不需要大量样本函数的集合，只需对一个样本函数求时间平均即可。甚至有时也不需要计算时的极限，况且也不可能。通常的做法是取一个有限的、计算系统能够承受的N求时间均值和时间方差。根据强调计算速度或精度的不同，可选择不同的算法。设随机数序列，一种计算均值的方法是直接计算下式式中，xn为随机数序列中的第n个随机数。另一种方法是利用递推算法，第n次迭代的均值也亦即前n个随机数的均值为迭代结束后，便得到随机数序列的均值递推算法的优点是可以实时计算均值，这种方法常用在实时获取数据的场合。当数据量较大时，为防止计算误差的积累，也可采用式中，m1是取一小部分随机数计算的均值。方差的计算计算方差也分为直接法和递推法。仿照均值的做法方差的递推算法需要同时递推均值和方差迭代结束后，得到随机数序列的方差为其它矩函数也可用类似的方法得到。统计随机数的概率密度直方图假定被统计的序列的最大值和最小值分别为a和b。将区间等分M（M应与被统计的序列的个数N相适应，否则统计效果不好。）份后的区间为， , ， , 。用，表示序列的值落在区间里的个数，统计序列的值在各个区间的个数，则就粗略地反映了随机序列的概率密度的情况。用图形方式显示出来就是随机数的概率密度直方图。实验过程和结果分析1. 均值和方差的检验 (1)均匀分布随机数 x=random(unif,3,6,1,1024) m=mean(x)m = 4.5064 d=var(x)d = 0.7523(2) 产生高斯分布,瑞利分布和分布的均值与方差 N=30000;g=-6:0.1:6;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);X2=G1.*G1+G2.*G2+G3.*G3+G4.*G4; m1=mean(G1)m1 = 0.0055 d1=var(G1)d1 = 1.0101 m2=mean(R)m2 = 1.2576 d2=var(R)d2 = 0.4320 m3=mean(X2)m3 = 4.0194 d3=var(X2)d3 = 8.16842. 概率密度直方图（1） 均匀分布随机数 x=random(unif,3,6,1,1024); subplot;hist(x,2:0.01:7);（2） 高斯分布,瑞利分布和分布 N=30000; g=-6:0.1:6; 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); X2=G1.*G1+G2.*G2+G3.*G3+G4.*G4; subplot(311);hist(G1,g); subplot(312);hist(R,0:0.05:6); subplot(313);hist(X2,0:0.02:30);实验结论1. 对实验一产生各种分布的均值和方差验证，结果如下（1）均匀分布m = 4.5064 d = 0.7523（2） 高斯分布,瑞利分布和分布的均值与方差m1 = 0.0055 d1 =1.0101m2 =1.2576 d2 = 0.4320m3 =4.0194 d3 =8.16842.概率密度直方图如图所示实验三 中心极限定理的验证实验目的 利用计算机产生均匀分布的随机数。对相互独立的均匀分布的随机变量做和，可以很直观看到均匀分布的随机变量的和，随着做和次数的增加分布情况的变化，通过实验对中心极限定理的进行验证。实验内容产生多组0,1区间上的均匀分布的随机数序列，各序列的对应元素做和，够成的和序列再进行随机数的概率密度直方图的统计，并作图显示。实验原理如果n个独立随机变量的分布是相同的，并且具有有限的数学期望和方差，当n无穷大时，它们之和的分布趋近于高斯分布。这就是中心极限定理中的一个定理。我们以均匀分布为例，来解释这个定理。若n个随机变量Xi (i=1,2,n)都为0,1区间上的均匀分布的随机变量，且互相独立，当n足够大时，其和的分布接近高斯分布。实验过程和结果分析 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); x=x1+x2+x3+x4+x5+x6+x7+x8; subplot;hist(x,0:0.05:8);实验结论若n个随机变量Xi (i=1,2,n)都为0,1区间上的均匀分布的随机变量，且互相独立，当n足够大时，其和的分布接近高斯分布。实验四 自相关函数的计算实验目的在随机信号理论中，自相关函数是非常重要的概念。在实际系统仿真中也会经常计算自相关函数。通过本试验学生可以亲自动手计算自相关函数，加深对概念的理解，并增强实际动手能力。实验内容用一个数学期望为零和非零，方差为某值的高斯分布随机数，作为样本序列求自相关函数的估值，并用图形显示。实验原理在实际应用中，我们可以把产生的随机数序列看成随机过程中的一个样本函数。如果平稳随机序列满足各态历经性，则统计自相关序列可用时间自相关序列代替。当数据的样本数有限时，也只能用有限个数据来估计时间自相关序列，统计自相关序列的估值。若各态历经序列X(n)的一个样本有N个数据，由于实序列自相关序列是对称的，自相关函数的估值为实验过程和结果分析 N=1024; xn1=random(norm,0,1,1,N); xn2=random(norm,1,1,1,N); Xk1=fft(xn1,2*N); Xk2=fft(xn2,2*N); Rx1=ifft(abs(Xk1).2)/N); Rx2=ifft(abs(Xk2).2)/N); m=-N:N-1; axis(-N N-1 -0.5 1.5); subplot(211); plot(m,fftshift(Rx1); subplot(212); plot(m,fftshift(Rx2); 实验结论用一个数学期望为零和非零，方差为某值的高斯分布随机数，作为样本序列求自相关函数的估值如上图所示。 实验五 功率谱密度实验目的在随机信号理论中，功率谱密度和自相关函数一样都是非常重要的概念。在实际系统仿真中也会经常计算。通过本试验学生可以亲自动手，加深对概念的理解，并增强实际动手能力。实验内容用实验四计算出的自相关函数的估值，作为样本序列求功率谱密度的估值，并用图形显示。实验原理一般把平稳随机序列的功率谱定义为自相关序列的傅里叶变换。如果自相关序列是周期序列，可仿照随机过程的情况，引人适当的d函数。平稳序列X(n)的功率谱与自相关序列的关系为与实平稳过程一样，实平稳序列的功率谱也是非负偶函数，即可以证明，功率谱还可表示为当X(n)为各态历经序列时，可去掉上式中的统计均值计算，将随机序列X(n)用它的一个样本序列x(n)代替。在实际应用中，由于一个样本序列的可用数据个数N有限，功率谱密度也只能是估计值式中，X(w)是x(n)的傅里叶变换。这是比较简单的一种估计方法，这种功率谱密度的估计方法称为周期图方法。如果直接利用数据样本做离散傅里叶变换，可得到X(w)的离散值。由于这种方法可借助FFT算法实现，所以得到了广泛的应用。实验过程和结果分析1.高斯分布样本序列求功率谱密度 N=1024; fs=1000; xn1=random(norm,0,1,1,N); xn2=random(norm,1,1,1,N); Sx1=abs(fft(xn1).2/N; Sx2=abs(fft(xn2).2/N; f=(-N/4+1:N/4)*fs/N; subplot(211);plot(f,10*log10( Sx1(1:N/2); subplot(212);plot(f,10*log10( Sx2(1:N/2);2. f1=-50HZ;f2=60HZ. N=1024; fs=1000; t=(-N/2+1:N/2)/fs; fai=random(unif,0,1,1,2)*2*pi; xn=cos(2*pi*-50*t+fai(1)+3*cos(1*pi*60*t+fai(2)+randn(1,N); Sx=abs(fft(xn).2/N; f=(-N/4+1:N/4)*fs/N; plot(f,10*log10( Sx(1:N/2);实验结论实验六 随机信号经过线性系统前后信号仿真实验目的系统仿真是信号仿真处理的一个重要部分，通过该实验要求学生掌握系统仿真的基本概念，并学会系统的仿真方法。实验内容仿真信号和加性噪声经过各种系统前后的自相关函数和功率谱密度并图示。实验原理需要先仿真一个指定系统，再根据需要仿真输入的随机信号，然后使这个随机信号通过指定的系统。通过对实际系统建模, 计算机可以对很多系统进行仿真。在信号处理中，一般将线性系统分解为一个全通放大器（或衰减器）和一个特定频率响应的滤波器。由于全通放大器可以用一个常数代替，因此线性系统的仿真往往只需设计一个数字滤波器。滤波器设计可采用MATLAB提供的函数，也可利用相应的方法自行设计。MATLAB提供了多个设计滤波器的函数，可以很方便地设计低通、带通、高通、多带通、带阻滤波器。实验过程和结果分析根据自己的仿真思路简要阐述，重点给出相应仿真图（结果），并进行分析！注：使用例：2.6-4中的随机信号作为输入信号：频率分量可以根据设计的系统任意选取。将输入信号分别通过设计的各种系统低通、带通、高通、多带通、带阻滤波器（要求做3种情况，多做可酌情加分！）1. 低通：N=1024;fs=1000;f1=30;f2=200;t=(0:N-1)/fs;fai=random(unif,0,1,1,2)*4*pi; %产生两个（0，4pi）的随机数xn=cos(2*pi*f1*t+fai(1)+3*cos(2*pi*f2*t+fai(2)+randn(1,N);Sx=abs(fft(xn).2)/N;h=fir1(50,0.2);H1=fft(h,N);Sxx=Sx.*(abs(H1).2);Rxx=fftshift(ifft(Sxx);w=(1:N)/N;t=(0:N-1)/N*(N/500);subplot(4,1,1);plot(w,10*log10(Sx(1:N);subplot(4,1,2);plot(w,abs(H1(1:N);subplot(4,1,3);plot(w,abs(Sxx(1:N);subplot(4,1,4);plot(t,Rxx);2.带通：N=1024;fs=1000;f1=30;f2=200;t=(0:N-1)/fs;fai=random(unif,0,1,1,2)*4*pi; %产生两个（0，4pi）的随机数xn=cos(2*pi*f1*t+fai(1)+3*cos(2*pi*f2*t+fai(2)+randn(1,N);Sx=abs(fft(xn,2*N).2)/N;h=fir1(101,0.2,0.6);H1=fft(h,2*N);Sxx=Sx.*(abs(H1).2)/(2*N);Rxx=fftshift(ifft(Sxx);w=(1:N)/N;t=(-N:N-1)/N*(N/500);subplot(4,1,1);plot(w,10*log10(Sx(1:N);subplot(4,1,2);plot(w,abs(H1(1:N);subplot(4,1,3);plot(w,abs(Sxx(1:N);subplot(4,1,4);plot(t,Rxx);3.多带通：N=1000;band=1:50:500/(1000/2);h3=fir1(200,band,DC-0);H3=fft(h3,N);f=(1:N/2)/(N/2);plot(f,20*log10(abs(H3(1:N/2);N=1024;fs=1000;f1=30;f2=200;t=(0:N-1)/fs;fai=random(unif,0,1,1,2)*4*pi; %产生两个（0，4pi）的随机数xn=cos(2*pi*f1*t+fai(1)+3*cos(2*pi*f2*t+fai(2)+randn(1,N);Sx=abs(fft(xn,2*N).2)/N;band=1:50:500/(1000/2);h3=fir1(200,band,DC-0);H3=fft(h3,2*N);Sxx=Sx.*(abs(H3).2)/(2*N);Rxx=fftshift(ifft(Sxx);w=(1:N)/N;t=(-N:N-1)/N*(N/500);subplot(4,1,1);plot(w,10*log10(Sx(1:N);subplot(4,1,2);plot(w,abs(H3(1:N);subplot(4,1,3);plot(w,abs(Sxx(1:N);subplot(4,1,4);plot(t,Rxx);实验结论实验结论1.当分别等于25Hz和200Hz时，画出其功率谱密度如下图2.（1）低通100Hz 归一化数字频率为0.2（2）经过低通滤波器后的结果3.（1）带通滤波器（2）加完带通滤波器后的结果4.多带通