作业1如图1.1所示一阶系统，系统传递函数为G(s)=1/(0.1s+1)，如果采用M序列作为输入信号进行系统辨识，采用5级移位寄存器产生M序列作为输入信号，取M序列的时钟脉冲△=15ms，a=2辨识该系统的脉冲响应。

并说明取5级移位寄存器合理与否。

图1.1 一阶RC系统答:1.解题步骤1.初始化参数，设置模型参数，设置产生M序列的各个关键参数；2.利用产生伪随机二进制序列信号的函数getPRBS产生M序列，并作为系统输入；3.通过系统模型，产生系统输出，并将输入输出画在同一图中；4.计算系统输入输出相关函数R xy；5.计算系统脉冲估计值ghat和系统真实脉冲输出g2.程序清单主程序clc; close all; clear all;%% InitializationR = 100e3; % system initialization resistance=100k ohmC = 1e-6; % capacitance=1uftc = R*C; % Time Constant% generate M-sequencen=5;a=2; % Level of the PRBSdel = 15e-3; % clock pulse periodN=2^n-1; % Period of M sequence total=2*N;% Generate m-sequence using the 'getPRBS' function Out = getPRBS(n,a,del,total);% Generate response y(t) of the systems = tf('s');G = 1/(tc*s+1);tf = total*del;tim = 0:del:tf-del;y = lsim(G, Out, tim);%plot input and output of the systemfigurestairs(tim,Out);axis([0 1.0 -2.5 2.5]);hold onplot(tim,y,'r');hold off% Compute Rxy(i*del)sum = 0.0;Rxy = [];iDel_vec=[];for i=1:Ntau=i-1;iDel_vec=[iDel_vec;tau*del];for j=1:Nsum=sum+sign(Out(j))*y(j+tau);endRxy_i = (a/N)*sum;sum=0.0;Rxy = [Rxy; tau Rxy_i];end% Compute ghat & gind = length(Rxy);C = -Rxy(ind, 2);S = (N+1)*a^2*del/N;Rxy_iDel = Rxy(:,2);ghat=(Rxy_iDel+ C )/S;ghat(1)=2*ghat(1);g = 10*exp(-10.*iDel_vec);Result = [Rxy ghat g];%list formdisp(' -------------------------------------------');disp(' i Rxy(iDel) ghat g');disp(' -------------------------------------------');disp(num2str(Result));disp(' -------------------------------------------');调用产生M序列的getPRBS函数function Out=getPRBS(n,a,del,total); %Get PRBS signal%parameters are n registers, a, altitude of m sequence, del, clock pulse, total, the length of m-sequence to be requiredOut = []; % Make Empty Out for storing binary sequence% Initialize n Registersfor i = 1:nR(i) = 1;endif (R(n)==1)Out(1) =-a;endif(R(n)==0)Out(1)=a;endfor i=2:totaltemp=R(1);R(1)= xor(R(n-2),R(n)); %modulo 2 adderj=2;while j<=n %registers shifttemp1=R(j);R(j)=temp;j=j+1;temp=temp1;endif (R(n)==1)Out(i) =-a;endif(R(n)==0)Out(i)=a;endend3. 运行结果图1.2 一阶RC 电路M 序列辨识响应图表1-I 输入输出相关函数、脉冲估计值、真实脉冲输出表-------------------------------------------------------------------------------------------i Rxy(iDel) ghat g--------------------------------------------------------------------------------------------0 -0.0319183 2.86143 10 1 0.529696 10.4984 8.60708 2 0.43794 9.01697 7.40818 3 0.358965 7.74186 6.37628 4 0.290991 6.64436 5.48812 5 0.232485 5.69973 4.72367 6 0.182129 4.88668 4.0657 7 0.138787 4.18689 3.49938 8 0.101482 3.58456 3.01194 9 0.0693728 3.06614 2.5924 10 0.0417366 2.61993 2.2313 11 0.0179498 2.23588 1.9205 12 -0.0025236 1.90531 1.65299--------------------------------------------------------------------------------------------timeA m p i t i t u d e一阶RC 系统输入输出曲线图续表1-I 输入输出相关函数、脉冲估计值、真实脉冲输出表-------------------------------------------------------------------------------------------i Rxy(iDel) ghat g-------------------------------------------------------------------------------------------13 -0.0201452 1.6208 1.4227414 -0.0353123 1.37591 1.2245615 -0.0483668 1.16514 1.0539916 -0.0596028 0.983723 0.9071817 -0.0692738 0.827577 0.78081718 -0.0775977 0.693181 0.67205519 -0.0847621 0.577505 0.57844320 -0.0909286 0.477942 0.49787121 -0.0962361 0.392248 0.42852122 -0.100804 0.31849 0.36883223 -0.104736 0.255006 0.31745624 -0.108121 0.200364 0.27323725 -0.111033 0.153334 0.23517726 -0.11354 0.112855 0.20241927 -0.115698 0.078014 0.17422428 -0.117556 0.0480262 0.14995629 -0.119154 0.0222155 0.12906830 -0.12053 0 0.11109-----------------------------------------------------------------------------选取5级移位寄存器作为输入信号合理。

