利用相关分析法辨识脉冲响应

利用相关分析法辨识脉冲响应自1205 刘彬 412511411 实验方案设计1.1 生成输入数据和噪声用M序列作为辨识的输入信号，噪声采用标准正态分布的白噪声。生成白噪声时，首先利用乘同余法生成U0,1均匀分布的随机数，再利用U0,1均匀分布的随机数生成标准正态分布的白噪声。1.2 过程仿真模拟过程传递函数，获得输出数据y(k)。采取串联传递函数仿真，用M序列作为辨识的输入信号。1.3 计算互相关函数其中r为周期数，表示计算互相关函数所用的数据是从第二个周期开始的，目的是等过程仿真数据进入平稳状态。1.4 计算脉冲响应估计值、脉冲响应理论值、脉冲响应估计误差脉冲响应估计值脉冲响应理论值脉冲响应估计误差1.5 计算噪信比信噪比2 编程说明M序列中，M序列循环周期取，时钟节拍=1Sec，幅度，特征多项式为。白噪声循环周期为。采样时间设为1Sec，3 源程序清单3.1 均匀分布随机数生成函数function sita=U(N)%生成N个0 1均匀分布随机数A=179; x0=11; M=215;for k=1:N x2=A*x0; x1=mod(x2,M); v1=x1/(M+1); v(:,k)=v1; x0=x1;endsita=v;end3.2 正态分布白噪声生成函数function v=noise(aipi)%生成正态分布N(0,sigma)sigma=1; %标准差for k=1:length(aipi) ksai=0; for i=1:12 temp=mod(i+k,length(aipi)+1; ksai=ksai+aipi(temp); end v(k)=sigma*(ksai-6);endend3.3 M序列生成函数function Np r M=createM(n,a)%生成长度为n的M序列,周期为Np，周期数为rx=1 1 1 1 1 1; %初始化初态for i=1:n y=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=2lenx-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)=0;y(1)=0;for k=2:n x(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:Np y(k)=0; for i=Np+1:(r+1)*Np y(k)=y(k)+u(i-k)*z(i); end y(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序列驱动，生成yz=y+v;R_Mz=RMz(Np,r,u,z); %计算相关函数% 计算脉冲响应估计值g_k=zeros(1,Np);for k=1:Npg_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_k)2/norm(Eg)2);ov=fangcha(v); %计算噪声方差oy=fangcha(y); %计算信号方差yita=sqrt(oy/ov); %计算信噪比End3.5 画图函数1%mainPlot.mfigure(1)for n=4:40 og yita=main(n); y1(n)=og;endy1=y1(4:40);plot(4:40,y1);xlabel(周期数);ylabel(脉冲响应估计误差); figure(2)for n=4:40 og yita=main(n); y2(n)=yita;endy2=y2(4:40);plot(4:40,y2);xlabel(周期数);ylabel(噪信比);3.5 画图函数2%mainPlot2.mN=252;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序列驱动，生成yz=y+v;R_Mz=RMz(Np,r,u,z); %计算相关函数% 计算脉冲响应估计值g_k=zeros(1,Np);for k=1:Npg_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 figure(1)plot(1:252,y,1:252,z);Legend(不含噪声的输出序列,含噪声的输出序列); figure(2)plot(1:63,g_k,1:63,Eg);Legend(脉冲响应估计值,脉冲响应理论值);4 数据记录表1 脉冲响应估计值与脉冲响应理论值的比较t1234567脉冲响应估计值0.79 0.92 1.02 1.04 1.05 1.01 0.92 脉冲响应理论值2.03 3.52 4.59 5.32 5.77 6.02 6.11 t891011121314脉冲响应估计值0.87 0.80 0.74 0.65 0.57 0.50 0.42 脉冲响应理论值6.07 5.94 5.74 5.49 5.21 4.91 4.60 t15161718192021脉冲响应估计值0.33 0.23 0.17 0.10 0.05 -0.01 -0.06 脉冲响应理论值4.29 3.99 3.69 3.40 3.12 2.86 2.62 t22232425262728脉冲响应估计值-0.10 -0.16 -0.19 -0.22 -0.25 -0.29 -0.28 脉冲响应理论值2.39 2.18 1.98 1.80 1.63 1.48 1.33 t29303132333435脉冲响应估计值-0.30 -0.31 -0.32 -0.36 -0.37 -0.39 -0.41 脉冲响应理论值1.20 1.09 0.98 0.88 0.79 0.71 0.64 t36373839404142脉冲响应估计值-0.44 -0.46 -0.47 -0.46 -0.49 -0.51 -0.52 脉冲响应理论值0.58 0.52 0.46 0.41 0.37 0.33 0.30 t43444546474849脉冲响应估计值-0.53 -0.54 -0.55 -0.55 -0.56 -0.54 -0.56 脉冲响应理论值0.27 0.24 0.21 0.19 0.17 0.15 0.13 t50515253545556脉冲响应估计值-0.57 -0.57 -0.56 -0.57 -0.57 -0.56 -0.55 脉冲响应理论值0.12 0.11 0.10 0.09 0.08 0.07 0.06 t57585960616263脉冲响应估计值-0.53 -0.52 -0.53 -0.52 -0.53 0.00 0.61脉冲响应理论值0.05 0.05 0.04 0.04 0.03 0.03 0.03 5 曲线打印图1 信噪比随着周期数增大的变化图2 脉冲响应计算误差随着周期数增大的变化图3 加入噪声前后的输出序列比较图4 脉冲响应理论值与估计值的比较6 结果分析6.1 信噪比脉冲响应计算误差随周期的变化随着周期数的增加，信噪比减小，说明噪声随着周期数的增加变得更强烈，而计算误差的减小表示周期数的增加使得不确定因素的影响减小，使得计算结果与实际更接近。6.2 加入噪声前后的输出序列比较加入噪声前后的变化并不大，说明噪声对输出序列影响不大，在第二个周期之后输出序列变得稳定，具有周期变化。6.3 脉冲响应理论值与估计值的比较随着时间的增加，脉冲响应理论值与估计值全部变小，且差值变小，与实验前的理论推导一致。7 实验体会通过本次试验，我不仅更深层次的了解了系统辨识的内容，而且也复习和运用了matlab，为以后的工作实践打了一些基础。- 9 -