随机信号的仿真与分析.doc

实验十八 随机信号的仿真与分析一实验目的利用计算机仿真随机信号，考察其数字特征，以此加深对满足各种分布的随机信号的理解。熟悉常用的信号处理仿真软件平台：matlab或c/c+语言.二实验原理 随机信号的产生和定义随机信号是随机变量在时间上推进产生的过程量，它同时具有过程性和不确定性。定义如下：给定参量集T与概率空间（, F, P），若对于每个，都有一个定义在（, F, P）上的实随机变量X(t)与之对应，就称依赖于参量t的随机变量族为一（实）随机过程或随机信号。 高斯分布随机信号统计分布是正态分布（高斯分布）的随机信号为高斯分布随机信号。高斯分布的随机变量概率密度函数满足下式： 均匀分布随机信号统计分布是均匀分布的随机信号为均匀分布随机信号。均匀分布的随机变量概率密度函数满足下式： 正弦随机信号给定具有某种概率分布的振幅随机变量A、角频率随机变量与相位随机变量，（具体概率分布与特性视应用而定），以（时间）参量t建立随机变量：。于是，相应于某个参量域T的随机变量族为正弦随机信号（或称为正弦随机过程）。5 贝努里随机信号贝努里随机变量X(s)基于一个掷币实验（s表示基本结果事件）：1表示s为正面，0表示s不为正面；s不为正面的概率为PX(s)=1=p，s为正面的概率为PX(s)=0=q，其中p+q=1。若无休止地在t=n (n=0, 1, 2, )时刻上，独立进行（相同的）掷币实验构成无限长的随机变量序列：,其中与n和s都有关，应记为X(n,s)，于是，而且有概率：其中, p+q=1。上述的随机变量序列：通常被称为随机序列（或随机过程），也被称为（离散）随机信号，即贝努里随机信号。正式的定义如下：给定某个序列随机实验，观测某事件B发生与否，建立事件B的指示函数，而且，序列随机实验间彼此统计独立并有相同的概率，于是，是一个（0，1）贝努里随机变量，相应的随机变量序列为（0，1）贝努里随机序列（或称随机信号，有时也称为随机过程）。三实验任务与要求 用matlab或c/c+语言编程并仿真。 生成满足几种概率分布的仿真随机信号，自己编写程序计算几种概率分布的仿真随机信号的特征。具体要求： 随机数的产生与测量：产生的10000个泊松分布随机数，计算它们的均值、方差与概率密度、频谱、功率谱密度，自相关函数，并绘出函数曲线。确定泊松过程是一个马尔可夫过程。 产生N(0,3)与N(2,3)的随机数，计算它们的概率密度、频谱、功率谱密度，自相关函数，并绘出函数曲线。确定它们是否属于白噪声。 统计分析：二维正态分布（X,Y），N(0，1；0，4；0.5)的联合概率密度函数为 其中1=0，1=1，2=0，2=4，=0.5。求 ；） ；）。 计算二维正态分布（X,Y），N(0，1；0，4；0.5) 的概率密度、频谱、功率谱密度，自相关函数，并绘出函数曲线。 按要求写实验报告。四实验仿真与分析1随机信号的仿真（1）高斯分布随机信号的仿真MATLAB程序：N=1000;M=0;SIGMA2=3;SIGMA=sqrt(SIGMA2);X=M+SIGMA.*randn(N,1);t=0:N-1plot(t,X);grid on title(高斯分布的随机信号);xlabel(t); 说明：M为均值，SIGMA2为其方差；用randn（N，1）函数产生一个N个样本的正态分布的列矢量；用X=M+SIGMA.*randn(N,1)仿真产生N（M,SIGMA2）的高斯分布的随机数；取t=0:N-1,用plot(t,X)对随机数进行仿真绘图。编程运行后得到的高斯分布随机信号的仿真图如下：(2)均匀分布随机信号的仿真程序：a=5;b=8;N=1000; x=b-(b-a).*rand(N,1);t=0:N-1;plot(t,x);title(均匀分布随机信号 );说明：用rand（N，1）函数产生一个N个样本的（0,1）均匀分布的列矢量；用x=b-(b-a).*rand(N,1)仿真产生（a,b）的均匀分布随机数；取t=0:N-1，用plot(t，x)对随机数进行仿真绘图；编程运行后得到的均匀分布随机信号的仿真图如下：（3）正弦随机信号的仿真程序：t=0:pi/50:2*pi;for i=1:100 A=rand(); w=rand()*2; B=rand()*2*pi; y(i)=A*sin(w*t(i)+B);endfigure(1);plot(y,-r);grid on;title(y=Asin(Omegat+Theta);xlabel(t);ylabel(y);编程运行后得到的均匀分布随机信号的仿真图如下：（4）贝努里随机信号的仿真：程序：p=0.5;N=50;ind=find(rand(N,1)p);z(1:N)=1;z(ind)=0;stairs(1:N,z(1:N);axis(0 N-1 -0.5 1.5);title(贝努里随机信号 );说明：P为某事件B发生的概率；用rand（N，1）函数产生一个N个样本的（0,1）均匀分布的列矢量；find(rand(N,1)p)找出满足大于p的（0,1）均匀分布的随机数所在的点所在的下标，并把那些点仿真为某事件B不发生，并用z(ind)=0将其置为0；用z(1:N)=1将仿真为某事件B发生的随机数的点置为1；用阶梯函数stairs(1:N,z(1:N)对其进行仿真绘图；用axis(0 N -0.5 1.5)对坐标进行控制，是图像更加美观。编程运行后得到的贝努里随机信号的仿真图如下：2、几种概率分布的仿真随机信号的特征（1）泊松分布 产生的10000个泊松分布随机数，计算均值和方差Matlab程序如下：y=poissrnd(1,1,10000);%产生泊松随机数10000个，bosong.mjunzhi=mean(y); %求均值disp(均值); %显示均值disp(junzhi);fangcha=var(y); %求方差disp(方差); %显示方差 disp(fangcha);figure(1)title(泊松分布随机数10000个);xlabel(x),ylabel(y)plot(y) 图像如下： 概率密度函数Matlab程序如下： y=poissrnd(1,1,10000);%产生泊松随机数10000个figure(1)subplot(2,1,1)f,yi=ksdensity(y); %概率密度plot(yi,f,-)grid ontitle(10000个泊松分布随机数的概率密度估计)z=poissfit(y) %泊松分布的参数估计y=poissrnd(1,1,100); %产生泊松随机数100个，将坐标放大figure(1)subplot(2,1,2)f,yi=ksdensity(y); %概率密度plot(yi,f,-)grid ontitle(100个泊松分布随机数的概率密度估计)命令行结果：均值 0.9993方差 1.0138程序运行后的图像： 频谱Matlab 程序如下：y=poissrnd(1,1,10000); %产生泊松随机数10000个pingpu=fft(y) %傅里叶变换求频谱figure(1)subplot(2,1,1)plot(abs(pingpu)title(10000个泊松分布随机数的频谱估计)grid onsubplot(2,1,2)y=poissrnd(1,1,100); %产生泊松随机数10000个pingpu=fft(y)plot(abs(pingpu)title(100个泊松分布随机数的频谱估计)grid on程序运行后的图像： 功率谱密度Matlab 程序如下：y=poissrnd(1,1,10000); %产生泊松随机数10000个figure(1)window=hann(10000)periodogram(y,window,512,10000)%功率谱密度title(功率谱密度曲线)程序运行后的图像： 自相关函数Matlab 程序如下：y=poissrnd(1,1,10000); %产生泊松随机数10000个R=xcorr(y)figure(1)subplot(2,1,1)plot(R)grid ontitle(10000个泊松随机数的自相关函数曲线) y=poissrnd(1,1,100); %产生泊松随机数10000个R=xcorr(y)subplot(2,1,2)plot(R)grid ontitle(100个泊松随机数的自相关函数曲线)程序运行后的图像： 确定泊松过程是马尔可夫过程证明： 马尔可夫过程： 如果第n次的取样试验只与第n-1次取样结果有关，而与小于n-1次取样结果无关。 设Xn,n=1,2,为一随机序列，其状态空间E=a1,a2,aN，若对于任意的n，满足 PXn=ain|Xn-1=ain-1,Xn-2=ain-2,X1=ai1=PXn=ain|Xn-1=ain-1 则称Xn为马尔可夫链。 泊松过程具有如下性质：1. PX（t0）=0=1；2. 增量N(tn)-N(tn-1)相互独立3. 增量均匀，即对任意s,t0 PN(s+t)-N(s)=n=PN(t)=n独立增量过程X1,X2, Xn彼此独立，并令X(t0)=X(0)=0Xn+1与所有X(ti),in独立PX(tn+1)=xn+1X(t0)=x0,X(t1)=x1,.,X(tn)=xn= Pxn+Xn+1=xn+1X(tn)=xn=PX(tn+1)=xn+1X(tn)=xn满足马尔可夫链，所以泊松过程是一个马尔可夫过程。（2）高斯分布AN(0,3)正态分布分布 概率密度函数Matlab 程序如下：N=1024;Fs=10000;m=0;%均值sigma2=3;%方差 sigma=sqrt(sigma2);u=randn(N,1);x=m+sigma.*u;%产生所需的高斯随机数%求高斯随机过程的概率密度 t=0:N-1;f,t=ksdensity(x)figure(1);subplot(1,2,1);plot(t,f);axis square;title(概率密度函数曲线);grid on y=normpdf(t,m,sigma);subplot(1,2,2);plot(t,y);axis square;title(理论概率密度函数曲线 );grid on 程序运行后的对应的图像： 频谱Matlab 程序如下：%求高斯随机过程的频谱figure(2);f=(0:N-1)/N*Fs; X=fft(x);%对高斯随机信号进行傅里叶变换XX=abs(X) %取信号的幅度 plot(f(1:N/2),XX(1:N/2);%画出频谱 title(频谱);xlabel(Frequency/Hz); grid on 程序运行后的对应的图像： 功率谱密度Matlab 程序如下：%利用periodogram函数计算功率谱 figure(3);nfft=1024;f=(0:N-1)/N*Fs; window=boxcar(length(R); Pxx,f=periodogram(R,window,nfft,Fs) %加窗的周期图谱估计法plot(f,Pxx); title(功率谱密度);Xlabel(Frequency/Hz) ;grid on 程序运行后的对应的图像： 自相关函数Matlab 程序如下：%求高斯随机过程的自相关函数figure(4);R,lags=xcorr(x,500,unbiased)plot(lags,R);title(自相关函数);grid on 程序运行后的对应的图像： 分析与确定是否是白噪声由图可知：其自相关函数只在某点上产生较大的冲击，冲击大小约为高斯随机过程的方差，在其他点都在零值的附近，近似为零值；在系统的带内，它的功率谱密度基本上是常数；均值近似为零；在广义上满足白噪声的条件，故其属于白噪声。BN(2,3)正态分布分布 概率密度函数Matlab程序：N=1024;Fs=10000;m=2;%均值sigma2=3;%方差 sigma=sqrt(sigma2);u=randn(N,1);x=m+sigma.*u;%产生所需的高斯随机数%求高斯随机过程的概率密度 t=0:N-1;f,t=ksdensity(x)figure(1);subplot(1,2,1);plot(t,f);axis square;title(概率密度函数曲线);grid on y=normpdf(t,m,sigma);subplot(1,2,2);plot(t,y);axis square;title(理论概率密度函数曲线 );grid on 程序运行后的对应的图像： 频谱Matlab 程序如下：%求高斯随机过程的频谱figure(2);f=(0:N-1)/N*Fs; X=fft(x);%对高斯随机信号进行傅里叶变换XX=abs(X) %取信号的幅度 plot(f(1:N/2),XX(1:N/2);%画出频谱 title(频谱);xlabel(Frequency/Hz); grid on 程序运行后的对应的图像： 功率谱密度Matlab 程序如下：%利用periodogram函数计算功率谱 figure(3);nfft=1024;f=(0:N-1)/N*Fs; window=boxcar(length(R); Pxx,f=periodogram(R,window,nfft,Fs) %加窗的周期图谱估计法plot(f,Pxx); title(功率谱密度);Xlabel(Frequency/Hz) ;grid on 程序运行后的对应的图像： 自相关函数Matlab 程序如下：%求高斯随机过程的自相关函数figure(4);R,lags=xcorr(x,500,unbiased)plot(lags,R);title(自相关函数);grid on 程序运行后的对应的图像： 分析与确定是否是白噪声由图可知：其功率谱为零，不满足为一个常数的条件；其自相关函数在零点只在某点上产生较大的冲击，但在其他点都不为零；均值不为零；故其不满足白噪声的条件，故其不属于白噪声。（3）二维正态分布相关计算 边缘密度函数Matlab程序：syms x y; f=exp(-2*(x*x-0.5*x*y+0.25*y*y)/3)/( 2*pi*sqrt(3) ); fy=int(f,x,-inf,inf)命令窗口显示的结果：fy = 281474976710656/6126469605667549*exp(-1/8*y2)*2(1/2)*3(1/2)*pi(1/2)由计算可知：=281474976710656/6126469605667549*2(1/2)*3(1/2)*pi(1/2)=0.1995因此，经过化简后得到， 的数学期望Matlab程序段：x_yf=abs(x+y)*f; Exy=int(int(x_yf,x,-inf,inf),y,-inf,inf)命令窗口显示的结果：Exy = 1125899906842624/6126469605667549*42(1/2)*pi(1/2)化简Exy结果值，得=0.2111 的概率密度分布Matlab程序段：P=int(int(f,x,-inf,0),y,-inf,0)命令窗口显示的结果：P = 140737488355328/18379408817002647*pi*6(1/2)*2(1/2)+70368744177664/6126469605667549*pi*6(1/2)*8(1/2)化简Exy结果值，得=0.33333、二维正态分布随机信号的特征 概率密度函数Matlab程序段：mu1=0,0; Sigma2=1 1; 1 4; % 输入均值向量和协方差矩阵X,Y=meshgrid(-4:0.1:4,-4:0.1:4);xy=X(:) Y(:); %产生网格数据p=mvnpdf(xy,mu1,Sigma2);P=reshape(p,size(X); %求取联合概率密度figure（1）,surf(X,Y,P)title(概率密度)程序运行后的对应的图像： 频谱Matlab程序段：signal=mvnrnd(mu1,Sigma2,100);signal_frequency=fft(signal);h=abs(signal_frequency);i=(0:length(signal_frequency)-1)*100/length(signal_frequency);figure（2）,plot(i,h),grid on;%绘制频谱曲线title(幅频响应)figure（3）,plot(i,angle(signal_frequency),grid on;title(相频响应)程序运行后的对应的图像： 功率谱