随机过程数学建模分析

随机过程数学建模分析

任何通信系统都有发送机和接收机，为了提高系统的可靠性，即输出信噪比，

通常在接收机的输入端接有一个带通滤波器，信道内的噪声构成了一个随机过

程，经过该带通滤波器之后，则变成了窄带随机过程，因此，讨论窄带随机过程

的规律是重要的。

一、窄带随机过程。

一个实平稳随机过程X（t）,若它的功率谱密度Sx（3）具有下述性质：

.=件（3），工|3|工9+

0others

中心频率为O）C，带宽为△3=230，当△3<<3c时，就可认为满足窄带条件。

若随机过程的功率谱满足该条件则称为窄带随机过程。若带通滤波器的传输函数

满足该条件则称为窄带滤波器。随机过程通过窄带滤波器传输之后变成窄带随机

过程。

图1为典型窄带随机过程的功率谱密度图。若用一示波器来观测次波形，则

可看到，它接近于一个正弦波，但此正弦波的幅度和相位都在缓慢地随机变化，

图2所示为窄带随机过程的一个样本函数。

图】典型窄带随机过程的功率谱密度图

图2窄带随机过程的•个样本函数

二、窄带随机过程的数学表示

1、用包络和相位的变化表示

由窄带条件可知，窄带过程是功率谱限制在。.附近的很窄范围内的一个随机

过程,从示波器观察（或由理论上可以推知）:这个过程中的一个样本函数（一个实现）

的波形是一个频率为人且幅度和相位都做缓慢变化的余弦波。

写成包络函数和随机相位函数的形式:

X(t)=A(t)*cos[(Dct+O(t)]

其中:A⑴称作X(t)的包络函数;①⑴称作X⑴的随机相位函数。包络随时间做

缓慢变化，看起来比较直观，相位的变化，则看不出来。

2^莱斯(Rice)表示式

任何一个实平稳随机过程X(t)都可以表示为：

X(t)=Ac(t)coso)ct-As(t)sino)ct

其中同相分量：

Ac(t)=X(t)cos(pt=X⑴cos(oct+^(Osincoct=LPfX(t)*2cos(oct]

正交分量：

As(t)=X(t)sin(pt=^(Ocoscoet—X(t)sincoct=LP[-X(t)*2sina)ct]

(欠(t)为X(£)的希尔伯特变换。LP[A]表示取A的低频部分)。4⑴和As⑴都

是实随机过程，均值为0,方差等于X⑴的方差。

三、窄带随机过程仿真建模要求

1、用Matlab编程仿真窄带随机信号：X⑴-(1+A(t))*cos(5t+(p)+n⑴。其中

包络A⑴频率为IKHz,幅值为IV。载波频率为：4KHz,幅值为IV,<p是一个

固定相位，n⑴为高斯白噪声，采样频率设为16KHz。实际上，这是一个带有载

波的双边带调制信号。

2、计算窄带随机信号的均值、均方值、方差、概率密度、频谱及功率谱密度、

相关函数，用图示法来表示。

3、窄带系统检测框图如图3所示。

图3窄带系统检测框图

4、低通滤波器设计：

低通滤波器技术要求：通带截止频率IKHz,阻带截止频率2KHz。过渡带:

IKHz,阻带衰减：＞35DB,通带衰减：＜1DB,采样频率：＜44.1KHz

5、计算a点、b点、Ac(t)、As(t)＞y(t)的均值、均方值、方差、频谱及功率谱

密度、相关函数，用图示法来表示。

四、建模仿真过程及结果(程序见附件)

1、根据要求得到X(t)的表达式：

x=(1+a).*cos(2*pi*4()(X)*t+2)+noisy/10;

其中：noisy为高斯白噪声，由wgn函数生成，

a=cos(2*pi*1000*t),

均值：Ex=mean(x),

方差：Dx=var(x),

计算可得：X⑴的均值为0.0019,

X(t)的方差为0.7590。

如图4所示，其中蓝色线为X⑴一个样本的时域波形，红色点连成的线为X(t)

的均值，绿色点连成的线为X(t)的方差。

图4窄带随机信号时域波形

2、求X(t)的概率密度，方法是将最大最小区间分成14等份，然后分别计算各个

区间的个数，如图02中柱形条所示，利用曲线拟合，得到合适的概率密度函数。

为了得到光滑的曲线，利用了多项式拟合，经过测试，9次拟合曲线比较符合要

求，获得的曲线如图5中曲线所示：

图5X⑴的概率分布密度函数

3、对X(t)进行频谱分析，在Matlab中，利用fft函数可以很方便得求得X(t)的频

谱，然后用abs和angle函数求得幅值和相位，画出图像如图6所示：

图6X(t)的频谱图

4、求X⑴的自相关函数，用xcorr函数求出自相关序列，得到X⑴自相关函数的

时域波形，如图7所示。

图7X⑴自相关函数的时域波形

5、对X(t)自相关函数进行fft变换，得到X(t)的功率谱密度，如图8所示。

图8X⑴的功率谱密度

6、建立滤波器，建立一个巴特沃思滤波器，对产生的x⑴进行检测。滤波器的幅

度谱和相位谱所示：

图9地通滤波器的幅度谱和相信诺

7、求Ac(t)的统计特性，Ac⑴为X⑴*2coso)ct通过低通滤波器的信号，

Ac(t)的均值Eh=-0.40754(带有直流分量)，

Ac⑴的均方值是E2h=0.2458

Ac(t)的方差Dh=0.0798

Ac(t)的波形如图10、图11所示：

图ioA⑴的时域波形图和频谱图

图11A,⑴的自相关函数的时域波形图和Ac⑴的功率谱密度

8、求As⑴的统计特性，As(。为X(t)*2cos3ct通过低通滤波器的信号，

As(t)的均值Eh=0.8972(带有直流分量)，

As⑴的均方值是E2h=1.1565

As(t)的方差Dh=0.3518

As⑴的波形如图13、图14所示：

图13As⑴的时域波形图和频谱图

图14As⑴的自相关函数的时域波形图和As⑴的功率谱密度

9^求出Y(t)的统计特性，Y(t)=Ac(t)coscoct-As(t)sino)ct,

其统计特性如下

输出信号Y⑴的均值Eh=-4.401le-004s

输出信号Y⑴的均方值E2h=3.0280

输出信号Y⑴的方差Dh=3.0303

Y⑴的仿真图形如图15、图16所示。

图15Y⑴的时域波形图和频谱图

图16Y⑴的自相关函数的时域波形图和Y(l)的功率谱密度

附件:

clc

fs=16000;%设定采样频率

N=1300;

n=0:N-l;%取的样本点数

t=n/fs;%获得以1/16000为时间间隔采样序列

noisy=wgn(I,N,0);%产生高斯白噪声

a=cos(2*pi*1000*t);%获取A⑴的采样点

x=(14-a).*cos(2*pi*4(X)0*t+2)+noisy/101%获取x⑴的采样点

%以t为横坐标画出x⑴的时域图型

figurs(l);subplot(2,l,l);plot(n,x);

axis([0140-33]),xlabelC采样点工ylabWX(l)/V);lill«窄带随机信号波形);giidon,

%求*⑴的统计特性并画出来

dispCX(t)的均值为);Ex=mcan(x);disp(Ex);%求X⑴均值

holdon：plot(n,Ex,'r.');

dispCX(t)的方差为);Dx=var(x);disp(Dx);%求x⑴方差

