实验2相关分析法辨识脉冲响应

实验2 相关分析法辨识脉冲响应

实验报告

专业：自动化

班级：

姓名：

日期：年月日

1．实验题目：

2．实验目的

通过仿真实验掌握利用相关分析法辨识脉冲响应的原理和方法。

6．程序代码

function [] = Lab_2()

%--------------------------------------------------------实验题目及初始化定义 disp('实验2 相关分析法辨识脉冲响应') disp(' ')

disp(' 本实验系统输入采用实验1中的M序列，输出受到白噪声污染.') ds = input(' 假设测量白噪声服从正态分布，均值为0，请输入方差 ds = '); disp(' ')

a0=65539;M=2147483647;x=123456;b=10000; U=[];X=[];Y=[];V=[];Z=[];T=[];

P=zeros(12000,6);Rmz=zeros(1,63);Rmz1=zeros(1,63); K=120;T1=8.3;T2=6.2;T0=1;K1=K/(T1*T2); E1=exp(-T0/T1);E2=exp(-T0/T2); Np=63;r=3;a=1;sum1=0;sum2=0;

%-------------------------------------------------生成M序列，获得输入数据u(k) for i=1:6 P(1,i)=1;

U(1)=-2*(P(1,6)-0.5); end

for j=2:12000 for i=2:6

P(j,i)=P(j-1,i-1); end

P(j,1)=mod(P(j-1,5)+P(j-1,6),2); U(j)=-2*(P(j,6)-0.5); end

%--------------------------------------------生成高斯白噪声，获得输入数据v(k) for i=1:12000

x=mod(a0*x+b,M); T(i)=x/M; end

aT=mean(T);vT=var(T);

for i=1:1000 tempsum=0; for j=1:12

tempsum=tempsum+T(12*i-j+1); end

V(i)=sqrt(ds)*(tempsum-12*aT)/sqrt(12*vT); end

aV=mean(V);vV=var(V);

%-------------------------------------- 模拟过程传递函数G（s）获得输出数据z(k) X(1)=0;Y(1)=0;Z(1)=Y(1)+V(1);

for j=2:1000

X(j)=E1*X(j-1)+T1*K1*(1-E1)*U(j-1)+T1*K1*(T1*(E1-1)+T0)*(U(j)-U(j-1))/T0; Y(j)=E2*Y(j-1)+T2*(1-E2)*X(j-1)+T2*(T2*(E2-1)+T0)*(X(j)-X(j-1))/T0; Z(j)=Y(j)+V(j);

end

%-------------------------------------------------------------绘制模拟结果图 disp('** 模拟过程传递函数G（s）获得：')

disp(' 系统白噪声v（k）图形如figure 1所示；') disp([' 该系统白噪声v（k）均值为 ' num2str(aV)]) disp([' 方差为 ' num2str(vV)]) disp(' 系统输出z（k）图形如figure 2所示。') disp(' ')

figure(1) plot(V)

title('系统白噪声v（k）图形')

figure(2) plot(Z)

title('系统输出z（k）图形')

%-------------------------------------------------------------计算互相关系数 for k=0:Np-1 for i=(Np+1):(r+1)*Np

Rmz(k+1)=Rmz(k+1)+1/(r*Np)*U(i-k)*Z(i); end end

disp('** 互相关函数计算时，本实验默认周期数r=3,并获得如figure 3所示Rmz图形。')

figure(3)

plot(1:Np,Rmz)

title('周期r=3时的互相关函数Rmz图形')

choice = input(' 是否查看r=1与r=3互相关函数仿真结果对比图？（按数字1查看，其他均忽略）');

if choice == 1

disp(' 周期r=1与r=3时的互相关函数Rmz图形对比如figure 4所示。') disp(' ')

for k=0:Np-1 for i=(Np+1):2*Np

Rmz1(k+1)=Rmz1(k+1)+1/(1*Np)*U(i-k)*Z(i); end end

figure(4)

plot(1:Np,Rmz,'r',1:Np,Rmz1,'g')

title('周期r=1与r=3时的互相关函数Rmz图形对比') legend('r=3','r=1','Location', 'Best') end

%----------------------------------------------绘制脉冲响应估计图形及仿真图形 c=-Rmz(Np-1);

for k=0:Np-1

g0(k+1)=K/(T1-T2)*(exp(-k*T0/T1)-exp(-k*T0/T2));

