SysIdentify2.rar_最小二乘估计_模型阶次_模型阶次辨识_系统辨识_阶次辨识

%%%%%%递推最小二乘估计和模型阶次辨识%%%%%% %%输入信号M序列的生成 %M序列生成器参数，取4个节点，特征多项式为F(s)=s^4+s+1 P = 4; Np = (2^P) - 1; %M序列循环周期 dt = 1; %时钟节拍 a = 1; %M序列幅度值 N = 534; %数据总长度 u = zeros(1,N); %定义输入数据向量 X = [1,zeros(1,P-1)]; %M序列初始状态（1,0,0,0） for k = 1:N X0 = X(4)+X(1); if X0 == 2 X0 = 0; end for i = P:-1:2 X(i) = X(i-1); end X(1) = X0; if X(1) == 0 u(k) = a; end if X(1) == 1 u(k) = -a; end u(k) = u(k) + a/Np; %去除直流分量 end %%白噪声的生成 %乘同余法的参数 M = 32768; %2的15次，周期为2的13次 A = 179; x0 = 11; %随机数生成的初值 %生成方差为SigmaV的白噪声 SigmaV = 1; %白噪声方差 v = zeros(1,N); %定义白噪声向量 for k = 1:N ksai = 0; for i = 1:12 x0 = A * x0; x0 = mod(x0,M); ksai = x0/M + ksai; end v(k) = SigmaV * (ksai-6); end %%系统输出数据的计算 %系统参数 lambda = 0.5; y = zeros(1,N); y(1) = 0;y(2) = u(1); for i = 3:N y(i) = 1.5*y(i-1)-0.7*y(i-2)+u(i-1)+0.5*u(i-2); end e = zeros(1,N); e(1) = 0;e(2) = 1.5*e(1)+v(2); for i = 3:N e(k) = 1.5*e(k-1)-0.7*e(k-2)+v(k); end z = y+e; %%系统模型辨识 %采用遗忘因子法，系统参数定义 L = 500; miu = 0.96; F500 = 3.01; Nbeg = 1;Nend = 10; nc = 0; for n = Nbeg:1:Nend sita = zeros(2*n,L+1); dsita = zeros(2*n,L+1); P = 10000*eye(2*n); J = zeros(1,L+1); for k = 1:L h = zeros(2*n,1); for j = 1:n if k-j+1>0 h(j) = -z(k-j+1); h(j+n) = u(k-j+1); end end K = P*h/(h'*P*h+miu); J(k+1) = miu*(J(k)+(z(k+1)-h'*sita(:,k))^2/(h'*P*h+miu)); sita(:,k+1) = sita(:,k)+K*(z(k+1)-h'*sita(:,k)); P = (eye(2*n)-K*h')*P/miu; if n == nc dsita(:,k+1) = [-1.5,0.7,1,0.5]'-sita(:,k+1); end end Lw = sqrt(J(L+1)/(L-2*n)); a = 0;b = 0; for i = 1:n a = a+sita(i); b = b+sita(n+i); end Kw = b/(1+a); if n == nc break; end F = (J(n)-J(n+1))/J(n+1)*(L-2*n-2)/2; if F <= F500 nc = n+1; end end %%计算系统噪信比 Y = 0; for k = 1:N Y = Y + y(k); end meanY = Y/N; %输出数据均值 Y0 = 0; for k = 1:N Y0 = Y0 + (y(k)-meanY)^2; end SigmaY = Y0/N; %输出数据方差 E = 0; for k = 1:N E = E + e(k); end meanE = E/N; %噪声均值 E0 = 0; for k = 1:N E0 = E0 + (e(k)-meanE)^2; end SigmaE = E0/N; %噪声方差 Yita = sqrt(SigmaE/SigmaY); %噪信比 sp = dsita(:,L+1)./[-1.5,0.7,1,0.5]'; s1=0;s2=0; for i = 1:2*nc s1 = s1+sp(i)^2; s2 = s2+dsita(i,L+1)^2; end s3 = 1.5^2+0.7^2+1^2+0.5^2; d1 = sqrt(s1); %参数估计平方相对偏差 d2 = sqrt(s2/s3); %参数估计平方根偏差 Kt = (1+0.5)/(1-1.5+0.7); dk = sqrt(Kw/Kt); %静态增益估计相对偏差 figure (1) %模型参数变化曲线 hold on title('模型参数变化曲线'); plot(sita(1,:),'k-.'); plot(sita(2,:),'b-'); plot(sita(3,:),'r--'); plot(sita(4,:),'g-*'); legend('a 理论值为-1.5','b 理论值为0.7','c 理论值为1','d 理论值为0.5'); ylabel('估计值'); grid on figure (2) %模型参数偏差变化图 hold on title('模型参数偏差变化曲线'); plot(dsita(1,:),'k-.'); plot(dsita(2,:),'b-'); plot(dsita(3,:),'r--'); plot(dsita(4,:),'g-.'); legend('a偏差','b偏差','c偏差','d偏差'); ylabel('估计值'); grid on fprintf('辨识的模型阶次为%d\n',nc); fprintf('噪信比为%f\n',Yita);