北理工随机信号分析实验

北理工随机信号分析实验
实验一随机序列的产生及数字特征估计一、实验目的
1、学习和掌握随机数的产生方法。

2、实现随机序列的数字特征估计。

二、实验原理
1、随机数的产生
随机数指的是各种不同分布随机变量的抽样序列(样本值序列)。

进行随机信号仿真分析时，需要模拟产生各种分布的随机数。

在计算机仿真时，通常利用数学方法产生随机数，这种随机数称为伪随机数。

伪随机数是按照一定的计算公式产生的，这个公式称为随机数发生器。

伪随机数本质上不是随机的，而且存在周期性，但是如果计算公式选择适当，所产生的数据看似随机的，与真正的随机数具有相近的统计特性，可以作为随机数使用。

(0，1)均匀分布随机数是最最基本、最简单的随机数。

(0，1)均匀分布指的是在[0，1]区间上的均匀分布，即 U(0，1)。

实际应用中有许多现成的随机数发生器可以用于产生(0，1)均匀分布随机数，通常采用的方法为线性同余法，公式如下: y,1,yky(modN)0nn,1
x,y/Nnn
x序列为产生的(0，1)均匀分布随机数。

,,n
下面给出了上式的3组常用参数:
710,，5101、，周期; N10,k7,,
83116,，5102、(IBM 随机数发生器)周期; N2,k23,,,，
9315,，2103、(ran0)周期; N21,k7,,,,
由均匀分布随机数，可以利用反函数构造出任意分布的随机数。

定理 1.1 若随机变量 X 具有连续分布函数F (x)，而R 为(0,1)均匀分布随
机变量，X
则有
,1X,F(R) x
由这一定理可知，分布函数为F (x)的随机数可以由(0,1)均匀分布随机数按上式进行变X
换得到。

2、MATLAB 中产生随机序列的函数
(1)(0，1)均匀分布的随机序列
函数:rand
用法:x = rand(m,n)
功能:产生m×n 的均匀分布随机数矩阵。

(2)正态分布的随机序列
函数:randn
用法:x = randn(m,n)
功能:产生m×n 的标准正态分布随机数矩阵。

2如果要产生服从分布的随机序列，则可以由标准正态随机序列产生。

N(,),,
(3)其他分布的随机序列
MATLAB 上还提供了其他多种分布的随机数的产生函数，下表列出了部分函
数。

MATLAB 中产生随机数的一些函数
3、随机序列的数字特征估计
对于遍历过程，可以通过随机序列的一条样本函数来获得该过程的统计特性。

这里我们
假定随机序列 X (n)为遍历过程，样本函数为x(n)，其中n=0,1,2，…，N-1。

那么，X (n)的
均值、方差和自相关函数的估计为
利用MATLAB 的统计分析函数可以分析随机序列的数字特征。

2
(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' 不做归一化处理。

3
三、实验内容
1、采用线性同余法产生均匀分布随机数1000 个，计算该序列均值和方差与理论
值之间的误差大小。

改变样本个数重新计算。

代码及结果:
>>Num=input('Num=');
>>N=2^31;
>>k=2^16+3;
>>Y=zeros(1,num); >>X=zeros(1,num); >>Y(1)=1;
>>for i=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
>>axis tight
1、 Num=1000 delta_=0.0110 delta_sigma=0.0011
4
0.9
0.8
0.7
0.6
0.5X(n)
0.4
0.3
0.2
0.1
1002003004005006007008009001000
n
2、 Num=5000 delta_m =2.6620e-04 delta_sigma =0.0020
0.9
0.8
0.7
0.6
0.5X(n)
0.4
0.3
0.2
0.1
500100015002000250030003500400045005000
n 可以看出，样本值取得越大，实际值与理论值越接近，相差越小。

5
2、参数为的指数分布的分布函数为 ,
,,x F,1,ex
利用反函数法产生参数为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)');axis tight; >>subplot(212);
>>plot(m,Rm,'k');xlabel('m');ylabel('R(m)');axis tight;
15
10
X(n)
5
1002003004005006007008009001000
n
6000
4000
R(m)
2000
-800-600-400-2000200400600800
m
DX =4.1201
理论上方差的值为1/(0.5^2)=4，实际值为4.1201，因为取样个数有限，导致存在一定
偏差。

但大体相近。

