实验一随机过程的模拟与特征估计.docx
《实验一随机过程的模拟与特征估计.docx》由会员分享，可在线阅读，更多相关《实验一随机过程的模拟与特征估计.docx（21页珍藏版）》请在冰点文库上搜索。
实验一随机过程的模拟与特征估计
实验一：
随机过程的模拟与特征估计
一、实验目的
了解随机过程特征估计的基本概念和方法，学会运用MATLAB软件产生各种随机过程，对随机过程的特征进行估计，并通过实验了解不同估计方法所估计出来的结果之间的差异。
二、实验原理
（1）高斯白噪声的产生
利用MATLAB函数randn产生
（2）自相关函数的估计
MATLAB自带的函数为xcorr（），阐述xcorr的用法
R=xcorr（x,y）或R=xcorr（x,y,’option’）用来求序列x（n）与y（n）的互相关函数
R=xcorr（x）或R=xcorr（x,’option’）用来求序列x（n）的自相关函数
option选项是:
‘biased’有偏估计，
‘unbiased’无偏估计，
‘coeff’m=0的相关函数值归一化为1
‘none’不作归一化处理
（3）功率谱的估计
利用周期图方法估计功率谱，
提示：
MATLAB自带的函数为periodogram（），阐述periodogram（）的用法;
阐述其它谱估计方法的用法。
[Pxx,w]=periodgram（x）
Pxx为对应频率w的功率谱密度值。
[Pxx,w]=periodgram（x,window）
window =boxcar（n）矩形窗（Rectangle Window）
window =triang（n）三角窗（Triangular Window）
window =hanning（n）汉宁窗（Hanning Window）
window =hamming（n）海明窗（Hamming Window）
window =blackman（n）布拉克曼窗（Blackman Window）
window=kaiser（n,beta）恺撒窗（Kaiser Window）
Window代表与x等长度的窗序列，对数据进行加窗。
其它谱估计方法：
相关函数法（BT法）
该方法先由序列x（n）估计出自相关函数R（n）,然后对R（n）进行傅立叶变换,便得到x（n）的功率谱估计。
r=xccor（x）;
R=fft（r）;
Pxx=abs（R）;
平均周期图法和平滑平均周期图法
对于周期图的功率谱估计,当数据长度N太大时,谱曲线起伏加剧,若N太小,谱的分辨率又不好,因此需要改进。
两种改进的估计法是平均周期图法和平滑平均周期图法。
Bartlett法:
Bartlett平均周期图的方法是将N点的有限长序列x（n）分段求周期图再平均。
Welch法:
Welch法对Bartlett法进行了两方面的修正,一是选择适当的窗函数w（n）,并在周期图计算前直接加进去,加窗的优点是无论什么样的窗函数均可使谱估计非负。
二是在分段时,可使各段之间有重叠,这样会使方差减小。
现代谱估计-非参数法
MTM法pmtm：
使用正交窗口来截取获得相互独立的改进周期图法功率谱估计，然后再把这些估计结果结合得到最终的估计。
随着NW的增大，窗的个数增多，会有更多的谱估计，从而谱估计的方差得到减小，但同时带来谱泄露的增大，而且正的谱估计的结果将会有更大的偏差。
MUSIC法pmusic：
基于矩阵特征分解的谱估计非参数方法，它把相关数据矩阵中的信息分类，把信息分配到信号的子空间或噪声的子空间。
它适合于普遍情况下的正弦信号参数估计的方法，是多信号分类法的简称。
特征向量法peig：
也是一种基于矩阵特征分解的谱估计非参数方法，它主要适用于混有噪声的正弦信号的功率谱估计，此方法利用相关矩阵的特征值来对MUSIC法公式中的求和进行加权得到的。
Matlab函数
（4）均值的估计
MATLAB自带的函数为mean（）
（5）方差的估计
MATLAB自带的函数为var（）
（6） AR
（1）模型的理论自相关函数和理论功率谱
对于AR
（1）模型
，自相关函数为
，其功率谱为
。
三、实验容（带*为选作）
1.相关高斯随机序列的产生
按如下模型产生一组随机序列
，其中
为均值为1，方差为4的正态分布白噪声序列。
（1）产生并画出a=0.8和a=0.2的x（n）的波形；
（2）估计x（n）的均值和方差；
（3）估计x（n）的自相关函数，并画出相关函数的图形
当a=0.8时
程序：
clc,clearall;
a=0.8;
w=1+sqrt（4）*randn（1,1000）;%w（n）均值为1方差为4的正态分布白噪声
x
（1）=w
（1）/（sqrt（1-a^2））;%x
（1）初值条件
forn=2:
1000
x（n）=a*x（n-1）+w（n）;
end
subplot（2,1,1）;
plot（x）;%画出x（n）的波形
title（'x（n）波形图'）;
axis（[01000-1020]）;%设定坐标区间
mean=mean（x）%估计x（n）的均值
var=var（x）%估计x（n）的方差
R=xcorr（x,'coeff'）;%估计x（n）自相关函数,归一化
subplot（2,1,2）;
plot（R）;%画出自相关函数
title（'x（n）自相关函数'）;
axis（[02000-0.31]）;
结果：
估计x（n）的均值及方差
mean=
4.6889
var=
10.8210
X（n）的波形及其估计自相关函数图形：
当a=0.2时
程序：
clc,clearall;
a=0.2;
w=1+sqrt（4）*randn（1,1000）;%w（n）均值为1方差为4的正态分布白噪声
x
（1）=w
（1）/（sqrt（1-a^2））;%x
（1）初值条件
forn=2:
1000
x（n）=a*x（n-1）+w（n）;
end
subplot（2,1,1）;
plot（x）;%画出x（n）的波形
title（'x（n）波形图'）;
axis（[01000-610]）;%设定坐标区间
mean=mean（x）%估计x（n）的均值
var=var（x）%估计x（n）的方差
R=xcorr（x,'coeff'）;%估计x（n）自相关函数,归一化
subplot（2,1,2）;
plot（R）;%画出自相关函数
title（'x（n）自相关函数'）;
axis（[02000-0.31]）
结果：
估计x（n）的均值及方差
mean=
1.2322
var=
3.9102
X（n）的波形及其估计自相关函数图形：
对比a=0.8和a=0.2时的情况，当a>0时，当a越大，x（n）变化幅度越大，变化越快，自相关函数下降越缓慢。
2.两个具有不同频率的正弦信号的识别
设信号为
，
，其中
为零均值正态白噪声，方差为
。
（1）假定
，针对
和
两种情况，使用周期图periodogram（）的方法估计功率谱。
取N=1000，由于本题中f<0.5,故周期图函数使用默认采样频率Fs=1Hz
程序：
clc,clearall;
f1=0.05;
f2=0.08;
w=sqrt
（1）*randn（1,1000）;%w（n）均值为0方差为1的正态分布白噪声
n=1:
1000;
x（n）=sin（2*pi*f1*n）+2*cos（2*pi*f2*n）+w（n）;
subplot（2,1,1）;
periodogram（x）;
title（'f1=0.05,f2=0.08时功率谱'）
f1=0.05;
f2=0.20;
y（n）=sin（2*pi*f1*n）+2*cos（2*pi*f2*n）+w（n）;
subplot（2,1,2）;
periodogram（y）;
title（'f1=0.05,f2=0.20时功率谱'）
f与图像上角频率对应关系
从两幅图中可以看出，在f1和f2对应的频率上功率能量明显较大，第一幅图w1=0.1pi，w2=0.16pi，第二幅图w1=0.1pi，w2=0.2pi，且w2对应能量幅度大于w1.
（2）假定
，针对
和
两种情况，用周期图periodogram（）的方法估计功率谱
程序
clc,clearall;
f1=0.05;
f2=0.08;
w1=sqrt
（1）*randn（1,1000）;%w（n）均值为0方差为1的正态分布白噪声
w2=sqrt（4）*randn（1,1000）;%w（n）均值为0方差为4的正态分布白噪声
n=1:
1000;
x（n）=sin（2*pi*f1*n）+2*cos（2*pi*f2*n）+w1（n）;
subplot（2,1,1）;
periodogram（x）;
title（'噪声方差为1时功率谱'）
y（n）=sin（2*pi*f1*n）+2*cos（2*pi*f2*n）+w2（n）;
subplot（2,1,2）;
periodogram（y）;
title（'噪声方差为4时功率谱'）
对比两幅图像，可以看出，噪声方差越大能量越大。
*（3）假定
选用不同的谱估计方法进行估计，并进行比较。
clc,clearall;
f1=0.05;
f2=0.08;
w=sqrt（4）*randn（1,1000）;%w（n）均值为0方差为4的正态分布白噪声
n=1:
1000;
x（n）=sin（2*pi*f1*n）+2*cos（2*pi*f2*n）+w（n）;
%自相关法
R=xcorr（x）;
subplot（5,1,1）;
Pw=fft（R/2000）;
f=2*（0:
length（Pw）-1）/length（Pw）;
plot（f,10*log10（abs（Pw）））;
axis（[01-5050]）;
title（'理论x（n）功率谱'）;gridon;
title（'自相关法'）
subplot（5,1,2）;
periodogram（x）;xlabel（''）;ylabel（''）;
title（'周期图法'）
subplot（5,1,3）;
window=hann（1000）;
periodogram（x,window）;xlabel（''）;
title（'汉宁周期图法'）
subplot（5,1,4）
pwelch（x,512）;xlabel（''）;ylabel（''）;
title（'平滑平均周期图法'）
subplot（5,1,5）
pyulear（x,999）;xlabel（''）;ylabel（''）;
title（'Yule--Walker方法'）
周期图法比自相关方产生图像更为平滑，汉宁窗周期图法产生的波形比无窗函数周期法减少了大量高次谐波，而平滑平均周期图法比汉宁窗周期图法产生波形更为平滑。
从图中看出，f1和f2对应频点能量明显高。
3.理论值与估计值的对比分析
设有AR
（1）模型，
，
W（n）是零均值正态白噪声，方差为4。
●用MATLAB模拟产生X（n）的500个样本，并估计它的均值和方差；
●画出X（n）的理论的自相关函数和功率谱；
●估计X（n）的自相关函数和功率谱。
clc,clearall;
holdoff;
a=-0.8;
w=1+sqrt（4）*randn（1,500）;%w（n）均值为1方差为4的正态分布白噪声
x
（1）=w
（1）/（sqrt（1-a^2））;%x
（1）初值条件
forn=2:
500
x（n）=a*x（n-1）+w（n）;
end
subplot（111）
plot（x）;%画出x（n）的波形
title（'x（n）波形图'）;
axis（[0500-2020]）;%设定坐标区间
mean=mean（x）%估计x（n）的均值
var=var（x）%估计x（n）的方差
fori=1:
500;
R（i）=4*（a^i）/（1-a^2）;
end;
forj=1:
999
ifj<500
R1（j）=R（500-j）;
elseifj>500
R1（j）=R（j-499）;
end
end
end
R2=xcorr（x）;%估计x（n）自相关函数,归一化
figure;
subplot（2,1,1）;
i=1:
999;
plot（i-500,R1（i））;%画出自相关函数
title（'理论x（n）自相关函数'）;
subplot（2,1,2）;
plot（R2）;%画出自相关函数
title（'估计x（n）自相关函数'）;
figure;
subplot（2,1,1）;
Pw=fft（R/50000）;
f=2*（0:
length（Pw）-1）/length（Pw）;
plot（f,10*log10（abs（Pw）））;
title（'理论x（n）功率谱'）;
subplot（2,1,2）;
periodogram（x,[],'twosided'）;
title（'估计x（n）功率谱'）;
axis（[02-4040]）
mean=
0.5657
var=
11.5339
可以看出估计的自相关函数与理论值大致吻合，但高频上有些噪声误差。
但随着样本数的增多，估计曲线接近理论值。
若采样点数为5000
4.随机信号通过线性系统分析
考虑图示系统
假定w为正态分布的随机序列和均匀分布的随机序列，分别估计输出的概率密度。
X[n]=w[n]+0.9w[n-1]-0.1w[n-2]
w为均匀分布的随机序列
clc,clearall;
w=rand（1,1000000）;
x
（1）=w
（1）;
x
（2）=w
（2）+0.9*w
（1）;
forn=3:
1000000
x（n）=w（n）+0.9*w（n-1）-0.1*w（n-2）;
end
subplot（2,1,1）;
ksdensity（x）;
title（'输出概率密度曲线'）;
subplot（212）;
m=-0.5:
0.1:
2
hist（x,m）;
title（'输出序列直方图'）;
w为正态分布的随机序列
可以看出，当输入随机信号服从正态分布，通过线性系统后，输出的概率密度曲线仍为正态分布；而输入随机信号服从均匀分布，输出概率密度曲线为三角形。
四、实验思考题
（1）自相关函数R（m）最大值应该在n=0，用MATLAB估计得到的结果与理论的结果相同吗？
为什么？
MATLAB中，自相关函数最大值在序列中点。
这是因为matlab中数组是从n=1开始定义的。
实际上与理论结果一样，matlab计算出的R（m）相当于实际的函数向右平移一半的序列点数。
（2）对两个正弦信号的识别，识别的效果如何？
信噪比对识别效果有何影响？
在第2问中，我们可以在估计功率谱看到两个明显凸起的频点，与题目中的正弦信号对应。
两个频点相距越近时，需要一定量的采样点才能确保识别。
信噪比越高，识别信号难度减低，所需采用点减少。
5、心得体会
通过此次实验，掌握了使用matlab进行随机序列数字特征估计的方法。
对之前的知识理解加深。
升级会员 