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

本科实验报告实验名称：窄带高斯随机过程的产生

一、实验目的

熟悉窄带随机过程的定义，了解窄带随机过程产生的原理与方法，最后估计实验产生的窄带随机过程的功率谱；掌握具有指定功率谱的随机过程产生方法，并以此产生窄带随机过程。

二、实验原理

（一）窄带随机过程的产生原理

窄带随机过程可以表示为下面的准正弦振荡的形式：

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 k

d

X e f T X t π∞

∞

=-=

=

∑

% 若k X 是零均值的高斯随机变量，那么()X t %也是零均值的高斯随机过程。若

{}k X 是两两正交的序列

()2

2

20()(())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 --+-L L ，其方差为2

2

()k k E X g =。2k g 应

与()0X G kf 成比例，即()20X k G g kf β=，则有

()()2

20

()M

M

M

B

X X

B

M

k k

k k M

M

k G f G kf df E X g β

=-=-=--=

=

=∑

∑

?

∑

即

()()

B

X B M

X

M

k G f G df kf β=--=?∑ 产生步骤： i. 根据要求的时长d T 确定01d f T = ，根据功率谱的带宽确定0

B M f = ii. 计算系数β

iii.

产生21M +个独立的高斯随机变量，

()0(0,)(,1,,0,,1~,)k X N k M X M G M M kf β=--+-L L

iv.

构建时域样本函数[]()02()M

f k M

j i t k k X X i i t X e π?=-=?=

∑

（三）用时域滤波法模拟任意随机过程

功率谱为1的白噪声通过线性系统，输出的是服从高斯分布的，且输出的功率谱为()()2

X G f H f =，因此要产生功率谱为()X G f 的有色高斯噪声，只需设计一个滤波器即可，该滤波器的传递函数应满足()()X H f G f =

三、实验内容

模拟产生一段时长5ms 的窄带高斯随机过程()X t 的样本函数。假定c A t ()与

s A t ()的功率谱密度均为4

1

()()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 ，其功率谱密度要求为()4

1

1(/)

X

f G f f =

+? （1）利用频域法产生()c A t 根据时长5d T ms =，可确定01

200d

f Hz T '=

= 1f kHz ?=是功率谱密度的3dB 带宽，取6B f =?，则030B

M f =

='

计算系数()()6000

604304

300001

1(/1000)199.981

1(/1000)

6200B

X

B M

M 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 β

=--+L

构建时域样本函数[](30

22003)0

j t k k i k X i X e π?=-=

∑

（2）利用时域滤波法产生()s A t

滤波器的传递函数()(

)4

1

1(/)X H f G f f f ==+?

有色高斯的功率谱密度为

()4

4/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 ππππ--?==+?-?-?-?-?g g g g ，

前两个极点与()H f 有关，后两个极点与()H f *有关，故滤波器的传递函数为

()2

/43/4

()()()

j j H f f f f e f f e ππ?=-?-?g g 做傅立叶反变换得系统的冲激响应为

()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 ωωω-=-

则()0

00

-2cos s

n

f n

n s

n H z e

z 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; %时长5ms

Ac=zeros(1,length(t));

for m=-M:M

Ac=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; % 时域长度5ms

fs = 1e5; % 采样频率10kHz

n = 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; %时长5ms

As=zeros(1,length(t));

for m=-M:M

As=As+z(m+M+1)*exp(I*2*pi*m*fo*t);

end;

figure

subplot 211

plot(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)'); %定义图名