holdon;plot(n,Dx,'g.');

%画出X⑴的概率分布函数

each=Iinspace(min(x),max(x),14);%将最大最小区间分成14等份，然后分别计算各个区间的个数

nr=hist(x,each);%计算各个区间的个数

nr=nr/length(x);%计算各个区间的个数归一化

subplot(2,l,2);p=polyfit(each,nr,9);%画出概率分布直方图

bar(each.nr)：%多项式拟合

holdon;plot(each.nr,'g')

eachi=-2:0.1:2;

nri=polyval(p,eachi);

plot(eachi,nri,'r')

axisiight;lille('X⑴概率密度分布');xlabel('X(t)');ykibel('P(x)');gridon;

%XtX⑴进行频谱分析

Fx=fTt(x,N);%对*⑴进行fft变换，在0〜16000区间内得到2N-1个频率值

magn=abs(Fx);%求x⑴幅值

xanglc=angle(Fx);%求X⑴相位

labelang=(O:length(x)-l)*16000/length(x);%在0〜16(X)0区间内求横坐标刻度

figu32);plot(labelang,magn*IO);%在0-16000区间内做频谱和相位图

axis([016000-0.5600]);xlabel('频率/Hz');ylabel('幅值');tille('X⑴频谱图'):gridon;

%求*(1)的自相关函数

[c,)ags]=xcorr(x,'cocff);%求出自相关序列

figurs(3);subplot(2,l,l);plot(lags/fs,c);%在时域内画自相关函数

axislight;xlabel(T);ylabel('Rx(T)，);出IR'X⑴的自相关函数');gridon;

%求*⑴的功率谱密度

long=lcngth(c);

Sx=fft(c,long);

labclx=(0:long-1)*2*pi;

p!ot_magn=10*log10(abs(Sx));

subplot(2,1,2);plot(labelx,plot_magn);%画功率谱密度

axisiighl;xkibel('w');ykibel('Sx(w)');liUerX⑴的功率谱密度');gridon;

%窄带系统检测

z1=2*cos(2*pi*4000*t);

z2=-2.*sin(2*pi*4000*t);

Ac=zl.*x;%滤波后生成Ac(t)

As=z2.*x;%滤波后生成As⑴

y=Ac.*cos(2*pi*4000*t)-As.*sin(2*pi*4000*t);

%滤波器设计

f_p=1000;Cs=1600:R_p=l;R_s=35;%设定滤波器参数；通、阻带截止频率，通、阻带衰减

Ws=2*Ls/fs;Wp=2*Lp/fs;%频率归一化

[nAVn]=buttord(Wp,Ws,R_p.R_s);%采用巴特沃思滤波器

[b,a]=butter(n,Wn);%求得滤波器传输函数的多项式系数

figure(4);

[H,W]=freqz(b,a);%求得滤波器传输函数的幅频特性

subp]ot(2J,l);plot(W*fs/(2*pi),abs(H));%在0〜2pi区间内作幅度谱

tille(低通滤波器幅度谱);gridon;

subplot(2,l,2);plot(W*fs/(2*pi),angle(H));%在。〜2pi区间内作相位谱

lillef低通滤波器相位谱');gridon;

%求Ac⑴滤波后的统计特性

mc=iilter(b,a,Ac);%上支路通过流波器Ac(t)

disp(，Ac(t)的均值);Eh=mean(mc)%求Ac⑴的均值

dispCAc⑴的均方值是');E2h=mc*mc7N%求Ac⑴的均方值

disp('Ac⑴的方差'):Dh=var(mc)%求Ac⑴的方差

%画Ac(l)的时域波形

figured);subplot(2,l,l);n=0:N-l;plot(n.mc);

