《Matlab仿真窄带随机过程.docx》由会员分享，可在线阅读，更多相关《Matlab仿真窄带随机过程.docx（8页珍藏版）》请在taowenge.com淘文阁网|工程机械CAD图纸|机械工程制图|CAD装配图下载|SolidWorks_CaTia_CAD_UG_PROE_设计图分享下载上搜索。
1、随机过程数学建模分析任何通信系统都有发送机和接收机，为了提高系统的可靠性，即输出信噪比， 通常在接收机的输入端接有一个带通滤波器，信道内的噪声构成了一个随机过 程，经过该带通滤波器之后，则变成了窄带随机过程，因此，讨论窄带随机过程 的规律是重要的。一、窄带随机过程。一个实平稳随机过程X（t）,若它的功率谱密度Sx（ 3）具有下述性质： JSx（3）3c，|3|屹+ 30others中心频率为（0c，带宽为3=230，当处时，就可认为满足窄带条件。 若随机过程的功率谱满足该条件则称为窄带随机过程。若带通滤波器的传输函数 满足该条件则称为窄带滤波器。随机过程通过窄带滤波器传输之后变成窄带随机 过程
2、。图1为典型窄带随机过程的功率谱密度图。若用一示波器来观测次波形，则 可看到，它接近于一个正弦波，但此正弦波的幅度和相位都在缓慢地随机变化， 图2所示为窄带随机过程的一个样本函数。图1典型窄带随机过程的功率谱密度图图2窄带随机过程的一个样本函数二、窄带随机过程的数学表示1、用包络和相位的变化表示由窄带条件可知，窄带过程是功率谱限制在g附近的很窄范围内的一个随机 过程,从示波器观察（或由理论上可以推知）:这个过程中的一个样本函数（一个实现） 的波形是一个频率为人且幅度和相位都做缓慢变化的余弦波。写成包络函数和随机相位函数的形式：X(t)=A(t)*coscoct+ 其中:A称作X(t)的包络函数
3、;称作X的随机相位函数。包络随时间做 缓慢变化，看起来比较直观，相位的变化，则看不出来。2、莱斯(Rice)表示式任何一个实平稳随机过程X(t)都可以表示为：X(t)=Ac(t) coscoct-As(t) sin3ct其中同相分量：Ac(t)= X(t) cos(pt= X(t) cosa)ct+ (Osin(oct=LPX(t) *2coscoct正交分量：As(t) = X(t)sin(pt=(O cos3ct X(t) sin3ct= LP|-X(t) *2sincoct(N(。为X()的希尔伯特变换。LPA表示取A的低频部分)。&和As都 是实随机过程，均值为0,方差等于X(t)的方
4、差。三、窄带随机过程仿真建模要求1、用 Matlab 编程仿真窄带随机信号：X(t)=(l+ A(t)*cos(a)ct+(p)+n(t)o 其中 包络A(t)频率为IKHz,幅值为IV。载波频率为：4KHz,幅值为1 V, p是一个 固定相位，n为高斯白噪声，采样频率设为16KHz。实际上，这是一个带有载 波的双边带调制信号。2、计算窄带随机信号的均值、均方值、方差、概率密度、频谱及功率谱密度、 相关函数，用图示法来表示。3、窄带系统检测框图如图3所示。图3窄带系统检测框图4、低通滤波器设计：低通滤波器技术要求：通带截止频率IKHz,阻带截止频率2KHz。过渡带： IKHz,阻带衰减：35D
5、B,通带衰减：1DB,采样频率：44.1 KHz5、计算a点、b点、Ac、As、y的均值、均方值、方差、频谱及功率谱 密度、相关函数，用图示法来表示。四、建模仿真过程及结果(程序见附件)1、根据要求得到X的表达式：x= (1+a) .*cos (2*pi*4000*t+2) +noisy/10;其中：noisy为高斯白噪声，由wgn函数生成，a=cos (2*pi*1000*t),均值：Ex=mean (x),方差：Dx=var(x),计算可得：X(t)的均值为0.0019,X(t)的方差为0.7590。如图4所示，其中蓝色线为X(t)一个样本的时域波形，红色点连成的线为X(t) 的均值，绿色
6、点连成的线为X(t)的方差。图4窄带随机信号时域波形2、求X的概率密度，方法是将最大最小区间分成14等份，然后分别计算各个 区间的个数，如图02中柱形条所示，利用曲线拟合，得到合适的概率密度函数。 为了得到光滑的曲线，利用了多项式拟合，经过测试，9次拟合曲线比较符合要 求，获得的曲线如图5中曲线所示：图5 X的概率分布密度函数3、对X进行频谱分析，在Matlab中，利用fft函数可以很方便得求得X(t)的频 谱，然后用abs和angle函数求得幅值和相位，画出图像如图6所示：图6 X的频谱图4、求X(t)的自相关函数，用xcorr函数求出自相关序列，得到X自相关函数的 时域波形，如图7所示。图
7、7 X自相关函数的时域波形5、对X自相关函数进行fft变换，得到X(t)的功率谱密度，如图8所示。图8 X的功率谱密度6、建立滤波器，建立一个巴特沃思滤波器，对产生的x(t)进行检测。滤波器的幅 度谱和相位谱所示：图9地通滤波器的幅度谱和相位谱7、求A*)的统计特性，Ac(。为X(t)*2cos3ct通过低通滤波器的信号，A。)的均值Eh = -0.4075 4 (带有直流分量)，A。)的均方值是E2h =0.2458Ac(t)的方差 Dh = 0.0798Ac的波形如图10、图11所示：图10 Ac的时域波形图和频谱图图11 Ac(。的自相关函数的时域波形图和Ac的功率谱密度8、求As(t)
8、的统计特性，As为X(t) *2cos(oct通过低通滤波器的信号，As(t)的均值Eh =0.8972 (带有直流分量)，As(。的均方值是E2h=1.1565As的方差Dh = 0.3518As(t)的波形如图13、图14所示：图13 As的时域波形图和频谱图图14 As的自相关函数的时域波形图和As的功率谱密度9、求出 Y(t)的统计特性，Y (t)=Ac(t) coso)ct-As(t) sin(oct,其统计特性如下输出信号Y的均值Eh=-4.401 le-004s输出信号Y的均方值E2h = 3.0280输出信号Y的方差Dh = 3.0303Y的仿真图形如图15、图16所示。图15
9、 丫的时域波形图和频谱图图16 丫的自相关函数的时域波形图和丫的功率谱密度附件:clcfs=l 6(X)0;%设定采样频率N=1300;n=0:N-l;%取的样本点数t=n/fs;%获得以1/16000为时间间隔采样序列noisy=wgn(l,N,0);%产生高斯白噪声a=cos(2*pi*1000*t);%获取A的采样点x=( 1 +a).*cos(2*pi*4000*t+2)+noisy/10;%获取 x(t)的采样点%以I为横坐标画出x的时域图型figure(l); subplot(2,l,l); plot(n,x);axis(l() 140 -3 引);xlabclC采样点);ylab
10、clCX(t)/V);titleC窄带随机信号波形);grid on;%求乂的统计特性并画出来disp(X的均值为);Ex=mean(x); disp(Ex);%求 X均值 hold on; plot(n,Ex,r.);disp(X的方差为);Dx=var(x); disp(Dx);%求 x方差 hold on; plot(n,Dx,g.);%画出X的概率分布函数each=linspace(min(x).max(x),14);%将最大最小区间分成14等份，然后分别计算各个区间的个数nr=hist(x,each);%计算各个区间的个数nr=nr/length(x);%计算各个区间的个数归一化su
11、bplot(2,l,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,nn/r)axis lighi;tiUe(X概率密度分布);xlabel(X);ylabel(P(x);grid on;%对X进行频谱分析Fx=fft(x,N);%对x进行fft变换，在0-16000区间内得到2N-1个频率值magn=abs(Fx);%求 x幅值xangle=angle(Fx); %求 X相位labcl
12、ang=(0:lcngth(x)-l)*16000/lcngth(x);%在 0 16000 区间内求横坐标刻度figurc(2); plot(labclang,magn* 10);%在 0 16000 区间内做频谱和相位图axis(0 16000 -0.5 600); xlabel(频率/Hz);ylabel(幅值频谱图);grid on; %求X的自相关函数 c,lags=xcorr(x,coeff);%求出自相关序列figure(3); subplot(2,l,l); plot(lags/fs,c);%在时域内画自相关函数axis tight; xlabel(T);ylabel(Rx(T
13、);Ele(X(t)的自相关函数)grid on;%求X的功率谱密度 long=Iength(c);Sx=ffl(c,long);lahelx=(0:long-1 )*2*pi;plot_magn= 10*log 10(abs(Sx);subplot(2,l,2); plot(labelx,plot_magn);% 画功率谱密度axis light;xlabel(w);ylabel(Sx(w);tiHe(X(t)的功率谱密度);grid on; %窄带系统检测7.1 =2.*cos(2*pi*4000*t);z2=-2.*sin(2*pi*4000%)；Ac=zl.*x; %滤波后生成AcAs
14、=z2.*x;%滤波后生成As(t)y=Ac.*cos(2*pi*4000*t)-As.*sin(2:1:pi*4000*t);%滤波器设计f_p=1000;f_s=1600;R_p=l;R_s=35;%设定滤波器参数；通、阻带截止频率，通、阻带衰减Ws=2*f_s/fs;Wp=2*Lp/fs;% 频率归一化ln,Wnj=buttord(Wp,Ws,R_p,R_s);% 采用巴特沃思滤波器lb,a=butter(n,Wn);%求得滤波器传输函数的多项式系数figure(4); H,W=freqz(b,a);%求得滤波器传输函数的幅频特性subplot(2J,l); plot(W*fs/(2*p
15、i),abs(H);%在 02pi 区间内作幅度谱lille(低通滤波器幅度谱grid on; subplot(2,l,2); plot(W*fs/(2*pi),angle(H);%在 02pi 区间内作相位谱title。低通滤波器相位谱);grid on;%求Ac滤波后的统计特性mc=filter(b,a,Ac); %上支路通过滤波器Ac(t)dispfAc的均值);Eh=mean(mc)%求Ac的均值disp(Ac的均方值是):E2h=mc*mc7N%求Ac的均方值dispCAc的方差);Dh=var(mc) %求Ac的方差%画Ac的时域波形figure(6); subplot(2,l,l
16、); n=0:N-l; plot(n,mc);axis(0 300-1 l);xlabel(，亲样点);ylabel(扁值XtitlefAc的时域波形);grid on;%画Ac的频谱图yc=fft(mc,length(mc);%对 Ac(。进行 fft 变换longc=length(yc);%求傅里叶变换后的序列长度labclx=(0: longe-1)*16000/longc; magnl=abs(yc);%求 Ac的幅值subplot(2,l,2); plot(labelx,magnl);%画 Ac的频谱图axis tight; xlabel(频率(Hz); ylabel(幅值)；tit
17、le(Ac(t)频谱图)；grid on;%求Ac的自相关函数c 1 Jags 1 =xcorr(mc?coeff);%求出 Ac的自相关序列figure(7); subplot(2,ld); plot(lagsi/fs,cl);%在时域内画 Ac的自相关函数xlabel(T);ylabel(Rx(T);axis tight;title(Ac(l)的自相关函数);grid on;%求Ac的双边功率谱Sac=fft(c 1 ,length(c 1);%对Ac的自相关函数进行傅里叶变换magnc=abs(Sac);%求Ac的双边功率谱幅值long=lcngth(Sac);%求傅里叶变换后的序列长度
18、labelc=(0: long-1)*16000/long;subplot(2,l,2); plot(labelc, 10*log 10(magnc);%画Ac的自相关函数频谱即为Ac的双边功率谱xlabelf频率(Hz);ylabel(，功率谱(dbW);axis tight;title(，Ac(t)的双边功率谱);grid on;%求得As的统计特性ms=filter(b,a,As);%对下支路信号进行滤波得As(t)disp(As的均值)；Eh=mean(ms) %求 As(t)的均值disp(As的均方值是)；E2h=ms*ms7N %求As的均方值dispCAs的方差)； Dh=va
19、r(ms) %求As的方差%作As的时域波形figure(8);subplot(2,1,1); n=O:N-l;plot(n,ms);%画出 As的时域波形axis(0 300 -0.5 2); xlabel(采样点);ylabel(幅值);title(As(t)的时域波形);grid on;%对As进行FFT变换并做频谱图ys=fft(ms,length(ms);%对 As进行 fft 变换longs=length(ys);%求傅里叶变换后的序列长度labe】x=(0: longs-1 )* 16000/longs;magn2=abs(ys);%求 As(t)的幅值subplot(2,l,2
20、); plot(labelx,magn2);%画出 As的频谱图axis tighl;xlabel(濒率(Hz);ylabel(幅值的频谱图);grid on;%求As的自相关函数c2,lags2=xcorr(ms,cocff);%求出 As的自相关序列flgure(9);subplot(2,1,1 );plot(lags2/fs,c2);%画出 As(t)自相关函数的时域波形xlabel(T);ylabel(Rx(T);axis lighl;litle(As的的自相关函数);grid on;%求As(t)的双边功率谱Sas=fft(c2,length(c2);%对As的自相关函数进行傅里叶变
21、换magnc=abs(Sac);%求As(t)的双边功率谱幅值long=length(Sas);%求傅里叶变换后的序列长度labels=(0:long-l )*16000/long;subplot(2,l,2); plot(labclc,10*log 10(magnc);%画 As的自相关函数频谱xlabel(濒率(Hz);ylabelC 功率谱(dbW);axis iighi:iiile(As 的双边功率谱)；%求y的统计特性disp(输出信号Y的均值);Eh=mean(y)%求输出信号Y的均值disp(输出信号Y的均方值上E2h=y*y，/N%求输出信号Y的均方值dispC输出信号Y的方差
22、);Dh=var(y)%求输出信号Y(t)的方差%作输出信号Y的时域波形figure( 10); subplot(2,1,1 );n=0:N-1 ;plot(n,y);axisdO 150-2 2);xlabelC采样点)ylabelf幅值);title(Y(t)的时域波形);grid on;%进行FFT变换并做频谱图yy=fft(y,length(y);%对相加后的信号进行fft变换longy=length(yy);% Y傅里叶变换后的序列长度labelx=(0:longy-1 )* 16000/longy;magn3=abs(yy);%豪 Y的幅值subplot(2,l,2); plot(
23、labelx,magn3);%做 Y的频谱图axis tight;x】abelC频率(Hz);ylabelC幅值);(。的频谱图);grid on;%求输出信号Y的自相关函数c3,lags31=xcoir(y,coeff);%求出 Y的自相关序列figured 1); subplot(2,l,l); plol(lags3/fs,c3);%画 Y自相关函数的时域波形xlabel(T);ylabel(Rx(T);axis tight;title(Y的的自相关函数上grid on;%求输出信号Y的双边功率谱Sy=fft(c3,lcngth(c3);%对Y的自相关函数进行傅里叶变换magny=abs(Sy);%求Y(l)双边功率谱幅值long=length(Sy);labely=(0:long-l)*l 6000/long;subplot(2,l,2); plot(labely, 10*log 10(magny); %*画 Y(t)的功率谱密度xlabel(濒率(Hz);ylabel(，功率谱(dbW);axis tight;tiUe(Y的双边功率谱);grid on;