matlab随机过程的非线性变换实验报告

随机过程的线性变换

姓名：徐延林学号：200904013026

专业：电子工程指导教师：谢晓霞

2012年4月17日

一、实验目的

了解随机过程线性变换的基本概念和方法，学会运用MATLAB 软件模拟各种随机过程的线性变换，对其结果进行仿真分析，并通过实验了解不同随机过程经过窄带系统的输出。

二、实验原理

（1）均匀分布白噪声序列利用MATLAB 函数rand 产生；

laplace 分布的白噪声表达式

()

()(0)2

c x m c f x e m --=

=白噪声 据此我们可以产生拉普拉斯白噪声序列。 （2）自相关函数的估计

||1

1

?()()()||N m x

n R m x n m x n N m --==+-∑

MATLAB 自带的函数为xcorr 。 （3）功率谱的估计

先估计自相关函数?()x

R m ，再利用维纳－辛钦定理，功率谱为自相关函数的傅立叶变换：1

(1)

()()N jm x x m N G R m e ωω+-=--=

∑

MATLAB 自带的函数为periodogram 、pyulear 或pburg 。 （4）均值的估计

1

1

1?()N x n m

x n N -==∑

MATLAB 自带的函数为mean 。 （5）方差的估计

1

221

1??[()]N x

x n x n m N σ

-==-∑

MATLAB 自带的函数为var 。

（6） ARMA 模型的理论自相关函数和理论功率谱

对于AR(1)模型()(1)()X n aX n W n =-+，其理论自相关函数和功率谱分别为

222

2()(0)1()(1)m

X X j a R m m a G ae ωσσω-?=≥?-?

??=

?-?

对于ARMA 模型

01201()(1)(2)()()(1)()

N M a X n a X n a X n a X n N b W n bW n b W n M +-+-+?+-=+-+?+- 其理论的功率谱密度为

220

()M

jkw

k k x N jkw

k

k b e

G w a e

σ-=-==∑∑

（7）白噪声过有限系统或宽带信号过窄带系统输出信号成正态分布。

三、实验内容及结果分析 1.PAM 信号的匹配滤波

假定信号为脉冲幅度调制（PAM ）信号，1

0()()M k s k s t A p t kt -==-∑，k A 等概率取

+1和-1两个值，1s t =，信号在信道中传输会受到加性高斯白噪声的污染，在接收端每一个脉冲要判断发射的是“1”还是“0”。 （1） 画出信号、信号加噪声的波形；

（2） 对匹配滤波器输出信号，每隔s t 秒进行取样（在每个脉冲的结尾时刻取

样），取样值与一门限（自行确定）进行比较，超过门限判“1”，低于

门限判“0”，画出匹配滤波器输出的波形，并标出取样值。

（3）产生10000个二进制数字（随机产生），统计输出端检测的误码率。

结果分析：

clear all;clc;

ts=1;s=0;h=0;out=0;wuma=0;c=1;t0=100;

a=2.*randint(1,100)-1; %随机等概率产生+1和-1

t=0:0.01:100;

for i=1:100

p1=rectpuls(t-0.5-i.*ts); %产生单个矩形脉冲

p2=rectpuls(-t-0.5-i.*ts+t0);

f1=a(i).*p1;

f11=a(i).*p1+0.8.*randn(1,10001);

f2=a(i).*p2;

h1=c.*p1; %产生单个矩形脉冲的匹配滤波器

y=conv(f11,h1); %求单个矩形脉冲过匹配滤波器输出波形 out=out+y;

s=s+f1;

h=h+f2;

end

t1=0:0.01:200;

figure;

plot(t1,out);title('非量化输出波形');

for i=1:10000;

if out(i.*2)>0 out(i.*2)=1;

else out(i.*2)=-1; %输出波形量化

end

output(i)=out(i.*2);

if output(i)~=s(i) wuma=wuma+1;

end

end

Pwuma=wuma./10000; %统计误码率

w=0.8.*randn(1,10001);

zs=s+w;

figure;