6
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)');axis tight; >>subplot(212);
>>plot(m,Rm,'k');xlabel('m');ylabel('R(m)');axis tight;
5
X(n)0
-5
1002003004005006007008009001000
n
4000
3000
R(m)2000
1000
-800-600-400-2000200400600800
m
Mx =0.9937 Dx = 3.8938 理论上的均值为1，方差为4。

而在实验中得到的均值为0.9937，方差为
3.8938。

考虑到取样点有限，误差可以接受，理论值和实验值基本相同。

7
实验二随机过程的模拟与数字特征
一、实验目的
1、学习利用MATLAB 模拟产生随机过程的方法。

2、熟悉和掌握特征估计的基本方法及其MATLAB 实现。

二、实验原理
1、正态分布白噪声序列的产生
MATLAB 提供了许多产生各种分布白噪声序列的函数，其中产生正态分布白噪
声序列的
函数为randn。

函数:randn
用法:x = randn(m,n)
功能:产生m×n 的标准正态分布随机数矩阵。

2如果要产生服从分布的随机序列，则可以由标准正态随机序列产生。

如果
N(,,,)
2N(0,1)，则。

,，,,,X~N(,)
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(...)
功能:实现周期图法的功率谱估计。

其中:
8
Pxx 为输出的功率谱估计值;
f 为频率向量;
w 为归一化的频率向量;
window 代表窗函数，这种用法对数据进行了加窗，对数据加窗是为了减少功
率谱估计
中因为数据截断产生的截断误差，下图列出了产生常用窗函数的MATLAB函
数。

nfft设定FFT算法的长度;
fs表示采样频率;
三、实验内容
1、按如下模型产生一组随机序列
x(n)0.8x(n1)(n),,，,
其中,(n)是均值为1，方差为 4的正态分布白噪声序列。

估计过程的自相
关函数和功率谱。

实验代码及结果:
>>y0=randn(1,500); %产生一长度为500的随机序列 >>y=1+2*y0;
>>x(1)=y(1);
>>n=500;
>>for i=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)的自相关函
数
9
>>plot(c);
>>title('R(n)');
>>p=periodogram(x); %用periodogram函数求功率谱密度 >>subplot(313); >>plot(p);
>>title('S(w)');
得到长度为500的样本序列分布、自相关函数及功率谱如下:
x(n)
20
-20050100150200250300350400450500
4R(n)x 102
1
001002003004005006007008009001000
S(w)
2000
1000
0050100150200250300
2、设信号为
w(n)其中，为正态分布白噪声序列，试在N =256和f,0.05,f,0.1212
N=1024点时，分别产生随机序列x(n)，画出x(n)的波形并估计x(n)的相关函数和
功率谱。

实验代码及结果:
(1)N=256时
>>N=256;
>>w=randn(1,N); %用randn函数产生一个长度为256的正态分布白噪声序列>>n=1:1:N;
10
>>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的样本序列分布、自相关函数及功率谱:
x(n)
5
-5050100150200250300
R(n)
1000
-10000100200300400500600
S(w)
100
50
0020406080100120140
(2)N=1024时
>>N=1024;
>>w=randn(1,N); %用randn函数产生一个长度为256的正态分布白噪声序列>>n=1:1:N;
11
>>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的样本序列、序列的自相关函数、序列的功率谱如下: x(n)
10
-10020040060080010001200
R(n)
5000
-500005001000150020002500
S(w)
400
200
00100200300400500600
12
实验三随机过程通过线性系统的分析一、实验目的
1、理解和分析白噪声通过线性系统后输出的特性。

2、学习和掌握随机过程通过线性系统后的特性，验证随机过程的正态化问题。

二、实验原理
1、白噪声通过线性系统
设连续线性系统的传递函数为H( )或H(s)，输入白噪声的功率谱密度为
SX( )=N0/2，那么系统输出的功率谱密度为
2S ()=|H()| (3.1) Y
输出自相关函数为
2R ()=H()| (3.2) Y
输出相关系数为
(3.3)
输出相关时间为
= (3.4) 0
输出平均功率为
2E=H()| (3.5)
上述式子表明，若输入端是具有均匀谱的白噪声，则输出端随机信号的功率谱主要由系统的幅频特性|H()|决定，不再是常数。

2、等效噪声带宽
在实际中，常常用一个理想系统等效代替实际系统的H()，因此引入了等效噪声带宽的概念，他被定义为理想系统的带宽。

等效的原则是，理想系统与实际系统在同一白噪声的激励下，两个系统的输出平均功率相等，理想系统的增益等于实际系统的最大增益。

