窄带高斯随机过程的产生.docx
《窄带高斯随机过程的产生.docx》由会员分享,可在线阅读,更多相关《窄带高斯随机过程的产生.docx(11页珍藏版)》请在冰豆网上搜索。
本科实验报告
实验名称:
窄带高斯随机过程的产生
一、实验目的
熟悉窄带随机过程的定义,了解窄带随机过程产生的原理与方法,最后估计实验产生的窄带随机过程的功率谱;掌握具有指定功率谱的随机过程产生方法,并以此产生窄带随机过程。
二、实验内容
本实验模拟产生一段时长为5ms的窄带高频随机过程X(t)的样本函数。
根据窄带随机过程的理论,X(t)可表示为
其中,Ac(t)和As(t)均为低频的高斯随机过程,因此,要模拟产生X(t),首先要产生两个相互独立的高斯随机过程Ac(t)和As(t),然后用两个正交载波cos2πf0t和sin2πf0t进行调制,如图所示。
假定Ac(t)和As(t)的功率谱密度均为,其中为功率谱密度的3dB带宽。
在3.7节中介绍了有色高斯随机过程的产生,请按照频域法或时域滤波器法分别产生时长5ms的低通过程Ac(t)和As(t),然后按图所示合成X(t),其中f0=1000/π,要求分别画出模拟产生的Ac(t)、As(t)、X(t)的波形。
三、实验原理
(一)、有色高斯随机过程的模拟——频域法
首先将X(t)进行周期延拓,得到一个周期信号,再对周期信号进行傅里叶级数展开,即
由于傅里叶级数是Xk的线性组合,所以,如果Xk是零均值的高斯随机变量,那么也是零均值高斯过程,如果{Xk}是两两正交的序列,则周期信号的功率谱为线谱,即
通过选择gk就可以得到期望的功率谱。
假定Gx(f)是带限的,即
(|f|>B)
那么,{gk2}只有有限项,即{},其中M=[B/f0],[·]表示取整,与此对应的傅里叶级数系数{Xk}也是2M+1项。
因此,只需产生2M+1个相互正交的零均值高斯随机变量{},其方差,并在1式中将时间限定为(0,Td)就可以得到模拟过程X(t)。
应与成比例,即,系数β的选择满足下式:
即
总结如下:
1.根据所需过程的时长Td确定频率f0,并确定傅里叶级数系数的长度M=[B/f0];
2.根据确定β;
3.产生2M+1个独立的高斯随机变量,即
4.构建时域样本函数
其中为任意小的时间间隔。
(二)、有色高斯随机过程的模拟——时域滤波法
功率谱为1的白噪声通过线性系统,输出的是服从高斯分布的,且输出的功率谱为,因此要产生功率谱为的有色高斯噪声,只需设计一个滤波器即可,该滤波器的传递函数应满足
图2:
时域滤波法产生有色高斯噪声的示意图
(三)、窄带随机过程的产生
X(t)=Ac(t)cos2pf0t-As(t)sin2pf0t
用相同估计方法产生两次窄带高斯序列,分别为Ac(t)和As(t),再带入上式与载波相乘并作变换,就得到了窄带随机过程。
四、实验过程及结果
(一)、有色高斯随机过程的模拟——频域法
1.因为Td=5ms,则f0=200Hz;由功率谱密度可知是功率谱密度的3dB带宽,严格来说,该过程带宽是无限的,但频率足够高时,功率谱密度已经很小,取。
故有M=30。
2.计算系数β:
3.产生2M+1个独立的高斯随机变量,即
构建时域样本函数
其中为任意小的时间间隔,这里取。
(二)、窄带随机过程的产生
用同样的方法产生两个独立的高斯随机信号As(t)和Ac(t),再用载波进行调制,即
可得到最终信号。
五、实验结论及分析
1.有色高斯随机过程的模拟——频域法
图3:
模拟产生的具有给定功率谱的高斯随机过程
2.有色高斯随机过程的模拟——时域滤波法
图4:
时域滤波法产生有色高斯噪声
3.窄带随机过程的产生
按图3、4所示方法产生Ac(t)和As(t),并进行载波调制,产生窄带高斯随机过程:
图5:
窄带高斯随机过程
六、心得体会
1.本实验锻炼了我的MATLAB编程能力,学到了随机信号模拟的基本函数;
2.本实验让我对有色高斯噪声有了更深入的认识,学会了模拟产生具有特定频率谱的高斯随机过程;
3.了解了频域法和时域滤波法的原理和思想;
4.锻炼了实践能力和自学能力。
七、代码附录
%窄带随机过程的产生
clc;clear;
%设置参数
fc=1000/pi;%信号的载波频率
dt=0.00001;%采样间隔
Td=0.005;%信号时长
df=1000;%3dB带宽
B=6*df;
fo=1/Td;%中心频率点
M=floor(B*Td);%傅里叶级数系数长度
m=[-M:
M];
I=sqrt(-1);%虚数i
x=0:
0.01:
10;
psd=1./(1+x.^4);%功率谱密度的函数表达式
symsfreal
power=vpa(int(1/(1+(f/1000)^4),-6000,6000),5);%功率绝对大小
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);
%作图显示
subplot211;
stem(m*fo,s/fo,'b');%点线图,横轴为频率,以fo为单位值,纵轴为功率谱相对值
holdon;
plot(f,psd0,'r');%连续的功率谱密度
axis([-8*df8*df01.2]);
xlabel('频率(Hz)');
ylabel('功率谱');
%生成时域信号对应的傅立叶变换
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;%时长5ms
X=zeros(1,length(t));
form=-M:
M
X=X+z(m+M+1)*exp(I*2*pi*m*fo*t);
end;
holdon;
subplot212;
plot(t*1000,real(X),'b');
xlabel('时间(ms)');
ylabel('X(t)');
T=0.005;%时域长度5ms
fs=100000;%采样频率10kHz
n=round(T*fs)+1;%采样点数
t=linspace(0,T,n);
W=randn(1,n);%高斯白噪声
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;%时长5ms
As=zeros(1,length(t));
form=-M:
M
As=As+z(m+M+1)*exp(I*2*pi*m*fo*t);
end;
figure
subplot211
plot(t,As)
xlabel('时间(s)');
ylabel('Z(t)');
%合成信号X(t)
t=0:
dt:
Td;
X=X.*cos(2*pi*fc*t*1000)-As.*sin(2*pi*fc*t*1000);
figure;
subplot211;
plot(t*1000,real(X),'b');
xlabel('时间(ms)');
ylabel('X(t)');
title('\fontsize{18}\sl合成信号X(t)');%定义图名