白噪声通过LTI的仿真.docx
《白噪声通过LTI的仿真.docx》由会员分享，可在线阅读，更多相关《白噪声通过LTI的仿真.docx（20页珍藏版）》请在冰豆网上搜索。
白噪声通过LTI的仿真
实验2白噪声通过LTI的仿真
1、实验目的
了解白噪声通过LTI系统的原理与处理方法，学会运用Matlab函数对随机过程进行均值、相关函数和功率谱的估计,并且通过实验分析理论分析与实验结果之间的差异。
2、实验原理
假定一具有单位方差的抽样序列｛X（n）｝的白噪声随机过程X（t）通过一脉冲响应为
的线性滤波器，绘出输入输出信号的均值、方差、相关函数及功率谱密度。
设系统冲激响应为h（n），传递函数
，
或者用Z变换，结果为
。
输入为X（n），输出为
，
均值关系：
，
若平稳有，
自相关函数关系，
，
当是平稳时候，有
题目中假设为白噪声，可以根据白噪声的性质进行理论计算。
白噪声的自相关函数，
这里，假设的是零均值和单位方差，于是
，而
对应的功率谱，
在这里，由于
，，a=0.95，可以算出输出信号的方差为，
可以用留数法简单计算出来。
下面对输入输出信号的均值、方差、相关函数及功率谱密度分别进行讨论。
均值变化
输入为白噪声，并且均值为0，按照理论公式，可得到
下面对实际值进行分析：
输入的随机序列，服从标准正态分布。
可以用下面的语句产生
x=randn（1,500）;%产生题设的随机序列,长度为500点
系统的冲激响应为
，可以用下面的语句产生这个冲激信号：
b=[1];a=[1,-0.5];%设置滤波器的参数，b为分子系数，a为分母系数
h=impz（b,a,20）;%得到这个系统的冲激响应，就是题设中的h（n）
输入信号通过线性系统，可以通过卷积的方法，或者用filter函数，
y1=filter（b,a,x）;%用滤波器的方法，点数为500点
y2=conv（x,h）;%通过卷积方法得到，点数为519点
实现的MATLAB代码如下：
clearall;
x=randn（1,500）;%产生题设的随机序列,长度为500点
b=[1];a=[1,-0.5];%设置滤波器的参数，b为分子系数，a为分母系数
h=impz（b,a,20）;%得到这个系统的冲激响应，就是题设中的h（n）
y1=filter（b,a,x）;%用滤波器的方法，点数为500点
subplot（2,1,1）;
plot（y1,'r'）;
Title（'邹先雄——用滤波器的方法，点数为500点'）;
x=randn（1,500）;
y2=conv（x,h）;%通过卷积方法得到，点数为519点
subplot（2,1,2）;
plot（y2,'b'）;
title（'邹先雄——通过卷积方法得到，点数为519点'）;
gridon;
下面画出两者得到波形的区别：
（为了保持一致，对y2的输出取前500点）
两者的输出波形近似一致，可以采用任意一个进行分析。
就采用y1进行讨论，
输出均值为：
y1_mean=mean（y1）;%进行时间平均，求均值
最终值为-0.0973，与理论的零值有一定误差，考虑到输入随机序列的均值不是0，
m_x=mean（x）=-0.0485,按照上面式子，得到m_y=m_xH（0）=2m_x=-0.0970理论值和实际值是非常吻合的。
附运行结果图：
*因为是随机序列，所以每次运行得到y1和m_x的值也是随机的，但是它们始终满足y1=2m_x
方差变化
输入信号方差的理论值就是1，按照公式，输出的功率谱为
下面对实际值进行分析，用y1_var=var（y1）;求得输出均值为1.3598，与理论值的1.3333有差距。
如图：
自相关函数的理论与实际值
理论值为：
在题设中，为白噪声，所以
所以，输出的自相关函数理论值为
可以得到，在零点的值就是1.3333，也就是输出信号的平均功率。
由MATLAb计算的结果为1.3608，这和计算结果非常接近，实际的自相关函数曲线为：
clearall;
x=randn（1,500）;%产生题设的随机序列,长度为500点
b=[1];a=[1,-0.5];%设置滤波器的参数，b为分子系数，a为分母系数
h=impz（b,a,20）;%得到这个系统的冲激响应，就是题设中的h（n）
y1=filter（b,a,x）;%用滤波器的方法，点数为500点
y2=conv（x,h）;%通过卷积方法得到，点数为519点
Y3=var（y1）
title（'自相关函数'）;
Ry=xcorr（y1,20,'coeff'）;%进行归一化的自相关函数估计，相关长度为20
n=-20:
1:
20;
stem（n,Ry,'MarkerFaceColor','red'）;
title（'邹先雄——实际的自相关函数曲线'）;
功率谱密度函数的理论与实际值
对于理论的功率谱密度，可以表示为，
而对于观测数据，可以用功率谱估计的方法得到功率谱密度。
首先，采用Welch法估计信号的功率谱。
它的原理是将数据分成等长度的小段，并且允
许数据的重叠，对每段进行估计，再进行平均，得到信号的功率谱。
在Matlab中有专用
函数pwelch，它的用法是：
[Px,f]=pwelch（X,WINDOW,NOVERLAP,NFFT,Fs,'onesided'）;%window是采用的数据
窗，NOVERLAP是重叠的数目，NFFT是做FFT的点数，Fs是采样频率，onesided是频率取值。
针对本例，可以用下面语句实现：
window=hamming（20）;%采用hanmming窗，长度为20
noverlap=10;%重叠的点数
Nfft=512;%做FFT的点数
Fs=1000;%采样频率，为1000Hz
x=randn（1,500）;%产生题设的随机序列,长度为500点
b=[1];a=[1,-0.5];%设置滤波器的参数，b为分子系数，a为分母系数
h=impz（b,a,20）;%得到这个系统的冲激响应，就是题设中的h（n）
y1=filter（b,a,x）;
y1_mean=mean（y1）;%进行时间平均，求均值
y1_var=var（y1）;%进行时间平均，求方差
Ry=xcorr（y1,20,'coeff'）;%进行归一化的自相关函数估计，相关长度为20
[Py,f]=pwelch（y1,window,noverlap,Nfft,Fs,'onesided'）;%估计功率谱密度
f=[-fliplr（f）f（1:
end）];%构造一个对称的频率，围是[-Fs/2,Fs/2]
Py=[-fliplr（Py）Py（1:
end）];%对称的功率谱
plot（f,10*log10（abs（Py））,'r'）;
title（'邹先雄——实际功率谱密度曲线'）;
gridon;
最后，得到的估计值为
根据上述值，可以计算出理论的功率,由于
可以用下面的语句实现：
w=2*pi*f/Fs;;%转化到数字域上面
H=（1+0.25-2*0.5*cos（w））.^（-1）;%系统函数模平方
Gy=H/max（H）;%归一化处理
Gy=10*log10（Gy）;%化成dB形式
plot（f,Gy,'b'）;
title（'邹先雄——理论功率谱密度曲线'）;
gridon;
画出的图形见下图：
这是理论的功率谱密度。
为了方便显示，将两幅图画在一起，便于比较。
window=hamming（20）;%采用hanmming窗，长度为20
noverlap=10;%重叠的点数
Nfft=512;%做FFT的点数
Fs=1000;%采样频率，为1000Hz
x=randn（1,500）;%产生题设的随机序列,长度为500点
b=[1];a=[1,-0.5];%设置滤波器的参数，b为分子系数，a为分母系数
h=impz（b,a,20）;%得到这个系统的冲激响应，就是题设中的h（n）
y1=filter（b,a,x）;
y1_mean=mean（y1）;%进行时间平均，求均值
y1_var=var（y1）;%进行时间平均，求方差
Ry=xcorr（y1,20,'coeff'）;%进行归一化的自相关函数估计，相关长度为20
[Py,f]=pwelch（y1,window,noverlap,Nfft,Fs,'onesided'）;%估计功率谱密度
f=[-fliplr（f）f（1:
end）];%构造一个对称的频率，围是[-Fs/2,Fs/2]
Py=[-fliplr（Py）Py（1:
end）];%对称的功率谱
Py=Py/max（Py）;%归一化处理
w=2*pi*f/Fs;;%转化到数字域上面
H=（1+0.25-2*0.5*cos（w））.^（-1）;%系统函数模平方
Gy=H/max（H）;%归一化处理
Gy=10*log10（Gy）;%化成dB形式
plot（f,10*log10（abs（Py））,'r',f,Gy,'b'）;
title（'邹先雄——实际功率谱和理论功率谱拟合'）;
legend（'','实际值','理论值'）;
gridon;
结果为：
从结果上可以看出来，两者存在着比较大的差距，这是由于输入随机序列的功率谱并不是常数的缘故，也就是输入不是严格的白噪声，所以会出现波动。
当随着数据值的增加，拟合的程度会有所改善。
3、实验容
假定一具有单位方差的抽样序列｛X（n）的白噪声随机过程X（t）通过一脉冲响应为
的线性滤波器，利用matlab工具绘出输入输出信号的均值、方差、相关函数及功率谱密度。
实现的MATLAB代码和结果如下:
clearall;
x=randn（1,500）;%产生题设的随机序列,长度为500点
b=[1];a=[1,-0.6];%设置滤波器的参数，b为分子系数，a为分母系数
h=impz（b,a,20）;%得到这个系统的冲激响应，就是题设中的h（n）
y1=filter（b,a,x）;%用滤波器的方法，点数为500点
figure
（1）
subplot（2,1,1）;
plot（y1,'r'）;
title（'邹先雄——用滤波器的方法，点数为500点'）;
y2=conv（x,h）;%通过卷积方法得到，点数为519点
subplot（2,1,2）;
plot（y2,'b'）;
title（'邹先雄——通过卷积方法得到，点数为519点'）;
Y2=mean（y1）%进行时间平均，求均值，为理论值
m_x=mean（x）%输出的实际值可以通过2m_x计算
Y3=var（y1）
figure
（2）
Ry=xcorr（y1,20,'coeff'）;%进行归一化的自相关函数估计，相关长度为20
n=-20:
1:
20;
stem（n,Ry,'MarkerFaceColor','red'）;
title（'邹先雄——实际的自相关函数曲线'）;
%实际功率谱密度
window=hamming（20）;%采用hanmming窗，长度为20
noverlap=10;%重叠的点数
Nfft=512;%做FFT的点数
Fs=1000;%采样频率，为1000Hz
x=randn（1,500）;%产生题设的随机序列,长度为500点
b=[1];a=[1,-0.5];%设置滤波器的参数，b为分子系数，a为分母系数
h=impz（b,a,20）;%得到这个系统的冲激响应，就是题设中的h（n）
y1=filter（b,a,x）;
y1_mean=mean（y1）;%进行时间平均，求均值
y1_var=var（y1）;%进行时间平均，求均值
Ry=xcorr（y1,20,'coeff'）;%进行归一化的自相关函数估计，相关长度为20
[Py,f]=pwelch（y1,window,noverlap,Nfft,Fs,'onesided'）;%估计功率谱密度
f=[-fliplr（f）f（1:
end）];%构造一个对称的频率，围是[-Fs/2,Fs/2]
Py=[-fliplr（Py）Py（1:
end）];%对称的功率谱
figure（3）
plot（f,10*log10（abs（Py））,'r'）;
title（'邹先雄——实际功率谱密度'）;
gridon;
%理论功率谱密度
w=2*pi*f/Fs;;%转化到数字域上面
H=（1+0.25-2*0.5*cos（w））.^（-1）;%系统函数模平方
Gy=H/max（H）;%归一化处理
Gy=10*log10（Gy）;%化成dB形式
figure（4）
plot（f,Gy,'b'）;
title（'邹先雄——理论功率谱密度'）;
gridon;
下图为随机产生的波形，以此来求相应的均值、方差、相关函数及功率谱密度。
自相关函数波形
功率谱密度波形
如下图截的部分程序及运行结果，其中y2为输出的均值，y3为方差
4.实验总结：
随机信号实验综合性强，运用到了随机信号分析，数字信号处理，概率论，matlab等课程知识，让我们尝试综合运用理论知识，这样，我们能够更深刻的理解理论知识和各个学科之间的联系。
实验过程中不懂的问题及时向老师提问，或者使用MATLAB自带的help来查看相关信息，对解决遇到的问题都非常有帮助。
实验中经常出现各种错误，往往因为各种细小错误的出现使得整个程序运行错误，比如全角和半角的切换，如果输入时没有考虑，那么编译极有可能不能通过。
在这次实验中我体会到：
实验就是一个发现错误并分析和改正错误的过程。
正因为有错误的出现才显示出实验的魅力。
如果我们只是将代码复制的MATLAB中运行一遍，不用自己改错，不用自己来尝试编写程序，我们的实验会很顺利结束，但是我们也不可能从这次实验中有所收获了。
升级会员 