Matlab仿真窄带随机过程

随机过程数学建模分析

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

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

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

的规律是重要的。

一、窄带随机过程。

一个实平稳随机过程X(t)，若它的功率谱密度具有下述性质：

中心频率为ω，带宽为△ω=2ω，当△ω<<ω 时，就可认为满足窄带条件。

c0c

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

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

过程。

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

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

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

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

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

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

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

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

c

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

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

c

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

X(t)=A(t)*cos[ωt+ Φ(t)]

c

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

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

2、莱斯（Rice）表示式

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

X(t)=A(t) cosωt-A(t) sinωt

ccSc

其中同相分量：

At＋ sinωt=LP[X(t) *2cosωt]

cccc

(t)= X(t) cosφt= X(t) cosω

正交分量：

A(t) = X(t)sinφt= cosωt— X(t) sinωt= LP[-X(t) *2sinωt]

Sccc

（LP[A]表示取A的低频部分）。A(t)和A(t)都

cS

是实随机过程，均值为0，方差等于X(t)的方差。

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

1、用Matlab 编程仿真窄带随机信号：X(t)=(1+ A(t))*cos(ωt+φ)+n(t)。 其中

c

包络A(t)频率为1KHz，幅值为l V。载波频率为：4KHz，幅值为l V，φ是一个

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

波的双边带调制信号。

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

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

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

图3 窄带系统检测框图

4、低通滤波器设计：

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

1KHz，阻带衰减：>35DB，通带衰减：<1DB，采样频率：≤44.1KHz

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

cS

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

四、建模仿真过程及结果（程序见附件）

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

x= (l+a) .*cos (2*pi*4000*t+2) +noisy/10;

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

a=cos (2*pi*l000*t)，

均值：Ex=mean (x)，

方差：Dx=var (x)，

计算可得：X(t)的均值为0.0019，

X(t)的方差为0.7590。

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

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

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

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

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

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

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

图5 X(t)的概率分布密度函数

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

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

图6 X(t)的频谱图

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

时域波形，如图7所示。

图7 X(t)自相关函数的时域波形

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

图8 X(t)的功率谱密度

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

度谱和相位谱所示：

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

7、求A(t)的统计特性，A(t)为X(t) *2cosωt通过低通滤波器的信号，

ccc

A(t)的均值Eh = -0.4075 4（带有直流分量），

c

A(t)的均方值是E2h =0.2458

c

A(t)的方差Dh = 0.0798

c

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

c

图10 A(t)的时域波形图和频谱图

c

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

c

8、求A(t)的统计特性，A(t)为X(t) *2cosωt通过低通滤波器的信号，

SSc

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

S

A(t)的均方值是E2h = 1.1565

S

A(t)的方差Dh = 0.3518

S

A(t)的波形如图13、图14所示：

S

图13 A(t)的时域波形图和频谱图

S

图14 A(t)的自相关函数的时域波形图和A(t)的功率谱密度

SS

9、求出Y(t)的统计特性，Y (t)=A(t) cosωt-At，

ccSc

(t) sinω

其统计特性如下

输出信号Y(t)的均值Eh = -4.4011e-004s

输出信号Y(t)的均方值E2h = 3.0280

输出信号Y(t)的方差Dh = 3.0303

Y(t)的仿真图形如图15、图16所示。

图15 Y(t)的时域波形图和频谱图

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

附件：

clc

fs=16000; %设定采样频率

N=1300;

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

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

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

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

x=(1+a).*cos(2*pi*4000*t+2)+noisy/10; %获取x(t)的采样点

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

figure(1); subplot(2,1,1); plot(n,x);

axis([0 140 -3 3]);xlabel('采样点');ylabel('X(t)/V');title('窄带随机信号波形');grid on;

%求X(t)的统计特性 并画出来

disp('X(t)的均值为'); Ex=mean(x); disp(Ex);%求X(t)均值

hold on; plot(n,Ex,'r.');

disp('X(t)的方差为');Dx=var(x); disp(Dx);%求x(t)方差

hold on; plot(n,Dx,'g.');

%画出X(t)的概率分布函数

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

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

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

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

bar(each,nr); %多项式拟合

hold on; plot(each,nr,'g')

eachi=-2:0.1:2;

nri=polyval(p,eachi);

plot(eachi,nri,'r')

axis tight;title('X(t)概率密度分布');xlabel('X(t)');ylabel('P(x)');grid on;

%对X(t)进行频谱分析

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

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

xangle=angle(Fx); %求X(t)相位

labelang=(0:length(x)-1)*16000/length(x); %在0~16000区间内求横坐标刻度

figure(2); plot(labelang,magn*10); %在0~16000区间内做频谱和相位图

axis([0 16000 -0.5 600]); xlabel('频率/Hz');ylabel('幅值');title('X(t)频谱图');grid on;

%求X(t)的自相关函数

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

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

axis tight; xlabel('T');ylabel('Rx(T)');title('X(t)的自相关函数');grid on;

%求X(t)的功率谱密度

long=length(c);

Sx=fft(c,long);

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

plot_magn=10*log10(abs(Sx));

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