原因是由5级移位寄存器产生的信号既拥有了较好的随机特性，其周期特性又能与本题中的系统相吻合。

由其产生的信号能够较好的激发出系统的输入输出特性。

作业2已知系统的差分方程为:() 1.5(1)0.7(2)(1)0.5(2)()(1)0.2(2)y k y k y k u k y k k k k εεε--+-=-+-+--+- 其中()k ε是均值为0，7.2σ=并服从正态分布的不相关随机噪声，()u k 采用4级移位寄存器产生的幅度为1的M 序列，1ms ∆=。

数据长度取N=240，请绘出Newton-Raphson 方法求参数的极大似然估计程序流程图并附上Matlab 程序。

1. 程序流程图图2.1 程序流程图2.程序清单clcclose allsigma=7.2;%均方差epsilon=0.001;%迭代终止条件total=16;N=240;V=sigma*randn(total,1); %噪声%M序列产生y1=1;y2=1;y3=1;y4=0;for i=1:15x1=xor(y3,y4);x2=y1;x3=y2;x4=y3;y(i)=y4;if y(i)>0.5,u(i)=-1;else u(i)=1;endy1=x1;y2=x2;y3=x3;y4=x4;end%最小二乘一般算法,产生初始估计值a1,a2,b1,b2;z=zeros(1,total);for k=3:totalz(k)=1.5*z(k-1)+0.7*z(k-2)+u(k-1)+0.5*u(k-2)+V(k);end%给样本系数矩阵H=[-z(2) -z(1) u(2) u(1);-z(3) -z(2) u(3) u(2);-z(4) -z(3) u(4) u(3);-z(5) -z(4) u(5) u(4);-z(6) -z(5) u(6) u(5);-z(7) -z(6) u(7) u(6);-z(8) -z(7) u(8) u(7);-z(9) -z(8) u(9) u(8);-z(10) -z(9) u(10) u(9);-z(11) -z(10) u(11) u(10);-z(12) -z(11) u(12) u(11);-z(13) -z(12) u(13) u(12);-z(14) -z(13) u(14) u(13);-z(15) -z(14) u(15) u(14)];Z=[z(3);z(4);z(5);z(6);z(7);z(8);z(9);z(10);z(11);z(12);z(13);z(14);z(15);z(16)];c=inv(H'*H)*H'*Z;a1=c(1);a2=c(2);b1=c(3);b2=c(4);d1=0.1;d2=0.1;theta=[a1,a2,b1,b2,d1,d2]';%参数估计初值v(1)=0;v(2)=0;d_theta1=zeros(6,1);d_theta2=zeros(6,1);v_da1(1)=0;v_da2(1)=0;v_db1(1)=0;v_db2(1)=0;v_dd1(1)=0;v_dd2(1)=0;v_da1(2)=0;v_da2(2)=0;v_db1(2)=0;v_db2(2)=0;v_dd1(2)=0;v_dd2(2)=0;j=1;bef=zeros(6,1);%执行算法while sum(abs(theta-bef))>epsilon%采集输入输出for i=1:N+3x1=xor(y3,y4);x2=y1;x3=y2;x4=y3;y(i)=y4;if y(i)>0.5,u(i)=-1;else u(i)=1;endy1=x1;y2=x2;y3=x3;y4=x4;endy(1)=0;y(2)=0;V=sigma*randn(N+3,1); %噪声y(1)=1;y(2)=0.01;for k=3:N+3y(k)=1.5*y(k-1)-0.7*y(k-2)+u(k-1)+0.5*u(k-2)+V(k)-V(k-1)+0.2*V(k-2);endJ_d=0;JJ_d=0;a1=theta(1);a2=theta(2);b1=theta(3);b2=theta(4);d1=theta(5);d2=theta(6);for k=3:N+3v(k)=y(k)+a1*y(k-1)+a2*y(k-2)-b1*u(k-1)-b2*u(k-2)-d1*v(k-1)-d2*v(k-2);%求取v(k) v_da1(k)=y(k-1)-d1*v_da1(k-1)-d2*v_da1(k-2);v_da2(k)=y(k-2)-d1*v_da2(k-1)-d2*v_da2(k-2);v_db1(k)=-u(k-1)-d1*v_db1(k-1)-d2*v_db1(k-2);v_db2(k)=-u(k-2)-d1*v_db2(k-1)-d2*v_db2(k-2);v_dd1(k)=-v(k-1)-d1*v_dd1(k-1)-d2*v_dd1(k-2);v_dd2(k)=-v(k-2)-d1*v_dd2(k-1)-d2*v_dd2(k-2);d_theta=[v_da1(k),v_da2(k),v_db1(k),v_db2(k),v_dd1(k),v_dd2(k)]';J_d=J_d+v(k)*d_theta;JJ_d=JJ_d+d_theta'*d_theta;endbef=theta;theta=theta-inv(JJ_d)*J_d;v(1)=v(N+1);v(2)=v(N+2);v_da1(1)=v_da1(N+1);v_da2(1)=v_da2(N+1);v_db1(1)=v_db1(N+1); v_db2(1)=v_db2(N+1);v_dd1(1)=v_dd1(N+1);v_dd2(1)=v_dd2(N+1); v_da1(2)=v_da1(N+2);v_da2(2)=v_da2(N+2);v_db1(2)=v_db1(N+2); v_db2(2)=v_db2(N+2);v_dd1(2)=v_dd1(N+2);v_dd2(2)=v_dd2(N+2); %求取误差error1(j)=-1.5-theta(1);error2(j)=0.7-theta(2);error3(j)=1-theta(3);error4(j)=0.5-theta(4);error5(j)=-1-theta(5);error6(j)=0.2-theta(6);v_error(j)=v(N+2);j=j+1;endtheta%输出估计参数%作图figure(1);plot(error1)hold onplot(error2)hold onplot(error3)hold onplot(error4)hold onplot(error5,'r')hold onplot(error6,'r')title('参数估计误差')xlabel('迭代次数')ylabel('误差')hold offfigure(2);plot(-1.5-error1,'b');hold onplot(0.7-error2,'c')hold onplot(1-error3,'g')hold onplot(0.5-error4,'y')hold onplot(-1-error5,'m')hold onplot(0.2-error6,'r')legend('a1','a2','b1','b2','d1','d2',-1) title('参数估计值变化') xlabel('迭代次数') ylabel('参数') hold off3. 程序运行结果表3-I 模型参数辨识结果被估参数 a1 a2 b1 b2 d1 d2 真实值 -1.5 0.7 1.01 0.5 -1.0 0.2 估计值-1.47110.68170.96870.5875-0.95260.191参数估计值的变化过程曲线如图3.2所示，参数误差估计曲线如图3.3所示。