白噪声通过LTI 的仿真

实验2 白噪声通过LTI 的仿真1、实验目的了解白噪声通过LTI 系统的原理与处理方法，学会运用Matlab 函数对随机过程进行均值、相关函数和功率谱的估计,并且通过实验分析理论分析与实验结果之间的差异。

2、实验原理假定一具有单位方差的抽样序列｛X(n)｝的白噪声随机过程X(t)通过一脉冲响应为的线性滤波器，绘出输入输出信号的均值、方差、相关函数及功率谱密度。

设系统冲激响应为h (n )，传递函数()()jnwn H w h n e∞-=-∞=∑，或者用Z 变换，结果为()()nn H w h n e∞-=-∞=∑。

输入为 X (n )，输出为()()()()()*k Y n h n X n h n k X k ∞=-∞==-∑，均值关系：()()()*Y X m n h n m n =，若平稳有， ()()0Y X m n m H =自相关函数关系， ()()()()121212,**,Y X R n n h n h n R n n =，当是平稳时候，有()()()()**Y X R m h m h m R m =-题目中假设为白噪声，可以根据白噪声的性质进行理论计算。

白噪声的自相关函数，()()2v v r m t σδ=这里，假设的是零均值和单位方差，于是()()x r m t δ=，而()1x R z =()()()1y R z H z H z -=对应的功率谱， ()()()221||12cos jw y x P w P w H e a w a ==-+在这里，由于()()1jw x x P w R e ==，，a=0.95，可以算出输出信号的方差为，()()2102jwyy yr R e dw ππσπ-==⎰可以用留数法简单计算出来。

下面对输入输出信号的均值、方差、相关函数及功率谱密度分别进行讨论。

均值变化输入为白噪声，并且均值为 0，按照理论公式，可得到()00Y X m m H ==下面对实际值进行分析：输入的随机序列，服从标准正态分布。

可以用下面的语句产生x = randn(1,500); % 产生题设的随机序列,长度为500 点系统的冲激响应为()()()0.5nh n u n =，可以用下面的语句产生这个冲激信号：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 代码如下：clear all;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 点'); grid on;下面画出两者得到波形的区别：（为了保持一致，对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，这和计算结果非常接近，实际的自相关函数曲线为：clear all;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'); % 进行归一化的自相关函数估计，相关长度为20n=-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窗，长度为20noverlap=10; % 重叠的点数Nfft=512; % 做FFT的点数Fs=1000; % 采样频率，为1000Hzx = 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('邹先雄——实际功率谱密度曲线');grid on;最后，得到的估计值为根据上述值，可以计算出理论的功率,由于可以用下面的语句实现：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('邹先雄——理论功率谱密度曲线');grid on;画出的图形见下图：这是理论的功率谱密度。

为了方便显示，将两幅图画在一起，便于比较。

window=hamming(20); % 采用hanmming窗，长度为20noverlap=10; % 重叠的点数Nfft=512; % 做FFT的点数Fs=1000; % 采样频率，为1000Hzx = 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(' ','实际值','理论值');grid on;结果为：从结果上可以看出来，两者存在着比较大的差距，这是由于输入随机序列的功率谱并不是常数的缘故，也就是输入不是严格的白噪声，所以会出现波动。

当随着数据值的增加，拟合的程度会有所改善。