subplot(2,1,1);plot(t,s);axis([0 100 -5

5]);legend('s(t)');title('random PAM signal')

subplot(2,1,2);plot(t,zs);axis([0 100 -5

5]);legend('s(t)+w(n)');title('random PAM + guass noise signal'); figure;

subplot(2,1,1);plot(0:0.01:99.99,output);axis([0 100 -5 5]);legend('output signal');title('量化后的输出波形');

sc=conv(zs,h); %求所有脉冲和信号的匹配输出 subplot(2,1,2);plot(t1,sc);legend('The output signal');title('经过相参累加器后输出波形');

经过运行程序得到相应结果如下（信噪比1:0.8）

102030405060708090100

-50

5

random PAM signal

102030405060708090100

-50

5

random PAM + guass noise signal

非量化输出波形

量化后输出波形

经过相参累加器后的输出

信噪比对误码率的影响(信号脉宽1s)

信噪1:0.1 1:0.5 1:0.8 1:1.5 1:2 1:3 1:100

从表格中的数据我们可以看到，信噪比越高，误码率越大，这从实际中也是很好理解的。为了研究误码率与信噪比的定量关系，我们对上述数据进行拟合得到如下结果：

由拟合结果看，误码率与信噪比是近似呈反指数关系的，关系式近似为

0.27

0.305x w

error p e

-=

根据我们自己的理解，当信噪比趋于零时，噪声已将信号完全“淹没”，此时的匹配滤波器就没有什么作用了，由于原信号是随机生成的二进制序列，那么此时的输出误码率应该趋近于50%，即理想的误码率与信噪比的关系应该为

0.5x w

error p e

α

-=

信号脉宽对误码率的影响(信噪比1:0.8)

从上面的结果来看，随着信号脉宽的加大，误码率逐渐降低。这是由于匹配滤波器总是在脉冲结束的时刻将信号的能量集中起来，信号的脉宽越大，集中的能量就越多，匹配滤波器的输出信噪比就越大，这就越容易“解码”。

单个脉冲信号的能量信号

2

Q a x t =?（x 表示脉冲信号的幅度，t ?表示脉宽）

所以匹配滤波器输出信号的信噪比正比于信号的脉宽，由上面的分析我们可知误码率是与信噪比呈反指数关系的，所以我们可以合理的猜测，误码率与脉冲信号的脉宽也是呈反指数关系的。其实验拟合效果如下：

由上面的拟合关系我们可以看到，误码率的确与信号脉宽呈反指数关系的，由于实验的随机性，上面的数据只能作为一种定性的描述，不能作为定量的关系式。

脉宽（*0.01s ） 1 2 3 4 5

2. 理论值与估计值的对比分析

2.1设有AR(1)模型，

=--+，

X n X n W n

()0.8(1)()

W(n)是零均值正态白噪声，方差为4。

（1）用MATLAB模拟产生X(n)的500个样本，并绘出波形；

（2）用产生的500个观测点估计()

X n的均值和方差；

（3）画出X(n)的理论的自相关函数和功率谱；

（4）估计X(n)的自相关函数和功率谱。

结果分析：

N=500;u=0;sigma2=4;a=-0.8;

W=u+sqrt(sigma2).*randn(1,N);

for n=1:N-1

x(1)=0;

x(n+1)=a.*x(n)+W(n+1);

end

figure;

stem(x,'.'); %画序列波形

mu=mean(x);

sigma=var(x);

for i=1:500

Rx0(i)=(sigma2.*(a.^i))./(1-(a.^2)); %计算理论自相关函数

end

for j=1:1000

if j<=500

Rx1(j)=Rx0(501-j);

else

Rx1(j)=Rx0(j-500);

end

end%将自相关函数处理成轴对称序列Rx2=xcorr(x)/500;

figure;

subplot(2,1,1),stem(-500:1:499,Rx1,'.'),title('理论自相关函数');

subplot(2,1,2),stem(-499:1:499,Rx2,'.'),title('估计自相关函数');

Pw=fft(Rx0/5000);

