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

本科实验报告实验名称：窄带高斯随机过程的产生一、实验目的熟悉窄带随机过程的定义，了解窄带随机过程产生的原理与方法，最后估计实验产生的窄带随机过程的功率谱；掌握具有指定功率谱的随机过程产生方法，并以此产生窄带随机过程。

二、实验原理（一）窄带随机过程的产生原理窄带随机过程可以表示为下面的准正弦振荡的形式：cos X t A t ωτϕτ0()=()[+()]或者表示为同相分量与正交分量的合成：00cos sin c s X t A t t A t t ωω()=()-()其中c A t ()与s A t ()均为低频变化的随机过程，可以通过模拟其分布及功率谱特性来实现窄带随机过程的产生。

（二）用频域法模拟任意随机过程模拟一个时长为d T 的高斯随机过程的一个样本函数()X t , 要求功率谱密度满足指定的形式，先将()X t 进行周期性延拓，并做DFS()0201()j k k f kdX e f T X t π∞∞=-==∑若k X 是零均值的高斯随机变量，那么()X t 也是零均值的高斯随机过程。

若{}k X 是两两正交的序列()2220()(())k k k k X g f k G f E X f g δ=-∞∞=-=∑即可以控制k g 得到期望的功率谱。

假定()(0)X G f B f =>，即()X G f 带限，则{}2k g 为有限项，对应的DFS系数{}k X 也为21M +项0()B M f ⎡⎤=⎢⎥⎣⎦，因此只需产生21M +个相互正交的零均值高斯随机变量{}101,,,,,,M M M M X X X X X --+-，其方差为22()k kE X g =。

2k g 应与()0X G kf 成比例，即()20X k G g kf β=，则有()()220()MMMBX XBMk kk k MMk G f G kf df E X g β=-=-=--===∑∑⎰∑即()()BX B MXMk G f G df kf β=--=⎰∑ 产生步骤： i. 根据要求的时长d T 确定01d f T = ，根据功率谱的带宽确定0B M f = ii. 计算系数βiii.产生21M +个独立的高斯随机变量，()0(0,)(,1,,0,,1~,)k X N k M X M G M M kf β=--+-iv.构建时域样本函数[]()02()Mf k Mj i t k k X X i i t X e π∆=-=∆=∑（三）用时域滤波法模拟任意随机过程功率谱为1的白噪声通过线性系统，输出的是服从高斯分布的，且输出的功率谱为()()2X G f H f =，因此要产生功率谱为()X G f 的有色高斯噪声，只需设计一个滤波器即可，该滤波器的传递函数应满足()()X H f G f =三、实验内容模拟产生一段时长5ms 的窄带高斯随机过程()X t 的样本函数。

假定c A t ()与s A t ()的功率谱密度均为41()()1(/)c s G f G f f f ==+∆，其中f ∆为功率谱密度的3dB 带宽。

按照频域法或时域滤波器法分别产生时长为5ms 的低通过程c A t ()和s A t ()，然后按图示合成()X t ，其中01000/()f Hz π=。

分别画出模拟产生的c A t ()、s A t ()和()X t 的波形。

1. 零均值高斯随机序列的产生产生一段5ms 的零均值高斯随机过程的两个样本()(),c s A t A t ，其功率谱密度要求为()411(/)Xf G f f =+∆ （1）利用频域法产生()c A t 根据时长5d T ms =，可确定01200df Hz T '== 1f kHz ∆=是功率谱密度的3dB 带宽，取6B f =∆，则030BM f =='计算系数()()60006043043000011(/1000)199.9811(/1000)6200BXB MM k X k G f df df f G kf k β=-=---+===+∑∑⎰⎰ 产生21M +个独立的高斯随机变量，4(0,)(30,29,,29,30)1(200/1~000)k X N k k β=--+构建时域样本函数[](3022003)0j t k k i k X i X e π∆=-=∑（2）利用时域滤波法产生()s A t滤波器的传递函数()()411(/)X H f G f f f ==+∆有色高斯的功率谱密度为()44/43/4/43/41()1(/)()()()()j j X j j f f f f f e f f e f f e f G f f e ππππ--∆==+∆-∆-∆-∆-∆，前两个极点与()H f 有关，后两个极点与()H f *有关，故滤波器的传递函数为()2/43/4()()()j j H f f f f e f f e ππ∆=-∆-∆ 做傅立叶反变换得系统的冲激响应为()0002cos t t e h t ωωω-=-式中02f ωπ=∆，输入为白噪声()W t ，则输出为有色高斯过程()()(t)X h W t t =*对()0002cos t t e h t ωωω-=-进行采样离散，采样频率为410s f Hz =，得离散序列()0/002cos /s n f s n n h e f ωωω-=-则()000-2cos snf nn sn H z ez f ωωω∞--=∞=-∑，这将涉及到IIR 滤波器的设计，过程比较繁琐，考虑到()(),c s A t A t 信号的一致性，这里仍然采取频域法2. 将()(),c s A t A t 调制后合成()X t根据题意()00()cos2()sin 2c s X A t f t t A t f t ππ=-，其中01000/f π= 故由仿真可得结果程序如下：%窄带随机过程的产生 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+x.^4); %功率谱密度的函数表达式power=2*df*sum(psd)*0.01; %功率绝对大小 ∑-B~B Gx() s=1./(1+((m*fo)/df).^4); %以fo 为单位，s 即为各个离散点处功率谱密度函数的值 beta=power/sum(s); %系数βs=beta*s; %s=∑Gx(kfo)，而所需的 ，故beta*s 即为所要的功率谱密度%原功率谱密度函数图 -8000Hz - 8000Hz f=[-8:0.01:8]*df;psd0=1./(1+(f/df).^4);%作图显示subplot 211;stem(m*fo,s/fo,'b'); %点线图，横轴为频率，以fo为单位值，纵轴为功率谱相对值hold on;plot(f,psd0,'r'); %连续的功率谱密度axis([-8*df 8*df 0 1.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=[zminus z0 zplus];%做反傅立叶变换，求出时域信号，即窄带随机过程频域法高斯有色信号X(t) Ac(t)t = 0:dt:Td; %时长5msAc=zeros(1,length(t));for m=-M:MAc=Ac+z(m+M+1)*exp(I*2*pi*m*fo*t);end;hold on;subplot 212;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=[zminus z0 zplus];%做反傅立叶变换，求出时域信号，即窄带随机过程频域法高斯有色信号X(t) Ac(t)t = 0:dt:Td; %时长5msAs=zeros(1,length(t));for m=-M:MAs=As+z(m+M+1)*exp(I*2*pi*m*fo*t);end;figuresubplot 211plot(t,As)xlabel('t (msec)');ylabel('As(t)');title('\fontsize{18}\sl时域滤波法产生的As(t)');%% Ac(t) As(t) 调制合成%合成信号X(t)t=0:dt:Td;X = Ac.*cos(2*pi*fc*t*1000) - As .*sin(2*pi*fc*t*1000 );figure;subplot 211;plot(t*1000,real(X),'b');xlabel('t (msec)');ylabel('X(t)');title('\fontsize{18}\sl 合成信号X(t)'); %定义图名。