利用相关分析法辨识脉冲响应.docx
采用相关分析法辨识脉冲响应1生成输入数据和噪声用M序列作为辨识的输入信号，噪声采纳标准正态分布的白噪声。生成白噪声时，首先采用乘同余法生成Uo1匀称分布的随机数，再采用U0,l匀称分布的随机数生成标准正态分布的白噪声。白噪声循环周期为215=32768O2过程仿真模拟过程传递函数G(s),获得输出数据y(k)。G(三)实行串联传递函数仿真，G(三)=。一!-!r,用M序列作为辨识的输入信号。G(三)采样时间”设为工T2ISec,K=120,=8.3Sec,7；=6.2SeC(1)惯性环节Mk)IX_Iy(k)5+1TA其中，T为惯性环节的时间常数，K为惯性环节的静态放大倍数。若采样时间记作4,则惯性环节的输出可写成：y(k)=e°fy(k-1)+TK(-e°f)u(k-1)+Tr(e-°,-1)+仅二D(2)传递函数G(三)仿真(串联)TiT25+1/T1s+lT2令&=/，则G(三)的表达框图为:Kl1湫)5+1/15+lT2y(k)=e°,y(k-1)+TK(I-e°,)u(k-1)+ TK(e°,-D +编程语句可写成：Kl=KQ*TD;F,=EXP(-T1);E2=EXP(-TqT2);x(0)=0;>(0)=0;for(k=;k<252;k+)x(k)=Ei*x(Jt-l)+*C,*(1-Ei)*m(-1)+Ti*7C1*,*(E1-l)+7；*WW-W(-l)/7；;y(k)=E2*y(k-1)+Tj*(I-E2)*x(k-1)+2*(E2-i)+*x)-x-i);u(k,-1)=u(k);x(k-)=x(k);y(A-l)=y(k);)；3、白噪声生成采用乘同余法生成U0,l匀称分布的随机数x/+1=Axi,(mod)W=2U0,lM其中,M=2匕=3276&(循环周期2=)A=179=3(mod8)XO=Il,(奇数)采用U0,l匀称分布的随机数生成正态分布的白噪声v(k)=v-6VV(0,l2)i=l)其中，标准差4分别取Ool,0.5。编程语句for(Z=1;A<252,Z+)ksai=O;for(f=l<12+)xi=A*xZ;H=MoD(就M);ksai=ksai+FLOAT(XiZ历);)；v(k)=Sigma*(ksai-6.0);；4、M序列生成用M序列作为辨识的输入信号，N序列的循环周期取Np=26-1=63,时钟节拍A=ISec,幅度。=1,规律"0”为a,规律“1”为a特征多项式自选，如F(三)=S6/1。生成M序列的结构图编程语句for(A:=k<252,无+)M(0)=M+M(2);ifM(O)=2then(0)=0;for(i=P;i>0;i-)M(z)=M(z-l);)；ifM(O)=Othenu(k)=a;ifM(O)=Ithenu(k)=-a;)；5、相互关函数的计算其中，r为周期数，i=Np+l表示计算相互关函数所用的数据是从其次个周期开头的，目的是等过程仿真数据进入平稳状态。6、c的补偿补偿量C应取一见:(NP-1),不能取一&於(Np)。由于尺虻伏)是周期函数，则有RAnNP)=R族(0)，故不能取一RAnN)。7、计算脉冲响应估量值 脉冲响应估量值虱k)=NP、R2)+c(/Vp÷l)6Tf 脉冲响应理论值go(Z)=-b""7i-叫12X脉冲响应估量误差UoW-u)2(g°(%)2k=8源程序清单8.1匀称分布随机数生成函数functionSita=U(N)%生成N个01匀称分布随机数A=179;xO=lhM=215;fork=l:Nx2=A*xO;xl=mod(x2,M);vl=xl(M+l);v(:,k)=vl;x=x1;endsita=v;end8.2 正态分布白噪声生成函数functionv=noise(aipi)%生成正态分布N(0,sigma)sigma=1;%标准差fork=klength(aipi)ksai=O;fori=l:12temp=mod(i+k,length(aipi)+l;ksai=ksai+aipi(temp);endv(k)=sigma*(ksai-6);endend8.3 M序列生成函数functionNprM=createM(n,a)%生成长度为n的M序列,周期为Np,周期数为rX=Ullll1;%初始化初态fori=l:ny=;x(2:6)=y(l:5);X=Xor(y(5),y(6);U(i)=y(6);endM=U*a;lenx=length(x);Np=2lenx-1;r=nNp;end8.4 过程仿真函数functiony=createy(u,K,T1,T2,T0)n=length(u);K1=K(T1*T2);El=exp(-TOTl);E2=exp(-T0T2);x(D=0;y(D=O;fork=2:nx(k)=El*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)TO;u(k-l)=u(k);x(k-l)=x(k);y(k-l)=y(k);endend8.5 相关函数计算函数functionR_Mz=RMZ(NP,r,u,z)r=r-l;y=zeros(l,Np);fork=l:Npy(k)=O;fori=Np+k(r+l)*Npy(k)=y(k)+u(i-k)*z(i);endy(k)=y(k)(r*N);endR_Mz=y;end8.6 主函数N=252;K=120;T1=8.3;T2=6.2;T0=1;a=1;Sita=U(N);%生成01匀称分布随机数v=noise(sita);%采用aipi生成正态分布白噪声Npru=createM(N,a);%生成长度为N的M序列y=createy(u,K,T1,T2,T0);%采用M序列驱动，生成yz=y+v;R_Mz=RMz(Np,r,u,z);%计算相关函数%计算脉冲响应估量值g_k=zeros(1,Np);fork=1:Npg_k(1,k)=(R_Mz(1,k)-R_Mz(Np-1)*Np(Np+1)*a*a*T0);end%计算脉冲响应理论值Eg=zeros(1,Np);fork=1:NpEg(1,k)=K(T1-T2)*(exp(-k*T0T1)-exp(-k*T0T2);endfigure(1)plot(1:252,y,1:252,z);Iegende不含噪声的输出序列？含噪声的输出序列)figure(2)plot(1:63,g_k,1:63,Eg);Iegend('脉冲响应估量值脉冲响应理论值')；9数据纪录表1脉冲响应估量值与脉冲响应理论值的比较t1234567脉冲响应估量值0.790.921.021.041.051.010.92脉冲响应理论值2.033.524.595.325.776.026.11t891011121314脉冲响应估量值0.870.800.740.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.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.06t57585960616263脉冲响应估量值-0.53-0.52-0.53-0.52-0.530.000.61脉冲响应理论值0.050.050.040.040.030.030.0310曲线打印图1加入噪声前后的输出序列比较图2脉冲响应理论值与估量值的比较作业：如图所示一阶系统，采纳5级移位寄存器产生M序列作为输入信号，辨识该系统的脉冲响应。M序列相关参数设置为：a=2,=15msR1I-IO(MQM)F-T要求：1、每个流程说明；2、流程图；3、源程序；4、仿真波形图;5、PPT和电子档
