系统建模与仿真之一种敷衍的阶跃响应面积法求传递函数

题目

给定函数模型G(s)=1/(6*s^3+12*s^2+3*s+1)，用step函数产生阶跃响应，在阶跃响应上加上高斯白噪声，以后利用面积法编程辨识模型的传递函数参数。

代码

clc;
clear;
close all;
%我采用的是直接计算c1、c2、c3的方法，因为仿真式中分子为1，只计算分母系数即可。
%如果a和b都计算因为采用的积分函数原因不知道为什么误差很大。
%仿真时间设置
delta_t=0.001;
time=0:delta_t:61;
varience=input('输入方差：');%方差输入
A=43;x0=41;M=2^18;Num=length(time);N=Num*30;%产生噪声
v=zeros(1,N);
for k = 1:N
    x2 = A*x0;
    x1 = mod(x2,M);
    v1 = x1/M;
    v(k) = v1;
    x0 = x1;
    v0 = v1;
end
B=reshape(v,Num,30);
noise=zeros(1,Num);
for i=1:Num-1
    noise(i)=(sqrt(varience)*(sum(B(i,:))-30/2))/(sqrt(30/12));
end
noise=[0 noise(1:1:Num-1)];
%理论
num=1;
den=[6 12 3 1];
G=tf(num,den);
Y = step(G, time);
Y_noise=Y+noise';
figure(1);
plot(time,Y_noise);
hold on
plot(time,Y,'LineWidth',2);
legend(["响应+噪声","原始响应曲线"]);
xlim([0 50]);
hold off
%积分算C1
inter1=cumsum(delta_t*(1-Y_noise));
c1=inter1(length(time));
den1=[c1 1];
G1=tf(num,den1);
Y1 = step(G1, time);
%积分算C2
inter2=cumsum(delta_t*(Y1-Y_noise));
inter2=cumsum(delta_t*inter2);
c2=inter2(length(time));
den2=[c2 c1 1];
G2=tf(num,den2);
Y2 = step(G2, time);
%积分算C3
inter3=cumsum(delta_t*(Y2-Y_noise));
inter3=cumsum(delta_t*inter3);
c3=trapz(time,inter3);
%3个结果都得到
den3=[c3 c2 c1 1];
G3=tf(num,den3);
Y3 = step(G3, time);
figure(2)
plot(time,Y);
hold on
plot(time,Y3);
xlim([0 50]);
legend(["原始响应曲线","辨识曲线"]);
%估计误差计算
error_sum=0;square_sum=0;
for count=1:length(time)
    error_sum=error_sum+(Y3(count)-Y(count))^2;
    square_sum=square_sum+Y(count)^2;
end
sigma_g=sqrt(error_sum/square_sum);
disp(['阶跃响应估计误差为：',num2str(sigma_g)]);

图像