g1(k+1)=Np/((Np+1)*a*a*T0)*(Rmz(k+1)+c); g2(k+1)=g0(k+1)-g1(k+1); end

disp('** 脉冲响应理论值g0、估计值g1，以及估计误差值g2如figure 5所示。') disp(' ')

figure(5)

plot(1:Np,g0,1:Np,g1,'r',1:Np,g2,'g') title('系统脉冲响应曲线')

legend('脉冲响应理论值','脉冲响应估计值','脉冲响应估计误差值','Location', 'Best')

%-------------------------------------------------------计算脉冲响应估计误差 for k=1:Np

sum1=sum1+g2(k)*g2(k); sum2=sum2+g0(k)*g0(k); end

dg=sqrt(sum1/sum2);

disp(['** 经计算，脉冲响应估计误差 dg = ' num2str(dg)]) disp(' ') end

7．实验结果及分析

图1 实验结果截图

在实验给定传递函数基础上，实验演示时，程序默认周期数r = 3，并选取白噪声方差ds为0.1，并具体得到以下实验结果。

图2 系统白噪声图形（ds=0.1）

系统输入采用M序列，输出受到服从正态分布的白噪声污染，且白噪声均值为0，方差为0.1，如图2所示。经验证，该白噪声序列均值为4.1529×10-15，方差为0.096993，与理论值几乎一致，符合要求。通过模拟过程传递函数G(s)，获得如图3所示系统输出图像。

图3 系统输出图形（ds=0.1）

计算互相关函数时，默认周期数r = 3，获得如图4所示图形中红色曲线所示。为比较所选周期数对于互相关系数的影响，另取r = 1，获得如图4中绿色曲线所示。进一步分析，在选定某一周期数之后，互相关系数的计算是从第二个周期开始的，目的是等仿真数据进入平稳状态。具体分析将在后文中展开。

图4 互相关函数图形（ds=0.1）

本实验中，补偿量c = 1.8996，仿真获得如图5所示脉冲响应曲线。如图所示，当取白噪声方差ds = 0.1时，脉冲估计值与脉冲理论住吻合的较好，此时，计算的脉冲响应估计误差为0.02937。进一步当ds = 0时，获得如图6所示的互相关函数图形对比图及图7所示的系统脉冲响应曲线图。

图5 系统脉冲响应曲线（ds=0.1）

图6 互相关函数图形（ds=0）

图7 系统脉冲响应曲线（ds=0）

为了进一步比较分析白噪声误差以及周期数选择对于系统相关分析法效果的影响，通过进一步的仿真，得到如表1所示的脉冲响应估计误差数据

表1 脉冲响应估计误差

白噪声方差r = 1 0.027047 0.0305280.0446320.0580940.0691240.078687 0.12152r = 3 0.027066 0.0293770.038930.0485850.0567080.063845 0.096366白噪声方差r = 1 0.17084 0.241230.381570.540080.764441.2098 1.7119 r = 3 0.1343 0.188730.297590.420730.595120.94142 1.3318

如表1数据所示，在选定的某一周期数中，系统脉冲响应估计误差随白噪声方差的增大而增大；而对于白噪声方差相同的情况下，白噪声方差较小时，r = 1的估计误差较小；随着白噪声方差的不断增大，r = 3的估计效果逐渐超越r = 1的情况，且效果越来越好。

实验过程中，遇到一个问题。Matlab中申请的数组下标从1开始，而计算互相关函数时，下标起始为0，周期为63。因而，若忽略这一事实，直接计算下标从1到63的互相关函数，绘制的图像会在最后一个数发生突变从而使图像上扬不光滑，而这一突变的值其实是起始点出的值。为解决这个问题，我将所有互相关函数下标后移一位，即Rmz（1）实为Rmz（0）的值，因而体现在图像上时，图像从k = 0开始绘制，终止于 k = 63，避免出现上扬，保证最终估计误差的准确性。

8．结论

1. 由背景知识可得，根据维纳-霍夫积分方程，只要记录x(t)、y(t)的值，并计算它

们的互相关函数，即可求得脉冲响应函数g(τ)。

2. 本实验利用相关分析法分析脉冲响应，得到脉冲响应的估计误差是随着输入白

噪声方差的增加而增大的，带有白噪声污染的输出z,在白噪声方差为0时与理

想输出y是重合的，即白噪声的方差越小对系统的输出干扰越小。

对于系统采用相关分析法估计，选择周期数越大，估计效果越好。