《利用相关分析法辨识脉冲响应教学总结》由会员分享，可在线阅读，更多相关《利用相关分析法辨识脉冲响应教学总结（12页珍藏版）》请在金锄头文库上搜索。
1、利用相关分析法辨识脉冲响应自 1205刘彬 412511411实验方案设计1.1生成输入数据和噪声用M序列作为辨识的输入信号，噪声采用标准正态分布的白噪声。生成白噪声时，首先利用乘同余法生成 U0,1均匀分布的随机数，再利用U0,1均匀分布的随机数生成标准正态分布的白噪声。1.2过程仿真模拟过程传递函数G(s)，获得输出数据y(k)。 G(s)采取串联传递函数仿真,G(s)盘立JT1T，用M序列作为辨识的输入信号。T21.3计算互相关函数1 (r 1)np兀NpU5其中r为周期数，i Np 1表示计算互相关函数所用的数据是从第二个周期开始的，目的是等过程仿真数据进入平稳状态。1.4计算脉冲响应
2、估计值、脉冲响应理论值、脉冲响应估计误差脉冲响应估计值g(k)2 十 RMz(k) Rmz(W 1)(Np 1)a t脉冲响应理论值g0(k)K k t/T1k t/T2e eT1T21.5计算噪信比脉冲响应估计误差ggo(k) 0(k)Npgo(k)k 12信噪比y(k) y v(k) v2编程说明M序列中，M序列循环周期取Np 26 1 63，时钟节拍t=isec,幅度a 1 , 特征多项式为F(s) s6 s5 1。白噪声循环周期为215 32768 0 G(s)采样时间T 设为 1Sec K 120, T1 8.3Sec, T2 6.2Sec3源程序清单3.1均匀分布随机数生成函数fu
3、n ctio n sita=U(N)%生成N个0 1均匀分布随机数A=179; x0=11; M=2A15;for k=1:Nx2=A*x0;x仁 mod(x2,M);v1=x1/(M+1);v(:,k)=v1;x0=x1;endsita=v;end3.2正态分布白噪声生成函数function v=no ise(aipi)%生成正态分布N(0,sigma)sigma=1; % 标准差for k=1:length(aipi)ksai=0;for i=1:12temp=mod(i+k,length(aipi)+1; ksai=ksai+aipi(temp);endv(k)=sigma*(ksai-
4、6);endend3.3 M 序列生成函数function Np r M=createM(n,a)%生成长度为n的M序列,周期为Np，周期数为rx=1 1 1 1 1 1; % 初始化初态for i=1:ny=x;x(2:6)=y(1:5);x(1)=xor(y(5),y(6);U(i)=y(6);endM=U*a;lenx=length(x);Np=2Ale nx-1; r=n/Np; end3.4 过程仿真函数function y=createy(u,K,T1,T2,T0)n=length(u);K1=K/(T1*T2);E1=exp(-T0/T1);E2=exp(-T0/T2);x(1)
5、=0;y(1)=0;for k=2:nx(k)=E1*x(k-1)+T1*K1*(1-E1)*u(k-1). +T1*K1*(T1*(E1-1)+T0)*(u(k)-u(k-1)/T0;y(k)=E2*y(k-1)+T2*(1-E2)*x(k-1).+T2*(T2*(E1-1)+T0)*(x(k)-x(k-1)/T0;u(k-1)=u(k);x(k-1)=x(k);y(k-1)=y(k);endend3.5 相关函数计算函数function R_Mz=RMz(Np,r,u,z)r=r-1;y=zeros(1,Np);for k=1:Npy(k)=0;for i=Np+1:(r+1)*Np y(
6、k)=y(k)+u(i-k)*z(i);endy(k)=y(k)/(r*Np);endR_Mz=y;end3.5 主函数function og yita=main(time)% 脉冲响应估计误差 og% 噪信比 yitaN=time*63;K=120; T1=8.3; T2=6.2; T0=1; a=1; sita=U(N); %生成0 1均匀分布随机数 v=noise(sita); %利用aipi生成正态分布白噪声Np r u=createM(N,a); % 生成长度为N的M序列 y=createy(u,K,T1,T2,T0); % 利用 M序列驱动，生成 y z=y+v;R_Mz=RMz(
7、Np,r,u,z); % 计算相关函数% 计算脉冲响应估计值g_k=zeros(1,Np);for k=1:Np g_k(1,k)=(R_Mz(1,k)-R_Mz(Np-1)*Np/(Np+1)*a*a*T0); end% 计算脉冲响应理论值Eg=zeros(1,Np);for k=1:Np Eg(1,k)=K/(T1-T2)*(exp(-k*T0/T1)-exp(-k*T0/T2);end% 计算脉冲响应估计误差og=sqrt( norm(Eg-g_kF2/norm(Eg)A2);ov=fangcha(v); % 计算噪声方差oy=fangcha(y); % 计算信号方差yita=sqrt(
8、oy/ov); % 计算信噪比 End3.5 画图函数 1%mainPlot.mfigure(1)for n=4:40og yita=main(n); y1(n)=og;end y1=y1(4:40); plot(4:40,y1);xlabel(周期数);ylabel(脉冲响应估计误差);figure(2)for n=4:40og yita=main(n);y2(n)=yita;end y2=y2(4:40); plot(4:40,y2);xlabel( 周期数 ); ylabel(噪信比);3.5 画图函数 2%mainPlot2.mN=252;K=120; T1=8.3; T2=6.2;
9、T0=1; a=1; sita=U(N);%生成0 1均匀分布随机数v=noise(sita); %利用aipi生成正态分布白噪声Np r u=createM(N,a); % 生成长度为N的M序列 y=createy(u,K,T1,T2,T0); % 利用 M序列驱动，生成 y z=y+v;R_Mz=RMz(Np,r,u,z); % 计算相关函数% 计算脉冲响应估计值g_k=zeros(1,Np);for k=1:Np g_k(1,k)=(R_Mz(1,k)-R_Mz(Np-1)*Np/(Np+1)*a*a*T0); end% 计算脉冲响应理论值Eg=zeros(1,Np);for k=1:N
10、pEg(1,k)=K/(T1-T2)*(exp(-k*T0/T1)-exp(-k*T0/T2);endfigure(1)plot(1:252,y,1:252,z);Legend(不含噪声的输出序列,含噪声的输出序列);figure(2)plot(1:63,g_k,1:63,Eg);Legend(脉冲响应估计值,脉冲响应理论值);4数据记录表1脉冲响应估计值与脉冲响应理论值的比较t1234567脉冲响应估计值0.790.921.021.041.051.010.92脉冲响应理论值2.033.524.595.325.776.026.11t891011121314脉冲响应估计值0.870.800.74
11、0.650.570.500.42脉冲响应理论值6.075.945.745.495.214.914.60t15161718192021脉冲响应估计值0.330.230.170.100.05-0.01-0.06脉冲响应理论值4.293.993.693.403.122.862.62t22232425262728脉冲响应估计值-0.10-0.16-0.19-0.22-0.25-0.29-0.28脉冲响应理论值2.392.181.981.801.631.481.33t29303132333435脉冲响应估计值-0.30-0.31-0.32-0.36-0.37-0.39-0.41脉冲响应理论值1.201.
12、090.980.880.790.710.64t36373839404142脉冲响应估计值-0.44-0.46-0.47-0.46-0.49-0.51-0.52脉冲响应理论值0.580.520.460.410.370.330.30t43444546474849脉冲响应估计值-0.53-0.54-0.55-0.55-0.56-0.54-0.56脉冲响应理论值0.270.240.210.190.170.150.13t50515253545556脉冲响应估计值-0.57-0.57-0.56-0.57-0.57-0.56-0.55脉冲响应理论值0.120.110.100.090.080.070.06t5
13、7585960616263脉冲响应估计值-0.53-0.52-0.53-0.52-0.530.000.61脉冲响应理论值0.050.050.040.040.030.030.035曲线打印图1信噪比随着周期数增大的变化刪H二垃一旦基盘洋图3加入噪声前后的输出序列比较图4脉冲响应理论值与估计值的比较6 结果分析6.1 信噪比脉冲响应计算误差随周期的变化随着周期数的增加，信噪比减小，说明噪声随着周期数的增加变得更强烈， 而计算误差的减小表示周期数的增加使得不确定因素的影响减小， 使得计算结果 与实际更接近。6.2 加入噪声前后的输出序列比较加入噪声前后的变化并不大， 说明噪声对输出序列影响不大， 在第二个周
