北理工随机信号分析实验

实验一随机序列的产生及数字特征估计实验目的1、学习和掌握随机数的产生方法。2、实现随机序列的数字特征估计。实验原理1、随机数的产生随机数指的是各种不同分布随机变量的抽样序列（样本值序列）。进行随机信号仿真分析时，需要模拟产生各种分布的随机数。在计算机仿真时，通常利用数学方法产生随机数，这种随机数称为伪随机数。伪随机数是按照一定的计算公式产生的，这个公式称为随机数发生器。伪随机数本质上不是随机的，而且存在周期性，但是如果计算公式选择适当，所产生的数据看似随机的，与真正的随机数具有相近的统计特性，可以作为随机数使用。(0，1)均匀分布随机数是最最基本、最简单的随机数。(0，1)均匀分布指的是在[0，1]区间上的均匀分布，即U(0，1)。实际应用中有许多现成的随机数发生器可以用于产生(0，1)均匀分布随机数，通常采用的方法为线性同余法，公式如下：序列为产生的(0，1)均匀分布随机数。下面给出了上式的3组常用参数：1、，周期；2、（IBM随机数发生器）周期；3、（ran0）周期；由均匀分布随机数，可以利用反函数构造出任意分布的随机数。定理1.1若随机变量X具有连续分布函数FX(x)，而R为(0,1)均匀分布随机变量，则有由这一定理可知，分布函数为FX(x)的随机数可以由(0,1)均匀分布随机数按上式进行变换得到。2、MATLAB中产生随机序列的函数（1）(0，1)均匀分布的随机序列函数：rand用法：x=rand(m,n)功能：产生m×n的均匀分布随机数矩阵。（2）正态分布的随机序列函数：randn用法：x=randn(m,n)功能：产生m×n的标准正态分布随机数矩阵。如果要产生服从分布的随机序列，则可以由标准正态随机序列产生。（3）其他分布的随机序列MATLAB上还提供了其他多种分布的随机数的产生函数，下表列出了部分函数。MATLAB中产生随机数的一些函数3、随机序列的数字特征估计对于遍历过程，可以通过随机序列的一条样本函数来获得该过程的统计特性。这里我们假定随机序列X(n)为遍历过程，样本函数为x(n)，其中n=0,1,2，…，N-1。那么，X(n)的均值、方差和自相关函数的估计为利用MATLAB的统计分析函数可以分析随机序列的数字特征。（1）均值函数函数：mean用法：m=mean(x)功能：返回按上面第一式估计X(n)的均值，其中x为样本序列x(n)。（2）方差函数函数：var用法：sigma2=var(x)功能：返回按上面第二式估计X(n)的方差，其中x为样本序列x(n)，这一估计为无偏估计。（3）互相关函数函数：xcorr用法：c=xcorr(x,y)c=xcorr(x)c=xcorr(x,y,'opition')c=xcorr(x,'opition')功能：xcorr(x,y)计算X(n)与Y(n)的互相关，xcorr(x)计算X(n)的自相关。option选项可以设定为：'biased'有偏估计，即'unbiased'无偏估计，即按上面第三式估计。'coeff'm=0时的相关函数值归一化为1。'none'不做归一化处理。实验内容1、采用线性同余法产生均匀分布随机数1000个，计算该序列均值和方差与理论值之间的误差大小。改变样本个数重新计算。代码及结果：>>Num=input('Num=');>>N=2^31;>>k=2^16+3;>>Y=zeros(1,num);>>X=zeros(1,num);>>Y(1)=1;>>fori=2:num>>Y(i)=mod(k*Y(i-1),N);>>end>>X=Y/N;>>a=0;>>b=1;>>m0=(a+b)/2;>>sigma0=(b-a)^2/12;>>m=mean(X);>>sigma=var(X);>>delta_m=abs(m-m0);>>delta_sigma=abs(sigma-sigma0);>>plot(X,'k');>>xlabel('n');>>ylabel('X(n)');>>delta_m>>delta_sigma>>axistightNum=1000delta_=0.0110delta_sigma=0.0011Num=5000delta_m=2.6620e-04delta_sigma=0.0020可以看出，样本值取得越大，实际值与理论值越接近，相差越小。2、参数为的指数分布的分布函数为利用反函数法产生参数为0.5的指数分布随机数1000个，测试其方差和相关函数。实验代码及结果：>>R=rand(1,1000);>>lambda=0.5;>>X=-log(1-R)/lambda;>>DX=var(X);>>[Rm,m]=xcorr(X);>>subplot(211);>>plot(X,'k');xlabel('n');ylabel('X(n)');axistight;>>subplot(212);>>plot(m,Rm,'k');xlabel('m');ylabel('R(m)');axistight;DX=4.1201理论上方差的值为1/(0.5^2)=4，实际值为4.1201，因为取样个数有限，导致存在一定偏差。但大体相近。3、产生一组N(1,4)分布的高斯随机数（1000个样本），估计该序列的均值、方差和相关函数。实验代码及结果：>>X=normrnd(1,2,[1,1000]);>>Mx=mean(X);Dx=var(X);>>[Rm,m]=xcorr(X);>>subplot(211);>>plot(X,'k');xlabel('n');ylabel('X(n)');axistight;>>subplot(212);>>plot(m,Rm,'k');xlabel('m');ylabel('R(m)');axistight;Mx=0.9937Dx=3.8938理论上的均值为1，方差为4。而在实验中得到的均值为0.9937，方差为3.8938。考虑到取样点有限，误差可以接受，理论值和实验值基本相同。实验二随机过程的模拟与数字特征一、实验目的1、学习利用MATLAB模拟产生随机过程的方法。2、熟悉和掌握特征估计的基本方法及其MATLAB实现。二、实验原理1、正态分布白噪声序列的产生MATLAB提供了许多产生各种分布白噪声序列的函数，其中产生正态分布白噪声序列的函数为randn。函数：randn用法：x=randn(m,n)功能：产生m×n的标准正态分布随机数矩阵。如果要产生服从分布的随机序列，则可以由标准正态随机序列产生。如果N(0,1)，则。2、相关函数估计MATLAB提供了函数xcorr用于自相关函数的估计。函数：xcorr用法：c=xcorr(x,y)c=xcorr(x)c=xcorr(x,y,'opition')c=xcorr(x,'opition')功能：xcorr(x,y)计算X(n)与Y(n)的互相关，xcorr(x)计算X(n)的自相关。option选项可以设定为：'biased'有偏估计。'unbiased'无偏估计。'coeff'm=0时的相关函数值归一化为1。'none'不做归一化处理。3、功率谱估计MATLAB函数periodogram实现了周期图法的功率谱估计。函数：periodogram用法：[Pxx,w]=periodogram(x)[Pxx,w]=periodogram(x,window)[Pxx,w]=periodogram(x,window,nfft)[Pxx,f]=periodogram(x,window,nfft,fs)periodogram(...)功能：实现周期图法的功率谱估计。其中：Pxx为输出的功率谱估计值；f为频率向量；w为归一化的频率向量；window代表窗函数，这种用法对数据进行了加窗，对数据加窗是为了减少功率谱估计中因为数据截断产生的截断误差，下图列出了产生常用窗函数的MATLAB函数。nfft设定FFT算法的长度；fs表示采样频率；三、实验内容1、按如下模型产生一组随机序列其中是均值为1，方差为4的正态分布白噪声序列。估计过程的自相关函数和功率谱。实验代码及结果：>>y0=randn(1,500);%产生一长度为500的随机序列>>y=1+2*y0;>>x(1)=y(1);>>n=500;>>fori=2:1:n>>x(i)=0.8*x(i-1)+y(i);%按题目要求产生随机序列>>x(n)=0.8x(n-1)+w(n)>>end>>subplot(311);>>plot(x);>>title('x(n)');>>subplot(312);>>c=xcorr(x);%用xcorr函数求x(n)的自相关函数>>plot(c);>>title('R(n)');>>p=periodogram(x);%用periodogram函数求功率谱密度>>subplot(313);>>plot(p);>>title('S(w)');得到长度为500的样本序列分布、自相关函数及功率谱如下：2、设信号为其中，为正态分布白噪声序列，试在N=256和N=1024点时，分别产生随机序列x(n)，画出x(n)的波形并估计x(n)的相关函数和功率谱。实验代码及结果：（1）N=256时>>N=256;>>w=randn(1,N);%用randn函数产生一个长度为256的正态分布白噪声序列>>n=1:1:N;>>f1=0.05;>>f2=0.12;>>x=sin(2*pi*f1*n)+2*cos(2*pi*f2*n)+w(n);%产生题目所给信号>>R=xcorr(x);%求x(n)的自相关函数>>p=periodogram(x);%求x的功率谱>>subplot(311);>>plot(x);title('x(n)');>>subplot(312);>>plot(R);title('R(n)');>>subplot(313);>>plot(p);title('S(w)');得到长度为256的样本序列分布、自相关函数及功率谱：（2）N=1024时>>N=1024;>>w=randn(1,N);%用randn函数产生一个长度为256的正态分布白噪声序列>>n=1:1:N;>>f1=0.05;>>f2=0.12;>>x=sin(2*pi*f1*n)+2*cos(2*pi*f2*n)+w(n);%产生题目所给信号>>R=xcorr(x);%求x(n)的自相关函数>>p=periodogram(x);%求x的功率谱>>subplot(311);>>plot(x);title('x(n)');>>subplot(312);>>plot(R);title('R(n)');>>subplot(313);>>plot(p);title('S(w)');得到长度为1024的样本序列、序列的自相关函数、序列的功率谱如下：实验三随机过程通过线性系统的分析一、实验目的1、理解和分析白噪声通过线性系统后输出的特性。2、学习和掌握随机过程通过线性系统后的特性，验证随机过程的正态化问题。二、实验原理1、白噪声通过线性系统设连续线性系统的传递函数为H()或H(s)，输入白噪声的功率谱密度为SX()=N0/2，那么系统输出的功率谱密度为SY(ω)=|H(ω)|2∙N02 输出自相关函数为RY(τ)=N04π-∞∞|H(ω)|2输出相关系数为γyτ=输出相关时间为τ0=0∞γy输出平均功率为EY2(τ)=N02π0∞|上述式子表明，若输入端是具有均匀谱的白噪声，则输出端随机信号的功率谱主要由系统的幅频特性|H(ω)|决定，不再是常数。2、等效噪声带宽在实际中，常常用一个理想系统等效代替实际系统的H(ω)，因此引入了等效噪声带宽的概念，他被定义为理想系统的带宽。等效的原则是，理想系统与实际系统在同一白噪声的激励下，两个系统的输出平均功率相等，理想系统的增益等于实际系统的最大增益。实际系统的等效噪声带宽为∆ωe=1|H(ω)|max20∞|H(ω)|或∆ωe=12j|H(ω)|max3、线性系统输出端随机过程的概率分布（1）正态随机过程通过线性系统若线性系统输入为正态过程，则该系统输出仍为正态过程。（2）随机过程的正态化随机过程的正态化指的是，非正态随机过程通过线性系统后变换为正态过程。任意分布的白噪声通过线性系统后输出是服从正态分布的；宽带噪声通过窄带系统，输出近似服从正态分布。三、实验内容1、仿真一个平均功率为1的白噪声带通系统，白噪声为高斯分布，带通系统的两个截止频率分别为3kHz和4kHz，估计输出的自相关函数和功率谱密度函数。（假设采样频率为10kHz）实验代码及结果：>>Fs=10000;%采样频率10KHz>>x=randn(1000,1);%产生随机序列>>figure(1);>>subplot(3,1,1);>>plot(x);gridon;%x的分布>>xlabel('t');>>subplot(3,1,2);>>x_corr=xcorr(x,'unbiased');%自相关函数>>plot(x_corr);>>gridon;>>subplot(3,1,3);>>[Pxx,w]=periodogram(x);%功率谱密度>>x_Px=Pxx;plot(x_Px);>>gridon;>>figure(2);>>subplot(2,1,1);>>xlabel('f/Hz');>>[x_pdf,x1]=ksdensity(x);%平滑密度分布函数估计>>plot(x1,x_pdf);>>gridon;>>subplot(2,1,2);>>f=(0:999)/1000*Fs;>>X=fft(x);%求DFT>>mag=abs(X);>>plot(f(1:1000/2),mag(1:1000/2));%幅频特性>>gridon;>>xlabel('f/Hz');>>figure(3);>>subplot(3,1,1);>>[b,a]=ellip(10,0.5,50,[3000,4000]*2/Fs);%构造带通滤波器>>[H,w]=freqz(b,a);>>plot(w*Fs/(2*pi),abs(H));%画出滤波器频率特性>>set(gcf,'color','white')>>xlabel('f/Hz');>>ylabel('H(w)');>>gridon;>>subplot(3,1,2);>>y=filter(b,a,x);%通过带通滤波器>>[y_pdf,y1]=ksdensity(y);%绘出通过后的概率分布>>plot(y1,y_pdf);>>gridon;>>y_corr=xcorr(y,'unbiased');%自相关函数>>subplot(3,1,3);>>plot(y_corr);>>gridon;>>figure(4);>>Y=fft(y);%DFT>>magY=abs(Y);>>subplot(2,1,1);>>plot(f(1:1000/2),magY(1:1000/2));%幅频特性>>gridon;>>xlabel('f/Hz');>>subplot(2,1,2);>>nfft=1024;>>index=0:round(nfft/2-1);%下标>>ky=index.*Fs./nfft;>>window=boxcar(length(y_corr));%矩形窗>>[Pyy,fy]=periodogram(y_corr,window,nfft,Fs);%求功率谱密度>>y_Py=Pyy(index+1);>>plot(ky,y_Py);>>gridon;得到：2、设白噪声通过下图所示的RC电路，分析输出的统计特性。（1）试推导系统输出的功率谱密度、相关函数、相关时间和系统的等效噪声带宽。（2）采用MATLAB模拟正态分布白噪声通过上述RC电路，观察输入和输出的噪声波形以及输出噪声的概率密度。（3）模拟产生均匀分布的白噪声通过上述RC电路，观察输入和输出的噪声波形以及输出噪声的概率密度。（4）改变RC电路的参数（电路的RC值），重做（2）和（3），与之前的结果进行比较。（1）由图中所示电路，根据电路分析的相关知识，可推导出输出功率谱密度为：相关函数为：相关时间为：等效噪声带宽为：（2）实验代码及结果：>>R=100;>>C=0.01;>>b=1/(R*C);>>n=1:1:500;>>h=b*exp(-n*b);>>x=randn(1,1000);>>y=conv(x,h);>>[fyy1]=ksdensity(y)>>subplot(3,1,1);>>plot(x);>>title('x(n)');>>subplot(3,1,2);>>plot(y);>>title('y(n)');>>subplot(3,1,3);>>plot(fy);>>title('fy');（3）实验代码及结果：>>R=100;>>C=0.01;>>b=1/(R*C);>>n=1:1:500;>>h=b*exp(-n*b);>>x=rand(1,1000);>>y=conv(x,h);>>[fyy1]=ksdensity(y);>>subplot(3,1,1);>>plot(x);>>title('x(n)');>>subplot(3,1,2);>>plot(y);>>title('y(n)');>>subplot(3,1,3);>>plot(fy);>>title('fy');实验结果：（4）改变RC值R=200，C=0.01;正态分布：均匀分布：R=10，C=0.01;正态分布:均匀分布:由图可得，系统相关时间与系统带宽成反比。另外，由输入输出波形可以看出，正态随机过程通过一个线性系统后，输出仍服从正态分布。而对于任意分布的白噪声，通过一个线性系统后，输出也服从正态分布。实验四窄带随机过程的产生及其性能测试一、实验目的1、基于随机过程的莱斯表达式产生窄带随机过程。2、掌握窄带随机过程的特性，包括均值（数学期望）、方差、相关函数及功率谱密度等。二、实验原理1.窄带随机过程的莱斯表达式任何一个实平稳窄带随机过程X(t)都可以表示为上式称为莱斯表达式，根据上式可以模拟产生窄带随机过程，具体过程下图所示。2.窄带随机过程包络与相位的概率密度包络的概率密度为fA相位的概率密度为fφφ3.窄带随机过程包络平方的概率密度包络平方的概率密度为fU三、实验内容1、按上图所示结构框图，基于随机过程的莱斯表达式，用MATLAB产生一满足条件的窄带随机过程。实验代码及结果：>>n=1:1:1000;>>h=exp(-n);>>c1=randn(1,1000);>>a=conv(c1,h);>>c2=randn(1,1000);>>b=conv(c2,h);>>fc=10000;>>x=zeros(1,1000);>>fori=1:1000>>x(i)=a(i)*cos(2*pi*fc*i)-b(i)*sin(2*pi*fc*i);>>end>>plot(x);>>title('窄带随机过程');实验结果：2、画出该随机过程的若干次实现，观察其形状。代码同上，得到图形：3、编写MATLAB程序计算该随机过程的均值函数