实际系统的等效噪声带宽为
2=H()| (3.6)
13
或
= (3.7) 3、线性系统输出端随机过程的概率分布
(1)正态随机过程通过线性系统
若线性系统输入为正态过程，则该系统输出仍为正态过程。

(2)随机过程的正态化
随机过程的正态化指的是，非正态随机过程通过线性系统后变换为正态过程。

任意分布
的白噪声通过线性系统后输出是服从正态分布的;宽带噪声通过窄带系统，输出近似服从正
态分布。

三、实验内容
1、仿真一个平均功率为1的白噪声带通系统，白噪声为高斯分布，带通系
统的两个截止频率分别为3kHz和4kHz，估计输出的自相关函数和功率谱密度函
数。

(假设采样频率为10kHz)
实验代码及结果:
>>Fs=10000; %采样频率10KHz >>x=randn(1000,1); %产生随机序列
>>figure(1);
>>subplot(3,1,1);
>>plot(x);grid on; %x的分布 >>xlabel('t');
>>subplot(3,1,2);
>>x_corr=xcorr(x,'unbiased'); %自相关函数 >>plot(x_corr);
>>grid on;
>>subplot(3,1,3);
>> [Pxx,w]=periodogram(x); %功率谱密度 >>x_Px=Pxx;plot(x_Px);
>>grid on;
>>figure(2);
>>subplot(2,1,1);
>>xlabel('f / Hz');
>> [x_pdf,x1]=ksdensity(x); %平滑密度分布函数估计 >>plot(x1,x_pdf);
>>grid on;
>>subplot(2,1,2);
>>f=(0:999)/1000*Fs;
>>X=fft(x); %求DFT >>mag=abs(X);
14
>>plot(f(1:1000/2),mag(1:1000/2)); %幅频特性 >>grid on;
>>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)');
>>grid on;
>>subplot(3,1,2);
>>y=filter(b,a,x); %通过带通滤波器 >> [y_pdf,y1]=ksdensity(y); %绘出通过后的概率分布 >>plot(y1,y_pdf);
>>grid on;
>>y_corr=xcorr(y,'unbiased'); %自相关函数 >>subplot(3,1,3);
>>plot(y_corr);
>>grid on;
>>figure(4);
>>Y=fft(y); %DFT >>magY=abs(Y);
>>subplot(2,1,1);
>>plot(f(1:1000/2),magY(1:1000/2)); %幅频特性 >>grid on;
>>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);
>>grid on;
得到:
15
高斯白噪声序列
5
-501002003004005006007008009001000
t
高斯白噪声自相关函数
2
-20200400600800100012001400160018002000
高斯白噪声功率谱密度
2
1
00100200300400500600
高斯白噪声一维概率密度函数
0.4
0.3
0.2
0.1
0-5-4-3-2-1012345
模拟高斯白噪声序列频谱
80
60
40
20
00500100015002000250030003500400045005000
f / Hz
16
带通滤波器
1
0.5H(w) 00500100015002000250030003500400045005000
f / Hz
带通滤波后一维概率密度函数
1
0.5
0-1.5-1-0.500.511.52
限带高斯白噪声自相关函数
0.2
-0.20200400600800100012001400160018002000
2、设白噪声通过下图所示的RC 电路，分析输出的统计特性。

(1)试推导系统输出的功率谱密度、相关函数、相关时间和系统的等效噪声带宽。

(2)采用MATLAB 模拟正态分布白噪声通过上述RC 电路，观察输入和输出的噪声波形以及输
出噪声的概率密度。

(3)模拟产生均匀分布的白噪声通过上述RC 电路，观察输入和输出的噪声波形以及输出噪声
的概率密度。

(4)改变RC 电路的参数(电路的RC 值)，重做(2)和(3)，与之前的结果进行比较。

17
(1)由图中所示电路，根据电路分析的相关知识，可推导出
N输出功率谱密度为: S(w),22222wCR，
,,NRC相关函数为: ,Re4RC
相关时间为: ,,RC
,等效噪声带宽为: B,2RC
(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);
>> [fy y1]=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');
x(n)5
-501002003004005006007008009001000
y(n)2
-2050010001500
fy4
2
00102030405060708090100
18
(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);
>> [fy y1]=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');
实验结果:
x(n)
1
0.5 001002003004005006007008009001000
y(n)
1
0.5
0050010001500
fy
4
2
00102030405060708090100
19
(4)改变RC值
R=200，C=0.01;
正态分布:
x(n)
5
-501002003004005006007008009001000 y(n)
2
-2050010001500
fy
4
2
00102030405060708090100
均匀分布:
x(n)
1
0.5 001002003004005006007008009001000 y(n)
1
0.5
0050010001500
fy
4
2
00102030405060708090100
20
R=10，C=0.01;
正态分布:
x(n)5
-501002003004005006007008009001000 y(n)1
-1050010001500
fy10
5
00102030405060708090100
均匀分布:
x(n)1
0.5
001002003004005006007008009001000
y(n)1
0.5
0050010001500
fy5
00102030405060708090100
由图可得，系统相关时间与系统带宽成反比。