axis tight;xlabel('w');ylabel('Sx(w)');title('X(t)的功率谱密度');grid on;

%窄带系统检测

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

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

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

As=z2.*x; %滤波后生成As(t)

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

%滤波器设计

f_p=1000;f_s=1600;R_p=1;R_s=35; %设定滤波器参数; 通、阻带截止频率,通、阻带衰减

Ws=2*f_s/fs;Wp=2*f_p/fs; %频率归一化

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

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

figure(4);

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

subplot(2,1,1); plot(W*fs/(2*pi),abs(H)); %在0～2pi区间内作幅度谱

title('低通滤波器幅度谱'); grid on;

subplot(2,1,2); plot(W*fs/(2*pi),angle(H)); %在0～2pi区间内作相位谱

title('低通滤波器相位谱'); grid on;

%求Ac(t)滤波后的统计特性

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

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

disp('Ac(t)的均方值是');E2h=mc*mc'/N %求Ac(t)的均方值

disp('Ac(t)的方差');Dh=var(mc) %求Ac(t)的方差

%画Ac(t)的时域波形

figure(6); subplot(2,1,1); n=0:N-1; plot(n,mc);

axis([0 300 -1 1]);xlabel('采样点');ylabel('幅值');title('Ac(t)的时域波形');grid on;

%画Ac(t)的频谱图

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

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

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

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

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

axis tight; xlabel('频率(Hz)'); ylabel('幅值'); title('Ac(t)频谱图'); grid on;

%求Ac(t)的自相关函数

[c1,lags1]=xcorr(mc,'coeff'); %求出Ac(t)的自相关序列

figure(7); subplot(2,1,1); plot(lags1/fs,c1); %在时域内画Ac(t)的自相关函数

xlabel('T');ylabel('Rx(T)');axis tight;

title('Ac(t)的自相关函数');

grid on;

%求Ac(t)的双边功率谱

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

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

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

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

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

xlabel('频率(Hz)');ylabel('功率谱(dbW)');axis tight;title('Ac(t)的双边功率谱');grid on;

%求得As(t)的统计特性

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

disp('As(t)的均值'); Eh=mean(ms) %求As(t)的均值

disp('As(t)的均方值是'); E2h=ms*ms'/N %求As(t)的均方值

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

%作As(t)的时域波形

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

axis([0 300 -0.5 2]); xlabel('采样点');ylabel('幅值');title('As(t)的时域波形');grid on;

%对As(t)进行FFT变换并做频谱图

ys=fft(ms,length(ms)); %对As(t)进行fft变换

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

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

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

subplot(2,1,2); plot(labelx,magn2); %画出As(t)的频谱图

axis tight;xlabel('频率(Hz)');ylabel('幅值');title('As(t)的频谱图');grid on;

%求As(t)的自相关函数

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

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

xlabel('T');ylabel('Rx(T)');axis tight;title('As(t)的的自相关函数');grid on;

%求As(t)的双边功率谱

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

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

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

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

subplot(2,1,2); plot(labelc,10*log10(magnc)); %画As(t)的自相关函数频谱

xlabel('频率(Hz)');ylabel('功率谱(dbW)');axis tight;title('As(t)的双边功率谱');

% 求y(t)的统计特性

disp('输出信号Y(t)的均值');Eh=mean(y) %求输出信号Y(t)的均值

disp('输出信号Y(t)的均方值');E2h=y*y'/N %求输出信号Y(t)的均方值

disp('输出信号Y(t)的方差');Dh=var(y) %求输出信号Y(t)的方差

%作输出信号Y(t)的时域波形

figure(10); subplot(2,1,1);n=0:N-1;plot(n,y);

axis([0 150 -2 2]);xlabel('采样点');ylabel('幅值');title('Y(t)的时域波形');grid on;

%进行FFT变换并做频谱图

yy=fft(y,length(y)); %对相加后的信号进行fft变换

longy=length(yy); %Y(t)傅里叶变换后的序列长度

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

magn3=abs(yy); %求Y(t)的幅值

subplot(2,1,2); plot(labelx,magn3); %做Y(t)的频谱图

axis tight;xlabel('频率(Hz)');ylabel('幅值');title('Y(t)的频谱图');grid on;

%求输出信号Y(t)的自相关函数

[c3,lags3]=xcorr(y,'coeff'); %求出Y(t)的自相关序列

figure(11); subplot(2,1,1); plot(lags3/fs,c3); %画Y(t)自相关函数的时域波形

xlabel('T');ylabel('Rx(T)');axis tight;title('Y(t)的的自相关函数');grid on;

%求输出信号Y(t)的双边功率谱

Sy=fft(c3,length(c3)); %对Y(t)的自相关函数进行傅里叶变换

magny=abs(Sy); %求Y(t)双边功率谱幅值

long=length(Sy);

labely=(0:long-1)*16000/long;

subplot(2,1,2); plot(labely,10*log10(magny)); %****画Y(t)的功率谱密度

xlabel('频率(Hz)');ylabel('功率谱(dbW)');axis tight;title('Y(t)的双边功率谱');grid on;