f=(0:length(Pw)-1)*1000/length(Pw);

figure;

subplot(2,1,1),plot(f,10*log10(abs(Pw))),title('理论功率谱'); subplot(2,1,2),periodogram(x,[],'twosided',512,1000),title('估计功率谱');

经过运行程序，得到x(n)的序列如下图所示，此时该序列的均值为-0.023，方差为9.952。并且该序列的理论和估计的自相关函数以及功率谱如下图示。

理论自相关函数

估计自相关函数理论功率谱

Frequency (Hz)

P o w e r /f r e q u e n c y (d B /H z )

估计功率谱

从上面的自相关函数序列中我们可以看到：估计的自相关函数基本上与理论的自相关函数轮廓符合，都是逐渐衰减的序列；所不同的是，随着时间差的增大理论的自相关函数序列是单调减小的，而估计的自相关函数不是单调减小的，它是一个波动的减小过程，这主要是因为信号噪声的影响。估计自相关函数的表达式为

||1

1?()()()||N m x

n R m x n m x n N m --==+-∑

它是与序列()

x n的表达式中含有噪声分量，所以这就导致x n的值相关的，而()

估计的相关函数必然是一个波动减小的过程，并且我们可以断定波动的程度是与噪声信号的方差呈正相关的。

从功率谱的理论和估计结果来看，我们同样可以发现，估计的双边功率谱与理论的功率谱基本吻合，但是估计的含有较大的波动性，其原因与上面分析的原因类似。

2.2设有ARMA（2,2）模型：

+---=+---X n X n X n W n W n W n

()0.3(1)0.2(2)()0.5(1)0.2(2)

W n是零均值正态白噪声，方差为4。

()

（1）用matlab模拟产生()

X n的500点样本函数，并绘出波形；

（2）用产生的500个观测点估计()

X n的均值和方差；

（3）画出理论的功率谱；

（4）估计()

X n的相关函数和功率谱。

结果分析：

clear all;clc;

N=500;u=0;sigma2=4;

W=u+sqrt(sigma2).*randn(1,N);

for n=1:N-2

x(1)=0;x(2)=0;

x(n+2)=W(n+2)+0.5.*W(n+1)-0.2.*W(n)-0.3.*x(n+1)+0.2.*x(n);

end

mu=mean(x);

sigma=var(x);

Rx=xcorr(x)/500;

figure;

subplot(2,1,1),stem(x,'.');title('x序列波形');

subplot(2,1,2),stem(-499:1:499,Rx,'.'),title('估计自相关函数');

f=0:5:1000;

middle1=1+0.5.*exp(-i.*f)-0.2.*exp(-i.*2.*pi.*f);

middle2=1+0.3.*exp(-i.*f)-0.2.*exp(-i.*2.*pi.*f);

g=sigma2.*(abs(middle1./middle2)).^2; G=-10.*log(g)./log(10)-20; figure;