axis([0300/l]);xlabel('采样点);ylabel('幅值然itle('Ac⑴的时域波形');gridon;

%画Ac⑴的频谱图

yc=fit(mc,length(nic));%对Ac⑴进行fft变换

longc=lcngth(yc);%求傅里「十变换后的序列长度

labelx=(0:longe-1)*16000/longc;

magnl=abs(yc);%求Ac⑴的幅值

subplot(2,1,2);plot(labelx.magnl);%画Ac⑴的频谱图

axislight;xlabel。频率(Hz)');ylabelf幡值')；tille('Ac(t)频谱图')；gridon;

%求A<:(。的自相关函数

[cIJags1]=xcorr(mc.,cocff);%求出Ac⑴的自相关序列

figur?⑺；subplot(2,l,l);plot(lagsl/fsxl);%在时域内画Ac(t)的自相关函数

xlabel(T);ylabel('Rx(T)');axislight;

iitle('Ac(t)的自相关困数')；

gridon;

%求Ac⑴的双边功率谱

Sac=fft(cl,length(cl));%对Ac⑴的自相关函数进行傅里叶变换

magnc=abs(Sac);%求Ac⑴的双边功率谱幅值

long=length(Sac);%求傅里叶变换后的序列长度

labelc=(0:long-l)*16000/Iong;

subplot(2,1,2);plot(labelc,1O*log10(magnc));%画Ac⑴的自相关函数频谱即为Ac⑴的双边功率谱

xlabel('频率(Hz)');ylabel('功率谱(dbW));axistighl;title('Ac⑴的双边功率谱);gridon;

%求得As⑴的统计特性

ms=filter(b,a,As);%对下支路信号进行滤波得As⑴

disp(As(t)的均值＞；Eh=mcan(ms)%求As(t)的均值

dispCAs⑴的均方值是')；E2h=ms*ms7N%求As⑴的均方值

disp('As(t)的方差')；Dh=var(ms)%求As⑴的方差

%作As⑴的时域波形

figure(8);subplot(2,l,1);n=0:N-l;plot(n,ms);%画出As⑴的时域波形

axis([O300-0.52]);xlabel('采样点');ylabel('幅值');tiUe('As⑴的时域波形');gridon;

%对As⑴进行FFT变换并做频谱图

ys=fft(ms,lcngih(ms));%对As⑴进行ff(变换

longs=length(ys);%求傅里叶变换后的序列长度

Iabclx=(0:longs-1)*16000/longs;

magn2=abs(ys);%求As(t)的幅值

subp]ot(2,l,2);plot(labelx,magn2);%画出As⑴的频谱图

axislight;xlabel('施率(Hz)');ylabel('幅值);title('As(t)的频谱N'):gridon;

%求As⑴的自相关函数

[c2,lags2]=xcorr(ms,'coeff);%求出As⑴的自相关序列

figur3(9);subplot(2,1,1);plot(lags2/fs,c2);%画出As(t)自相关函数的时域波形

xlabel(T);ylabel('Rx(T)');axistight;title('As⑴的的自相关函数以gridon;

%求As⑴的双边功率谱

Sas=fft(c2,length(c2));%对As⑴的自相关函数进行傅里叶变换

magnc=abs(Sac);%求As⑴的双边功率谱幅值

long=length(Sas);%求傅里叶变换后的序列长度

labels=(0:long-l)*16000/long;

subp]ot(2.1.2):plotdabelc.10*logIO(magnc)):%画As⑴的自相关函数频谱

xlabd(濒率(Hz/);ylabel(‘功率谱(dbW));axistighl;title('As⑴的双边功率谱);

%求y⑴的统计特性

dispC输出信号Y⑴的均值):Eh=mcan(y)%求输出信号Y⑴的均值

disp(输出信号Y⑴的均方值);E2h=y*y7N%求输出信号Y⑴的均方值

disp(输出信号Y(t)的方差RDh=var(y)%求输出信号Y(l)的方差

%作输出信号Y⑴的时域波形

figuc(10);subplot(2.1.1