白噪声通过LTI 的仿真

1、实验2 白噪声通过LTI 的仿真1、实验目的了解白噪声通过LTI 系统的原理与处理方法，学会运用Matlab 函数对随机过程进行均值、相关函数和功率谱的估计,并且通过实验分析理论分析与实验结果之间的差异。2、实验原理假定一具有单位方差的抽样序列X(n)的白噪声随机过程X(t)通过一脉冲响应为 的线性滤波器，绘出输入输出信号的均值、方差、相关函数及功率谱密度。设系统冲激响应为h(n)，传递函数 ，或者用Z 变换，结果为 。输入为 X (n)，输出为 ，均值关系： ，若平稳有， 自相关函数关系， ，当是平稳时候，有 题目中假设为白噪声，可以根据白噪声的性质进行理论计算。白噪声的自相关函数， 这里，

2、假设的是零均值和单位方差，于是，而 对应的功率谱， 在这里，由于，a=0.95，可以算出输出信号的方差为， 可以用留数法简单计算出来。下面对输入输出信号的均值、方差、相关函数及功率谱密度分别进行讨论。均值变化输入为白噪声，并且均值为 0，按照理论公式，可得到 下面对实际值进行分析：输入的随机序列，服从标准正态分布。可以用下面的语句产生x = randn(1,500); % 产生题设的随机序列,长度为500 点系统的冲激响应为，可以用下面的语句产生这个冲激信号：b=1;a=1,-0.5; % 设置滤波器的参数，b 为分子系数，a 为分母系数h=impz(b,a,20); % 得到这个系统的冲激响

3、应，就是题设中的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,

4、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); % 进行时间

5、平均，求均值最终值为-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 有差距。如图：自相关函数的理论与实际值理论值为： 在题设中，为白噪声，所以 所以，输出的自相关

6、函数理论值为 可以得到，在零点的值就是 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('自相关

7、函数'); Ry=xcorr(y1,20,'coeff'); % 进行归一化的自相关函数估计，相关长度为20n=-20:1:20;stem(n,Ry,'MarkerFaceColor','red');title('邹先雄实际的自相关函数曲线');功率谱密度函数的理论与实际值对于理论的功率谱密度，可以表示为， 而对于观测数据，可以用功率谱估计的方法得到功率谱密度。首先，采用 Welch 法估计信号的功率谱。它的原理是将数据分成等长度的小段，并且允许数据的重叠，对每段进行估计，再进行平均，得到信号的功率谱。在Matlab 中有

8、专用函数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); % 产生题设的随机序列,长度为

9、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'); % 进行归一化的自相关函数估计，相关长度为20Py,f=pwelch(y1,window,noverlap,Nfft,Fs, 'onesided'); % 估计功率谱密度f=-fliplr(f) f(1:

10、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,&#

11、39;b');title('邹先雄理论功率谱密度曲线');grid on;画出的图形见下图：这是理论的功率谱密度。为了方便显示，将两幅图画在一起，便于比较。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); % 得到这个系统的冲激响应，就是

12、题设中的h（n）y1=filter(b,a,x);y1_mean=mean(y1); % 进行时间平均，求均值y1_var=var(y1); % 进行时间平均，求方差Ry=xcorr(y1,20,'coeff'); % 进行归一化的自相关函数估计，相关长度为20Py,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/m

13、ax(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;结果为：从结果上可以看出来，两者存在着比较大的差距，这是由于输入随机序列

14、的功率谱并不是常数的缘故，也就是输入不是严格的白噪声，所以会出现波动。当随着数据值的增加，拟合的程度会有所改善。3、实验内容假定一具有单位方差的抽样序列X(n)的白噪声随机过程X(t)通过一脉冲响应为 的线性滤波器，利用 matlab 工具绘出输入输出信号的均值、方差、相关函数及功率谱密度。实现的MATLAB代码和结果如下:clear all;x = randn(1,500); % 产生题设的随机序列,长度为500点b=1;a=1,-0.6; % 设置滤波器的参数，b为分子系数，a为分母系数h=impz(b,a,20); % 得到这个系统的冲激响应，就是题设中的h（n）y1=filter(b,

15、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

16、(y1,20,'coeff'); % 进行归一化的自相关函数估计，相关长度为20n=-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.

17、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'); % 进行归一化的自相关函数估计，相关长度为20Py,f=pwelch(y1,window,noverlap,Nfft,Fs, 'onesided'); % 估计功率谱密度f=-fliplr(f) f(1:end); % 构造一个对称的

18、频率，范围是-Fs/2, Fs/2 Py=-fliplr(Py) Py(1:end); % 对称的功率谱figure(3)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形式figure(4)plot(f,Gy,'b');title('邹先雄理论功率谱密度');grid on;下图为随机产生的波形，以此来求相应的均值、方差、相关函数及功率谱密度。自相关函数波形功率谱密度波形如下图截的部分程序及运行结果，其中y2为输出的均值，y3为方差4.实验总结： 随机信号实验综合性强，运用到了随机信号分析，数字信号处理，概率论，matlab等课程知识，让我们尝试综合运