subplot(2,1,1),plot(f,G),title(理论功率谱');grid on ;

subplot(2,1,2),periodogram(x,[],'twosided',512,1000),title('估计功率谱');

经过运行程序，得到x(n)的序列如下图所示，此时该序列的均值为0.061，方差为3.927。并且该序列的理论和估计的自相关函数以及功率谱如下图示。

x 序列波形图

估计自相关函数

理论功率谱

Frequency (Hz)

P o w e r /f r e q u e n c y (d B /H z )

估计功率谱

3.随机过程过线性系统分析

假定滤波器为课本图3.7给出的RC 电路（低通滤波器）。 （1） 将低通滤波器转化成数字低通滤波器；

（2） 产生一组均匀分布的白噪声序列，让这组白噪声序列通过数字低通滤

波器，画出输出序列的直方图，并与输出的理论分布进行比较；

（3） 产生一组拉普拉斯分布的白噪声序列，让这组白噪声序列通过数字低

通滤波器，画出输出序列的直方图，并与输出的理论分布进行比较；

（4） 改变滤波器的参数（电路RC 值），重做（1）~（3），并与前一次的结

果进行比较。

结果分析：

clear all ;clc;

R=1000;C=10.^(-6);w0=1./(R.*C); w=0:1:1000;

H=1./(1+i.*w.*R.*C);

%均匀分布白噪声过低通系统

w1=rand(1,1001); W1=fft(w1,1001); OUT1=W1.*H; out1=ifft(OUT1);

[f1,xi1]=ksdensity(out1); k1=0:0.01:1; figure;

plot(w,abs(H));grid on ;title('一阶低通数字滤波器'); figure;

subplot(2,1,1);plot(xi1,f1);grid on ;title('均匀分布白噪声输出概率分布'); subplot(2,1,2);hist(out1,k1);grid on ;title('均匀分布白噪声输出分布直方图');

%拉普拉斯分布白噪声过低通系统 s=1/2*(2*rand(1,1001)-1);

w2=-0.20*sign(s).*log(1-2*abs(s)); W2=fft(w2,1001); OUT2=W2.*H; out2=ifft(W2);

[f2,xi2]=ksdensity(out2); k2=-1.5:0.01:1.5; figure;

subplot(2,1,1);plot(xi2,f2);grid on ;title('拉普拉斯分布白噪声输出概率分布'); subplot(2,1,2);hist(out2,k2);grid on ;title('拉普拉斯分布白噪声输出分布直方图');

该题中的低通滤波器的截止角频率

1c w RC

=

即系统带宽c BW w =，改变,R C 的值就可以改变系统带宽，下面我们研究不同带宽下系统的输出。

1）1,1R k C uF =Ω=

0.7

0.75

0.8

0.85

0.9

0.95

1

一阶数字滤波器

均匀分布白噪声输出概率分布

均匀分布白噪声输出分布直方图

拉普拉斯白噪声输出概率分布

拉普拉斯白噪声输出分布直方图

2）10,1R k C uF =Ω=

00.1

0.20.30.40.50.60.70.80.91一阶数字滤波器

均匀分布白噪声输出概率分布

均匀分布白噪声输出分布直方图

拉普拉斯白噪声输出概率分

布

拉普拉斯白噪声输出分布直方图

3）10,100R k C uF =Ω=

00.1

0.20.30.40.50.60.70.80.91一阶

数字滤波器

均匀分布白噪声输出概率分布

均匀分布白噪声输出分布直方图

拉普拉斯白噪声输出概率分布

拉普拉斯白噪声输出分布直方图

从上面的结果中我们可以看到：随着R,C 值的增大，滤波器的截止频率降低，系统带宽越来越小，输出的信号分布越来越接近于高斯分布。

理论上，白噪声过限带系统输出应该为高斯分布。在试验中由于低通滤波器只是一阶的，滤波效果不是特别好，这就导致了输出波形只是近似高斯分布，而且输出结果又一定的随机性，上图只是挑选了效果比较好的几幅图；然而，从另一方面看，随着系统带宽的减少，输出越来越接近高斯分布的事实，也从另一个侧面验证了白噪声过限带系统输出接近高斯分布的结论。

四、心得体会

总体上来说，感觉这次实验的挑战性比较大。其一是因为一些实验结论在事先我们是不知道的，探究性比较强，比如说PAM 信号的误码率问题；其二是感觉时间有点紧，因为现在各种事情好多，真有点分身乏术的感觉。

俗话说有挑战才有动力，也正因为有挑战所以当自己完成任务的时候成就感就比较强。

这次实验一方面加深了我对所学知识的理解，另一方面也培养了我对matlab 的学习兴趣。一个人为了一个问题翻看各种书籍其实回想起来也是一种乐趣，会有问题得到解答时的兴奋，也会有找不到答案的沮丧，各种滋味也只有自己体会过了才会有收获吧。

最后，我觉得探究性的题目做起来更能调动我的积极性，得到的结论也许是错的，但只要是自己未知的，那种感觉也是很好的。我们的时间和精力是有限的，如果我们有与时间和精力赛跑的心情，我相信我们的效率会大大提高的。