另外，由输入输出波形可以看出，正态随
机过程通过一个线性系统后，输出仍服从正态分布。

而对于任意分布的白噪声，通过一个线
性系统后，输出也服从正态分布。

21
实验四窄带随机过程的产生及其性能测试一、实验目的
1、基于随机过程的莱斯表达式产生窄带随机过程。

2、掌握窄带随机过程的特性，包括均值(数学期望)、方差、相关函数及功率谱密度等。

二、实验原理
1.窄带随机过程的莱斯表达式
任何一个实平稳窄带随机过程X (t)都可以表示为
上式称为莱斯表达式，根据上式可以模拟产生窄带随机过程，具体过程下图所示。

2.窄带随机过程包络与相位的概率密度
包络的概率密度为，服从瑞利分布。

相位的概率密度为，呈均匀分布。

3.窄带随机过程包络平方的概率密度
包络平方的概率密度为0，为指数概率密度函数。

22
三、实验内容
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);
>>for i=1:1000
>> x(i)=a(i)*cos(2*pi*fc*i)-b(i)*sin(2*pi*fc*i);
>>end
>>plot(x);
>>title('窄带随机过程');
实验结果:
窄带随机过程
2
1.5
1
0.5
-0.5
-1
-1.501002003004005006007008009001000
23
2、画出该随机过程的若干次实现，观察其形状。

代码同上，得到图形: 窄带随机过程1.5
1
0.5
-0.5
-1
-1.501002003004005006007008009001000
窄带随机过程1.5
1
0.5
-0.5
-1
-1.5
-201002003004005006007008009001000
24
窄带随机过程
1.5
1
0.5
-0.5
-1
-1.501002003004005006007008009001000
、编写MATLAB程序计算该随机过程的均值函数、自相关函数、功率谱、包3 络、包络平方及相位的一维概率密度，画出相应的图形并给出解释。

实验代码及结果:
>>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);
>>for i=1:1000
>> x(i)=a(i)*cos(2*pi*fc*i)-b(i)*sin(2*pi*fc*i); >>end
>>m=mean(x)
>>figure(1)
>>plot(m);
>>title('均值')
>>R=xcorr(x);
25
>>figure(2)
>>plot(R);
>>title('自相关函数')
>> [S,w]=periodogram(x); >>figure(3)
>>plot(S);
>>title('功率谱密度')
>>B=zeros(1,1000);
>>for i=1:1000
>> B(i)=sqrt(a(i)^2+b(i)^2); >>end
>> [fB2 j]=ksdensity(B); >>figure(4)
>>plot(fB2);
>>title('包络概率密度')
>>B=zeros(1,1000);
>>for i=1:1000
>> B(i)=(a(i)^2+b(i)^2); >>end
>> [fB2 j]=ksdensity(B); >>figure(5)
>>plot(fB2);
>>title('包络平方概率密度') >>for i=1:1000 >> fai(i)=atan(b(i)/a(i)); >>end
>> [fp j]=ksdensity(fai); >>figure(6); >>plot(fp);
>>title('相位一维概率密度函数') 实验结果: 26
均值
1.5
1
0.5
-0.5
-100.20.40.60.811.21.41.61.82
均值m = 0.0038
自相关函数
200
150
100
50
-500200400600800100012001400160018002000 27
功率谱密度0.9
0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
00100200300400500600
包络一维概率密度函数1.5
1
0.5
00102030405060708090100
28
包络平方一维概率密度函数
2.5
2
1.5
1
0.5
00102030405060708090100
相位一维概率密度函数
0.35
0.3
0.25
0.2
0.15
0.1
0.05
00102030405060708090100分析:生成的两个高斯白噪声，分别通过低通滤波器得到a(t)和b(t)。

用莱斯表达式的原理产生一个窄带随机过程。

由实验结果可以看出，窄带随机过程的均值为零，包络服从瑞利分布，相位服从均匀分布，而包络的平方呈指数型分布。

29。