实验二：窄带高斯随机过程的产生

本科实验报告实验名称：窄带高斯随机过程的产生一、实验目的熟悉窄带随机过程的定义，了解窄带随机过程产生的原理与方法，最后估计实验产生的窄带随机过程的功率谱；掌握具有指定功率谱的随机过程产生方法，并以此产生窄带随机过程。二、实验原理窄带随机过程的产生原理窄带随机过程可以表示为下面的准正弦振荡的形式：X(t)=A(t)cos[®T+9(T)]0或者表示为同相分量与正交分量的合成：X(t)=A(t)coswt-A(t)sinwtc0s0其中A(t)与A(t)均为低频变化的随机过程，可以通过模拟其分布及功率谱cs特性来实现窄带随机过程的产生。用频域法模拟任意随机过程模拟一个时长为T的高斯随机过程的一个样本函数X(t),要求功率谱密d度满足指定的形式，先将X(t)进行周期性延拓，并做DFSX(t)=兰Xej (f=丄)k oTk=—g d若X是零均值的高斯随机变量，那么X(t)也是零均值的高斯随机过程。若k｛X｝是两两正交的序列kG(f)=兰g28(f-kf) (g2=E(|X2))X k 0 k kk=-x即可以控制g得到期望的功率谱。k假定G(f)=0 (|f|>B)，即G(f)带限，则12｝为有限项，对应的DFSTOC\o"1-5"\h\zX X k系数(X｝也为2M+1项(M=—)，因此只需产生2M+1个相互正交的零均值k fL丿°」高斯随机变量(X,X ，…,X，…,X ,X｝，其方差为E(|X|2)=g2。g2应-M -M+1 0 M-1M IkI k k与G(f)成比例，即g2=pG(kf)，则有X0 k X0

JBG(f)df=区E(XF)=送g2=卩区G(f)_BX 1k k X0— k=_M k=_M k=_M即JBG(f)df卩=—B-迓G(kf)X0k=_M产生步骤：根据要求的时长T确定f=丄，根据功率谱的带宽确定M=—TOC\o"1-5"\h\zd 0T fd “0计算系数卩产生2M+1个独立的高斯随机变量，X〜N(0,卩G(kf)) (k=_M，-M+1,…,0,…,M_1,M)k X0构建时域样本函数X[i]=X(iAt)=区kk=_M用时域滤波法模拟任意随机过程功率谱为1的白噪声通过线性系统，输出的是服从高斯分布的，且输出的功率谱为G(f)=H(f)2，因此要产生功率谱为G(f)的有色高斯噪声，只需设X 1 X计一个滤波器即可，该滤波器的传递函数应满足H(f)「G(f)有色高斯噪声三、实验内容模拟产生一段时长5ms的窄带高斯随机过程X(t)的样本函数。假定A(t)与11+11+(f/Af)4，其中Af为功率谱密度的A(t)的功率谱密度均为G(f)=G(f)=s c s3dB带宽。按照频域法或时域滤波器法分别产生时长为5ms的低通过程A(t)和cA(t)，然后按图示合成X(t)，其中f=1000/兀(Hz)。分别画出模拟产生的A(t)、s 0 cA(t)和X(t)的波形。s1.零均值高斯随机序列的产生产生一段5ms的零均值高斯随机过程的两个样本A(tA(t)和X(t)的波形。s1.零均值高斯随机序列的产生产生一段5ms的零均值高斯随机过程的两个样本A(t),A(t)，其功率谱密c s度要求为GX(f)-(1)利用频域法产生A(t)c1根据时长T=5ms，可确定f= =200Hzd 0TdBAf=1kHz是功率谱密度的3dB带宽，取B=6Af，则M=—=30丿0计算系数卩丿化°)〃—'6000 df刀八 —-60001+(f/1000)4 —[99.986G(f)耳k—-M1+(200k/1000)4k=-30独立的高斯随机变量，Xk〜N(0，1+(200k/1000)4)(k=-30,-29,…，29,30)构建时域样本函数x[/]= Xej2兀200k(zAt)kk=-30频谡法产生的加⑴M込I(2)利用时域滤波法产生A(t)s有色高斯的功率谱密度为GX(f)=1+(f/Af)4=_ (Af)4(f-Af•ej^/4)(f-Af»ej3兀/4)(f-Af»e-加/4)(f-Af»e-j3兀/4)前两个极点与H(f)有关，后两个极点与H心)有关，故滤波器的传递函数为滤波器的传递函数H滤波器的传递函数H(f)=、G(f)二\：1+(f/Af)4H(f)(AfH(f)(f一fej/4)(f-Af•ej3兀/4)做傅立叶反变换得系统的冲激响应为h(t)=-2®e-%tcoswt0 0式中w=吊Af，输入为白噪声W(t)，则输出为有色高斯过程0X(t)=h(t)*W(t)对h对h(t)=-2we-w°tcoswt进行采样离散,00米样频率为f=104Hz，得离散序s列h列h(n)=-2w0e-w0n/fscoswn0则H则H(z)=^-2we-w0:01n=-gCOSW0fZ-n'这将涉及到IIR滤波器的设计’过程s比较繁琐，考虑到A(比较繁琐，考虑到A(t),A(t)信号的一致性，这里仍然采取频域法c s对域滤济法产主的AQW2.将A(t),A(t)调制后合成X(t)c s根据题意X(t)=A(t)cos2兀ft—A(t)sin2兀ft，其中f二1000/兀c 0 s 0 0故由仿真可得结果程序如下：%窄带随机过程的产生clc;clear;%设置参数fc=1000/pi;dt=1e-5;Td=5e-3;df=1e3;%信号的载波频率%米样间隔%信号时长%3dB带宽B=6*df;fo=1/Td;%中心频率点M=floor(B*Td);%傅里叶级数系数长度m=[-M:M];I=sqrt(-1);%虚数i%%频域法的功率谱密度图Ac(t)x=0:0.01:10;psd=1./(1+x44);%功率谱密度的函数表达式power=2*df*sum(psd)*0.01;%功率绝对大小E-B~BGx()s=1./(1+((m*fo)/df).A4);%以fo为单位，s即为各个离散点处功率谱密度函数的值beta=power/sum(s);%系数Bs=beta*s;%s=EGx(kfo),而所需的，故beta*s即为所要的功率谱密度%原功率谱密度函数图-8000Hz-8000Hzf=[-8:0.01:8]*df;psd0=1./(1+(f/df).A4);%作图显示subplot211;stem(m*fo,s/fo,'b');%点线图，横轴为频率，以fo为单位值，纵轴为功率谱相对值holdon;plot(f,psd0,'r');%连续的功率谱密度axis([-8*df8*df01.2]);xlabel('frequency(Hz)');ylabel('PSD');title('\fontsize{18}\sl频域法离散采样后的功率谱密度与原功率谱密度');legend('频域法的功率谱密度','原功率谱密度')；%生成时域信号对应的傅立叶变换z0=randn(1)；z0=z0*sqrt(s(M+1))；zplus=sqrt(s(M+2:2*M+1)/2).*(randn(1,M)+I*randn(1,M))；zminus=conj(fliplr(zplus))；z=[zminusz0zplus]；%做反傅立叶变换，求出时域信号，即窄带随机过程频域法高斯有色信号X(t)Ac(t)t=0:dt:Td；%时长5msAc=zeros(1,length(t))；form=-M:MAc=Ac+z(m+M+1)*exp(I*2*pi*m*fo*t)；end；holdon；subplot212；plot(t*1000,real(Ac),'b')；xlabel('t(msec)')；ylabel('Ac(t)')；title('\fontsize{18}\sl频域法产生的Ac(t)')；%%时域滤波法的功率谱密度图As(t)T=0.005； %时域长度5msfs=1e5； %采样频率10kHzn=round(T*fs)+1；%采样点数t=linspace(0,T,n)；W=randn(1,n)；%高斯白噪声%plot(W)；w0=sqrt(2)*pi*df；h=-2*w0*exp(-w0*t).*cos(w0*t);Y=conv(W,h);As=T*Y(1:n);%生成时域信号对应的傅立叶变换z0=randn(1);z0=z0*sqrt(s(M+1));zplus=sqrt(s(M+2:2*M+1)/2).*(randn(1,M)+I*randn(1,M));zminus=conj(fliplr(zplus));z=[zminusz0zplus];%做反傅立叶变换，求出时域信号，即窄带随机过程频域法高斯有色信号X(t)Ac(t)t=0:dt:Td;%时长5msAs=zeros(1,length(t));form=-M:MAs=As+z(m+M+1)*exp(I*2*pi*m*fo*t);end;figuresubplot211plot(t,As)xlabel('t(msec)');ylabel('As(t)');title('\fontsize{18}\sl时域滤波法产生的As(